Testing of the PAW QM/MM Coupling: Molecular Dynamics

In order for the QM/MM electrostatic coupling scheme to be used for molecular dynamics simulations, it must allow for energy conserving dynamics to be performed. In this section we examine this aspect of the PAW QM/MM implementation. As a preliminary test we have examined the energy conservation associated with the QM/MM dynamics of an ethane molecule in ethane at 200 K (liquid). In this simulation 17 ethane molecules have been simulated with periodic boundary conditions in a 25 Å cubic cell where one of the molecules was treated at the density functional level. Jorgensen's OPLS-AA¹⁵⁹ molecular mechanics force field parameters were used. The OPLS-AA charges assigned to the MM carbon and hydrogen atoms are -0.18 and +0.06 e, respectively. A molecular dynamics time step of 3.0 a.u. was used with the standard Verlet propagator. Masses of 12.0 amu and 1.5 amu were used for carbon and hydrogen atoms, respectively. The system was slowly heated to 200 K and equilibrated for 500 fs. No thermostats were applied to any of the three subsystems (wavefunction, MM nuclei, QM nuclei). The details of the Car-Parrinello simulation of the QM subsystem are essentially same as those presented in Section 5.3, the only notable difference being that of the simulation temperature and the polarizable electrostatic coupling.

The electrostatic coupling energy is plotted in Figure 7.7b. More specifically, this is the interaction energy between the MM point charges and the model charge density of the QM ethane molecule as defined in Equation 7-6. For most of the simulation the electrostatic coupling energy is small, fluctuating around 2 kcal/mol. However, at a simulation time of approximately 800 fs, a close contact with a MM ethane molecule increases the interaction energy to a value of roughly 10 kcal/mol. Notably, the enhanced electrostatic interaction does not effect either the drift of the total conserved energy or the amplitude of its oscillations.

One component of the combined Car-Parrinello Lagrangian that might be sensitive to the electrostatic coupling is the dynamics of the wave function. The fluctuating electrostatic potential due to the MM charges may cause the wave function to become 'hot' such that there is a degenerative deviation from the Born-Oppenheimer potential surface. Worse still, intermolecular collisions may cause more violent perturbations in the wave function that can potentially dislodge the wave function from the Born-Oppenheimer potential surface. Although this is likely to be exhibited in the energy conservation, the kinetic energy of the wave function is most sensitive to it. Plotted in Figure 7.7c is the kinetic energy of the wave function for both the polarizable electrostatic coupling run and the mechanically coupled QM/MM run. We notice that there is a small, but equal rate of drift in both cases. Application of a wave function thermostat as prescribed by Parrinello and Blöchl¹⁴⁹ will minimize the deviation of the wave function from the Born-Oppenheimer state and counter the drift in both cases. Importantly, at the 800 fs mark where there is an enhanced electrostatic coupling energy, there is no significant increase in the wave function kinetic energy and therefore no evidence of the wave function dislodging from the Born-Oppenheimer surface.

In this section we have only demonstrated the energy conservation of the coupling in the PAW QM/MM implantation. Future studies should include a more thorough examination of the potential consequences of the QM/MM electrostatic coupling on the Car-Parrinello coupled dynamics. In this case, such studies will have ramifications to other Car-Parrinello QM/MM implementations where the electrostatic coupling allows for polarization of the wave function.

The two component QM/MM interaction potential involving a van der Waals component and an electrostatic component (including polarization of the wave function) has been shown by others to represent the true quantum mechanical interaction potential well enough to provide quantitative results for solvation effects. The electrostatic coupling scheme used in the PAW QM/MM implementation involves a multipolar expansion of the true density that is strictly valid only for long range interactions. In this section we will examine if this coupling scheme can be used to correctly model the true quantum mechanical interaction potential at critical short range distances.

