Importance of Three-Body Interactions in Molecular Dynamics

Article pubs.acs.org/JCTC

Importance of Three-Body Interactions in Molecular Dynamics Simulations of Water Demonstrated with the Fragment Molecular Orbital Method Spencer R. Pruitt,† Hiroya Nakata,‡ Takeshi Nagata,§ Maricris Mayes,∥ Yuri Alexeev,† Graham Fletcher,† Dmitri G. Fedorov,*,§ Kazuo Kitaura,⊥ and Mark S. Gordon*,# †

Argonne Leadership Computing Facility, Argonne National Laboratory, 9700 S. Cass Avenue, Lemont, Illinois 60439, United States Department of Fundamental Technology Research, R&D Center Kagoshima, Kyocera Corporation, 1-4 Kokubu Yamashita-cho, Kirishima-shi, Kagoshima 899-4312, Japan § Nanosystem Research Institute, National Institute of Advanced Industrial Science and Technology, 1-1-1 Umenzono, Tsukuba, Ibaraki 305-8568, Japan ∥ Department of Chemistry and Biochemistry, University of Massachusetts Dartmouth, 285 Old Westport Road, Dartmouth, Massachusetts 02747-2300, United States ⊥ Graduate School of System Informatics, Kobe University, 1-1 Rokkodai-cho, Nada-ku, Kobe 657-8501, Japan # Department of Chemistry and Ames Laboratory, Iowa State University, 201 Spedding Hall, Ames, Iowa 50011, United States ‡

S Supporting Information *

ABSTRACT: The analytic first derivative with respect to nuclear coordinates is formulated and implemented in the framework of the three-body fragment molecular orbital (FMO) method. The gradient has been derived and implemented for restricted second-order Møller−Plesset perturbation theory, as well as for both restricted and unrestricted Hartree−Fock and density functional theory. The importance of the three-body fully analytic gradient is illustrated through the failure of the two-body FMO method during molecular dynamics simulations of a small water cluster. The parallel implementation of the fragment molecular orbital method, its parallel efficiency, and its scalability on the Blue Gene/ Q architecture up to 262 144 CPU cores are also discussed.

1. INTRODUCTION The ability to simulate water with accurate methods is important for the computational study of biological and chemical systems, for which an accurate description of the hydrogen-bonding network in liquid water1 is crucial. When several layers of solvent are added, even for small inorganic solutes, the size of the system becomes too large to treat with conventional ab initio electronic structure methods. Biochemical systems, such as proteins, pose an even bigger problem in terms of computational cost. In general, ab initio simulations of solvated systems are too time-consuming to perform practical simulations with traditional approaches. A number of fragment-based methods have been proposed to improve the computational efficiency for large systems,2−20 and most of them involve calculations of agglomerates of fragments in one way or another to include many-body interactions. A recent application of one such method, the fragment molecular orbital (FMO) method, to water clusters21 has highlighted the importance of three-body terms in the relative energetics of water clusters. The requirement for methods to provide accurate three-body effects in water energetics emphasizes the © 2016 American Chemical Society

need for three-body analytic gradients as a prerequisite for practical applications, such as dynamics simulations. In the FMO method,22−25 one divides a large molecular system into fragments and performs ab initio calculations for individual fragments (monomers) and fragment pairs (dimers) in the presence of an electrostatic potential (ESP) due to the remaining fragments in order to reduce the overall computational scaling. The total two-body FMO (FMO2) energy is obtained by summing the monomer energies and the pair interaction energies (obtained from dimers). In an effort to increase the accuracy of the method, this summation was extended to include three-body (trimer) contributions (FMO3).26,27 It has been demonstrated that three-body effects play an integral role in systems dominated by quantum manybody effects, in particular in water clusters.21,28−36 The importance of three-body and higher order effects in cluster sizes as small as six water molecules has been shown by Xantheas,28−34 as well as Bates and Tschumper,37 to account Received: December 21, 2015 Published: February 25, 2016 1423

DOI: 10.1021/acs.jctc.5b01208 J. Chem. Theory Comput. 2016, 12, 1423−1435

Article

Journal of Chemical Theory and Computation for up to ∼30% of the total cluster interaction energy. These studies highlight the importance of accurately describing manybody interactions in aqueous systems during dynamical simulations. The FMO method has recently been used to perform molecular dynamics (MD) simulations.38−50 However, in many of these applications, an approximate gradient that lacked response terms was used. These response terms arise because of the use of a fixed monomer ESP during dimer calculations, neglecting the mutual self-consistency between monomers and dimers. Without fully analytic gradients, energy conservation in MD was shown to be poor.50 A complete implementation of the fully analytic FMO2 gradient that includes these response terms was recently developed at the Hartree−Fock (HF), density functional theory (DFT),51,52 and second-order Møller−Plesset perturbation theory (MP2) levels.53−56 An approximate FMO3 gradient26,49 has also been used for MD simulations.57,58 There is clearly a need for a fully analytic FMO3 gradient, and the present study formulates and implements the FMO3 gradient at the restricted HF59 (RHF), MP2,60,61 and DFT as well as restricted open-shell HF (ROHF), unrestricted HF (UHF), and DFT (UDFT) levels. Additionally, the importance of an efficient implementation cannot be overstated. Performing simulations efficiently on massively parallel computers is essential in order to simulate systems with more than a few atoms on a realistic time scale. The scalability of the FMO method on the Blue Gene/P (BG/P) supercomputer (Intrepid) hosted by the Argonne Leadership Computing Facility (ALCF) has been demonstrated.62 The present study builds upon the previous work and describes additional algorithmic strategies necessary for FMO-MD simulations to scale efficiently on the current BG/Q supercomputer (Mira) at the ALCF.63 The manuscript is organized as follows: Section 2 provides an overview of the FMO methodology, and the mathematical formalism and implementation details of the fully analytic FMO3 gradient. Section 3 gives an outline of the computational details of the investigation of three-body effects during dynamics simulations of water clusters, as well as the calculations used to test and validate the fully analytic FMO3 gradient. Section 4 presents the results of FMO2 MD simulations on a small water cluster, the accuracy of the fully analytic FMO3 gradient, and parallel efficiency and scalability data. Section 5 concludes.

∑

(μμKv + vμKv)

(4)

K ≠X

Throughout this study, the italic lower case Roman (ijkl...) and Greek (μνρσ...) indices denote the MOs and AOs, respectively. The one-electron and two-electron contributions in VIJμν are, respectively, K uμν =

∑ A∈K

K νμν =

∑

μ

−ZA ν |r − RA|

(5)

K Dλσ (μν|λσ )

(6)

λσ ∈ K

DKλσ

where is the density matrix element of monomer fragment K and (μν|λσ) is a two-electron integral in the atomic orbital (AO) basis. It is important to realize that FMO2 is not a purely two-body method. The polarization of all N fragments is done selfconsistently, taking into account N-body effects. A detailed explanation of how the many-body polarization effects are treated in FMO2 is given elsewhere.64 The dimer and trimer density matrix differences, ΔDIJ and ΔDIJK in eq 1, are defined by ΔDIJ = DIJ − (DI ⊕ D J )

