By means of a variational procedure spin-spin coupling constants between nuclei in molecules are evaluated. The Hamiltonian is divided into a part depending on electron spin operators and a part involving only electron orbital angular momenta. A variational function for the first part is constructed by expansion of the wave function in terms of electron spin operators, the first two terms of which are the electron spin densities at the nuclei of interest and succeeding terms contain electric moments of the unperturbed electronic charge distribution of the ground state of the molecule. Retaining terms up to and including the electric dipole moment, a two parameter variational function is used for minimization and an expression for the coupling constant is derived. This result is applied to the hydrogen molecule using valence bond, molecular orbital, and James-Coolidge type wave functions. C13, N14, O17, and F19 coupling constants with protons in CH4, NH4+, H2O, and HF, respectively, are computed with valence bond and molecular orbital functions. The relation of these results to the conventional perturbation theory approach is discussed, particularly in regard to effective ``average energy denominators'' which can be computed by the present procedure. A variational method for electron orbital contributions is briefly considered as well as the calculation of ``long range'' coupling constants.