Article pubs.acs.org/JCTC
Benchmark Relative Energies for Large Water Clusters with the Generalized Energy-Based Fragmentation Method Dandan Yuan,† Yunzhi Li,† Zhigang Ni,† Peter Pulay,‡ Wei Li,*,† and Shuhua Li*,† †
School of Chemistry and Chemical Engineering, Key Laboratory of Mesoscopic Chemistry of Ministry of Education, Institute of Theoretical and Computational Chemistry, Nanjing University, Nanjing 210023, P. R. China ‡ Department of Chemistry and Biochemistry, University of Arkansas, Fayetteville, Arkansas 72701, United States S Supporting Information *
ABSTRACT: The generalized energy-based fragmentation (GEBF) method has been applied to investigate relative energies of large water clusters (H2O)n (n = 32, 64) with the coupled-cluster singles and doubles with noniterative triple excitations (CCSD(T)) and second-order Møller−Plesset perturbation theory (MP2) at the complete basis set (CBS) limit. Here large water clusters are chosen to be representative structures sampled from molecular dynamics (MD) simulations of liquid water. Our calculations show that the GEBF method is capable of providing highly accurate relative energies for these water clusters in a cost-effective way. We demonstrate that the relative energies from GEBFMP2/CBS are in excellent agreement with those from GEBF-CCSD(T)/ CBS for these water clusters. With the GEBF-CCSD(T)/CBS relative energies as the benchmark results, we have assessed the performance of several theoretical methods widely used for ab initio MD simulations of liquids and aqueous solutions. These methods include density functional theory (DFT) with a number of different functionals, MP2, and density functional tightbinding (the third generation, DFTB3 in short). We find that MP2/aug-cc-pVDZ and several DFT methods (such as LC-ωPBED3 and ωB97XD) with the aug-cc-pVTZ basis set can provide satisfactory descriptions for these water clusters. Some widely used functionals (such as B3LYP, PBE0) and DFTB3 are not accurate enough for describing the relative energies of large water clusters. Although the basis set dependence of DFT is less than that of ab initio electron correlation methods, we recommend the combination of a few best functionals and large basis sets (at least aug-cc-pVTZ) in theoretical studies on water clusters or aqueous solutions. doubles with noniterative triple excitations (CCSD(T))15 at the complete basis set (CBS) limit (denoted as CCSD(T)/ CBS) are well-known to provide highly accurate descriptions for various molecular clusters (including water clusters).16 However, due to the high computational cost, this method has only been used to study small water clusters.17−20 In order to extend quantum chemistry calculations to large systems, theoretical chemists have developed various linear scaling quantum chemistry methods. Among these methods, energy-based fragmentation methods received much attention in recent years.21−27 The main idea of these fragmentation methods is to divide a large target molecule or molecular cluster into a series of small subsystems and obtain the ground-state energy (or properties) of this system as the linear combination of the ground-state energy (or properties) of small subsystems, which can be calculated with the existing quantum chemistry packages. Popular energy-based fragmentation methods include, for example, the generalized energy-based fragmentation (GEBF) method,25,28 the systematic molecular fragmentation (SMF) method, 26,29 the molecule tailoring approach
1. INTRODUCTION Water has been a highly active area of both experimental and computational research. Over the last several decades, numerous theoretical methods1−5 have been applied to study water clusters and aqueous solutions, with the purpose of understanding the structure and bonding character of water clusters and liquid water, and the chemical process in aqueous solutions. Among various theoretical methods, density functional theory (DFT)3,4 is the most popular method. However, the DFT results strongly depend on the exchange and correlation functionals.6,7 For example, liquid water from BLYP-based Car−Parrinello molecular dynamics (CPMD)8 simulations is overstructured and has too low self-diffusion coefficient, compared to the observed experimental facts.7 The simple and inexpensive empirical DFTB3 method (the third generation of density functional tight binding)9,10 has also been used for MD simulations of liquid water.11−13 This method is applicable to very large systems but cannot be expected to yield accurate results. On the other hand, ab initio electron correlation methods such as second-order Møller−Plesset perturbation theory (MP2)14 are substantially more accurate for water clusters or liquid water. Higher level electron correlation methods such as coupled cluster singles and © XXXX American Chemical Society
Received: March 17, 2017 Published: May 7, 2017 A
DOI: 10.1021/acs.jctc.7b00284 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX
Article
Journal of Chemical Theory and Computation (MTA),27,30 and so on. They mainly differ in the way of constructing subsystems. A number of high-level ab initio calculations on water clusters have been performed with fragment-based methods.31−34 For example, for 10 low-lying (H2O)20 clusters, we have applied the GEBF method within the explicitly correlated CCSD(T)-F12a method35 at the aug-ccpVDZ basis set to obtain CCSD(T)/CBS relative energies.31 With MTA, CCSD(T) and MP2 relative energies have been obtained for (H2O)16 and (H2O)17 clusters.32 The EE-3B (electronically embedded three-body) method,36 combined with CCSD(T)-F12,37 has also been demonstrated to reproduce CCSD(T)/CBS energies well for (H2O)26 clusters.34 For large (H2O)64 clusters, relative energies obtained with the conventional MP2 method and fragment-based methods at the 6-311++G** basis set have been reported.38 However, CCSD(T) investigations on large water clusters such as (H2O)32 and (H2O)64 with large basis sets have not been performed. Obviously, highly accurate ab initio investigations of large water clusters are very important for us to understand chemical and physical properties of liquid water, because the bonding nature of larger water clusters is much closer to that of bulk water. In this work, the GEBF method is used to study relative energies of large water clusters (H2O)n (n = 32,64) at the CCSD(T)/CBS level. The GEBF method developed by us has been applied to compute energies, geometries, and molecular properties for a broad range of molecular clusters and macromolecules at different theory levels.25,28,31,39−45 For each water cluster, we have taken 10 representative structures from classic MD simulations of liquid water. With these highlevel ab initio relative energies as the reference data, we assess the performance of several theoretical methods including DFT, MP2, and DFTB3 methods that are usually used in simulations of liquid water and aqueous solutions. We hope that the evaluation of these theoretical methods for larger water clusters provides useful information on the accuracy of these methods for liquid water simulations. The paper is organized as follows. In section 2, we first give a brief description of the GEBF method, followed by the computational details. In section 3, we first validate the accuracy of GEBF calculations for two water clusters and then report their CCSD(T)/CBS relative energies. Next, the performance of several theoretical methods is evaluated, with GEBF-CCSD(T) relative energies as the reference data. Finally, a summary is given in section 4.
system. A GEBF calculation generally consists of the following steps: (1) A target system is divided into many fragments. For water clusters under study, each water molecule is naturally defined as a fragment. (2) For each fragment, a primitive subsystem that consists of this fragment (central fragment) and several spatially close fragments within a distance threshold (ξ) is built. In building subsystems, additional hydrogen atoms should be added for valence saturation if a fragment has dangling bonds. For water clusters, this step is not needed. To control the size of primitive subsystems, we also define a parameter η to be the maximum number of fragments in primitive subsystems. If the number of fragments in a subsystem (according to the distance threshold) is more than η, then only (η-1) nearest fragments from the central fragment are retained so that the subsystem contains only η fragments. In addition, if two fragments separated by less than 2ξ do not occur simultaneously in any primitive subsystem, then a primitive subsystem consisting of these two fragments is constructed. (3) Once all primitive subsystems are available, derivative subsystems have to be constructed with the inclusion-exclusion principle, to cancel the overlapping of primitive subsystems. The algorithm of constructing the derivative subsystems has been described in our previous works.28 (4) Each subsystem is placed into the presence of background point charges generated by all atoms outside this subsystem. The atomic point charges on all atoms are obtained by performing conventional calculations on all primitive subsystems with an iterative procedure.28 (5) Perform conventional ab initio calculations with the existing quantum chemistry program on all embedded subsystems and obtain the total energy of the target system through eq 1. For (H2O)32 or (H2O)64 clusters, we have taken 10 representative structures from classic molecular dynamics (MD) simulations for liquid water (NPT ensemble, 300 K, 500 water molecules in periodic cubic box, TIP3P water model46) for study. For each cluster, one representative structure is shown in Figure 1. The coordinates of all structures
2. METHODOLOGY AND COMPUTATIONAL DETAILS In the GEBF approach, the ground-state energy of a large system can be obtained from ground-state energies of a series of electrostatically embedded subsystems28 M
Etotal
M
= ∑ CmEm̃ − (∑ Cm − 1) ∑ m
m
A
∑ A