(7)

ΔDIJK = DIJK − (DI ⊕ D J ⊕ DK )

(8)

The direct sums in eqs 7 and 8 mean blockwise additive insertions of monomer matrices into dimer and trimer matrices, respectively. The internal fragment RHF energies in eq 1 may be written in the form EX′ =

∑

X X Dμν hμν +

μν ∈ X

∑

+

1 2

∑ μνλσ ∈ X

⎡ X X 1 X X⎤ D D − Dμλ Dνσ ⎥(μν|λσ ) ⎣⎢ μν λσ ⎦ 2

X Dμν Pμν X + EXNR

μν ∈ X

(9)

hXμν

where is the X-mer one-electron Hamiltonian and the nuclear repulsion energy is

′ − ΔEIJ′ − ΔEIK ′ − ΔE′JK ) (ΔEIJK

I>J>K

∑

(3)

N

EXNR =

N

+

′ = EIJK ′ − EI′ − E′J − EK′ ΔEIJK

VμXv =

N

∑

(2)

where E′X is the internal fragment energy for fragment X (X = I, IJ, or IJK for monomers, dimers, or trimers, respectively). The internal energy for X differs from the RHF, MP2, or DFT energy of X by the value of the ESP contribution Tr(DXVX), subtracted in an a posteriori fashion after the self-consistent field (SCF) converges. VX is the matrix form of the ESP for the given X-mer due to the electron densities and nuclei of the remaining fragments, i.e.,

2. METHODOLOGY 2.1. Fragment Molecular Orbital Method. In the following equations, the FMO energy expression22−25 and the corresponding gradient equations are briefly summarized. The general expression for the FMO3 total energy27 applicable to RHF, MP2, and DFT is given by EFMO3 = EFMO2 +

ΔEIJ′ = EIJ′ − EI′ − E′J

IJK

IJK

IJ

IJ

[Tr(ΔD V ) − Tr(ΔD V )

∑ ∑ B ∈ X A(∈ X ) > B

ZAZB RAB

(10)

I>J>K

For MP2 there is an extra term that adds the correlation energy. For DFT there is an exchange−correlation term (not shown) and two-electron terms in eq 9 that may be modified depending on the functional. For fragmentation across a covalent bond, the hybrid orbital projection (HOP) contribution

− Tr(ΔDIK V IK ) − Tr(ΔD JK V JK )] N

EFMO2 =

N

N

∑ EI′+ ∑ ΔEIJ′ + ∑ Tr(ΔDIJ V IJ ) I

I>J

I>J

(1) 1424


Article

Journal of Chemical Theory and Computation occ

∑ 2⟨i|P ̂

X

∑

|i ⟩ =

i∈X

∂CμXi

X X Dμν Pμν

∂a

(11)

μν ∈ X

X

∑ Bk |θk⟩⟨θk|

(12)

k∈X

where |θk⟩ is a hybrid orbital and the universal constant Bk is usually set to 106 au. 2.2. Fully Analytic Three-Body FMO-RHF Gradient. To obtain the FMO3 gradient, one takes the analytic derivative of the energy in eq 1. The FMO3 gradient is thus divided into the FMO2 gradient and additional three-body corrections. The derivations for FMO2 can be found elsewhere.53 The derivation of the three-body corrections includes (a) terms involving internal energies E′X and (b) terms involving the ESP energies. The differentiation of EX′ with respect to nuclear coordinate a is ∂EX′ = ∂a

∑

X Dμν

1 + ∂a 2

+

X Dμν

∑

X ∂Pμν

∂a

μν ∈ X

∂ Tr(ΔDIJK V IJK ) ∂a

∑ μνλσ ∈ X occ

∑

−2

i ,j∈X

occ

Sjia , X F′jiX

−4∑

∑

i∈X r∈X

Uria , XVriX

∂E NR + X ∂a

∑ [2(ij|kk) − (ik|jk)] + PijX k∈X

U̅ a , M , IJK − 2

λσ ∈ L

∂(μν|λσ ) ⎤⎥ ∂a ⎥⎦

M ∂Sμν

∂a

+ U̅ a , IJK , IJK

∑ ∑

L(IJK ) a , L Sμν ΔXμν

occ

∑ ∑ ∑ ΔDμνIJK Uria ,L(μν|ri)

∑

+4

L Dλσ

L ≠ I , J , K μν ∈ L vir

(19)

L ≠ I , J , K μν ∈ IJK r ∈ L i ∈ L

where X , IJK Wμν =

(14)

1 4

X IJK X Dμλ Vλσ Dσν

∑

(20)

λσ ∈ X

and L(IJK ) ΔXμν =

1 4

⎡ ⎤ L⎢ IJK L ΔDζη Dμλ (ζη|λσ )⎥Dσν ∑ ⎢⎣ ⎥⎦ λσ ∈ L ζη ∈ IJK

∑

(21)

Similar to FMO2,53 in FMO3 it is also possible to simplify the response terms that arise in eq 13 and those coming from Tr(ΔDIJVIJ) and Tr(ΔDIJKVIJK), by summing the U̅ a , X , X contributions according to eq 1. The sum of all such contributions is

occ

N

(15)

N

N

U̅ a = − ∑ U̅ a , I , I − ∑ ΔU̅ a , IJ − I

where the MO-based projection operator matrix PijX is introduced for convenience

N

+

I>J

ΔU̅ a , IJK

∑ I>J>K

N

∑ ΔŨ a ,IJ + ∑ I>J

∑

M , IJK Wμν

M ∈ {I , J , K }

vir

The quantities C are the MO expansion coefficients. The selfconsistent charge process ensures the self-consistency of monomer densities with respect to each other but not with respect to dimer or trimer densities. When dimers and trimers are computed, the monomer densities used in the ESP are frozen. Therefore, the response terms (Ua,X ri ) appear in the gradient. The internal fragment Fock matrix element Fij′X is given in terms of MOs:

PijX =

∑

∑

−

∑

∂a

∑

+2

⎡ ∂u L ⎢ μν + ⎢⎣ ∂a

IJK ∂Sμν IJK , IJK Wμν

μν ∈ IJK

⎡ X X 1 X X ⎤ ∂(μν|λσ ) ⎢⎣DμνDλσ − DμλDνσ ⎥⎦ ∂a 2

∂a

Fij′ X = hijX +

L≠I ,J ,K

∑

−2

X ∂Sμν CμXi* CvjX

μν ∈ X

∑

μν ∈ IJK

where the overlap derivative matrix element is defined by

∑

IJK ΔDμν

∑

=

(13)

Sija , X =

(18)

m∈X

M ∈ {I , J , K } μν ∈ M

X ∂hμν

μν ∈ X

a,X X Umi Cμm

∑

To obtain the occupied-virtual orbital response Ua,X ri , one must solve the coupled-perturbed HF (CPHF) equations. This procedure is discussed in subsection 2.3. There are two kinds of differentiation of ESP-related terms in eq 1: (a) dimer terms of the form Tr(ΔDIJVIJ), whose derivative is given elsewhere,53 and (b) a trimer-specific term,

must be considered in eq 9. The sum on the left-hand side of eq 11 is over occupied (occ) molecular orbitals. The HOP operator is defined by P̂ =

occ + vir

=

(ΔŨ

a , IJK

− ΔŨ

a , IJ

− ΔŨ

a , IK

− ΔŨ

a , JK

X X CμXi*Pμν Cvj

(22) (16)

μν ∈ X

where

In order to facilitate treatment of the response terms, the following equation defined in a previous study56 is introduced: occ

U̅ a , X , Y = 4 ∑

ΔU̅ a , IJ = U̅ a , IJ , IJ − U̅ a , I , I − U̅ a , J , J a , IJ ΔŨ = U̅ a , IJ , IJ − U̅ a , I , IJ − U̅ a , J , IJ

vir

∑ Uria ,XVriY

i∈X r∈X

)

I>J>K

ΔU̅ a , IJK = U̅ a , IJK , IJK − U̅ a , I , I − U̅ a , J , J − U̅ a , K , K − ΔU̅ a , IJ − ΔU̅ a , IK − ΔU̅ a , JK

(17)

Y The response terms Ua,X ri associated with the ESP Vri in the MO basis arise from the expansion of the MO coefficient derivatives.65 Equation 17 sums the Ua,X ri terms multiplied by the VYri terms in order to simplify the gradient derivations. The derivatives of the MO coefficients can be written as

a , IJK ΔŨ = U̅ a , IJK , IJK − U̅ a , I , IJK − U̅ a , J , IJK − U̅ a , K , IJK

(23)

In eq 22, the dimer- and trimer-related orbital response terms (U̅ a , IJ , IJ , U̅ a , IJK , IJK ) need not be evaluated because they cancel 1425


Article

Journal of Chemical Theory and Computation N

out, and the only terms that require the solution of CPHF equations come from the monomer responses Ua,I ri . Similar to FMO2,56 U̅ a is zero when no approximations are applied to the ESP or when all the ESP contributions are computed with the ESP point charge (ESP-PC) approximation.56 Otherwise, U̅ a is not equal to zero and must be computed in order to obtain a fully analytic gradient expression. U̅ a may be regarded as a compensation term that is created by defining the internal fragment energies and the ESP contribution separately in eq 1. Without the ESP approximations, the contribution of the occupied-virtual orbital responses to the FMO3 gradient is obtained by collecting the quadruple sum in eq 19 N

ℜa =

ℜa , X = 4

I>J

M≠X

×

occ

∑ ∑ ∑ ∑ ΔDμνX Uria ,M(μν|ri) (25)

Equation 24 can be further simplified as vir

ℜa =

occ

∑ ∑ LriM Uria ,M

(26)

r∈M i∈M

where N

LriM =

∑ ∑

IJ 4ΔDμν (μν|ri)

I > J ≠ M μν ∈ IJ N

+

∑

[

I>J>K ≠M

−

∑

∑

IJK 4ΔDμν (μν|ri)

μν ∈ IJK

IJ 4ΔDμν (μν|ri) −

μν ∈ IJ

−

∑ μν ∈ JK

∑

M≠X

r ∈ M i ∈ M A ∈ M μν ∈ X

X a,M ΔDμν Uri μ

1 ν |r − RA|

∑∑

M M M (CrMλ CiM σ + Ciλ Crσ )Sλσ

(28)

where RMX is the distance between the fragments M and X and RESP‑PC is the predefined threshold for switching from the use of Coulomb integrals to the Mulliken point charge approximation. With the ESP-PC approximation, one cannot use the simplification of eq 24 into eq 26; instead, one must use eq 24 with the responses given in eq 28. Although not explicitly shown in the derivations, the FMO3 gradient can be used with approximations in which sums over all dimers and trimers are replaced by the sums over closely separated dimers and trimers. Then, one must also consider the electrostatic (ES) interactions of far-separated dimers. Note that there are no contributions from ES trimers, because of the additivity of ES interactions. The treatment of ES dimers and their analytic derivatives is the same as in FMO2.68 In summary, without the point charge approximation to the ESP, the contributions of the orbital response terms to the gradient are only expressed as in eq 24. On the other hand, when applying the ESP-PC approximation, one must explicitly consider all of the contributions shown in eqs 22 and 28. Because the point charge approximation dramatically reduces the computational cost, this study focuses on the FMO3 analytic gradient with the point charge approximation, and the importance of the response terms is discussed subsequently. In this work, the analytic FMO3 gradient was developed for the restricted closed and open shell, as well as unrestricted open shell Hartree−Fock methods, RHF, ROHF, and UHF, respectively.55,69−72 2.4. Fully Analytic Three-Body FMO-MP2 and FMODFT Gradients. In this work, the fully analytic FMO3-MP2 gradient is developed, with the restriction that the ESP approximation ESP-PC may not be used, as was also the case for the analytic FMO2-MP2 gradient.73 To attain a fully analytic gradient, appropriate response terms are added to the approximate FMO3-MP2 gradient,66 following the FMO2-MP2 derivation and extending the FMO3-RHF gradient discussed above. The FMO3-DFT analytic gradient has also been implemented as part of the effort that is discussed in this work, extending the previously implemented FMO2-DFT analytic gradient.74 The DFT gradient formulation can also be obtained by combining the FMO2-DFT analytic gradient with the FMO3-RHF development. In this work, the analytic FMO3 gradient was developed for both restricted and unrestricted DFT methods, RDFT and UDFT, respectively.71 Details on the formulations of the analytic FMO3-DFT and FMO3-MP2 gradients can be found in the Supporting Information. 2.5. Algorithmic Advances. As mentioned in the Introduction, an efficient implementation of the FMO3 gradient is necessary in order to harness available petascale resources toward the simulation of systems that contain thousands of atoms and that require the explicit treatment of three-body interactions. The FMO3 method has been parallelized using the generalized distributed data interface

where the response contribution of fragment X, ℜa , X , is represented in terms of the responses Ua,M ri of ESP fragment M

M ≠ X μν ∈ X r ∈ M i ∈ M

occ

∑∑ ∑ ∑

λ∈A σ∈M

(ℜa , IJK − ℜa , IJ − ℜa , IK − ℜa , JK )

vir

vir

∑ RMX ≥ RESP ‐ PC

(24)

N

μν ∈ X r ∈ M i ∈ M

N

+2

I>J>K

ℜa , X = 4

occ

RMX < RESP ‐ PC

N

∑ ℜa ,IJ + ∑

vir

∑ ∑ ∑ ΔDμνX Uria ,M(μν|ri)

∑

IK 4ΔDμν (μν|ri)

μν ∈ IK JK 4ΔDμν (μν|ri)]

(27)

In order to obtain the fully analytic FMO3 gradient, one must add (a) the total response contribution in eq 24 to the approximate gradient66 that neglects all response contributions in the derivatives above and (b) the ESP penalty response term in eq 24, unless it is zero as stated above. 2.3. Point Charge Approximation to the Electrostatic Potential. For FMO2, the derivations of the analytic gradient for approximate ESPs can be found elsewhere.67 For FMO3, with the ESP point charge approximation (ESP-PC), the contributions of the occupied-virtual orbital responses to the FMO3 gradient shown in eq 25 for far separated fragments are replaced with those calculated using Mulliken point charges instead of the electron density: 1426


Article

Journal of Chemical Theory and Computation (GDDI)75 in GAMESS. When GDDI is used, all nodes are divided into groups, and a single group performs an individual fragment calculation. The number of groups can be very large; for example, one can assign a single trimer per GDDI group. Because of the use of approximations, only compact trimers are calculated, so that their number increases linearly with the number of fragments N if the system is sufficiently large.25 However, the prefactor can be on the order of dozens. For a calculation of 4096 water molecules, the actual number of trimers to be computed is 473 315. Then, if a trimer is computed by a GDDI group, one must open 1 419 945 files simultaneously (3 files per trimer calculation). The files opened are (a) dictionary files, containing various information, such as one-electron integrals, MO coeffients, Fock matrix, etc., (b) second order self-consistent field (SOSCF) converger data files containing Fock matrix derivatives, and (c) DFT grid data. The latter only appears in DFT. This is a very heavy load for modern supercomputers which often have no local hard disks attached. Another important advance is made in this work by implementing a “file-less” execution of GAMESS. GAMESS has several SCF convergers, all of which use different sets of scratch files. The research reported here implemented a fully file-less execution of GAMESS for the SOSCF algorithm.76 In this context “file-less” only applies to scratch files. There is still an input file and an output file for each GDDI group. However, the amount of information written to the output file was considerably reduced. To achieve a file-less FMO-MP2 calculation, the distributed-data-based MP2 gradient of Fletcher et al.77 can be used. The differences between the new and old SOSCF algorithm are as follows. The old algorithm uses a file to store Fock matrix derivatives and other related information, which are necessary to obtain MO coefficients in the SOSCF algorithm. In the new algorithm these data are stored in dynamically allocated memory, distributed among CPU cores. There is a small difference in terms of parallelization: in the old algorithm a loop running through the SOSCF data accumulated over SCF iterations runs redundantly on all cores, whereas in the new algorithm the loop is parallelized over available cores. However, the cost of this loop is small, and overall there is no significant difference in parallel performance except for the fact that files are not used.

calculated by double differencing, with coordinate shifts of 0.0005 Å. Both the energy and the gradient were calculated with slightly more accurate integrals than the default setting, using the Rys quadrature method for all energy and gradient integrals. The FMO method includes approximations that can be used to reduce the number of explicit QM dimer and trimer calculations. The point charge approximation to the ESP (ESPPC) relevant to the fully analytic gradient formalism has been discussed in detail in Section 2.3. In the separated dimer approximation (ES-DIM), the energy difference in eq 2 is computed as the electrostatic interaction between the charge distributions of two fragments. In the COR-SD approximation applicable to FMO-MP2 only, the correlation energy contribution to the interaction in eq 2 is neglected for far separated dimers. Trimer interactions for far separated trimers are neglected in the I-TRIM approximation. An application of these approximations to water clusters can be found elsewhere,21 with the formalism derived in a recent review.2 The values of unitless23 thresholds, denoted by R, were set to RES‑DIM = 4.0, RESP‑PC = 4.0, and RI‑TRIM = 2,2,2,2. For MP2, RCOR‑SD was set to the same value as RES‑DIM. The ability of the FMO3 MD method to perform an extended MD simulation on water is presented in Section 4.3. A simulation using the new fully analytic FMO3 gradient at the MP2/6-31G(d,p) level of theory was performed on cluster sizes of (H2O)8 and (H2O)16 to demonstrate energy conservation. Both clusters were equilibrated, using the NVT ensemble and a Nosé−Hoover thermostat78 to control the temperature, for 2 ps with a step size of 0.5 fs. After equilibration, production simulations on the (H2O)8 and (H2O)16 clusters using the NVE ensemble were performed using a step size of 0.5 fs for 40 and 10 ps, respectively. A spherical boundary potential with the form V=

1 2

∑ Sforce·(RA − Roff )2 A

(29)

was applied to prevent evaporation of water molecules. The value of RA in eq 29 is taken to be the distance between atom A and the center of mass, and Roff is the user-defined radius of the spherical potential. The user-defined value of Sforce was set to 3 kcal/(mol Å2) to contain all of the atoms within Rof f from the center of the defined sphere. The potential V was applied to an atom when RA > Roff. Values for FMO-specific approximations used in Section 4.3 for the MD simulation were RESP‑PC = 0.0, RES‑DIM = 3.25, RCOR‑SD = 3.25, and RI‑TRIM = 1.25, −1, 2.0, 2.0. Finally, the parallel scalability of the FMO3 method using the new fully analytic gradient is discussed in Section 4.4. Benchmark calculations on large water clusters ranging in size from 512 to 4096 water molecules were performed on the Blue Gene/Q computer Mira at the ALCF. Single-point FMO2 and FMO3 fully analytic gradient calculations at the RHF-D/augcc-pVDZ level of theory were performed on up to 262 144 cores to demonstrate scalability of the FMO method to petascale core counts. The practical application to MD simulations of water exploiting the scalability of the FMO3 fully analytic gradient at the MP2/6-31G(d,p) level of theory is also shown.

3. COMPUTATIONAL DETAILS Initial simulations were performed on a small test system of eight water molecules to assess the performance of previously implemented versions of the FMO method. Section 4.1 presents the analysis of an FMO2-MP2/6-31+G(d) MD simulation of 770 fs performed using a 0.25 fs step size and an NVT ensemble. A Nosé−Hoover thermostat with chains78 was used to maintain the system temperature at 300 K. The potential energy surface (PES) of the last 48 fs of this FMO2MP2 MD simulation was examined in greater detail using RHFD, with the 2010 version (D3) of the Grimme empirical dispersion correction,79 and MP2 with the 6-31+G(d) basis set. The accuracy of the FMO3 analytic gradient is compared in Section 4.2 to that of the numerical gradient using a cluster of 64 water molecules calculated at the RHF/6-31G(d), MP2/631G(d), and DFT/6-31G(d) levels of theory. For open-shell methods, tetramethylpiperidine N-oxide (TEMPO) was solvated in 8 ethanol molecules and divided into 6 fragments, following an earlier study.71 The numerical energy gradient is

4. RESULTS 4.1. Failure of FMO2-MD for Water. Since the implementation of the fully analytic gradient for FMO2,53 it 1427


Article

Journal of Chemical Theory and Computation

Figure 1. FMO2/MP2/6-31+G(d) MD/failure for a cluster of eight water molecules. A proton (highlighted in green) is transferred from one water molecule to another after 770 fs of simulation time. Dashed black lines indicate hydrogen bonds, and solid red/white lines indicate covalent bonds.

Figure 2. Potential energy surface of a cluster of eight water molecules, calculated with MP2/6-31+G(d)-based methods.

for biological simulations, a more rigorous solution is to more accurately describe the many-body interactions through improved theoretical methods. To illustrate the problem that can occur when one uses FMO2 to simulate a system (like water) that has significant many-body effects, a single FMO2/MP2 MD simulation was performed on an 8-water cluster, using a fragmentation scheme of one water molecule per fragment. Figure 1 shows the transfer of a proton (highlighted in green) from one water molecule to an adjacent molecule after 770 fs of simulation time. Note that, at 300 K, the concentration of H+ in pure water is on the order of 1.0 × 10−7 mol/L,80 so dissociation is a rare event. The PES of the last 48 fs of the simulation was investigated to determine the cause of the proton transfer that is not expected in such simulations. The geometry at each femtosecond was extracted, and single-point energy calculations were performed using FMO2 and FMO3 at both the RHF and MP2 levels of theory. The behavior of the PES for MP2- and RHF-based methods is shown in Figures 2 and 3, respectively. It is apparent from Figure 2 that FMO2 fails to properly reproduce the repulsive

has become feasible to perform FMO2 MD molecular dynamic simulations. However, for systems such as water where intermolecular interactions play an important role beyond the two-body level, FMO2 MD is prone to failure. This section illustrates this failure for water clusters as small as eight water molecules and shows that the failure manifests as a proton transfer between two water molecules. Although the accuracy of FMO2 is sufficient in most applications, such as geometry optimizations, the thermal motion appears to add an additional driving force, leading to molecules coming close to each other, at which point the neglected three-body interactions become important. Indeed, the problem of insufficient accuracy of FMO2 was pointed out previously by Pruitt et al.21 and by Komeiji et al.38 Pruitt et al. demonstrated that FMO3, but not FMO2, reproduces the relative MP2 energies of water clusters, due to the omission of explicit many-body effects in FMO2. Komeiji et al.43 dealt with the spurious proton transfer between water molecules through “dynamic refragmentation” in order to prevent failure of the simulations. While this approach allows simulations to be performed and is an important development 1428


Article


Figure 3. Potential energy surface of a cluster of eight water molecules, calculated with RHF/6-31+G(d)-based methods.

wall of the PES as the proton begins to transfer. Instead, FMO2 provides a favorable path to proton transfer, producing a progressive lowering of the total energy as the simulation proceeds. This is clearly incorrect, with full MP2 results showing increasingly disfavored geometries as the water molecules move closer together. In contrast to FMO2, FMO3 reproduces the MP2 PES nearly exactly. A similar situation can be seen in Figure 3 for RHF-based methods, with FMO2 producing a relatively flat PES as the simulation progresses, while full RHF and FMO3 show more disfavored geometries. The addition of an empirical dispersion correction to the RHF methods (RHF-D) does little to help FMO2 produce better results, only translating the PES down in energy for all methods. The effect of the flatter FMO2-RHF PES can be seen when the total length of the simulation is compared to FMO2-MP2 MD simulations. While FMO2-MP2/6-31+G(d) fails after only 770 fs, FMO2-RHF/6-31+G(d) simulations are capable of running for picoseconds. Electron correlation effects, neglected by RHF, increase the binding and accelerate the spurious proton transfer. Eventually, FMO2-RHF suffers from the same problem as FMO2-MP2; however, other factors can cause failure more quickly, including an increase in either basis set or water cluster size. The incorrect decrease on the FMO2 potential energy surface is the main source of the inaccurate PES. Since full RHF and MP2 calculations, as well as FMO3, are capable of accurately reproducing the PES, the origin of the problem is a “fragmentation collapse” due to truncation of quantum manybody effects at the two-body level in FMO2. The failure of FMO2 to capture the many-body effects when using a fragmentation scheme of one water molecule per fragment is a possible source of the fragmentation collapse. Although the use of different fragmentation schemes such as two or more water molecules per fragment might result in

improved behavior of the potential energy surface, such a scheme is not amenable to dynamics simulations, in which the motion of water molecules away from each other is commonplace. Similarly, it has been suggested that dynamic refragmentation might fix the simulation failures;43 however, the PES produced from such simulations is suspect due to potential discontinuities. Therefore, the development of fully analytic energy gradients for the FMO3 method is the most promising approach to obtain accurate FMO MD simulations of aqueous chemical systems. 4.2. FMO3 Fully Analytic Gradient. The accuracy of the analytic FMO3 gradient for RHF, MP2, and DFT is demonstrated by a direct comparison with the numerical gradient. For this purpose, the analytic and numerical gradient elements for a cluster of 64 water molecules were calculated at the three levels of theory, and TEMPO solvated in ethanol was calculated for open-shell methods. Figure 4 shows a summary of the accuracy that is attained when the newly implemented response terms are included, compared to the previous gradient implementation. The horizontal axis in Figure 4 is the number of Cartesian coordinate elements, and the vertical axis is the error in the analytic gradient relative to the numerical gradient. Neglecting the response terms in eq 24 and eq 28 (blue dashed line in Figure 4) typically produces a maximum error of 0.003567 au/bohr and an RMSD of 0.000650 au/bohr. Inclusion of the response terms in eq 24 and eq 28 (red solid line in Figure 4) dramatically improves the accuracy of the gradient, producing a maximum error compared to the numerical gradient of 0.000031 au/bohr, while the RMSD falls to 0.000003 au/bohr. The magnitude of these errors is sufficiently small to enable accurate MD simulations, as well as geometry optimizations. For the open-shell methods the total error in the gradient is small for the particular system chosen as a test. However, 1429


Article


Figure 4. Errors in the new fully analytic gradients for (a) RHF, (H2O)64, (b) MP2, (H2O)64, (c) DFT, (H2O)64, (d) ROHF, TEMPO+(C2H5OH)8, (e) UHF, TEMPO+(C2H5OH)8, and (f) UDFT, TEMPO+(C2H5OH)8 using the 6-31G(d) basis set compared to fully numeric gradient elements. The gradient errors without the inclusion of response terms (no SCZV) are plotted in dashed blue, while the gradient errors using the new fully analytic gradient including response terms (SCZV) are plotted in solid red. The abscissa represents the Cartesian coordinates of the atoms. Errors are in Hartree/Bohr.

neglecting response terms creates an error too large even for geometry optimizations, and even more so for MD. 4.3. FMO3-MP2MD Energy Conservation. To measure the ability of the new fully analytic gradients to facilitate MD simulations, an FMO3/MP2/6-31G(d,p) MD simulation was

run on a system of (H2O)8 (for comparison with Section 4.1) for 40 ps with a step size of 0.5 fs, using the NVE ensemble (Figure 5a). Figure 5a shows that the addition of explicit, fully ab initio three-body effects to the simulation of water removes the “fragmentation collapse” exhibited by the FMO2 MD 1430


Article


Figure 5. Energy conservation for a FMO3-MP2MD simulation on (a) (H2O)8 and (b) (H2O)16 using the 6-31G(d,p) basis set over 40 and 10 ps, respectively.

4.4. Parallel Scalability. The FMO method can take advantage of multilevel parallelism81 to reduce the time-tosolution of a simulation. To test the efficacy of the distributed memory-based file storage for FMO calculations, single-point FMO2 and FMO3 gradient calculations were performed with both the disk-based and memory-based implementations on a system of 1024 water molecules. The file-less execution mode introduced in this work on the BG/Q architecture yields a reduction in wall time of ∼13% for FMO2 compared to diskbased implementations and a wall time reduction of ∼20% for FMO3. Consider the best wall time of 273 s for an FMO2 energy+gradient calculation using the aug-cc-pVDZ basis set on a cluster of 512 water molecules with 32 768 cores. A 13% reduction in wall time per MD time-step (41 s per step) would accumulate to a savings of ∼342 h of wall time, or ∼11 million core hours, over a 30 ps simulation. Such large core hour requirements illustrate the challenges in performing MD

simulations performed on (H2O)8, discussed in Section 4.1. The accuracy of the fully analytic three-body gradient demonstrated in Section 4.2 is also shown to produce good energy conservation. As shown in Figure 5, the standard deviation of the total energy for the simulation is calculated to be 0.11 kcal/mol, with a median absolute deviation (MAD) of 0.09 kcal/mol. A second FMO3/MP2/6-31G(d,p) MD simulation was run on a system of (H2O)16 for 10 ps with a step size of 0.5 fs, using the NVE ensemble (Figure 5b). The accuracy of the fully analytic three-body gradient demonstrated in Section 4.2 shows good energy conservation for this larger system as well. The standard deviation of the total energy for the simulation is calculated to be 0.09 kcal/mol, with a MAD of 0.06 kcal/mol. Future enhancements to the FMO-MD method (such as constraint algorithms) are expected to improve energy conservation and increase the usable time-step size, thereby reducing the simulation wall time. 1431


Article


Figure 6. Water clusters (512, 1024, 2048, and 4096 water molecules from left to right) used for FMO3-RHF-D benchmarking on Mira.

Figure 7. FMO2-RHF-D3 (left) and FMO3-RHF-D3 (right) MD step time using the aug-cc-pVDZ basis set. Calculations performed on Mira in c16 mode.

of the FMO3 MD method to perform practical MD simulations at the MP2 level of theory, the parallel scalability of the FMO3MP2 fully analytic gradient was tested on a system of 64 water molecules at the FMO3-MP2/6-31G(d,p) level of theory using cutoff values shown to provide good energy conservation for MD simulations for the electrostatic dimer, point charge ESP approximations, and trimer approximations. A cluster size of 64 water molecules was chosen, since this is the generally accepted minimum system size required for MD simulations of water using periodic boundary conditions (no periodic boundary conditions were used for these benchmarks). The parallel efficiency for one MD time step on Mira, using 64 nodes (1 node per fragment) as the baseline, is shown in Table 1. The parallel efficiency shows the percentage of achieving the maximum theoretical performance, so that 90.3% for 128 nodes means that the computed speedup was 0.903 × (128/64) = 1.8 compared to the baseline. The parallel

simulations using ab initio methods, even when employing fragmentation methods such as FMO. To investigate the scalability of the FMO method on the BG/Q architecture, spherical water clusters of 512, 1024, 2048, and 4096 molecules (Figure 6) were generated using the program SOLVATE.82 For each cluster size, a single fully analytic gradient calculation was performed with both FMO2 and FMO3 at the RHF-D3 level of theory. Since only the parallel scalability is tested here, default cutoff values for the electrostatic dimer, point charge ESP approximations, and trimer approximations were used. The aug-cc-pVDZ basis set was chosen to illustrate computational requirements for a moderate size basis set. Scalability of both FMO2 and FMO3 is shown up to 262 144 cores in Figure 7. Wall times presented are for execution of the FMO program within the GAMESS program package. To assess the scaling for both FMO2 and FMO3 gradient calculations with respect to system size, the wall times for each cluster size are compared (Figure 7). For FMO2 calculations, doubling the size of the system increases the wall clock time by a factor between 2 and 2.5, achieving agreement with previously published results.62 Doubling of the system size for FMO3 calculations increases wall clock times by a factor of ∼2 in the best case and a factor of ∼3.4 in the worst case. Although the parallel scalability of the FMO-RHF-D3 method demonstrates a possible application of the FMO method on massively parallel computers with large system sizes, practical application of FMO MD simulations requires a wall time per MD time step of under 60 seconds. Additionally, the inclusion of explicit electron correlation is imperative for the proper description of aqueous systems. To illustrate the ability

Table 1. Time to Solution (in seconds) and Parallel Efficiency for an FMO3-MP2/6-31G(d,p) Fully Analytic Gradient Calculation on (H2O)64a number of nodes number of cores time to solution parallel efficiency a

1432

128

256

512

1024

2048

4096

2048

4096

8192

16384

32768

65536

880

507

261

136

77

55

90.3%

78.4%

76.1%

73.1%

64.5%

45.2%

Calculations performed on Mira in c16 mode. DOI: 10.1021/acs.jctc.5b01208 J. Chem. Theory Comput. 2016, 12, 1423−1435

Article

Journal of Chemical Theory and Computation efficiency could be improved with a fine-tuning of the number of GDDI groups. The group division per calculation is shown in Table 2, using a node division scheme of one node per GDDI

increase the number of GDDI groups and improve the efficiency. Figure 8 shows the wall time required and strong scaling (scaling on an increasing number of cores with a fixed system size) for one MD time step on Mira. As the number of processors increases from 2048 to 65 536, the wall time for a single FMO3-MP2MD time step drops from 880 to 58 s (shown in Figure 8). While this wall time is still too long to run simulations longer than 50−100 ps, it is an important step toward fully ab initio MP2 dynamics. Efforts to reduce this wall time are already underway, for example, through the use of constraint algorithms such as RATTLE83 in conjunction with the FMO method84 to increase the time-step size above 0.5 fs.

Table 2. GDDI Group Division on Mira for All Water Cluster Sizes Benchmarked number of molecules

number of nodes

GDDI monomer groups

GDDI dimer groups

GDDI trimer groups

512

128 256 512 1024 2048 128 256 512 1024 2048 4096 128 256 512 1024 2048 4096 8192 128 256 512 1024 2048 4096 8192 16384

128 256 512 512 512 128 256 512 1024 1024 1024 128 256 512 1024 2048 2048 2048 128 256 512 1024 2048 4096 4096 4096

128 256 512 1024 2048 128 256 512 1024 2048 4096 128 256 512 1024 2048 4096 8192 128 256 512 1024 2048 4096 8192 16384

128 256 512 1024 2048 128 256 512 1024 2048 4096 128 256 512 1024 2048 4096 8192 128 256 512 1024 2048 4096 8192 16384

1024

2048

4096

5. CONCLUSIONS The work presented here has demonstrated the need for threebody interactions when simulating aqueous systems using the FMO method, especially in molecular dynamics simulations at room temperature. The failure of a two-body approximation for a simulation of water clusters has been analyzed in detail. The results exemplify the need for higher-order many-body interactions during such simulations. To enable FMO MD simulations that include water (or other polar protic solvents), the fully analytic three-body gradient has been derived and implemented in the GAMESS program package. The energy conservation in dynamics simulations has been tested and shown to be capable of accurate simulations using a time step of 0.5 fs. The file-less execution mode developed in this work is an important step in using GAMESS on large multinode computer systems. The ability of the FMO method to efficiently use petascale processor counts has been demonstrated. The dispersion interaction is important for water properties,36 and it should be included in simulations, for example, by using RHF-D, DFT-D, or MP2. The massively parallel FMO3/MD simulation engine reported here can be used in future applications of FMO, in particular to those in aqueous solution. The failure of FMO2 for water clusters is an important finding and should be taken into consideration when performing MD simulations with any method that omits many-body interactions.

group when the number of fragment calculations allows for it. For example, a 4096-node calculation can use a maximum of 4096 GDDI groups to increase efficiency by reducing the group size. In this work, the definition of a “node” in GDDI was taken to be a physical node. However, it is possible to define logical nodes, by dividing one physical node into multiple logical nodes. On most systems, such an approach can be used to

Figure 8. Scalability (left) and strong scaling curves (right) for a FMO3-MP2MD step time on a cluster of 64 water molecules using the 6-31G(d,p) basis set. Blue: Actual scaling. Red: Ideal scaling. Calculations performed on Mira in c16 mode. 1433


Article


■

(19) He, X.; Sode, O.; Xantheas, S. S.; Hirata, S. J. Chem. Phys. 2012, 137, 204505. (20) Hirata, S.; Valiev, M.; Dupuis, M.; Xantheas, S. S.; Sugiki, S.; Sekino, H. Mol. Phys. 2005, 103, 2255. (21) Pruitt, S. R.; Addicoat, M. A.; Collins, M. A.; Gordon, M. S. Phys. Chem. Chem. Phys. 2012, 14, 7752. (22) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Chem. Phys. Lett. 1999, 313, 701. (23) Nakano, T.; Kaminuma, T.; Sato, T.; Fukuzawa, K.; Akiyama, Y.; Uebayasi, M.; Kitaura, K. Chem. Phys. Lett. 2002, 351, 475. (24) Fedorov, D. G.; Kitaura, K. J. Phys. Chem. A 2007, 111, 6904. (25) The Fragment Molecular Orbital Method: Practical Applications to Large Molecular Systems; CRC Press: Boca Raton, FL, 2009. (26) Fedorov, D. G.; Kitaura, K. J. Chem. Phys. 2004, 120, 6832. (27) Fedorov, D. G.; Kitaura, K. Chem. Phys. Lett. 2006, 433, 182. (28) Bulusu, S.; Yoo, S.; Apra, E.; Xantheas, S. S.; Zeng, X. C. J. Phys. Chem. A 2006, 110, 11781. (29) Burnham, C. J.; Xantheas, S. S. J. Chem. Phys. 2002, 116, 1479. (30) Fanourgakis, G. S.; Schenter, G. K.; Xantheas, S. S. J. Chem. Phys. 2006, 125, 141102. (31) Hodges, M. P.; Stone, A. J.; Xantheas, S. S. J. Phys. Chem. A 1997, 101, 9163. (32) Xantheas, S. S. J. Chem. Phys. 1994, 100, 7523. (33) Xantheas, S. S.; Burnham, C. J.; Harrison, R. J. J. Chem. Phys. 2002, 116, 1493. (34) Xantheas, S. S.; Dunning, T. H. J. Chem. Phys. 1993, 98, 8037. (35) Pruitt, S. R.; Steinmann, C.; Jensen, J. H.; Gordon, M. S. J. Chem. Theory Comput. 2013, 9, 2235. (36) Pruitt, S. R.; Leang, S. S.; Xu, P.; Fedorov, D. G.; Gordon, M. S. Comput. Theor. Chem. 2013, 1021, 70. (37) Bates, D.; Tschumper, G. S. J. Phys. Chem. A 2009, 113, 3555. (38) Komeiji, Y.; Nakano, T.; Fukuzawa, K.; Ueno, Y.; Inadomi, Y.; Nemoto, T.; Uebayasi, M.; Fedorov, D. G.; Kitaura, K. Chem. Phys. Lett. 2003, 372, 342. (39) Ishimoto, T.; Tokiwa, H.; Teramae, H.; Nagashima, U. Chem. Phys. Lett. 2004, 387, 460. (40) Ishimoto, T.; Tokiwa, H.; Teramae, H.; Nagashima, U. J. Chem. Phys. 2005, 122, 094905. (41) Mochizuki, Y.; Komeiji, Y.; Ishikawa, T.; Nakano, T.; Yamataka, H. Chem. Phys. Lett. 2007, 437, 66. (42) Sato, M.; Yamataka, H.; Komeiji, Y.; Mochizuki, Y.; Ishikawa, T.; Nakano, T. J. Am. Chem. Soc. 2008, 130, 2396. (43) Komeiji, Y.; Ishikawa, T.; Mochizuki, Y.; Yamataka, H.; Nakano, T. J. Comput. Chem. 2009, 30, 40. (44) Komeiji, Y.; Mochizuki, Y.; Nakano, T.; Fedorov, D. G. J. Mol. Struct.: THEOCHEM 2009, 898, 2. (45) Fujita, T.; Watanabe, H.; Tanaka, S. J. Phys. Soc. Jpn. 2009, 78, 104723. (46) Fujiwara, T.; Mochizuki, Y.; Komeiji, Y.; Okiyama, Y.; Mori, H.; Nakano, T.; Miyoshi, E. Chem. Phys. Lett. 2010, 490, 41. (47) Fujiwara, T.; Mori, H.; Mochizuki, Y.; Tatewaki, H.; Miyoshi, E. J. Mol. Struct.: THEOCHEM 2010, 949, 28. (48) Sato, M.; Yamataka, H.; Komeiji, Y.; Mochizuki, Y.; Nakano, T. Chem. - Eur. J. 2010, 16, 6430. (49) Komeiji, Y.; Mochizuki, Y.; Nakano, T. Chem. Phys. Lett. 2010, 484, 380. (50) Brorsen, K. R.; Minezawa, N.; Xu, F.; Windus, T. L.; Gordon, M. S. J. Chem. Theory Comput. 2012, 8, 5008. (51) Hohenberg, P.; Kohn, W. Phys. Rev. 1964, 136, B864. (52) Kohn, W.; Sham, L. J. Phys. Rev. 1965, 140, A1133. (53) Nagata, T.; Brorsen, K.; Fedorov, D. G.; Kitaura, K.; Gordon, M. S. J. Chem. Phys. 2011, 134, 124115. (54) Nagata, T.; Fedorov, D. G.; Ishimura, K.; Kitaura, K. J. Chem. Phys. 2011, 135, 044110. (55) Nakata, H.; Schmidt, M. W.; Fedorov, D. G.; Kitaura, K.; Nakamura, S.; Gordon, M. S. J. Phys. Chem. A 2014, 118, 9762. (56) Nagata, T.; Fedorov, D. G.; Kitaura, K. Chem. Phys. Lett. 2009, 475, 124.

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.5b01208. Derivations for the analytic gradient of FMO3-DFT and FMO3-MP2 (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*(D.G.F.) E-mail: [email protected]. *(M.S.G.) E-mail: [email protected]. Funding

This work was supported by the Office of Science, U.S. Department of Energy, under Contract DE-AC02-06CH11357, and by a U.S. National Science Foundation Software Infrastructure (SI2) grant, ACI-1450217. D.G.F. was supported by the Next Generation Super Computing Project, Nanoscience Program (MEXT, Japan), and Computational Materials Science Initiative (CMSI, Japan). Some of the computations reported here were performed on the Iowa State University Cyence cluster, obtained via a National Science Foundation MRI grant, at Iowa State University. An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS The authors also acknowledge Dr. Kurt Brorsen and Dr. Federico Zahariev for helpful discussions. REFERENCES

(1) Ball, P. Nature 2008, 452, 291. (2) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Chem. Rev. 2012, 112, 632. (3) Otto, P.; Ladik, J. Chem. Phys. 1975, 8, 192. (4) Gao, J. L. J. Phys. Chem. B 1997, 101, 657. (5) Tong, Y.; Mei, Y.; Zhang, J. Z. H.; Duan, L. L.; Zhang, Q. G. J. Theor. Comput. Chem. 2009, 8, 1265. (6) Gao, J.; Wang, Y. J. Chem. Phys. 2012, 136, 071101. (7) Kobayashi, M.; Nakai, H. J. Chem. Phys. 2013, 138, 044102. (8) Collins, M. A. Phys. Chem. Chem. Phys. 2012, 14, 7744. (9) Kobayashi, M.; Yoshikawa, T.; Nakai, H. Chem. Phys. Lett. 2010, 500, 172. (10) He, X.; Merz, K. M., Jr. J. Chem. Theory Comput. 2010, 6, 405. (11) Söderhjelm, P. r.; Kongsted, J.; Ryde, U. J. Chem. Theory Comput. 2010, 6, 1726. (12) Frank, A.; Möller, H. M.; Exner, T. E. J. Chem. Theory Comput. 2012, 8, 1480. (13) Kiewisch, K.; Jacob, C. R.; Visscher, L. J. Chem. Theory Comput. 2013, 9, 2425. (14) Gordon, M. S.; Smith, Q. A.; Xu, P.; Slipchenko, L. V. Annu. Rev. Phys. Chem. 2013, 64, 553. (15) Friedrich, J.; Yu, H.; Leverentz, H. R.; Bai, P.; Siepmann, J. I.; Truhlar, D. G. J. Phys. Chem. Lett. 2014, 5, 666. (16) Hirata, S.; Gilliard, K.; He, X.; Li, J.; Sode, O. Acc. Chem. Res. 2014, 47, 2721. (17) Gilliard, K.; Sode, O.; Hirata, S. J. Chem. Phys. 2014, 140, 174507. (18) Willow, S. Y.; Salim, M.; Kim, K. S.; Hirata, S. Sci. Rep. 2015, 5, 14358. 1434


Article

Journal of Chemical Theory and Computation (57) Mori, H.; Hirayama, N.; Komeiji, Y.; Mochizuki, Y. Comput. Theor. Chem. 2012, 986, 30. (58) Matsuda, A.; Mori, H. J. Solution Chem. 2014, 43, 1669. (59) Roothaan, C. Rev. Mod. Phys. 1951, 23, 69. (60) Pople, J. A.; Krishnan, R.; Schlegel, H. B.; Binkley, J. S. Int. J. Quantum Chem. 1979, 16, 225. (61) Møller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618. (62) Fletcher, G. D.; Fedorov, D. G.; Pruitt, S. R.; Windus, T. L.; Gordon, M. S. J. Chem. Theory Comput. 2012, 8, 75. (63) Argonne Leadership Computing Facility. http://www.alcf.anl. gov. (64) Fedorov, D. G.; Kitaura, K. J. Comput. Chem. 2007, 28, 222. (65) Yamaguchi, Y.; Schaefer, H. F., III; Osamura, Y.; Goddard, J. A New Dimension to Quantum Chemistry: Analytical Derivative Methods in Ab Initio Molecular Electronic Structure Theory; Oxford University Press: New York, 1994. (66) Mochizuki, Y.; Nakano, T.; Komeiji, Y.; Yamashita, K.; Okiyama, Y.; Yoshikawa, H.; Yamataka, H. Chem. Phys. Lett. 2011, 504, 95. (67) Nagata, T.; Fedorov, D. G.; Kitaura, K. Chem. Phys. Lett. 2012, 544, 87. (68) Nakata, H.; Nagata, T.; Fedorov, D. G.; Yokojima, S.; Kitaura, K.; Nakamura, S. J. Chem. Phys. 2013, 138, 164103. (69) Pruitt, S. R.; Fedorov, D. G.; Kitaura, K.; Gordon, M. S. J. Chem. Theory Comput. 2010, 6, 1. (70) Pruitt, S. R.; Fedorov, D. G.; Gordon, M. S. J. Phys. Chem. A 2012, 116, 4965. (71) Nakata, H.; Fedorov, D. G.; Nagata, T.; Yokojima, S.; Ogata, K.; Kitaura, K.; Nakamura, S. J. Chem. Phys. 2012, 137, 044110. (72) Nakata, H.; Fedorov, D. G.; Yokojima, S.; Kitaura, K.; Nakamura, S. Theor. Chem. Acc. 2014, 133, 1477. (73) Nagata, T.; Fedorov, D. G.; Li, H.; Kitaura, K. J. Chem. Phys. 2012, 136, 204112. (74) Brorsen, K. R.; Zahariev, F.; Nakata, H.; Fedorov, D. G.; Gordon, M. S. J. Chem. Theory Comput. 2014, 10, 5297. (75) Fedorov, D. G.; Olson, R. M.; Kitaura, K.; Gordon, M. S.; Koseki, S. J. Comput. Chem. 2004, 25, 872. (76) Chaban, G.; Schmidt, M. W.; Gordon, M. S. Theor. Chem. Acc. 1997, 97, 88. (77) Fletcher, G. D.; Schmidt, M. W.; Gordon, M. S. In Advances in Chemical Physics; Prigogine, I., Rice, S. A., Eds.; John Wiley & Sons: Hoboken, NJ, U.S.A., 1999; Vol. 110, p 267. (78) Harvey, S. C.; Tan, R.; Cheatham, T. E., III J. Comput. Chem. 1998, 19, 726. (79) Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. J. Chem. Phys. 2010, 132, 154104. (80) Stillinger, F. H. In Theoretical Chemistry, Advances and Perspectives, Eyring, H., Henderson, D., Eds.; Academic Press: New York, 1978; Vol. 3, pp 177−234. (81) Fedorov, D. G.; Olson, R. M.; Kitaura, K.; Gordon, M. S.; Koseki, S. J. Comput. Chem. 2004, 25, 872. (82) Grubmüller, H. Solvate, version 1.0.1; http://www.mpibpc.mpg. de/grubmueller/solvate. (83) Andersen, H. C. J. Comput. Phys. 1983, 52, 24. (84) Pruitt, S. R.; Brorsen, K. B.; Gordon, M. S. Phys. Chem. Chem. Phys. 2015, 17, 27027.

1435


Importance of Three-Body Interactions in Molecular Dynamics

Recommend Documents