III. Semiempirical Molecular Orbital Methods
Quantum mechanical methods for the study of molecules can be divided into two categories: ab initio and semiempirical models. Ab initio methods refer to quantum chemical methods in which all the integrals are exactly evaluated in the course of a calculation. Ab initio methods include Hartree-Fock (HF) or molecular orbital (MO) theory, configuration interaction (CI) theory, perturbation theory (PT), and density functional theory (DFT). Ab initio methods that include correlation can have an accuracy comparable with experiment in structure and energy predictions. However, a drawback is that ab initio calculations are extremely demanding in computer resources, especially for large molecular systems. Semiempirical quantum chemical methods lie between ab initio and molecular mechanics (MM). Like MM, they use experimentally derived parameters to strive for accuracy; like ab initio methods, they are quantum-mechanical in nature. Semiempirical methods are computationally fast because many of the difficult integrals are neglected. The error introduced is compensated through the use of parameters. Thus, semiempirical procedures can often produce greater accuracy than ab initio calculations at a similar level.
1. Hartree-Fock Theory
We start with the stationary state Schrodinger equation:
The nonrelativistic, time independent, fixed nuclei Hamiltonian is:
where A and B designate nuclei, and I and j electrons, Z are the atomic numbers, and atomic units are used in eq (1). The many-electron wave-function Ψ is approximated as a product of one-electron functions φ(I) – the orbital approximation.
where A is the antisymmetrizer, ensuring the wavefunction obeys the Pauli exclusion principle, and O(s) is a spin projection operator that ensures that the wavefunction remains an eigenfunction of the spin-squared operator, S2.
Each molecular orbital is then expanded as a linear combination of atom orbitals, or basis functions, χu:
Utilizing the variational principle <Ψ| H |Ψ>/<Ψ|Ψ> ³ E(exp), Ψ is varied with respect to C, and an eigenvalue equation that yields the molecular orbitals and energies is obtained.
where f is an effective one-electron Fock operator with matrix elements given by
The one-electron matrix elements are give by
and two-electron integrals by
P is the first-order Fock-Dirac (or one-particle) density matrix
where ni is occupation number, 2, 1, or 0.
The Fock equations can be solved by matrix diagonalization of the matrix equations
where S is the overlap matrix and C is a square matrix with the ith column being the MO coefficients of the ith molecular orbital.
The following steps are common to all MO procedures:
1. Calculate the integrals needed to form the Fock matrix F.
2. Calculate the overlap matrix S.
3. Diagonalize S. W+SW = D
4. Form S-1/2 = WD-1/2W+.
5. Form the F matrix from equation (6).
6. Form F’ = S-1/2FS-1/2.
7. Diagonalize F’ for the MO eigenvalues E, V+F’V = E.
8. Back transform V to obtain the MO coefficients C, C = S-1/2V.
9. Form the density matrix P.
- Check P for convergence. If not converged, repeat steps 5-10 until self-consistent electronic field (SCF) is obtained.
2. Approximate MO theories.
Pople and coworkers introduced a series of zero-differential overlap approximations in 1965. The hierarchy of integral approximations and the effect on the Fock matrix formation are given below.
Complete Neglect of Differential Overlap (CNDO) method.
In this approximation, all integrals involving different atomic orbitals, χu, are ignored; <uv|st> = δ(u,v) δ(s,t)<uv|st> = <uu|ss>. Thus, the overlap matrix becomes the unit matrix, S = 1, and the Fock matrix elements:
Parameterization and implementation scheme of the CNDO method was also proposed by Pople:
* All two-center two-electron integrals between a pair of atoms are set equal: <uv|st> (= <uu|ss>) = γAB, where γAB is a function of atoms A and B, depending on the interatomic distance RAB.
* The off-diagonal one-electron, or resonance integrals, are set proportional to S:
where βA is a parameter that depends only on the nature of atom A.
* Electron-core attraction interactions for a given pair of atoms are set equal: <u| VB|v> = δ(u,v) VAB, with VAB = <uA| ZB/R | vA >.
* The one-center core integral Uuu is approximated by assuming that the same atomic orbital for the atom are appropriate for the positive ion: Uuu = -Iu – (ZA -1)γAA. This represents the energy required to remove an electron from atomic orbital χu in the fully ionized atom. An alternative would have been to derive Uuu from the electron affinity – the energy gained when a fully ionized atom captured an electron. CNDO/2 utilized the average of the two.
* A minimum basis set of valence orbitals was used chosen, using Slater type orbitals (STO).
where Y(θ,φ) is the real spherical harmonics.
These approximations reduce the Fock equations from the full Roothaan form to
However, it was soon realized that the electron-core attraction and the electron-electron repulsion are imbalanced in the CNDO ZDO approximation, and the appropriate two-center two-electron integral was used to approximate the electron-core attraction.
The total electronic energy was calculated from
Intermediate Neglect of Differential Overlap (INDO).
In INDO, the constraint in CNDO that monocentric two-electron integrals be equal was removed. Consequently, for atoms with s and p orbitals, there five unique two-electron one-center integrals: <ss|ss> = gss, <ss|pp> = gsp, <pp|pp> = gpp, <pp|p’p’> = gp2, p ¹ p’, and <sp|sp> = hsp. The diagonal Fock matrix element in INDO is
in which eq (15) was used in place of the electron-core attraction integral.
Since INDO and CNDO execute on a computer at about the same speeds and INDO contains some important integrals neglected in the CNDO method, INDO performs much better than CNDO especially in prediction of molecular spectral properties. One most successful program is that developed by Ridley and Zerner in 1973. The model, CNDO/S, was parameterized by carrying out CI-singles (CIS) calculations, and based one atomic spectroscopic data. In general, the INDO/S model reproduces the excitation energies of transitions below 40,000 cm-1 within 2000 cm-1. Another well-known program at this level is Dewar’s MINDO/3 model (1975). Here, all quantities that entered the Fock matrix and the energy expression were treated as free parameters. As a result, the MINDO/3 model achieved an impressive predictive power.
Neglect of Diatomic Differential Overlap (NDDO).
The NDDO model can be derived from the replacement:
This leads to the following Fock matrix elements:
Here, all two-electron two-center integrals involving charge clouds arising from pairs of orbitals on an atom were retained. This model allowed lone-pair lone-pair repulsions to be represented. Since there are four valence orbitals (for an sp atom), there are 10 unique pairs, giving rise to 100 integrals. For (spd) basis, 2025 two-electron integrals are needed compared to INDO’s 4.
The NDDO model is the only model so far that really relates to an actual basis set. The first practical NDDO model was introduced by Dewar and Thiel in 1977. Thiel has continued recently to expand the original MNDO model to include d orbitals. The model was parameterized on molecular geometries, heats of formation, dipole moments, and ionization potentials. A feature in the MNDO model is that the core-core repulsion term was made a function of the electron-electron repulsion integrals:
A major fatal problem in the MNDO model is its inability to reproduce hydrogen bonding interactions due to a spurious repulsion at just outside chemical bonding distances. The solution to this problem was to assign a number of spherical Gaussians to mimic the correlation effects. This is the Austin Model 1 (AM1):
This increased the number of parameters from the original 7 to 13-16 per atom. The AM1 model is undoubtedly a much improvement over the MNDO model for a wide range of properties. Stewart’s Parametric Method Number 3 (PM3) model treat all quantities that enter the Fock matrix as free parameters, while MNDO and AM1 derive the one-center two-electron integrals from atomic spectroscopy.
3. Integral evaluation
While in CNDO and INDO, most electron repulsion integrals were neglected, and those that are left had a common value (γAB) for a given pair of atoms. On the other hand, in NDDO, there are 22 non-vanishing integrals between first-row atoms. The NDDO electron repulsion integrals are determined in semiempirical programs in terms of multipole-multipole interactions, a method introduced by Dewar and Thiel in 1977. The atomic orbitals used in this treatment are the Slater-Zener orbitals, which are products of a radial function Rnl(r) and a normalized real spherical harmonic Ylm(θ,φ), with quantum numbers n, l, and m.
where ξ is the orbital exponent of the Slater AO, and Pl|m|(cos θ) is an associated Legendre function. Within MNDO, AM1, and PM3, only s and p basis AO were used (l = 0, 1). The multipole moments Mlm of a charge distribution ρ(r, θ, φ) is defined by
The notation of multipole moments may be compared with the common notations: M00 = q;
M10 = μz; M11 = μx; M1-1 = μy; etc.
The semiempirical formalism for the NDDO repulsion integrals were derived by first assuming that the two interacting charge distributions do not overlap. The solution to this problem is to introduce normalized real spherical harmonics and a bipolar expansion for 1/r12.
After substituting the 1/r12 expansion into the two-electron integral expression and comparison with the definition of multipoles, the two electron repulsion integrals can be expressed by
where the coefficients depends on l, m and dlm. Eq (24) has the correct asymptotic behavior for RAB ® ¥. For short distances, the assumption of zero diatomic overlap is no longer valid. Therefore, eq (23) must be modified semiempirically to ensure the behavior at short distances. The modification was provided by applying the Klopman formula for monopole-monopole interactions between charge distributions uu and vv:
The nonvanishing multipoles Mlm in the NDDO integrals are represented by configurations of 2l point charges of magnitude e/2l. The interactions between the multipoles are calculated by applying the Klopman formula to each point charge pair.
There are two values needed to be determined in this treatment, the distance separating the point charges for each multipole configuration and the parameters ρ0. The first was determined analytically by comparing the exact expression for the multipoles defined by eq (23), and the results from the point charges. These distances depend on the orbital exponents. ρ0 were chosen to yield the correct one-center limit for the monopole-monopole, dipole-dipole, and quadrupole-quadrupole interactions.
Consider, for example, the integral <spz|spz>. The only nonvanishing multipole for a charge distribution of sp is the dipole term μz. Therefore, <spz|spz> = [ μz, μz ].
In hybrid QM/MM application, the one-electron integral between a charge distribution uv in the QM region and a point charge qmm in the MM region is approximated, consistent with the MNDO/AM1/PM3 procedure, by the two-electron integrals: <u|qmm/Rmm|v> = <uv|smmsmm>. The relevant integrals, thus, include <ss|ss>, <pπ pπ|ss>, <pzpz|ss>, and <spz|ss> types. In ab initio/DFT hybrid QM/MM models, the one-electron integrals are, of course, computed analytically as do all other integrals (Freindorf, J. Comput. Chem. 17, 386 (1996)).
4. Computer programs.
The MOPAC package developed by J. J. P. Stewart was a public domain program up to MOPAC 7, and MOPAC 93. It contains the MINDO/3, MNDO, AM1, and PM3 models, which allows a variety of properties to be computed. There is also a commercial version of MOPAC as well as AMPAC. The latter contains the newest SAM1 model. The most familiar ab initio package is the Gaussian series of programs. More than a dozen of other similar programs are available. A public domain ab initio program, the GAMESS program, is also available. For spectroscopy calculations, Zerner’s INDO/S would be an excellent choice. The most sophisticated ab initio package on spectroscopy would be the MOLCAS from Roos.