A constitutive equation for liquid‐crystal polymer solutions is developed from phase space kinetic theory by modeling the polymer solution as a collection of rigid dumbbells suspended in a Newtonian solvent and subjected to anisotropic hydrodynamic drag. The isotropic to nematic phase transition is predicted by including polymer–polymer interaction through the Maier–Saupe mean‐field interaction potential. The diffusion equation for the configurational distribution function is the same as that used by Doi in his theory for concentrated solutions of rod‐like polymers. However, the expression obtained for the stress tensor contains an additional term involving the velocity gradient. The theory predicts the transitions to two nematic phases with prolate and oblate orientations, respectively, that evolve from the isotropic state as a transcritical bifurcation in the polymer concentration. The effects of a variety of deformations on the equilibrium phase behavior are presented. Multiple steady‐state nematic phases are predicted for both steady‐state shear and elongational flows; distinctly different states are identified for uniaxial and biaxial elongational flows. Material functions in shear and shearfree flows are determined, and the effect of anisotropic hydrodynamic drag on the rods is studied. Asymptotic behavior of the viscosity at high shear rates shows the importance of the additional term in the stress tensor expression.