Ab Initio Path Integral Molecular Dynamics and Monte Carlo

Chapter 15

Downloaded by STANFORD UNIV GREEN LIBR on May 16, 2012 | http://pubs.acs.org Publication Date (Web): February 6, 2012 | doi: 10.1021/bk-2012-1094.ch015

Ab Initio Path Integral Molecular Dynamics and Monte Carlo Simulations for Water Trimer and Oligopeptide Takatoshi Fujita,*,1 Masa-Aki Kusa,2 Takayuki Fujiwara,2 Yuji Mochizuki,2,3 and Shigenori Tanaka*,1 1Graduate

School of System Informatics, Kobe University, 1-1 Rokkodai, Nada, Kobe 657-8501, Japan 2Department of Chemistry and Research Center for Smart Molecules, Faculty of Science, Rikkyo University, 3-34-1 Nishi-ikebukuro, Toshima-ku, Tokyo 171-8501, Japan 3Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505 Japan *E-mails: [email protected] (T. Fujita); [email protected] (S. Tanaka)

We have performed ab initio path integral molecular dynamics and Monte Carlo simulations for water trimer and oligopeptide. In the first part, we illustrate the path integral molecular dynamics method based on fragment molecular orbital (FMO-PIMD) method. The FMO-PIMD method is applied to water trimer and glycine pentamer to investigate nuclear quantum effects on the structure and molecular interactions. In the second part, we employ the Møller-Plesset perturbation theory and discuss interplay of nuclear quantum effects and electron correlations.

© 2012 American Chemical Society In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.


Introduction Hydrogen-bonded systems are important in various fields of chemistry and biology. Theoretical treatments of hydrogen-bonded systems require accurate calculations of potential energy surface and consideration of quantum effects of nuclei. Earlier studies have shown that nuclear quantum effects (NQEs) of hydrogen atoms influence the structural and dynamical properties of these systems (1–9). Ab initio path integral simulation is one of the useful simulation methods suitable to treat hydrogen-bonded systems. Marx and Parrinello (10) and Tuckerman et al. (11) originally formulated ab initio path integral molecular dynamics (PIMD) method based on the density functional theory and the Car-Parrinello molecular dynamics (MD) method. This approach has been used in the studies of various hydrogen-bonded systems (12–16). Shiga, Tachikawa, and co-workers (17–19) later developed an ab initio PIMD method employing molecular orbital theory. Semiempirical molecular orbital theories have also been combined with the PIMD methods as an alternative approach to treat large systems (20). While most path integral simulations have employed the MD method to sample the nuclear configuration space, a path integral Monte Carlo (PIMC) (21) have also been employed by several authors (22–24). Pierleoni have formulated the PIMC method with quantum Monte Carlo method and applied it to the study of metallic hydrogen (24). In the first part of this article, we report the PIMD method to treat large hydrogen-bonded systems with high accuracy (25). Our approach is based on the fragment molecular orbital (FMO) method. The FMO method (26–33) is a highly parallelizable quantum-chemical method suitable for the accurate electronic structure calculations of large systems. The FMO-based molecular dynamics (FMO-MD) method have also been developed (34) and applied to hydrated molecular systems (35). The incorporation of three-body term (30, 31) in the FMO method improves the accuracy of the FMO-MD (36, 37). Nagata et al. (38) have developed fully analytical gradients in the FMO method, in which the self-consistent Z-vector method is derived for solving the response term due to the external potentials. Mochizuki et al. (39) have implemented the energy gradient of the second-order Møller–Plesset perturbation theory in conjunction with the FMO method. Fujita et al. (40) have recently introduced the periodic boundary condition on the FMO method and calculated the radial distribution functions of liquid water. With these developments of the FMO-MD methods in mind, FMO-based PIMD method is expected to be a useful simulation method to treat large hydrogen-bonded systems with inclusion of NQE. In the present study, we develop and implement the ab initio PIMD using the FMO method (FMO-PIMD). This FMO-PIMD method will provide an ab initio tool to simulate large hydrogen-bonded systems with high accuracy. Water trimer and glycine pentamer are studied using the FMO-PIMD method, where the electronic structure is first considered at the Hartree-Fock (HF) level. Analyzing the NQEs on the structure and the molecular interactions in these systems, we discuss the usefulness of the FMO-PIMD method. 188 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.


In the second part of this article, the PIMC calculation together with Møller-Plesset (MP) perturbation theory is carried out and compared with the previous HF-based PIMD results so as to examine the dependence of the NQEs on the approximation level of molecular orbital theories. This ab initio comparison is analogous to that by Habershon et al. (7). They have analyzed the NQEs on the dynamics of liquid water using several force fields for water and shown that the PIMD simulations that employ rigid water models and flexible water models where intramolecular interactions are described by simple harmonic functions overestimate the NQEs on the diffusivities in liquid water. They have also remarked that the NQEs on the hydrogen bonds (HBs) consist of two competing effects: Intermolecular zero-point energies and tunneling effects destabilize the HBs, which is canceled by intramolecular zero-point energies that lead to the increase in the O-H bond lengths and result in the stabilization of the HBs. This picture has recently been supported by Li et al. (9). They have analyzed the NQEs on the HBs in several hydrogen-bonded systems and suggested that the NQEs weaken the weak HBs but strengthen the relatively strong HBs. The electronic structure was considered at the HF level in the first part, while the HF method overestimates the O-H harmonic frequencies (41) and underestimates the hydrogen-bond interactions (42). Consequently, the PIMD methods with the HF may overestimate the NQE that destabilizes the HBs. The tendency similar to that by Habershon (7) may also be observed in the case of molecular orbital theories. To analyze this issue, the PIMC method with the MP perturbation theory was applied to water trimer and compared with previous PIMD-HF results. Since we are interested in interplay between electron correlations and the NQEs, the FMO method is not employed in this study. We then discuss how the NQEs depend on the level of molecular orbital theories.

