Here we correct, extend, and clarify results concerning the spin Hamiltonian ℋS used to describe the ground manifold of Hubbard models for magnetic insulators in the presence of spin-orbit interactions. Most of our explicit results are for a tetragonal lattice as applied to some of the copper oxide lamellar systems and are obtained within the approximation that ℋS consists of a sum of nearest-neighbor bond Hamiltonians. We consider both a "generic" model in which hopping takes place from one copper ion to another and a "real" model in which holes can hop from a copper ion to an intervening oxygen 2p band. Both models include orbitally dependent direct and exchange Coulomb interactions involving two orbitals. Our analytic results have been confirmed by numerical diagonalizations for two holes occupying any of the 3d states and, if applicable, the oxygen 2p states. An extension of the perturbative scheme used by Moriya is used to obtain analytic results for ℋS up to order t2 (t is the matrix of hopping coefficients) for arbitrary crystal symmetry for both the "generic" and "real" models. With only direct orbitally independent Coulomb interactions, our results reduce to Moriya’s apart from some minor modifications. For the tetragonal case, we show to all orders in t and λ, the spin-orbit coupling constant, that ℋS is isotropic in the absence of Coulomb exchange terms and assuming only nearest-neighbor hopping. In the presence of Coulomb exchange, scaled by K, the anisotropy in ℋS is biaxial and is shown to be of order Kt2λ2. Even when K=0, for systems of sufficiently low symmetry, the anisotropy in ℋS is proportional to t6λ2 when the direct on-site Coulomb interaction U is independent of the orbitals involved and of order t2λ2 otherwise. These latter results apply to the orthorhombic phase of La2CuO4.