Matrix elements of a spin-correlated crystal field (SCCF) operator have been evaluated using exact calculations for two-body double tensor operators and compared with earlier approximate calculations for Gd3+ in lanthanum ethylsulphate. In this latter case reductions to products of matrix elements of one-body operators were used. These approximations are shown to be accurate for the determination of zero-field splittings (ZFS) only to within 40%. Nonstandard rank combinations resulting from relativistic effects have also been incorporated which are shown to yield negligible contributions to the ZFS. The two-body form of the SCCF has been made relativistic as well in order to investigate the importance of new tensorial terms for the calculations of the ZFS. Again, their contribution has been shown to be very small.