Chapter 1 Introduction 1 General Introduction

Using the principle that interaction energies between separated densities can be expressed entirely through their electrostatic multipoles, and making the assumption this is valid for our solvent-solute model, we use the model charge density to express the electrostatic coupling between the QM and MM regions. The total electrostatic interaction energy between the fitted Gaussian charge distribution of the QM region and a point charges of the MM system, q_j, is expressed in equation 7-6 and represented pictorially in Figure 7.6.

Up to this point, the electrostatic coupling appears to be the same as the simple mechanical coupling where no polarization of the QM wave function occurs. How then is the polarization achieved in this scheme? In the Car-Parrinello molecular dynamics framework the polarization of the wave function can be expressed through the force on the expansion coefficients in the coupled equations of motion.

We first calculated the potential of the Gaussian charge distribution due to the MM point charges which is given in Equation 7-8 for each Gaussian.

The potential acting on the wave function is then obtained using the chain rule and the partial derivatives of the fitted charges with respect to the expansion coefficients.

Thus, the potential acting on the Gaussian charge distribution (Eqn. 7-8) is transferred to the force on the wavefunction via Equation 7-9. It is this bridging relationship that mediates the polarization of the wave function by the MM charge distribution through the model charge distribution.

For the forces on the QM atoms, there is an additional force resulting from the sensitivity of the fitted charges on the atomic positions. The so called Pulay forces²²⁵ are given by equation 7-11.

The Pulay term expressed in equation 7-11 ensures that the total energy is consistent with the gradients, a necessary ingredient for performing energy conserving dynamics. The exact expressions for the partial derivatives of the fitted charges with respect to the expansion coefficients used in equation 7-9 and with respect to the nuclear position used in equation 7-11 are detailed elsewhere since they were originally derived for the charge isolation scheme in PAW.¹⁴²

Although the electrostatic coupling approach used in the PAW QM/MM implementation is exact, the multipolar expansion of the true density is strictly only valid for long-range interactions. Thus, the model charge density may not properly represent the electrostatic potential acting on the MM charges for close range interactions between the solvent and QM solute molecules. Since the total QM/MM interaction potential is empirical in nature, it is possible that optimization of the adjustable parameters within the potential will be able to compensate for these defects (if any) in the PAW QM/MM electrostatic coupling scheme. In the following sections this and other aspects of the electrostatic coupling will be tested.