Methods Path Integral Molecular Dynamics Based on Fragment Molecular Orbital Method In the FMO method, the calculated system is first divided into a lot of fragments, and molecular orbital calculations for fragment monomers and dimers are performed to approximate the total energy and other properties (26–29). The Born-Oppenheimer (BO) energy of the system can be calculated as follows:

where Nf is the number of fragments; EI and EIJ denote the energies of fragment monomer and dimer, respectively, which are calculated under the electrostatic potentials from other fragments. The atomic forces acting on the nuclei are obtained by taking the derivatives of the BO energy (28, 38). We have then developed the PIMD methods based on the FMO method. Detailed 189 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.


implementation of our method is found elsewhere (25). ABINIT-MPX program package was employed for FMO calculations (27, 29, 32, 33). FMO-MD and FMO-PIMD methods are applied to water trimer and glycine pentamer in the present study. In both simulations, electronic structure calculations were carried out at the HF level, and temperature was set at 300 K. In the calculations of water trimer, 6-31G** basis set was employed, and the imaginary time slice of 32 was used in the PIMD calculation. For glycine pentamer, STO-3G basis set was employed, and the imaginary time slices of 16 was applied. Other details of the calculations are found in Ref . (25). Path Integral Monte Carlo with Møller-Plesset Perturbation Theory In the second part of this article, we show the results of the PIMC simulations with MP perturbation theory for water trimer. Here, we employ the PIMC methods because the gradient calculation of the MP perturbation theory is expensive. The second-order MP perturbation (MP2) theory was used to take account of electron correlations at the least computational cost. In addition to the MP2, we employ third-order MP perturbation (MP3) theory. The O-O distances of the water trimer optimized at the MP2 level are shorter than those optimized at the post MP2 level (41, 43–46). We discuss the effects of higher-order electron correlations on the structure in connection with the NQEs. We have employed the normal-mode sampling of PIMC (22, 23) which relies on the following effective potential:

where P is imaginary time slices, N is the number of atoms, and s and I refer to the indices of P and N, respectively. wP=P1/2/kBT, where kB is the Boltzmann constant, and T is temperature. q are normal mode coordinates and related to atomic Cartesian coordinates R through a linear transformation, where zero-frequency mode (s=1) correspond to centroid coordinates which describe the center of mass motions of the polymer. The normal mode masses MI(s) is related to the atomic masses MI through MI(s)=MIλ(s), where λ(s) is given by

E is the Born-Oppenheimer potential energy calculated by quantum mechanical methods. Following procedure was employed as a trial displacement. First, centroid coordinates (s=1) or other coordinates (s≠1) are chosen with each probability of 0.5. In the movements of centroid coordinates, we attempt to move the centroid coordinates of one atom. In the movements of other coordinates, we attempt to simultaneously move several normal mode coordinates other than the centroid coordinates. We performed this trial move in such a way that the acceptance ratio becomes 0.4-0.5 in the system of water trimer. 190 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.

Water trimer was employed also in this study. MC and PIMC simulations using MP2 and MP3 were performed at the temperature of 300K, where 6-31G** basis sets was employed. Simulations were carried out for 300000 steps; equilibration runs were carried out for 100000 steps, and following 200000 steps were used for calculations of physical properties. The imaginary time slices of 16 was employed for the PIMC simulations. We compare the MP2 and the MP3 results with that of the former HF and discuss the dependence of the NQEs on the level of molecular orbital theories.


Results and Discussion FMO-PIMD for Water Trimer We first show the optimized structure of the water trimer, before illustrating the results of the MD and the PIMD simulations. Figure 1 represents the global minimum (GM) structure of the water trimer optimized by the FMO-HF/6-31G** method. It is well established (47) that the GM structure of the water trimer is a cyclic structure where each water monomer simultaneously acts as both hydrogen-bond donor and acceptor. The O-O distances, O-O-O angles, and the interfragment interaction energies (IFIEs) of the FMO method are also shown in Fig. 1. Since each water molecule was assigned as one fragment, the IFIEs describe the intermolecular interactions between the water molecules.

Figure 1. Optimized structure of water trimer calculated by FMO-HF/6-31G** method. The O-O distances (ROO), the O-O-O angles (θOOO), and the interfragment interaction energies (ΔE) are shown in units of Å, degrees, and kcal/mol, respectively. In order to examine the NQEs on the structure of water trimer, probability distributions of the O-O distances and the O-O-O angles were calculated and are shown in Figs. 2(a) and 2(b), respectively. In the quantum system, the extended separation of the O-O distance was found up to the range longer than 4.5Å, which is not observed in the classical system. The differences between the quantum and the classical systems for nuclear degrees of freedom are also found in the distribution of the O-O-O angles. The maximum peaks are located at around 60 degrees in both the MD and the PIMD 191 In Advances in Quantum Monte Carlo; Tanaka, S., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 2012.


simulations, hence the cyclic structures are dominant. The distribution of the OO-O angles in the MD simulations extends within the range of 60±20 degrees, indicating that the cyclic structure was maintained throughout the simulation. On the other hand, in the PIMD simulation, the distribution is also observed in the range of 20< θOOO

Ab Initio Path Integral Molecular Dynamics and Monte Carlo

Recommend Documents