The equivalent operator method has been extended to anisotropic intersite electron-electron interactions between a central paramagnetic ion and its ligand neighbours. The relativistic character of both the interaction operator, in Breit-Pauli form, and the electronic wavefunction has been fully incorporated. The tensorial form of the resultant effective hamiltonian has been analysed and selection rules for its matrix elements have been determined. The importance of relativistic crystal fields is illustrated by examples in which the zero-field splitting is calculated for both S-state ions ( e.g. Mn 2+, Gd 3+) and non-S-state ions ( e.g. Sm 2+, Eu 3+, Sm 3+, Er 3+) in various host lattices. In the former group, a much more pronounced relativistic effect on the zero-field splitting is exhibited.