Generalized Switch Functions in the Multilevel Many-Body Expansion

Article pubs.acs.org/JCTC

Generalized Switch Functions in the Multilevel Many-Body Expansion Method and Its Application to Water Clusters Guo Dong Chen, Jingwei Weng, Guoliang Song,* and Zhen Hua Li* Collaborative Innovation Center of Chemistry for Energy Material, Shanghai Key Laboratory of Molecular Catalysis & Innovative Materials, Department of Chemistry, Fudan University, Shanghai 200433, China S Supporting Information *

ABSTRACT: The many-body expansion (MBE) method is the basis of many fragment-based methods and is widely applied to the computation of large molecular systems. To reach linear-scaling computation, a cutoff must be used to discard those subsystems with long interfragment distances. However, this leads to a discontinuous potential energy surface (PES) that would cause various problems in geometry optimizations and molecular dynamics simulations. To solve this problem, we present a generalized-switch-function (GSF) approach to smooth the PES computed by the MBE method with the use of a cutoff distance. The GSFs are naturally normalized and are permutation invariant. This approach can be applied to adaptively computing any order of many-body correction energies with multilevel computational methods and is a dynamic subsystem approach. We have applied the two versions of our method, GSF-MBE(m)/L1 and GSF-MBE(m)/(L1:L2:L3), to water clusters. Thorough tests show that our method can indeed give smooth potential-energy surface and is linear scaling but without losing much accuracy for very large water clusters with appropriately chosen cutoff distances.

1. INTRODUCTION Quantum mechanics methods provide a systematic way to approach the best energy of a molecule system.1 However, these methods become impracticable when the molecule contains a large number of atoms because of the unacceptable scaling of the computational cost with system size. Although density functional theory (DFT) has made electronic structure computations relatively inexpensive, its O(N2∼3) scaling still prevents it from tackling very large systems such as water nanoparticles, proteins, or condensed materials. For this bottleneck to be overcome, a variety of fragmentbased methods have been developed in the last two decades. This technology allows one to apply ab initio computation on increasingly larger molecules and chemical systems.2−26 In 1994, Gadre and co-workers proposed a method called the molecular tailoring approach (MTA).14 In 1999, Kitaura and co-workers proposed the fragment molecular orbital (FMO) theory.15 Then, many fragment-based methods have been developed such as the molecular fractionation with conjugated caps (MFCC) scheme,16,17 generalized-energy-based fragmentation (GEBF),18−21 systematic molecular fragmentation (SMF),22,23 combined fragmentation method (CFM),24 kernel energy method (KEM), 25 and generalized many-body expansion (GMBE).26 Among these fragment-based methods, many-body expansion (MBE) is both conceptually and practically a simple approach15−27 and is widely used in the calculation of large molecular systems.7−12 In fact, most fragment-based methods are somehow built on MBE.28 © 2017 American Chemical Society

In the MBE method, if an entire system is partitioned into M fragments, the final energy of the entire system is expressed as the sum of energies of the M noninteracting fragments plus the many-body correction energies between pair fragments (2body), three fragments (3-body), and so on.27 These fragments are usually called subsystems. The expansion is formally exact if the expansion is full, i.e., the many-body correction energies are summed to M-body correction energy. However, the computational cost does not decrease when the expansion is full. The key to reduce computational cost is to truncate the high-order terms, which are usually very small and can be ignored or approximately computed by low-level methods.29−31 The MBE method truncated after m-body correction is abbreviated as MBE(m) and its computational scaling is O(Mm). For some systems with large many-body effects, for example, water systems using one water as a fragment to compute the MBE energy with chemical accuracy, four- or even five-body correction shall be considered, which leads to O(M4) or O(M5) scaling.31,32 Therefore, the computational costs for large systems would quickly become tremendous. Many efforts have been dedicated to improving the performance of the MBE method. Truhlar’s research group proposed electrostatically embedded many-body expansion (EE-MBE)3 to accelerate the convergence of the expansion. Podeszwa and co-workers proposed a stratified approximation many-body approach (SAMBA),31 applying low-level methods Received: February 10, 2017 Published: April 19, 2017 2010

DOI: 10.1021/acs.jctc.7b00144 J. Chem. Theory Comput. 2017, 13, 2010−2020

Article

Journal of Chemical Theory and Computation

the exploration of structural and binding properties of water clusters provide a key for understanding bulk water in its liquid and solid phase,42 water clusters have received extensive scientific and technological interest. Many benchmarks based on the accurate yet expensive CCSD(T)43 method have been applied to the computation of water clusters, but its seventh-power scaling with system size prevents it from application to large water clusters. Xantheas and co-workers have reported CCSD(T)/aug-cc-pVTZ energies for (H2O)16 and (H2O)17,44 each of which took roughly 3.33 h utilizing 120,000 AMD Opteron cores. The largest water cluster fully computed at the CCSD(T)/aug-cc-pVTZ level contains only 20 water molecules.45 Some benchmark calculations employed MBE-based methods and can be applied to large water clusters with larger basis sets. For example, relatively large water clusters have been examined with the nbody:many-body QM:QM technique.46,47 MTA and GEBF have been used for the benchmark computation of (H2O)20.48,49 There are also numerous research works on water systems using cheaper methods, for example, MP2 and DFT methods.50−53 With our GSF scheme, we can push the limit of high-accuracy calculations to larger water clusters and perform molecular dynamics simulations for large water systems using accurate first-principle methods.

to compute high-order terms to curb computational cost. Saha and Raghavachari developed a dimer-based MBE method, dimers of dimers (DOD),33 for the study of water clusters. Song and co-workers proposed an extended energy divide-andconquer (E-EDC) method34 to approach an accuracy of a higher-order expansion method by extrapolating energies computed at a lower-order expansion method. However, the rapidly growing number of high order terms in these methods is still a problem. For linear-scaling computation for the MBE(m) method to be reached, a distance cutoff can be applied when computing the n-body (2 ≤ n ≤ m) correction energies.35 Ř ezác and Salahub presented a multilevel fragment-based approach (MFBA)8 to decrease the number of long-range many-body correction energies to be computed. Truhlar and co-workers applied a cutoff on correlation energy in EE-MBE,4 trying to reduce the cost of the computation. However, discarding the terms by a hard distance cutoff leads to discontinuous potential-energy surface (PES). Geometry optimization and molecular dynamics simulation thus may encounter various problems. This is a common problem among fragment-based methods. To smooth the potential energy surface, one can design appropriate switch functions to smoothly switch the n-body correction energies from its real value to zero. For the 2-body correction energies or any terms that depends only on one variable, which is usually the distance between the fragments, the switch function is very simple and has been used by some MBE-based methods. For example, Bowman and co-workers used a switch function to adapt shortand long-range 2-body potential in their many-body force field for water.36,37 However, as those terms depend on more than one distance, such as 3-body terms, this is not an easy task. Some groups are working on it and have made some progresses. Paesani and co-workers applied a simple switch function on the 3-body terms in their MB-pol potential.38 The switch function of Paesani and co-workers is not a generalized one and cannot be easily extended to those terms beyond 3body terms. Pinski and Csáyi designed a delicate algorithm that can dynamically determine which atoms belong to a fragment according to the bonding connections between the atoms in their reactive many-body expansion (RMBE) method.39 However, without using a cutoff between fragments, RMBE “expansion to order n would scale with the n-th power of the system size”. In the present work, we try to design a generalized switchfunction (GSF) approach to smooth the PES of the MBE(m) method with a cutoff. This scheme should be a generalized one and can be applied to adaptively compute any n-body correction energies. To reach a balance between accuracy and efficiency, this scheme should allow the computation of n-body subsystems with large n-body correction energies by a highlevel method and those with small correction energies by a lowlevel method, i.e., adaptively computing n-body correction energies with multilevel methods. In addition, the switch functions should be normalized and permutation invariant. Herein, we apply our scheme to water, which may be one of the most important substances in the world. It is a major component of our bodies, and it covers more than two-thirds of the earth’s surface. However, because of the complex hydrogenbonding networks formed through highly directional noncovalent interactions, water is notoriously difficult to model.40 Because the properties of water clusters can be used to test the models constructed for molecular simulation of water41, and

