Chapter 1 Introduction 1 General Introduction

The other broad approach to introducing solvent effects into the quantum mechanical potential surface is to treat the solvent molecules explicitly. With standard electronic structure calculations, this is achieved by surrounding the solute molecule with solvent molecules - all of which are treated quantum mechanically. For ab initio level calculations, generally only a few solvent molecules can be included. Although the interactions between the solvent molecules and solute molecule are treated rigorously - a few solvent molecules does not simulate the bulk solvent. Therefore, even qualitative conclusions can be dubious in nature. Recently, the Car-Parrinello¹³⁰ molecular dynamics method has allowed for simulations of bulk liquid to be carried out at the density functional level.^211-213 Although the approach is promising, it is still impractical for anything but the smallest solute and solvent molecules (e.g. water in water and methanol in water). The combined QM/MM approach seems well suited for performing explicit solvation simulations. A natural partition exists such that the solute is treated with quantum mechanics, while the explicit solvent molecules are handled by a much more efficient MM force field. Pioneering the approach, Gao and co-workers^{10,14,27,208,214,215} have demonstrated that the QM/MM method can be used to predict solvent polarization effects at the quantitative level for such properties as reaction barriers, equilibrium constants and solvation free energies.

In this chapter our efforts towards introducing explicit solvation effects into the Car-Parrinello framework with our PAW QM/MM implementation are presented. First we will briefly describe how the bulk liquid is to be simulated with periodic boundary conditions. Next we will provide a description of the QM/MM intermolecular potential that will be used to describe the interaction between the QM solute and MM solvent. Following this general discussion, details of the PAW QM/MM electrostatic coupling will then be given. Then the PAW QM/MM model will be evaluated for its ability to represent the true intermolecular interaction potential and to perform energy conserving dynamics. This work has the ultimate goal of allowing us to introduce solvent effects into our dynamical simulations of olefin polymerization catalysts. Although this has not yet been achieved, this chapter provides a foundation for future work.

7.2 Periodic Boundary Conditions

One of the most fundamental problems with applying the explicit solvation model involves the proper simulation of the bulk solvent. Consider the simulation of 1000 water molecules. At a density of 1.0 g/mL, these 1000 molecules would only fill a 30x30x30 Å container. Furthermore, the surface tension would force the system to assume a spherical shape. In this microscopic drop of water, nearly half of the solvent molecules in our simulation would be at the surface of the drop. This will create many undesired effects, since it is known that water molecules at the liquid-vacuum surface behave quite differently than in the bulk. Thus, even a simulation of a 1000 water molecules, which is quite substantial in computational terms, is not a good model of the bulk liquid.

The most common way to realistically simulate bulk liquid and to minimize the "surface" effects is to apply the technique of periodic boundary conditions.^◊ Conceptually, the system is placed in a simulation box that is surrounded by an infinite number of images of itself as depicted in Figure 7.2. The particles in the central cell interact with the images such that there is no liquid-vacuum surface. Furthermore, when a molecule leaves the central cell, it enters the opposite face with the same velocity as sketched in Figure 7.2. In this way, the bulk liquid is approximated. Since the central cell is exactly replicated in all directions, a periodic crystal is actually simulated. Thus care must be taken to insure that the periodic system is a good approximation to the macroscopic bulk liquid. Generally this is achieved by using a large simulation cell so as to minimize the artificial effects imposed by the symmetry and periodicity of the cell. The size of the cell necessary for a particular simulation is primarily governed by the range of the intermolecular potentials involved and the range of the phenomenon under investigation.³¹ For determining equilibrium behaviour of liquids away from any phase transitions, periodic boundary conditions have been shown to be an excellent approximation. A detailed discussion of the periodic boundary approximation and its limitations can be found in the book of Allen and Tidlesly.³¹

The circle surrounding atom 1 with a "cut-off" radius R_c is used to show which pair potentials are included and which are discarded with the minimum image convention. For example, particle 2 of the central cell lies within the radius, R_c, and therefore this interaction is included in determining the force on atom 1. All other interactions with the periodic images of particle 2 are then discarded. Particle 4 of central cell lies outside the cut-off radius and therefore this interaction is not included. However, its periodic image in the cell labeled 'B' does lie within the cut-off radius and, therefore, it is only this interaction that contributes to the total force on atom 1.

The minimum image convention can be enforced by imposing a cut-off function to all non-bonded potentials that all interactions involving distances greater than the cut-off radius, R_c, are set to zero. If the cut-off radius is too long, the minimum image convention will be violated. For the example shown in Figure 7.3, if the cut-off radius must be less than L/2, where L is the dimension of the simulation box. In other words, with R_c > L/2, two images of the same particle may lie within the cut-off radius. The truncation of the infinite sums in the application of the periodic boundary conditions offers substantial computational savings. The minimum image convention is a natural truncation scheme and it is one that allows for easy vectorization³¹ of the computer code, a valuable computational asset. For these reasons, the application of periodic boundary conditions with the minimum image convention is the most widely used method for simulating liquids at the atomic level.

For the simulation of some systems, such as an ionic liquid, the neglect of the long range electrostatic interactions which can effect the calculation of structural and thermodynamic properties.³¹ In this case, special methods to approximate the infinite summations have to be implemented. One such technique used to evaluate the interaction of a charge with all of its periodic images is the Ewald summation method.¹⁶⁹ At this point, however, these have not been implemented into our QM/MM code.

In this section we begin to discuss QM/MM solute-solvent intermolecular potential. For the time being let us confine ourselves to the situation where the QM/MM boundary does not cross any covalent bonds. This would be the case if we were to treat our solute molecules entirely at the quantum mechanical level and the solvent molecules at the molecular mechanics level. With this condition, the QM and MM molecules interact only through non-bonded potentials. These non-bonded interactions can be conveniently separated into two parts, an electrostatic component and a so called van der Waals component. Here, we provide a general introduction to how these non-bonded interactions are represented within a combined QM/MM potential with an emphasis on how the hybrid potential can be calibrated to properly represent 'true' intermolecular potential.

In order to model the polarization of the solute by the solvent, the charges assigned to the MM atoms should act to distort the QM wave function. This type of electrostatic coupling involves optimizing the wave function in the electrostatic field created by the MM charges. In a conventional Hartree-Fock or DFT QM/MM calculation this is most often achieved by the addition a MM charge/QM-electron term to the one-electron matrix of the Fock equations.^{8,13,162,163,216} Equation 7-1 expresses the Kohn-Sham Fock matrix element, where the last term in the bracket introduces the effect of the MM partial charges.

In the presence of the MM charge distribution, the QM wave function will be polarized as to enhance the favourable electrostatic interactions. It is important to note that there is an energy penalty for the distortion of the wave function from its gas phase state which is sometimes termed the distortion energy, ∆E_dist. From classical linear response theory¹⁴ the distortion energy is one-half of the gain in the interaction energy from polarizing the wavefunction.^◊