Lecture Notes

In this section we discuss some practical issues with regard to the implementation of a molecular mechanics (MM) or classical force field into quantum mechanical (QM) calculations. Although, the procedure appears to be general for both ab initio and semiempirical QM methods, our focus will be on the latter because from all practical considerations it is the only viable method at the present time and in the near future that can be applied to molecules of biological interest. Perhaps ab initio/DFT methods are only applicable to solute molecules of less than a few (3-5) non-hydrogen atoms with a modest 3-21G basis set in molecular dynamics and Monte Carlo simulations. An attractive alternative to semiempirical methods is the empirical valence bond (EVB) approach that has been extensively used by Warshel and his coworkers. Here, I will only try to briefly describe it since Dr. J.-K. Hwang of the National Tsing Hua University is an expert in this field.

The key in a hybrid QM/MM potential is to introduce the solvent electric field into quantum mechanical calculations of the solute molecule on the fly during a computer simulation such that the molecular wavefunction of the solute will be polarized by the dynamic change of the surrounding solvent molecules. The method was first described by Warshel in 1976 and a detailed prescription for MO calculations was presented by Field, Bash and Karplus in 1990. The idea was to divide a condensed phase system into a QM region and a MM region, plus an appropriate boundary treatment to mimic the bulk effects. Consequently the effective Hamiltonian of the system is written as follows:

(1)

The van der Waals term is to ensure that the qm and mm systems will not get too close because of a lack of electronic structural description of the solvent mm system. It turns out that the most significant part is to account for the short range electron repulsions; however, for convenience as well as for inclusion some dispersion interactions, a Lennard-Jones term is typically used for H^vdW.

(2)

where S and M are the number of solvent (MM) and solute (QM) atoms, and σ and ε are empirical parameters. It should be emphasized that these parameters must be carefully examined for a given QM (model, basis set, and level of theory) and MM (force fields) combination. However, on the other hand, it gives us the opportunity to optimize the performance of the hybrid QM/MM potential. Therefore, even though the semiempirical methods are generally not very good in describing hydrogen bonding interactions, a hybrid QM/MM model can yield good results thanks to the parameter optimization. There are two publications describing the vdW parameterization, one for semiempirical and one for ab initio HF/3-21G calculations.

In the hybrid QM/MM potential, solute-solvent (or qm/mm) interactions is given by

(3)

where Ψ is the wavefunction of the solute molecule in solution, which minimizes the energy of the Hamiltonian of eq (1). The HF-SCF calculations were carried out using the modified Fock matrix:

(4)

where F^o is the Fock matrix of the isolated "solute" molecule (in vacuo) and the one-electron integral due to the presence of the MM electrostatic potential - resulting from point charges is:

(5)

In ab initio method, eq (5) will be evaluated analytically for all basis functions, whereas in semiempirical method with NDDO approximation, the last term of eq (5) is used and χ_u and χ_v are orbitals on the same atom.

With the modification of the Fock matrix, eq (4), standard SCF calculations are carried out, yield the total energy of the solute (E_qm) plus solute-solvent electrostatic interaction energies (E^el_qm/mm). The solute-solvent van der Waals interaction (E^vdW_qm/mm) and solvent energies (E_mm) are determined separately by force field, or empirical potentials. In Monte Carlo simulations, the total energy (E_tot = E_qm + E_mm + E^el_qm/mm + E^vdW_qm/mm) is all that is needed to carry out the Metropolis procedure. In molecular dynamics simulations, the first derivative of the E_tot is needed to obtain forces to integrate Newton’s equations of motion.

(6)

Evaluation of the MM forces was straightforward, while procedures for computing the QM gradients are available in QM packages.

A procedure for carrying out Monte Carlo simulations with the hybrid QM/MM potential is given in J. Phys. Chem. 96, 537 (1992), molecular dynamics in J. Comput. Chem. 11, 700 (1990), empirical valence bond potential in Chem. Rev. 93, 2523 (1993), and density functional in J. Phys. Chem. 97, 8050 (1993); 97, 11868 (1993).

A straightforward extension of the hybrid HF-QM/MM model is to use the method of configuration interaction (CI). This allows the correlation energy to be obtained in ab initio calculations, and solvent effects on excited states to be investigated in fluid simulations. As in gas phase CI calculations, we represent the exact wavefunction as a linear combination of N-electron trial functions and use the linear variational method. Suppose for an even number of electron system, the molecule is adequately represented by a closed-shell restricted HF determinant, |Ψ₀>. If we have solved Roothaan’s equations in a finite basis set and obtained 2K spin orbitals {φ_i}, the determinant formed from the N lowest energy spin orbitals is |Ψ₀>. In addition to |Ψ₀>, a large number of N-electron determinants from the 2K spin orbitals can be formed. These other determinants together form a basis to expand the exact many-electron wavefunction |Φ₀>:

(7)

where the second and third terms are the singly excited and doubly excited determinants, respectively. By forming a full CI matrix and finding its eigenvalues, we obtain an upper bound to the ground and excited states energies of the system. Since the effective Hamiltonian of equation (1) was used, the solvent effects on the excited states are naturally included in the CI calculation. We can write down the total ground state energy of the solution system in a hybrid QM-CI/MM treatment by

(8)

Here I have emphasized that the total energy can be decomposed into three terms, the energy of the solute molecule in solution, a solute-solvent interaction energy, and solvent-solvent interaction energies. It also needs to be pointed out that the sum of the first two terms in eq (8) is the first root (lowest eigenvalue) of the CI matrix. Solute-solvent interaction energy in the CI representation is given as follows:

(9)

where P^CI,g is the ground state one-electron density matrix. A direct application of the hybrid QM-CI/MM potential in fluid simulations leads to an estimate of absorption spectroscopy of chromophores in solution, and a prediction of the solvatochromic shifts.

(10)

In addition, the difference in free energy of solvation between the ground and excited states of molecules can also be estimated via free energy perturbation method.

In practice, it is not possible to do full CI calculations, even for very small molecules, in fluid simulations. One is forced to use a computationally viable scheme by truncating the full CI matrix, or equivalently the CI expansion for the wavefunction. One possibility is to consider configurations that differ from the ground state determinant by a certain number of spin orbitals, say 4. Then, the CI would contain single, double, triple, and quadruple excitations from the selected trial orbitals. Alternatively, one may truncate the full CI to include only certain types of excitations. The simplest for estimating correlation energy is to include all single and double excitations from |Ψ₀>. For spectroscopy, the CIS method is also often used.

Key references: Luzhkov & Warshel, J. Am. Chem. Soc. 113, 4491 (1991); Gao, ibid. 116, 9324 (1994); Thompson & Schenter, J. Phys. Chem. 99, 6374 (1995).

At this point, we have implicitly assumed that only point charges in the MM region contribute to the modified the Fock matrix in the hybrid QM/MM method. This allows the wavefunction and charge distribution of the solute molecule to be polarized (for details of polarization energy calculation, see Science, 258, 631 (1992)). However, the solvent charges are kept fixed, without responding to the solute charge reorganization as well as the solvent dipole reorientations. Thus, it is important to treat both the solute and solvent polarizations on the same footing. Now, the MM solvent molecules are represented by a set of van der Waals parameters {σ_s, ε_s), a set of charges {q_s}, and a set of atomic polarizabilities {α_s}. Therefore, an additional term is added to the system effective Hamiltonian:

(10)

where the third term represents the interaction between the QM region and the induced dipoles of the solvent MM region. It is given below.

(11)

where μ_s are MM induced dipoles determined by

(12)

The solvent polarization energy is given as follows:

(13)

In these equations, F_s^qm(Ψ) is the electric field generated by the QM molecule at solvent position s. Since equation (10) depends on {μ_s} via eq (11), yielding the wavefunction Ψ, and {μ_s} depend on Ψ via eq (12). The solution of the HF equations and the convergence of MM induced dipoles are couple, and must be solved iteratively. The procedure may be summarized in the following scheme.

It turns out that a significant amount of computer time is spend in the computation of the MM induced dipoles. Using semiempirical QM models, it spends more computing time in determining the MM polarization energy than does the QM energy.

For excited state calculations, the fast dielectric response to the solute charge distribution up on excitation can be approximated by recomputing the solvent induced dipoles in the excited solute field, F^CI,excit(Ψ^CI,excit):

(14)

(15)

Note that the change of sign for the second term in eq (15) is due to the fact the CI calculations were performed on the basis of the ground state MM induced dipoles. The correction to the interaction between the excited state wavefunction and re-relaxed MM induced dipoles is estimated classically and twice the energy cost for creating these dipoles.

It appears that the next correction to the hybrid QM/MM potential is to include exchange repulsion interactions between the QM and MM region in QM calculations, and to handle charge transfer between the two regions. In this regard, a sort of "quantum mechanical" description of the "MM" region seems to be necessary.

Key references: Thompson & Shenter, J. Phys. Chem. 99, 6374 (1995); Gao, J. Comput. Chem. submitted for publication.

The EVB method has been effectively used by Warshel and his coworkers in hybrid QM/MM calculations. With the S_N2 reaction of Cl^- + CH₃Cl in water as an example, and consider the three configurations:

The potential surface may be simply represent by two electronic states, describing the reactant and product wave functions:

(16)

where R describes the distances of the forming and breaking bonds.

The matrix elements are given formally by H₁₁⁰ = < ψ₁|H|ψ₁>, H₂₂⁰ = <ψ₂|H|ψ₂>, and H₁₂⁰ = <ψ₁| H |ψ₂>. Although these matrix elements may be evaluated by quantum mechanical methods, the assumptions made in this approach made it to have large errors. So, instead, these matrix elements are represented empirically (force field) by fitting the gas-phase information. Thus, the gas phase energies ε_i^o are represented by

(17)

A similar equation is written for ε₂⁰, the product state gas phase energy. And, the off-diagonal element is given by

(18)

To include the solvent effect, the diagonal matrix elements are modified and assume H₁₂ = H₁₂^0.

(19)

By diagonalizing the Hamiltonian matrix, the analytical ground state potential energy surface is obtained:

(20)

Key references: Warshel, "Computer Modeling of Chemical Reactions in Enzymes and in Solutions", Wiley, 1991.; Hwang, et al. J. Am. Chem. Soc. 110, 5297 (1988).

1. Some relevant references may be consulted with regard to implementing DFT method into a hybrid potential. Wesolowski & Warshel, J. Phys. Chem. 97, 8050 (1993). Stanton et al. J. Phys. Chem. 97, 11868 (1993). C. Alhambra, K. Byun, J. Gao, "The Geometry of Water in Liquid Water from Hybrid Ab Initio-Monte Carlo and Density Functional-Molecular Dynamics Simulations." In Combined Quantum Mechanical and Molecular Mechanical Methods, J. Gao and M. A. Thompson, Eds., ACS Symposium Series 712, American Chemical Society: Washington, D.C., 1998, pp. 35-49.