We have devised a Gaussian constraint algorithm that makes the angular velocity of the director of a liquid crystal, consisting of uniaxial molecules, a constant of motion. By setting the angular velocity equal to zero, a director based coordinate system becomes an inertial frame. This also prevents the director reorientation from interfering with the tails of the time correlation functions. The constraint algorithm consequently makes it possible to correctly evaluate phase functions, time correlation functions, and transport coefficients relative to a director based coordinate system. We have applied the constraint algorithm combined with both equilibrium and nonequilibrium molecular dynamics methods to calculate the thermal conductivity of two nematic liquid crystals consisting of prolate and oblate soft ellipsoid fluids, respectively. In the prolate fluid, the thermal conductivity parallel to the director lambda parallel to parallel to is greater than the thermal conductivity perpendicular to the director lambda perpendicular to perpendicular to. In the oblate fluid, the reverse is true lambda perpendicular to perpendicular to > lambda parallel to parallel to. The constraint algorithm has also been used to calculate the torque exerted by the temperature gradient on the molecules. The prolate ellipsoids are twisted toward the perpendicular orientation relative to the temperature gradient. The oblate ellipsoids are twisted toward the parallel orientation. This phenomenom can be explained by postulating a quadratic coupling between the symmetric traceless order tenser and the temperature gradient. One should also note that in both systems, the molecules orient in such a way that the entropy production is minimized.