2. METHODOLOGY 2.1. Many-Body Expansion. In the MBE method, an entire system is divided into M interacting fragments, named monomers. The energy of the entire system is given by the expressions M

E=

M

∑ Ei1 + ∑ dE(n) i1

M

dE(2) =

(1)

n=2

M

∑ ΔEi1i2 = ∑ (Ei1i2 − Ei1 − Ei2) i1< i2

i1< i2

(2)

M

dE(3) =

∑

ΔEi1i2i3

i1< i2 < i3 M

=

∑

(Ei1i2i3 − Ei1i2 − Ei1i3 − Ei2i3 + Ei1 + Ei2 + Ei3)

i1< i2 < i3

(3) M

dE(4) =

∑

ΔEi1i2i3i4

i1< i2 < i3 < i4 M

=

∑

(Ei1i2i3i4 − Ei1i2i3 − Ei1i2i4 − Ei2i3i4 + Ei1i2 + Ei1i3

i1< i2 < i3 < i4

+ Ei1i4 + Ei2i3 + Ei2i4 + Ei3i4 − Ei1 − Ei2 − Ei3 − Ei4) (4) ···

where dE(n) (n = 2 ∼ M) is the total n-body correction energies for n-body subsystems, ΔEi1i2⋯in is the n-body correction energy of an n-body subsystem, (i1i2⋯in), formed by the i1th, i2th, i3th, ···, and inth monomers, Ei1 is the energy of monomer i1, and Ei1i2⋯in is the energy of the subsystem (i1i2···in). Here, an n-body subsystem is a part of the system formed by n monomers selected out of M monomers. Usually, the expansion is 2011


Article

Journal of Chemical Theory and Computation truncated after 4-body in the common application.27 One can turn to ref 27 for the general formulas. 2.2. Generalized Switch Functions. An n-body subsystem can be viewed as a polyhedron with n vertices. The polyhedron has C2n pairs of vertices. For each pair of vertices i1 and i2 (1 ≤ i1 < i2 ≤ n), a weight can be assigned to it. Here, we define two weight functions: wi1i2 and its reverse function 1 − wi1i2. The weight wi1i2 is a function of the distance between i1 and i2, and it can be any smooth function that at least its first derivative is continuous. Here, we choose the popular fifth-order polynomial smooth function54

wi1i2

⎧ 1, xi1i2 < 0 ⎪ ⎪ = ⎨−6xi51i2 + 15xi41i2 − 10xi31i2 + 1, 0 ≤ xi1i2 < ⎪ ⎪ 0, xi1i2 ≥ 1 ⎩

Cn2 + 1

∑

sg(n) = 1 (13)

g=1

2.3. Application of Generalized Switch Functions to MBE. For a subsystem, we can draw two different types of lines between each pair of monomers according to the distance between the two monomers: a solid line represents a short distance whereas a dashed line represents a long distance. A connectivity, which can be used to represent an n-body subsystem, is then defined by a set of these C2n lines. An n-body 2 subsystem then has a total of 2Cn possible connectivities. For the 2- and 3-body subsystems, the two and eight possible connectivities, respectively, are graphically shown in Figure 1.

⎫ ⎪ ⎪ 1⎬ ⎪ ⎪ ⎭ (5)

R i1i2 − R C

where xi1i2 = ΔR and Ri1i2 is the distance between vertices i1 and i2. RC and ΔR are predefined parameters. ΔR is the width of the buffer zone that determines how fast wi1i2 decays from one to zero. Both RC and ΔR should be carefully determined for the systems of interest to find a compromise between accuracy and efficiency. To simplify the derivation, we replace the two-dimensional index (i1, i2) by a one-dimensional index p: w1 = w12, w2 = w13, ..., wC2n = wn−1,n. The index p can be mapped to (i1, i2) by the equation 1 1 p = − i12 − i1 + (i1 − 1)Cn2 + i2 2 2

Figure 1. Connectivities and connectivity groups for the 2- and 3-body subsystems where the solid line stands for close distance between two monomers.

(6)

