Chapter 1 Introduction 1 General Introduction

atom typea	description	Ro (Å)	Do (kcal/mol)
CT	sp3 carbon	3.816	0.1094
CA	aromatic sp2 carbon	3.816	0.0860
HC	hydrogen bound to CT	2.918	0.0150
N3	sp3 nitrogen	3.750	0.1700
P	phosphate phosphorus	4.200	0.2000
O	sp2 oxygen	3.322	0.2100
OS	sp3 oxygen in ethers	3.367	0.1700
F	fluoride	3.500	0.0610
Li	lithium cation	2.274	0.0183

Molecular mechanics simulations using the simple van der Waals and fixed point charge intermolecular potential have a rich history for simulating liquids.³¹ With the adjustment of the van der Waals parameters and the partial atomic charges, these simple potentials have been successfully applied to simulate liquids as complex as water.^32,217,218 A good example of this is the TIP3P³² force field developed by Jorgenson and co-workers to simulate liquid water. The TIP3P model employs partial charges of +0.417e on the hydrogen atoms and a charge of -0.834e on the oxygen atom of the water molecules. For the van der Waals interaction, only a single Lennard-Jones potential is employed between the oxygen atoms of the water molecules. Despite the simplicity, Monte Carlo simulations of liquid water at 25°C and 1 atm have yielded excellent results for both thermodynamic and structural properties of water. The density, heat of vaporization and the heat capacity at constant pressure where found to be within 2%, 0.5% and 6% of the experimental values, respectively.³² For other liquids, mostly organic solvents, simple intermolecular potentials of similar quality have also been developed by Jorgensen and coworkers.^159,219

The pioneering work by Jorgensen and others acts as a good starting point for incorporating solvent effects via the QM/MM method. Since the solvent-solvent interactions have been carefully parameterized, the application of the QM/MM method can focus on the development of realistic intermolecular interactions between the QM and MM molecules. This is achieved by the optimization of the adjustable parameters in our hybrid potential, namely the van der Waals parameters involving the QM and MM atoms. Gao and coworkers were the first to do this successfully, where van der Waals parameters were optimized for their combined semi-empirical AM1/TIP3P potential. The potential was parameterized to reproduce a comprehensive set of gas-phase solute-water interaction energies obtained from experiment and high level ab initio calculations. Using this parameterization, Gao and co-workers were able to estimate free energies of hydration in good agreement with experiment.²⁰⁸ Parameterization specific to each kind of QM/MM implementation is necessary, even if the same MM water model is used. This is because wave functions of different levels of theory react differently to the same perturbation, in this case the MM point charges. Furthermore, the approximations of the QM/MM model and differences in the implementation can be remedied by the empirical refinement of the potential.

In this and the previous section, we have reviewed how solvent simulations are generally performed with computational methods, particularly, in the QM/MM approach. We now turn to a specific discussion of our efforts towards incorporating solvent effects into our PAW QM/MM molecular dynamics simulations. First, we will discuss the nature of the QM-MM electrostatic interactions used in our implementation. This will be followed by our parameterizations of the PAW QM/MM hybrid potential to match the 'true' potential.

QM/MM electrostatic coupling refers to the Coulombic interaction of the MM charge distribution with that in the QM system. The most simple coupling method, aside from its absence, does not involve the QM wave function directly. Instead, the interaction of the QM system with the MM charges is calculated using atomic point charges (or equivalent) assigned to the QM atoms.^13,164 The point charges are generally extracted from the QM wave function using some sort of charge density partitioning scheme^220,221 such as a Mulliken population analysis. Alternatively, they are adjusted to reproduce the electrostatic potential of the wave function in what is called an electrostatic potential fit (ESP fit).^222-224 Although the QM charge density can be adequately represented by the point charge model, this coupling method (sometimes called mechanical coupling)¹³ does not allow for the distortion or polarization of the wave function by the MM solvent molecules. In general, this is a necessity for properly simulating the effect of the solvent. Charge density polarization effects are introduced 'quantum mechanically' where the wave function is optimized in the electrostatic field due to the MM charge density. As outlined in the previous section, this is generally achieved by adding an one electron integral term into the Fock matrix elements.

The charges are chosen such that the electrostatic multipole moments of the true density within each simulation cell are reproduced. In other words a multipolar expansion of the true density is performed such that the long-range electrostatic potential is reproduced. This seemly complex fit can be accomplished at a relatively small expense in the plane wave formalism.¹⁴² Following a Fourier transform of the true density, the requirements can be expressed in a linear system of equations which can solved for the Q's. Since the Fourier transform of the density is ordinarily performed in the Car-Parrinello method, the additional cost of the multipolar expansion of the true density is negligible and essentially involves the inversion of an N_gxN_g matrix where N_g is the total number of Gaussians. Blöchl's scheme provides an accurate multipolar expansion of the true continuous charge density. The fitted charges are found to reproduce point charges^◊ obtained from experimental multipole moments to within 0.01e and dipole moments to within 0.2 D for a set of small organic and inorganic molecules.