Molecular dynamics simulations of planar elongational flow in a nematic liquid crystal model system based on the Gay–Berne fluid were undertaken by applying the SLLOD equations of motion with an elongational velocity field or strain rate. In order to facilitate the simulation, Kraynik–Reinelt periodic boundary conditions allowing arbitrarily long simulations were used. A Lagrangian constraint algorithm was utilized to fix the director at different angles relative to the elongation direction, so that the various pressure tensor elements could be calculated as a function of this angle. This made it possible to obtain accurate values of the shear viscosities which were found to agree with results previously obtained by shear flow simulations. The torque needed to fix the director at various angles relative to the elongation direction was evaluated in order to determine the stable orientation of the director, where this torque is equal to zero. This orientation was found to be parallel to the elongation direction. It was also noted that the irreversible entropy production was minimal when the director attained this orientation. Since the simulated system was rather large and fairly long simulation runs were undertaken it was also possible to study the cross coupling between the strain rate and the order tensor. It turned out to be very weak at low strain rates but at higher strain rates it could lead to break down of the liquid crystalline order.