These connectivities can be classified into different groups according to the number of their dashed lines in an ascending order: a 2-body subsystem has only two groups of connectivities, a 3-body subsystem has four groups, and an nbody subsystem has C2n + 1 groups. In principle, for each group we can use a different computational method, for example, groups with large many-body effects computed with a higherlevel method so that subsystems with short intermonomer distances can be treated with more accurate and expensive methods and subsystems with long intermonomer distances can be treated with less accurate and cheaper methods or even not to be computed at all. For a subsystem (i1i2···in), its n-body correction energy computed at computational level Lg is (Lg) (n) g) ΔE(L i1i2···in. For each ΔEi1i2···in, we give it a weight sg . Then, we obtain a general equation for the adaptive computation of dE(n)

Then, we define a reversing operator R̂ p that transfers the weight wp to its reverse weight 1 − wp R̂ pwp = 1 − wp

(7)

For each pair of vertices, we can assign either wp or R̂ pwp to it. Then, we can define switch functions that are the products of C2n weight functions. These switch functions can be classified by the number of reverse weights in it. In this paper, the switch functions with the same number of reverse weights are defined as the same group g and g = nrw + 1, where nrw (0 ≤ nrw ≤ Cn2) is the number of reverse weights. The switch function of the g group, s(n) g , is expressed as s1(n) = w1w2 ··· wCn2 s2(n) =

(8)

Cn2

M

∑ R̂ p(w1w2 ··· wCn2)

(9)

p=1

s3(n) =

Cn2 − 1

Cn2

∑ ∑

R̂ pR̂ q(w1w2 ··· wCn2)

p=1 q=p+1

s4(n) =

Cn2 − 2 Cn2 − 1

Cn2

∑ ∑ ∑

dE

R̂ pR̂ qR̂ r(w1w2 ··· wCn2)

p=1 q=p+1 r=q+1

(n)

=

Cn2 + 1

∑ ∑

(L )

(sg(n)ΔEi1i2g··· in) (14)

i1< i 2 ···< in g = 1

C2n

Applying eq 14, we seemingly need to compute + 1 times g) to obtain all ΔE(L i1i2···ins and this does not reduce computation costs at all. However, we actually do not need to compute C2n + 1 times for a given subsystem. For example, for a 2-body subsystem (i1i2) with an intermonomer distance shorter than (L2) RC, s(2) 2 = 0, and we do not need to compute ΔEi1i2 . If the distance is longer than RC + ΔR, s(2) 1 = 0 and we do not need to 1) compute ΔE(L i1i2 . Only when the distance is between RC and RC (L2) 1) + ΔR do we need to compute both ΔE(L i1i2 and ΔEi1i2 . Therefore, for a subsystem (i1i2···in), we only need to compute (n) g) those ΔE(L i1i2···ins with nonzero sg .

(10)

(11)

sC(n2)+ 1 = R1̂ R̂ 2 ··· R̂ Cn2(w1w2 ··· wCn2) = (1 − w1)(1 − w2)···(1 − wCn2) n

(12)

s(n) g s,

These switch functions, are called GSFs. They are permutation invariant and naturally normalized 2012


Article


all connected through solid lines in all four connectivities. Furthermore, one can choose to compute the contribution from the one to seven connectivities by the same level of computational method and discard the contribution of the eighth connectivity. It should also be noted that, usually, it would be expected that the connectivity with fewer dashed lines has stronger interaction energy so that it has larger contribution to dE(n) because the more dashed lines the connectivity has means the subsystem has more scattered monomers. However, some specific geometries with fewer dashed lines show strong many-body effects also, for example, the subsystem with a linear arrangement of monomers in water clusters.55 Therefore, it is important to compute those connectivities with large manybody effects with an accurate computational method. This can be realized by carefully selecting a proper grouping scheme. For example, for 3-body subsystems, grouping one to four connectivities into a big group would include those subsystems with linear arrangement of monomers into one group, whereas for 4-body subsystems, grouping those connectivities with more than three solid lines into one big group would certainly do for this purpose. Obviously, systematic tests are needed to design a scheme with less computational costs yet with high accuracy. In the present study when applying eq 19 to water clusters, unless otherwise noted when computing the n-body correction energies with n > 2, we chose to group those connectivities with zero and only one dashed line into a big group and compute them by the same level of computational method L1. For 2body subsystems, there are only two connectivities, and thus, the connectivity with a solid line, i.e., with a short distance between two water molecules, was computed, whereas the one with a dashed line, i.e., with a long distance between two water molecules, was not computed. Thus, we have

The force acting on an atom can be adaptively computed using a similar equation by simply differentiating eq 14. Other properties such as atomic charges can also be computed adaptively in a similar way, and the equations are not presented here. Applying eq 14, in principle, we can treat each group with a corresponding computational method. Therefore, for n-body subsystems there are a total of (C2n + 1)! combinations given C2n + 1 computational methods. We can also further group the C2n + 1 groups into bigger groups to reduce the levels of the computational methods. Then, we obtain eq 15 if we use G to represent these bigger groups M

dE(n) =

∑ ∑ (SG(n)ΔEi(Li ···) i ) G 12

n

(15)

i1< i 2 ···< in G

with SG(n) =

∑ sg(n) (16)

g∈G

In eq 15, when only two levels, L1 and L2, of computational methods are used, or say to group all connectivities into just two big groups, the expression of dE(n) and the computational implementation of eq 15 are very simple due to the fact that GSFs are normalized M

dE(n) =

∑

1) 2) [S1(n)ΔEi(L + (1 − S1(n))ΔEi(L ] 1i 2 ··· in 1i 2 ··· in

i1< i 2 ···< in

(17)

Because s(n) g is a function of RC, a more general expression for eq 14 is Cn2 + 1

M

dE

(n)

=

∑ ∑∑ i1< i 2 ···< in

k

(L (R (k)))

(sg(n)(R C(k))ΔEi1i2g··· iCn )

g=1

S1(n)

(18)

Lg (R(k) C )

where is the level of the computational method for the (1) (2) (3) group g for the kth cutoff distance R(k) C with RC < RC < RC ···. To further reduce computational costs and reach linear scaling, we can use only one cutoff distance and set the L2 method to be none, i.e., discarding the contribution of those subsystems with long intermonomer distances. Then, eq 17 can be simplified as

∑ i1< i 2 ···< in

(20)

(n) where s(n) 1 and s2 are given by eqs 8 and 9, respectively. This simplified scheme is denoted as GSF-MBE(m)/L1, where L1 is the level of the computational method, which is usually composed of a method and a basis set. For water clusters, a water molecule is chosen as a monomer, and the distance between two monomers is defined as the distance between the two oxygen atoms. Testing using this simplified GSF scheme can help us to determine the essential components of the error induced by using a cutoff.

M

dE(n) =

⎧ s (n) , n=2 ⎪ 1 =⎨ ⎪ s (n) + s (n) , n > 2 ⎩1 2

1) S1(n)ΔEi(L 1i 2 ··· in

(19)

3. COMPUTATIONAL DETAILS We have written a Java program to drive the Gaussian 09 software56 implementing MBE with GFSs (GSF-MBE). The program generates Gaussian input files and then reads output files to compute n-body correction energies and gradients. It should be noted that the energy values in the summary block of Gaussian outputs do not provide adequate numerical precision and one needs to grep the data from the lines that print them in a higher precision such as “SCF Done”. The formatted checkpoint file is the best place to obtain energies, gradients, and other data. To verify the continuity of the PES calculated by the GSFMBE method, we first tested it by a model potential and then an NVE MD simulation of 5 ps at 400 K on the water decamer whose geometry was taken from ref 57. In the MD simulation of the water decamer, the integrator used was velocity Verlet

It is worth noting here that our GSF scheme is a dynamic subsystem scheme: When all the monomers in a subsystem are very far away from each other, we can discard it according to eq 19. If they are close to each other again during a geometry optimization or a molecular dynamics simulation, the subsystem can be computed again. We must emphasize that our GSF scheme, especially eq 18, is a very general scheme. Applying even its simplest version, i.e., eq 19, there are still many choices of how to group all the connectivities into two big groups. For the 2-body subsystems, the choice is very simple because there are only two connectivities. For the 3-body subsystems, one can choose to compute the contribution from the first connectivity and discard the contribution of the other seven connectivities. Alternatively, one can choose to group the one to four connectivities into one group because the three monomers are 2013


Article


scheme or using a larger basis set such as aug-cc-pVTZ.67,68 To examine the effect of geometry on the convergence behavior of dE(n) with respect to RC, we computed the W25-A−W25-D clusters using the ωB97X-D/cc-pVDZ method. We have also tested the convergence behavior of dE(n) with respect to RC for the EE-MBE method proposed by Truhlar and co-workers. The EE-MBE method with our GSF approach is abbreviated as EEGSF-MBE(m). Fixed point charges of −0.778 and 0.389 au were used for the oxygen and hydrogen atoms, respectively, to be the background charges when computing the energies of the monomers and n-body (2 ≤ n ≤ m) subsystems.4 Truhlar and co-workers found that the error of using fixed background charges is negligible.3−5 Using fixed charges also simplifies the computation of energy gradients; otherwise, there are response terms that turn out to be important for energy conservation in MD simulations.69 To examine the convergence behavior of the correlation energy component in dE(n) with respect to RC, W25-A−W25-D were computed at the MP2/aug-cc-pVTZ level of theory. W25-1−W25-4 clusters were computed using the ωB97X-D/cc-pVDZ method to investigate the convergence behavior of dE(n) for the water clusters carved from liquid bulk water. In all these tests, 0.5 Å was used for ΔR throughout. Last, we applied the GSF-MBE(4)/ωB97X-D/cc-pVDZ method with RC = 8.0 Å and ΔR = 0.5 Å to (H2O)60 and (H2O)100 clusters, whose geometries were carved from the liquid bulk water and then optimized using the ωB97X-D/ccpVDZ method. In the above tests on medium and large water clusters, our focus is the convergence behavior of the GSF-MBE(m) method with respect to RC, that is to say, the error introduced by the use of a cutoff. The truncation error of the MBE(m) method itself is not the focus of our present study, which has been examined by other research works, for example in ref 32.

with a step size of 0.5 fs. The initial momenta of the atoms were set randomly according to the Maxwell−Boltzmann distribution at 400 K. At each step of the MD simulation, the energy and gradients were calculated by the GSF-MBE(4) method with three levels of computational methods, denoted as GSFMBE(4)/(L 1:L2 :L3 ). The details of the GSF-MBE(4)/ (L1:L2:L3) method are presented in Table 1. The cutoff distance used is RC = 3.0 Å with ΔR = 1.0 Å. Table 1. Details of the GSF-MBE(4)/(L1:L2:L3) Method Used for the MD Simulation of the Water Decamer nbody

number of reverse weights in s(n) g

connectivity group g

computational level

2

0 1 0 1 2 3 0 1 2 3 4 5 6

1 2 1 2 3 4 1 2 3 4 5 6 7

MP2/cc-pVDZ none MP2/cc-pVDZ

3

4

HF/cc-pVDZ none MP2/cc-pVDZ HF/cc-pVDZ none

Then, we tested the effect of the choice of RC on the convergence behavior of dE(n) computed by the GSF-MBE(4)/ L1 method for the (H2O)25 clusters, which are W25-I, W25-A, W25-B, W25-C, W25-D, W25-1, W25-2, W25-3, and W25-4. W25-I is the global minimum at the MTA-MP2/aug-cc-pVDZ level reported by Sahu and co-workers.48 W25-A−W25-D are four low-energy isomers reported by Furtado and co-workers.58 W25-1−W25-4 are four isomers carved from a snapshot of an equilibrated simulation (T = 298 K and p = 1 atm) with 8000 rigid water molecules using the TIP4P59 force field with periodic boundary conditions and then optimized by the M052X/cc-PVDZ method. The MD simulation was performed with the velocity Verlet integrator and a time step of 2.0 fs. To investigate the effect of different computational methods on the convergence behavior of dE(n) with respect to RC, we computed the W25-I cluster with the HF, MP2,60 ωB97X-D,61 and M052X62 methods in combination with the cc-pVDZ63 basis set. The MP2 method is believed to give results almost as accurate as the expensive CCSD(T) method for water clusters.40 ωB97X-D and M05-2X methods are reported by Truhlar and co-workers to be the two accurate methods for water clusters.64 The use of a moderate cc-pVDZ basis set is a low-cost choice, but it brings large basis set superposition error (BSSE) to many-body interaction energy by stabilizing the subsystems more than those of separated monomers.65,66 However, because of the fact that higher-order correction terms neglected in the MBE(m) method tend to also be stabilizing,66 a fortuitous cancellation of errors due to the use of a small basis set and the truncation of higher-order terms may improve the accuracy of a low-order MBE(m) method. Anyhow, in the present study we focus mainly on the errors induced by cutoff. BSSE is not taken into consideration, and therefore, a moderate basis set was used. For the actual application of MBE(m) method to water clusters, however, one needs to reduce or eliminate BSSE by either employing a counter-poise correction

4. RESULTS AND DISCUSSION 4.1. Smoothness of the PES. We designed a test to verify that GSFs can adapt multilevel methods as well as keep the PES smooth. Here, we tested the smoothness of the 3-body correction energy of a subsystem (i1i2i3), ΔEi1i2i3, computed employing eq 15. We grouped the eight connectivities (Figure 1) into three groups: G = 1:1−4, G = 2:5−7, and G = 3:8. Assuming the 3-body correction energy of the subsystem 1) (i1i2i3) computed at three different levels of method are ΔE(L i1i2i3 (L2) (L3) = 1, ΔEi1i2i3 = 2, and ΔEi1i2i3 = 3 of arbitrary units, respectively, then ΔEi1i2i3 = S1(3) + 2S2(3) + 3S3(3) = s1(3) + s2(3) + 2s3(3) + 3s4(3) (22) (3) s(3) 1 −s4

where are given by eqs 8−11 with n = 3. For 3-body subsystems, there are three intermonomer distances: r1, r2, and r3. Setting RC = 5 and ΔR = 3 of arbitrary units and keeping r3 to be 10 of arbitrary units, we plotted ΔEi1i2i3 as a function of the other two distances in Figure 2. From Figure 2 we can see that our GSFs can smoothly adapt the three energy levels. Therefore, when a subsystem changes its geometry, our GSF scheme can smoothly switch from one computational method to another so that subsystems can be adaptively treated by multilevel methods. It shall be noted that no extra basins or peaks appear on the PES so that our GSFs will not cause problems for geometry optimizations and molecular dynamics simulations. 2014


Article


determines how fast the many-body correction energy computed is switched from one level of computational method to another (for the GSF-MBE(m)/L1 method, the many-body correction energy is switched from nonzero to zero). Moreover, we can see from Table 2 that the σ and MAD values at ΔR = 0.5 Å are very close to the values at ΔR = 1.0 Å. This indicates that ΔR = 0.5 Å is suitable for water clusters. Therefore, ΔR is set to 0.5 Å for all the other tests in the present study because a smaller ΔR results in lower computational cost. 4.2. Convergence Behavior of dE(n) with Respect to RC. For determining the essential components of the error induced by using a cutoff and a reasonable cutoff distance RC for the computation of water clusters, we examine the error of dE(n) (Δ[dE(n)]) as a function of RC. Here, Δ[dE(n)] is defined as the difference between the dE(n) computed by the GSF-MBE(m)/ L1 method with a cutoff distance of RC and that without the use of cutoff, or that is to say, with an infinitively long RC

Figure 2. Imaginary PES for 3-body correction energy given by eq 22 with r3 = 10, RC = 5, and ΔR = 3 of arbitrary units.

Furthermore, another test was carried out by performing an NVE MD simulation for (H2O)10 with the GSF-MBE(4)/ (L1:L2:L3) method with L1 = MP2/cc-pVDZ, L2 = HF/ccpVDZ, and L3 = none (see Table 1 for details of the method). As shown in Figure 3, the total energy of the water cluster is

Δ[dE(n)] = dE(n)(R C) − dE(n)(R C = ∞)

To facilitate the comparison among water clusters of different sizes, Δ[dE(n)] values were further divided by the number of the waters in the cluster, and thus all Δ[dE(n)] values reported are per water molecule. First, we applied the GSF-MBE(4)/L1 method (L1 = HF/ccpVDZ, MP2/cc-pVDZ, ωB97X-D/cc-pVDZ, and M05-2X/ccpVDZ) to the W25-I cluster, which is the global minimum of (H2O)25.48 The error of dE(n) (Δ[dE(n)], n = 2, 3, 4) and the overall error of the GSF-MBE(4)/L1 method (4BE, = Δ[dE(2)] + Δ[dE(3)] + Δ[dE(4)]) as a function of RC are plotted in Figure 4. As shown in the figure, both wave function theory (WFT) methods and DFT methods show similar behavior with respect to RC. The errors dramatically decrease when RC is beyond 5.5 Å. Because of the strong many-body interaction in water,42 the errors converge very slowly with RC and show an oscillating behavior: The error plots show a valley at around 5.5 Å, a peak at around 6.5 Å, and a valley at around 7.2 Å. At 8.0 Å, which is almost the diameter of the cluster, the error begins to reduce to a negligible level. Comparing with Δ[dE(3)] and Δ[dE(4)], Δ[dE(2)] is relatively larger. Second, we applied the GSF-MBE(4)/ωB97X-D/cc-pVDZ method to the W25-A−W25-D clusters to investigate the effect of the water-cluster geometry on the convergence behavior of Δ[dE(n)]. The error profiles are shown in Figure 5. Although similar oscillating behavior of Δ[dE(n)] and 4BE as that on the W25-I cluster is observed on the W25-A−W25-D clusters, the geometry of the cluster does change the convergence behavior. For example, the curve for W25-A and W25-B dE(n) shows a faster convergence than those for W25-I, W25-C, and W25-D. The oscillating behavior in Δ[dE(n)] vs RC in water clusters is primarily due to the fact that water clusters have short-rangeordered structure so that the number of n-body subsystems is not a monotonous function of RC. This is also the reason that different clusters show different convergence behaviors because they have different radial distribution functions. For the four clusters, the mean unsigned error (MUE) of the GSF-MBE(4)/ ωB97X-D/cc-pVDZ method, i.e., mean unsigned 4BE, is 0.11 kcal mol−1 H2O−1 at RC = 6.0 Å, 0.030 kcal mol−1 H2O−1 at RC = 7.0 Å, 0.017 kcal mol−1 H2O−1 at RC = 8.0 Å, and from 8.0 Å, the MUE begins to vanish. The EE-MBE method proposed by Dahlke and Truhlar takes long-range electrostatic interaction and higher order polarization effects into consideration by using background charges

Figure 3. Energy conservation during an MD simulation for (H2O)10 with the GSF-MBE(4)/(L1:L2:L3) method with L1 = MP2/cc-pVDZ, L2 = HF/cc-pVDZ, and L3 = none (see Table 1 for details of the method).

well-conserved during the 5 ps simulation time. The standard deviation (σ) of the total energy is 0.029 kcal/mol with a mean absolute deviation (MAD) of 0.024 kcal/mol. We also performed three NVE MD simulations at 400 K for (H2O)8 with the GSF-MBE(4)/HF/cc-pVDZ method (RC = 3 Å) using three different ΔR values with a time step of 0.5 fs and a total simulation time of 5 ps. The plots of the total energies vs simulation time are presented in the Supporting Information. and the σ and MAD values of the total energies are listed in Table 2. The results show that the GSF-MBE(4)/HF/6-31G method still conserves the total energy well but that the fluctuation in the total energies becomes larger for smaller ΔR. Therefore, for smaller ΔR a smaller time step should be used. This is not surprising because ΔR is a parameter that Table 2. Standard Deviation (σ) and Mean Absolute Deviation (MAD) of the Total Energies of (H2O)8 in the NVE MD Simulation ΔR (Å)

σ (kcal/mol)

MAD (kcal/mol)

1.0 0.5 0.2

0.042 0.036 0.076

0.033 0.038 0.061

(23)

2015


Article


Figure 4. Error of dE(n) (n = 2, 3, 4) and the overall error for the GSF-MBE(4)/L1 method (4BE) for the W25-I cluster with respect to RC: (a) L1 = HF/cc-pVDZ, (b) L1 = MP2/cc-pVDZ, (c) L1 = ωB97X-D/cc-pVDZ, and (d) L1 = M05-2X/cc-pVDZ.

Figure 5. Error of dE(n) (n = 2, 3, 4) and the overall error for the GSF-MBE(4)/ωB97X-D/cc-pVDZ method (4BE) for the four low-energy (H2O)25 isomers: (a) W25-A, (b) W25-B, (c) W25-C, and (d) W25-D.

when computing the energy of a subsystem.3 Thus, an improvement in the cutoff error of the EE-GSF-MBE(m)/L1 method shall be expected. Figure 6 illustrates the errors as a function of RC for the EE-GSF-MBE(m)/ωB97X-D/cc-pVDZ method for the water clusters from W25-A to W25-D. Similar oscillating behavior as those in Figures 4 and 5 is also observed. In contrast to the conventional GSF-MBE(m)/L1 method, Δ[dE(3)] and Δ[dE(4)] converge much faster than Δ[dE(2)] for the EE-GSF-MBE(m)/L1 method. Δ[dE(2)] contributes most to 4BE, and thus, 4BE shows a very similar convergence behavior as that of Δ[dE(2)]. This phenomenon indicates that the use of background charges to compute the energies of subsystems indeed have included the majority of many-body interactions into 2-body terms. For the four clusters, the MUE of the EE-GSF-MBE(4)/ωB97X-D/cc-pVDZ method is 0.016 kcal mol−1 H2O−1 at RC = 6.0 Å, 0.0078 kcal mol−1 H2O−1 at RC = 7.0 Å, and 0.0037 kcal mol−1 H2O−1 at RC = 8.0 Å. It is worth noting here that because the computational bottleneck of

the MBE(m) method is the computation of m-body subsystems, it is possible to further reduce the computational costs while maintaining the accuracy of the EE-GSF-MBE(m)/ L1 method by using a smaller RC for n-body subsystems with n > 2 and a larger RC for 2-body subsystems. This is very easy to realize with our GSF scheme. However, there are still debates on the use of background charges to improve the convergence behavior of the low-order EE-MBE methods.32,70 According to the research work of Lao et al.,70 the error per monomer of the EE-MBE method is in fact a size-intensive quantity and depends on the choice of embedding charges. Further investigation on the convergence behavior of the EE-MBE methods is still necessary. Because correlation energy is typically short ranged as compared to HF energy, it may be reasonable to apply a cutoff only on correlation energy while keeping HF energy untouched.4,46,47,49 Here, we have computed the W25-A− W25-D clusters with the GSF-MBE(4)/MP2/aug-cc-pVTZ 2016


Article


Figure 6. Error of dE(n) (n = 2, 3, 4) and the overall error (4BE) for the EE-GSF-MBE(4)/ωB97X-D/cc-pVDZ method for the four low-energy (H2O)25 isomers: (a) W25-A, (b) W25-B, (c) W25-C, and (d) W25-D.

Figure 7. Correlation-energy components of the errors in dE(n) (Δ[dE(n)(E2)], n = 2, 3, 4) and the overall error (4BE(E2)) for the GSF-MBE(4)/ MP2/aug-cc-pVTZ method for the four low-energy (H2O)25 isomers: (a) W25-A, (b) W25-B, (c) W25-C, and (d) W25-D.

MBE(4)/MP2/cc-pVDZ method (Figure 4b) for the W25-I cluster. Although Δ[dE(n)(E2)] does converge slightly faster than the HF component, our tests indicate that it is not fast enough so that we can use a much shorter RC for correlation energy than that for HF energy to save computational costs. For the four clusters, the mean unsigned 4BE(E2) is 0.14 kcal mol−1 H2O−1 at RC = 6.0 Å, 0.032 kcal mol−1 H2O−1 at RC = 7.0 Å, and 0.0053 kcal mol−1 H2O−1 at RC = 8.0 Å. The reason for the slow convergence of Δ[dE(n)(E2)] is that the pure correlation energy of an interaction cannot be extracted by simply subtracting post-HF methods with the HF method. For example, a 2-body pure electrostatic interaction at the HF/ aug-cc-pVTZ level is different from a 2-body pure electrostatic interaction at the MP2/aug-cc-pVTZ level because the monomers are so much better described at the post-HF level that the better permanent multipoles and polarizabilities

method to investigate the convergence behavior of the secondorder correlation energy (E2) with respect to RC. The correlation-energy components of the errors in dE(n) (Δ[dE(n)(E2)], n = 2, 3, 4) and the overall error (4BE(E2)) for the GSF-MBE(4)/MP2/aug-cc-pVTZ method, which are equivalent to the differences of the errors computed by the GSF-MBE(4)/MP2/aug-cc-pVTZ and GSF-MBE(4)/HF/augcc-pVTZ methods, for the four clusters are plotted as a function of RC in Figure 7. The results indicate that the correlationenergy component in the many-body-correction energies computed by the MBE methods indeed converge slightly faster than the HF component. The lack of oscillating behavior in the plots in Figure 7 indicates that the oscillating behavior observed for the MBE method is mainly caused by the oscillation of the HF component in the many-body-correction energies. This is also confirmed by the very similar convergence behavior of the GSF-MBE(4)/HF/cc-pVDZ method (Figure 4a) and the GSF2017


Article


Figure 8. Error of dE(n) (n = 2, 3, 4) and the overall error (4BE) for the GSF-MBE(4)/ωB97X-D/cc-pVDZ method for the four (H2O)25 isomers carved from bulk: (a) W25-1, (b) W25-2, (c) W25-3, and (d) W25-4.

Table 3. Mean Errors (in kcal mol−1 H2O−1) of the GSFMBE(4)/ωB97X-D/cc-pVDZ Method at Different Cutoff Distances (RC) for the Nine (H2O)25 Clusters

improve the interaction such that a simple subtraction extracts both the correlation energy and the difference of the electrostatic interaction. However, it is worth mentioning here that basis-set size has a large effect on the convergence behavior of Δ[dE(n)(E2)], where results using a smaller basis set indicate that a much shorter RC can be used for correlation energy (see Figure S4 for the results for W25-I using the ccpVDZ basis set). The tests we have presented are all for the low-energy clusters that are less representative of the structures found in the bulk.71 To circumvent this drawback, we computed the other four (H2O)25 clusters, W25-1−W25-4, carved from a snapshot of an MD simulation using the TIP4P potential with the GSF-MBE(4)/ωB97X-D/cc-pVDZ method. Figure 8 illustrates the profiles of Δ[dE(n)] (n = 2, 3, 4) and 4BE with respect to RC. Different from the low-energy clusters, most of the four quantities are positive. This is in agreement with ref 71 that optimized clusters do not have the same structure as the bulk water. In addition, they show an attenuated oscillating behavior compared with those low-energy structures W25-I (Figure 4) and W25-A−W25-C (Figure 5). There are two possible reasons for this phenomenon. One is that the W25-1− W25-4 clusters have a less-ordered structure than those fully optimized low-energy structures. The other is that the electric fields in optimized clusters tend to align inducing larger dipoles, which interact more strongly with the fields, whereas in clusters carved from liquid water, they tend to cancel so as to reduce the induction and many-body interaction. Indeed, W25-1, W25-3, and W25-4 show a slightly faster convergence behavior than those low-energy structures: For the four clusters, the MUE of the GSF-MBE(4)/ωB97X-D/cc-pVDZ method is 0.029 kcal mol−1 H2O−1 at RC = 7.0 Å and 0.013 kcal mol−1 H2O−1 at RC = 8.0 Å. Both values are slightly smaller than those corresponding values for the W25-A−W25-D clusters. Table 3 summarizes the mean errors of the GSF-MBE(4)/ ωB97X-D/cc-pVDZ method for the nine (H2O)25 clusters at different RC values. The errors show an obvious oscillating behavior, but the trend is that longer RCs have smaller error.

a

RC (Å)

5.0

6.0

7.0

8.0

9.0

MSEa MUEb RMSEc

0.21 0.25 0.28

0.037 0.10 0.12

0.011 0.041 0.054

−0.0015 0.017 0.019

−0.0028 0.0060 0.0071

Mean signed error. bMean unsigned error. cRoot-mean-square error.

However, a longer RC also implies higher computational costs because it needs to compute more m-body subsystems. The error of a GSF-MBE(m) method has two sources: the one from RC so that only a part of n-body subsystems (2 ≤ n ≤ m) are computed and the one that is the truncation error of the MBE(m) method itself. Therefore, we need to squeeze the error originated from RC to an acceptable level to say 0.02 kcal mol−1 H2O−1 for MUE, which corresponds to a total MUE of 0.5 kcal mol−1 H2O−1 for the (H2O)25 clusters. For the GSFMBE(4)/ωB97X-D/cc-pVDZ method, an MUE of 0.02 kcal mol−1 H2O−1 corresponds to an RC of ∼8.0 Å. On the other hand, if one is satisfied with an MUE of 0.04 kcal mol−1 H2O−1, i.e., a total of 1.0 kcal/mol for the (H2O)25 clusters, the corresponding RC is 7.0 Å. When RC goes from 7.0 to 8.0 Å, the root-mean-square error (RMSE) decreases 64% while the computational cost roughly increases just 24%. Therefore, an RC of 8.0 Å may be a cost-efficiency choice for the GSFMBE(4)/ωB97X-D/cc-pVDZ method and probably for the other GSF-MBE(4)/L1 method as well because they have similar convergence behavior (Figure 4). We also have performed additional tests on larger water clusters, (H2O)60 and (H2O)100 using the GSF-MBE(4)/ ωB97X-D/cc-pVDZ method with RC = 8.0 Å and ΔR = 0.5 Å. For the (H2O)60 cluster, the error is 0.048 kcal mol−1 H2O−1, whereas for the (H2O)100 cluster, the error is 0.021 kcal mol−1 H2O−1. However, only 43% of the 4-body subsystems are computed for the (H2O)60 cluster, and only 16% of the 4-body subsystems are computed for the (H2O)100 cluster. The results 2018


Article


To investigate the convergence behavior of dE(n) with respect to RC, we applied two versions of our GSF-MBE(m) method, GSF-MBE(m)/(L1:L2:L3) (eq 15 with three computational methods) and GSF-MBE(m)/L1 (eqs 19 and 20), to water clusters. We found that the error induced by RC is not monotonous with respect to RC and shows an oscillating behavior. The error profile with respect to RC strongly depends on the structural conformation of water clusters where optimized clusters show stronger oscillating behavior and clusters carved from liquid bulk show attenuated oscillating behavior. We also found that, because of the much faster convergence behavior of the n-body correction energy with n > 2 for the EE-GSF-MBE(m)/L1 method, a short RC can be used for the n-body correction energy with n > 2. For a series of (H2O)25 clusters, we found that, for the GSF-MBE(4)/ωB97XD/cc-pVDZ method, an RC of 8.0 Å yields an average mean unsigned error of 0.017 kcal mol−1 H2O−1. For (H2O)100, using this cutoff saves ∼84% computational resources with an error of just 0.021 kcal mol−1 H2O−1. Because our GSF-MBE(m) is very general and there are many possible combinations of the choice of computational methods, cutoff distances, and connectivity groups, our future studies will focus on designing the best GSF-MBE(m) methods for the computation of a large variety of large systems. It is also worthwhile to extend our GSF scheme to combine the dynamic monomer scheme of the reactive many-body expansion (RMBE) method39 of Pinski and Csáyi and the dynamic subsystem scheme of the current work into a unified scheme.

indeed indicate that a dramatic amount of computational sources could be saved for larger water clusters without losing much accuracy. To illustrate that the GSF-MBE(m) method can be linear scaling, in Figure 9 we plotted the wall times for computing the

Figure 9. Wall times for computing the 2-, 3-, and 4-body subsystems of water clusters carved from bulk liquid using the GSF-MBE(4)/ MP2/cc-pVDZ method with RC = 8.0 Å and ΔR = 0.5 Å.

many-body subsystems using the GSF-MBE(4)/MP2/ccpVDZ method with RC = 8.0 Å and ΔR = 0.5 Å as a function of water-cluster size. The wall times were estimated from the average computation times of the n-body subsystems in the W25-I cluster using the GSF-MBE(4)/MP2/cc-pVDZ method, where the calculations were performed using a single processor of the E5-2685(v3) Intel CPU. The water clusters were carved from bulk liquid from an MD simulation for which the details are given in the Computational Details. It can be seen that the wall times increase linearly with respect to the number of water molecules in the cluster. From the plots we also notice that the slope for the 4-body subsystems is much steeper than those for the 3- and 2-body subsystems. Therefore, the number of 4-body subsystems are still too many for very large water clusters. The computational costs can be further reduced by using the EEGSF-MBE(m)/L1 method with a shorter RC for n-body subsystems with 2 < n ≤ m.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00144. Energy conservation for the NVE MD simulations of the (H2O)8 cluster and Cartesian coordinates of the water clusters (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected].

5. CONCLUSIONS For the many-body expansion (MBE) method, to reach linearscaling computation, a cutoff distance (RC) must be used. However, the use of a hard cutoff distance will cause discontinuity of the potential energy surface (PES). Herein, we present a generalized-switch-function (GSF) approach to solve this problem. The construction of GSFs is to view an nbody subsystem as a polyhedron with n vertices connected by solid lines and dashed lines where the solid line corresponds to a short distance and a dashed line corresponds to a long distance. A set of these C2n lines then defines a connectivity. These connectivities can be grouped into C2n + 1 groups according to the number of dashed lines in the connectivities in an ascending order. For each group g, we can assign a computational method Lg to it. Thus, for an n-body subsystem (i1i2···in), we compute its n-body correction energy at level Lg, (n) g) ΔE(L i1i2···in, and give it a weight as expressed by the GSF sg (eqs (Lg) 2 8−12). Summation of s(n) g ΔEi1i2···in from g = 1 to Cn + 1 gives the contribution of the subsystem (i1i2⋯in) to the total n-body correction energy dE(n) (eq 14). This is the foundation of our GSF approach to the MBE method.

ORCID

Zhen Hua Li: 0000-0002-5636-9865 Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS This work was supported by the National Natural Science Foundation of China (21573044). REFERENCES

(1) Szabo, A.; Ostlund, N. S. Modern Quantum Chemistry: Introduction to Advanced Electronic Structure Theory; Dover Publications: Mineola, NY, 1996. (2) Dahlke, E. E.; Truhlar, D. G. J. Phys. Chem. B 2006, 110, 10595. (3) Dahlke, E. E.; Truhlar, D. G. J. Chem. Theory Comput. 2007, 3, 46. (4) Dahlke, E. E.; Truhlar, D. G. J. Chem. Theory Comput. 2007, 3, 1342. (5) Dahlke, E. E.; Leverentz, H. R.; Truhlar, D. G. J. Chem. Theory Comput. 2008, 4, 33. (6) Fedorov, D. G.; Kitaura, K. J. Phys. Chem. A 2007, 111, 6904. (7) Addicoat, M. A.; Collins, M. A. J. Chem. Phys. 2009, 131, 104103. 2019


Article

Journal of Chemical Theory and Computation (8) Ř ezác, J.; Salahub, D. R. J. Chem. Theory Comput. 2010, 6, 91. (9) Weiss, S. N.; Huang, L.; Massa, L. J. Comput. Chem. 2010, 31, 2889. (10) Beran, G. J. O.; Nanda, K. J. Phys. Chem. Lett. 2010, 1, 3480. (11) Wang, X.; Liu, J.; Zhang, J. Z. H.; He, X. J. Phys. Chem. A 2013, 117, 7149. (12) Tan, H. J.; Bettens, R. P. A. Phys. Chem. Chem. Phys. 2013, 15, 7541. (13) Fujimoto, H.; Koga, N.; Fukui, K. J. Am. Chem. Soc. 1981, 103, 7452. (14) Gadre, S. R.; Shirsat, R. N.; Limaye, A. C. J. Phys. Chem. 1994, 98, 9165. (15) Kitaura, K.; Ikeo, E.; Asada, T.; Nakano, T.; Uebayasi, M. Chem. Phys. Lett. 1999, 313, 701. (16) Zhang, D. W.; Zhang, J. Z. H. J. Chem. Phys. 2003, 119, 3599. (17) Zhang, D. W.; Chen, X. H.; Zhang, J. Z. H. J. Comput. Chem. 2003, 24, 1846. (18) Li, W.; Li, S.; Jiang, Y. J. Phys. Chem. A 2007, 111, 2193. (19) Hua, S.; Hua, W.; Li, S. J. Phys. Chem. A 2010, 114, 8126. (20) Hua, S.; Li, W.; Li, S. ChemPhysChem 2013, 14, 108. (21) Li, S.; Li, W.; Ma, J. Acc. Chem. Res. 2014, 47, 2712. (22) Deev, V.; Collins, M. A. J. Chem. Phys. 2005, 122, 154102. (23) Collins, M. A.; Cvitkovic, M. W.; Bettens, R. P. A. Acc. Chem. Res. 2014, 47, 2776. (24) Le, H. A.; Tan, H. J.; Ouyang, J. F.; Bettens, R. P. A. J. Chem. Theory Comput. 2012, 8, 469. (25) Huang, L.; Massa, L.; Karle, J. Int. J. Quantum Chem. 2005, 103, 808. (26) Richard, R. M.; Herbert, J. M. J. Chem. Phys. 2012, 137, 064113. (27) Richard, R. M.; Lao, K. U.; Herbert, J. M. J. Chem. Phys. 2014, 141, 014108. (28) Gordon, M. S.; Fedorov, D. G.; Pruitt, S. R.; Slipchenko, L. V. Chem. Rev. 2012, 112, 632. (29) Beran, G. J. O. J. Chem. Phys. 2009, 130, 164115. (30) Hopkins, B. W.; Tschumper, G. S. Mol. Phys. 2005, 103, 309. (31) Góra, U.; Podeszwa, R.; Cencek, W.; Szalewicz, K. J. Chem. Phys. 2011, 135, 224102. (32) Yuan, D.; Shen, X.; Li, W.; Li, S. Phys. Chem. Chem. Phys. 2016, 18, 16491. (33) Saha, A.; Raghavachari, K. J. Chem. Theory Comput. 2014, 10, 58. (34) Song, G.-L.; Li, Z. H.; Fan, K. J. Chem. Theory Comput. 2013, 9, 1992. (35) Collins, M. A.; Bettens, R. P. A. Chem. Rev. 2015, 115, 5607. (36) Wang, Y.; Huang, X.; Shepler, B. C.; Braams, B. J.; Bowman, J. M. J. Chem. Phys. 2011, 134, 094509. (37) Yu, Q.; Bowman, J. M. J. Chem. Theory Comput. 2016, 12, 5284. (38) Babin, V.; Medders, G. R.; Paesani, F. J. Chem. Theory Comput. 2014, 10, 1599. (39) Pinski, P.; Csányi, G. J. Chem. Theory Comput. 2014, 10, 68. (40) Howard, J. C.; Tschumper, G. S. WIREs Comput. Mol. Sci. 2014, 4, 199. (41) Baranyai, A.; Kiss, P. T. J. Chem. Phys. 2009, 131, 204310. (42) Ludwig, R. Angew. Chem., Int. Ed. 2001, 40, 1808. (43) Bartlett, R. J. WIREs Comput. Mol. Sci. 2012, 2, 126. (44) Yoo, S.; Aprà, E.; Zeng, X. C.; Xantheas, S. S. J. Phys. Chem. Lett. 2010, 1, 3122. (45) Aprà, E.; Rendell, A. P.; Harrison, R. J.; Tipparaju, V.; deJong, W. A.; Xantheas, S. S. Liquid water: obtaining the right answer for the right reasons. In Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis, 2009; No. 66. (46) Bates, D. M.; Smith, J. R.; Janowski, T.; Tschumper, G. S. J. Chem. Phys. 2011, 135, 044123. (47) Tschumper, G. S. Chem. Phys. Lett. 2006, 427, 185. (48) Sahu, N.; Gadre, S. R.; Rakshit, A.; Bandyopadhyay, P.; Miliordos, E.; Xantheas, S. S. J. Chem. Phys. 2014, 141, 164304. (49) Wang, K.; Li, W.; Li, S. J. Chem. Theory Comput. 2014, 10, 1546. (50) Pruitt, S. R.; Nakata, H.; Nagata, T.; Mayes, M.; Alexeev, Y.; Fletcher, G.; Fedorov, D. G.; Kitaura, K.; Gordon, M. S. J. Chem. Theory Comput. 2016, 12, 1423.

(51) Anacker, T.; Friedrich, J. J. Comput. Chem. 2014, 35, 634. (52) Gillan, M. J.; Alfè, D.; Michaelides, A. J. Chem. Phys. 2016, 144, 130901. (53) Xantheas, S. S. J. Chem. Phys. 1995, 102, 4505. (54) Heyden, A.; Truhlar, D. G. J. Chem. Theory Comput. 2008, 4, 217. (55) Ouyang, J. F.; Bettens, R. P. A. J. Chem. Theory Comput. 2016, 12, 5860. (56) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision C.1; Gaussian, Inc.: Wallingford, CT, 2009. (57) Temelso, B.; Archer, K. A.; Shields, G. C. J. Phys. Chem. A 2011, 115, 12034. (58) Furtado, J. P.; Rahalkar, A. P.; Shanker, S.; Bandyopadhyay, P.; Gadre, S. R. J. Phys. Chem. Lett. 2012, 3, 2253. (59) Jorgensen, W.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (60) Møller, C.; Plesset, M. S. Phys. Rev. 1934, 46, 618. (61) Chai, J.-D.; Head-Gordon, M. Phys. Chem. Chem. Phys. 2008, 10, 6615. (62) Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Theory Comput. 2006, 2, 364. (63) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (64) Leverentz, H. R.; Qi, H. W.; Truhlar, D. G. J. Chem. Theory Comput. 2013, 9, 995. (65) Ouyang, J. F.; Cvitkovic, M. W.; Bettens, R. P. A. J. Chem. Theory Comput. 2014, 10, 3699. (66) Richard, R. M.; Lao, K. U.; Herbert, J. M. J. Chem. Phys. 2013, 139, 224102. (67) Ouyang, J. F.; Bettens, R. P. A. J. Chem. Theory Comput. 2015, 11, 5132. (68) Richard, R. M.; Lao, K. U.; Herbert, J. M. J. Phys. Chem. Lett. 2013, 4, 2674. (69) Brorsen, K. R.; Minezawa, N.; Xu, F.; Windus, T. L.; Gordon, M. S. J. Chem. Theory Comput. 2012, 8, 5008. (70) Lao, K. U.; Liu, K.-Y.; Richard, R. M.; Herbert, J. M. J. Chem. Phys. 2016, 144, 164105. (71) Friedrich, J.; Yu, H.; Leverentz, H. R.; Bai, P.; Siepmann, J. I.; Truhlar, D. G. J. Phys. Chem. Lett. 2014, 5, 666.

2020


Generalized Switch Functions in the Multilevel Many-Body Expansion

Recommend Documents