For this purpose we will explore the interaction between a water molecule and a Li⁺ cation. The interaction potential will be determined by the PAW QM/MM method where the water molecule is treated in the QM region and the Li⁺ ion is partitioned to the MM region as illustrated in Figure 7.8. The reference potential used to judge the quality of the PAW QM/MM method will be the 'pure QM' potential calculated with gradient corrected density functional theory. The details of the QM component of the PAW QM/MM calculations are as follows. The wave function of the water molecule was expanded in plane waves up to an energy cutoff of 30 Ry within a 18.4 Å cubic cell. The Becke¹⁰³-Perdew86^104,105 gradient corrected exchange-correlation functional was used in conjunction with the local density parameterization of Perdew and Zunger.¹¹⁸ For the pure QM calculation of the reference potential energy surface, the ADF program (version 2.3) was employed. Here the standard ADF triple- STO basis with polarization functions were employed. The Becke-Perdew86^103-105 exchange-correlation functional was also used in the ADF calculations. For both the ADF and PAW calculations, a [He] frozen core was utilized for oxygen. For all potential surfaces shown, the geometry of the water molecule was fixed (R_OH = 0.974 Å and HOH = 103.58°) and a C_2V symmetry was maintained such that the oxygen atom is oriented towards the Li⁺ ion as shown in Figure 7.8. The pure QM reference potential was calculated using ADF instead of PAW because the memory requirements for calculating the long range Li-water interactions with PAW were not available. (For long separations a large PAW simulation cell is required).

Li-O distance	a	PAW/AMBER interaction energyb (kcal/mol)					PAW fitted charges (e)
(Å)	(kcal/mol)				∆Edist			oxygen	hydrogenc
1.00	261.48	18914.68	-130.57	19010.4	34.80			-1.42	0.71
1.50	-10.17	68.28	-102.02	141.44	28.86			-1.42	0.71
1.60	-20.16	0.71	-86.59	63.98	23.32			-1.34	0.66
1.65	-23.32	-15.22	-79.88	43.69	20.96			-1.30	0.65
1.70	-25.58	-24.84	-73.79	30.11	18.84			-1.27	0.63
1.75	-27.13	-30.41	-68.34	20.91	17.01			-1.23	0.62
1.80	-28.12	-33.40	-63.38	14.62	15.35			-1.20	0.60
1.85	-28.66	-34.73	-58.96	10.28	13.94			-1.18	0.59
1.90	-28.87	-35.01	-54.91	7.26	12.62			-1.15	0.57
1.95	-28.81	-34.63	-51.29	5.15	11.49			-1.12	0.56
2.00	-28.52	-33.83	-47.96	3.66	10.45			-1.10	0.55
2.10	-27.53	-31.61	-42.19	1.85	8.72			-1.06	0.53
2.20	-26.20	-29.13	-37.41	0.92	7.36			-1.02	0.51
2.30	-24.70	-26.69	-33.38	0.44	6.24			-0.99	0.49
2.50	-21.61	-22.40	-27.03	0.06	4.56			-0.94	0.47
2.70	-18.70	-18.96	-22.36	-0.03	3.43			-0.89	0.45
3.00	-15.09	-15.03	-17.30	-0.05	2.32			-0.84	0.42
3.50	-10.85	-10.76	-12.04	-0.02	1.30			-0.79	0.39
4.00	-8.14	-8.09	-8.85	-0.014	0.77			-0.75	0.37
5.00	-5.10	-5.07	-5.39	-0.004	0.32			-0.70	0.35
6.00	-3.50	-3.48	-3.63	-0.001	0.15			-0.68	0.34
7.00	-2.55	-2.54	-2.62	-0.001	0.07			-0.66	0.33
8.00	-1.96	-1.95	-1.98	-0.000	0.02			-0.65	0.32
10.00	-1.25	-1.24	-1.25	-0.000	0.00			-0.64	0.32
12.00	-0.89	-0.85	-0.86	0.000	0.00			-0.63	0.31

The long-range part of the PAW/AMBER potential compares well with the reference potential. For distances of 3.0 Å and greater, the deviation of the hybrid potential from the reference is no more than 0.15 kcal/mol. For this range, the electrostatic contribution to the total interaction energy dominates, with the van der Waals energy accounting for less than 1% of the total. The agreement in this region shows that the model charge density used in the PAW QM/MM coupling represents the true density exceptionally well for long distance interactions.