The application of the Hartree-Fock method of obtaining approximate
molecular wavefunctions for this choice of Hamiltonian proceeds in a manner
similar to that used to obtain atomic HF wavefunctions.
Once again, the total electronic wavefunction is represented by an
antisymmetrized product of single particle wavefunctions, {
}
For molecular systems, the separation of radial and angular portions of the
wavefunction is no longer possible, and so a numerical representation of
would require that each 
be represented on an extensive
3-dimensional grid of points. The evaluation of two-center integrals over
this grid becomes impractically expensive for systems of even modest size.
Instead, a different sort of finite basis set
must be chosen to represent the single particle wavefunctions.
Traditionally, the view of molecular electronic structure has been built
upon the better-understood framework of atomic electronic structure.
Therefore, it is not surprising that the most popular basis sets
employed have been closely related to the hydrogenic solutions which form
the framework for the understanding of many-electron atoms.
For solving the non-relativistic wave equation, it
would appear that perhaps a good basis set would be comprised of
functions 
which are centered on the constituent atoms of the
molecular system of interest with a functional form given by
where the index a denotes the atomic center to which the spherical
coordinates refer, and the exponential scale factor, 
, may be
optimized for the atomic system of interest.
This choice of basis functions, known as Slater type functions (STF),
makes solution of the multi-center electron repulsion integrals which arise in
the HF equations extremely difficult to solve accurately.
For this reason,
the radial portion of the basis functions is more commonly represented by
Gaussian type functions (GTF) to yield the basis functions
where, again, the 
's may be optimized for the atomic system of
interest.
This choice of basis set does have significant drawbacks, however, as is
evidenced in the ability of this sort of basis set to represent the ground
state of the hydrogen atom. Replacing the exact wavefunction of the
1s electron of ground state H with a single normalized GTF with a
variationally optimized exponent, 
, yields an electronic energy
which is in error by 15%. The use of a linear combination of n normalized
GTF's
with variationally optimized 
and 
yields smaller errors
for higher values of n: 2.8% for n = 2, 0.15% for n = 4, 0.016%
for n = 8[7].
The reason for the inadequacy of the GTF's is directly
associated with their asymptotic behavior, both at the nucleus, and as
. Although some attempts have been made to implement
techniques which employ more compact basis sets constructed solely of STF's,
the extra work associated with optimizing the various parameters associated
with the GTF basis has proven to be far more manageable than the difficulties
which arise with solving electron repulsion integrals involving STF's.
The form of the atomic solutions of the DHF equations points toward the
selection of a four-component basis which is analogous to 2.62.
Because two-electron integrals analogous to the Coulomb and exchange
integrals of the HF equations also arise in the DHF equations, it is
necessary to instead utilize a basis set analogous to 2.63.
Because the large and small components of the DHF wavefunction exhibit
distinct roles in the description of molecular wavefunctions, it is useful
to describe each with a separate basis set.
Furthermore, it will prove to be useful to construct the basis functions
such that the two component large and small basis spinors are all
eigenfunctions of 
, so that the large and small
basis sets are given by
where the two component spinors, {
} are the product of a
function of the spatial coordinates, 
, and one of the
eigenvectors of 
:
In a molecular calculation, linear combinations of Cartesian Gaussian type functions (CGTF) centered on the a particular nucleus are typically employed as spatial basis functions in lieu of the products of radial functions and spherical harmonics which the atomic solutions would appear to suggest as a logical basis. A general CGTF is given by
This choice of basis greatly simplifies the analytic evaluation of integrals
involving these functions via the Gaussian product theorem and various
recursive relationships which exist among CGTFs.[8]
Functions of this type with i+j+k=q may be classified as having angular
momentum q, since linear combinations of each group of 
CGTFs with
the same q value give rise to 2q+1 functions which are directly related to
the solid spherical harmonics of angular momentum q: 
. The transformation to this set of ``pure angular
momentum'' basis functions does not affect the evaluation of integrals
involving these functions, though it can affect the size of the basis set,
since, for q > 1, the number of CGTFs is greater than the number of pure
angular momentum functions.
The molecular single particle 4-spinors,
, may be approximated by a
linear combination of the large and small component basis functions
centered on the molecule's constituent atoms:
where 
and 
are the number of large and small component
basis functions, respectively.
Adopting a shorthand notation, which will allow simplification of
expressions involving the evaluation of operators in the basis, the 4-spinor
becomes j while the two- and four-spinor basis functions
and 
become 
and 
,
respectively, where X is either L or S, and s is either 
or
. The value of the DHF
electronic energy, first presented in ( 2.60), may be expressed in
this new shorthand as
The next task is to derive a set of Fock-like equations from this Energy
expression. If we express 
as a sum of the most general sort of
4-spinor basis function, 
, it becomes possible to
express the DHF energy in terms of these orbitals and their expansion
coefficients, 
, so that
The DHF energy expression may be more compactly expressed through the use of the intermediates
to yield
where the repeated indices are assumed to be summed. This energy expression, along with the imposition of the orthonormality condition
Allows for the construction of the Lagrangian
This Lagrangian may subsequently be subjected to the condition that the
first derivative of 
with respect to the orbital expansion
coefficients, 
, be equal to zero. This method of
variationally solving for the optimal Dirac-Hartree-Fock wavefunction is
completely analogous to the method put forth by Roothaan for the
non-relativistic case[9].
The imposition of the stationary condition gives rise to a set of
linear equations which may be expressed as
where
and the matrix elements of the Dirac-Fock-Roothaan operator 
in the GTO basis are given by
The DC subscript on the Dirac-Fock operator acknowledges that this
operator is based upon the choice of the Dirac-Coulomb Hamiltonian.
For these equations to be useful, however, they must be cast in terms of
integrals involving only the scalar contracted CGTF basis functions,
{ 
}.
The matrix elements of the one-particle Dirac Hamiltonian in the basis are
given by
A significant number of these terms integrate to zero once the explicit form of the basis functions is inserted, and the matrix element are evaluated. For the V integrals, only those matrix elements involving two large component basis functions with the same spin symmetry or two small component basis functions with the same spin symmetry are non-zero, so that
where
For the 
integrals, on the other hand, only the matrix elements between
a large and small basis function will have a non-zero value. Integration
over the spin variables leaves elements which are non-zero for any
combination of spins on the two basis functions:
The overlap integrals may be reduced in a similar fashion:
Finally, reduction of a general four-component two-electron integral to integrals over the scalar spatial basis functions yields
This reduction of the pertinent integrals makes it possible to take advantage of the lore which has accumulated regarding the efficient evaluation of non-relativistic electronic integrals. Since the basic quantities which must be evaluated are the same as for the non-relativistic case, highly optimized non-relativistic integral codes may be adapted to evaluate relativistic integrals.
