CBS

Article pubs.acs.org/JPCA

Toward an Accurate and Inexpensive Estimation of CCSD(T)/CBS Binding Energies of Large Water Clusters Nityananda Sahu, Gurmeet Singh, Apurba Nandi, and Shridhar R. Gadre* Department of Chemistry, Indian Institute of Technology Kanpur, Kanpur 208 016, India S Supporting Information *

ABSTRACT: Owing to the steep scaling behavior, highly accurate CCSD(T) calculations, the contemporary gold standard of quantum chemistry, are prohibitively difficult for moderate- and large-sized water clusters even with the highend hardware. The molecular tailoring approach (MTA), a fragmentation-based technique is found to be useful for enabling such high-level ab initio calculations. The present work reports the CCSD(T) level binding energies of many low-lying isomers of large (H2O)n (n = 16, 17, and 25) clusters employing aug-cc-pVDZ and aug-cc-pVTZ basis sets within the MTA framework. Accurate estimation of the CCSD(T) level binding energies [within 0.3 kcal/mol of the respective full calculation (FC) results] is achieved after effecting the grafting procedure, a protocol for minimizing the errors in the MTA-derived energies arising due to the approximate nature of MTA. The CCSD(T) level grafting procedure presented here hinges upon the well-known fact that the MP2 method, which scales as O(N5), can be a suitable starting point for approximating to the highly accurate CCSD(T) [that scale as O(N7)] energies. On account of the requirement of only an MP2-level FC on the entire cluster, the current methodology ultimately leads to a cost-effective solution for the CCSD(T) level accurate binding energies of large-sized water clusters even at the complete basis set limit utilizing off-the-shelf hardware.

■

theory is required.14−17 In particular, the MP2 method has been shown to give reliable geometries and energies for hydrogenbonded systems. With the recent advances in computer technology in the past decade or so, the application of MP2 theory to larger-sized water clusters has become feasible. Many investigations, especially from the group of Xantheas18−25 have reported MP2 and CCSD(T) calculations for small to large water clusters using Dunning’s correlated basis sets,21−23 such as aug-cc-pVDZ (aVDZ) and aug-cc-pVTZ (aVTZ) ones. Some large-scale calculations carried out on leadership class computer architectures are listed below. The most timeconsuming (T) part of a CCSD(T)/aVTZ single-point calculation22 for the cage (H2O)24 required ∼3 h on all 223 200 processors of Oak Ridge National Laboratory’s “jaguar” supercomputer with a total of 150 TeraBytes (TB) of memory. Recently, Xantheas and co-workers23 carried out ab initio calculations for (H2O)16 and (H2O)17 at MP2 and CCSD(T) levels of theory. The single-point energies at the CCSD(T)/aVTZ//MP2/aVTZ level of theory for five lowlying isomers of (H2O)16 and two low-lying isomers of (H2O)17 cluster were reported.23 The single-point energy calculation for (H2O)16 at the CCSD(T)/aVTZ level of theory using 120 000

INTRODUCTION The study of molecular clusters bonded by weak intermolecular interactions1−3 such as CH···π, π···π, hydrogen bonding and van der Waals (vdW) forces, is an active research area4 in the contemporary science. Among all the weakly bonded molecular clusters, understanding the structures and energetics of the water clusters represents one of the most fascinating fields in the current scientific endeavors. Several theoretical and experimental investigations have been performed in the last three decades for probing the nature and strength of the interactions present in small- to large-water clusters.5−15 In earlier studies, many research groups have explored small- to large-water clusters using ab initio methods ranging from Hartree−Fock (HF) to second-order Møller−Plesset (MP2) or higher theories. These include the work of Burnham et al.8 utilizing the model potential-based parametrization to study the water clusters (n = 2−21), using first-principles ab initio data for the water dimer. Maheshwary et al.9 reported the structures and stability of (H2O)n clusters for n = 8−20 at HF level of theory. Due to limited computational resources before the 2000, earlier studies on water cluster were restricted mainly to either small-sized ones or HF/DFT levels of theory.10−13 For the purpose of a quantitative description of water clusters, the use of the wave function-based methods that include electron correlation such as MP2 and coupled-cluster singles and doubles with perturbative triples [CCSD(T)] © 2016 American Chemical Society

Received: May 4, 2016 Revised: June 25, 2016 Published: June 28, 2016 5706

DOI: 10.1021/acs.jpca.6b04519 J. Phys. Chem. A 2016, 120, 5706−5714

Article

The Journal of Physical Chemistry A

estimating binding energies. They have reported46 only the BSSE-corrected binding energies for (H2O)16 cluster isomers, although the full calculations of Xantheas and co-workers23 are not corrected for BSSE. Another noteworthy work for the “water community” is the recent study by Medders et al.47 wherein they have reported a potential (MB-pol) for estimating accurate binding energies for water clusters through representation of many-body interactions. However, this method is so far applied to water clusters (H2O)n only for n = 4, 5, and 6 and the results are compared with BSSE-corrected CCSD(T)-F12 full calculations. It would be interesting to assess the accuracy of the method for larger (H2O)n clusters. Very recently, Herbert and co-workers48 have discussed the effect of integral cutoff (τints) and SCF convergence parameter (τSCF) on the accuracy of the binding energies estimated by many-body expansion. The number of tetrameric fragments for many-body expansion for (H2O)n clusters for n = 16, 17, and 25 would turn out to be ∼103 to 104. On the contrary, the numbers of fragments within MTA for (H2O)n clusters for n = 16, 17, and 25 are respectively 94, 82, and 123. In the current work, all the fragment calculations are done with the default parameter τints = 10−10, and τSCF = 10−8 used by the Gaussian09 package. Making these parameters tighter by a factor of 100 each is seen to change the fragment energies only after seven decimal digits. Hence the MTA-based energies are expected to be stable to at least 6 decimal digits. The current work reports binding energies of low-lying isomers of (H2O)16 and (H2O)17 clusters at MP2 and CCSD(T) levels of theory employing MTA in combination with the grafting methodology. As the highlight of the work, in the grafting correction to CCSD(T) level, the MP2 level energies have been used as an starting point for extrapolation to CCSD(T). The CCSD(T) level grafting procedure presented here hinges upon the well-known fact that the MP2 method [which scales as O(N5)] captures most of the electron correlation energy for systems such as water clusters and can be a suitable starting point for approximating the highly accurate CCSD(T) energies [scaling requirement being O(N7)].29,49 In a manner similar to the energy correction, the CCSD(T)/CBS involves the first extrapolation of the MP2 energy to MP2/CBS limit followed by the addition of a CCSD(T) correction using a small basis set. The ultimate goal of the exercise is to produce the CBS limit to binding energies (employing limited off-the-shelf hardware) that agrees with their FC counterparts to within 1 kcal/mol (the currently accepted definition of chemical accuracy).

processors took 3.33 h with a memory usage of 820 MB on each processor.23 However, such calculations have two major drawbacks: (i) such large computational resources are not widely available and (ii) the computational codes should be highly parallelized for the effective use of several thousand processors. Hence the introduction of cost-effective approximate fragmentation methods based on divide-and-conquer strategy is warranted.26−30 To facilitate such high-level calculations, especially the MP2-level ab initio calculations, Gadre and co-workers31−33 developed a fragment-based technique, viz. molecular tailoring approach (MTA), for treating large systems using off-the-shelf high-end hardware, for which the full calculations (FC) on the entire system are either impossible or formidable. Some of the notable investigations on large water clusters include the MP2/ aVDZ level geometry optimization and vibrational frequency calculation of (H2O)n clusters for n = 16, 20, 25, 32 and the CCSD(T)/aVDZ//MP2/aVDZ single-point energy evaluation of various isomers of (H2O)16 and (H2O)17 clusters.34−39 The MTA methodology in combination with the newly proposed grafting procedure has shown its potential for estimating the electronic energies of a variety of molecular clusters to within ∼0.4 mH of their FC counterparts.34−39 As noted by Shields and co-workers,40 MP2 describes hydrogen-bonded systems well, yet a higher-order electron correlation correction in the form of a CCSD(T) calculation is usually necessary to achieve benchmark-quality energies. Besides this, the performance of MP2 method varies for different topology of the hydrogen bonds, pointing to the use of a better method.41 Further, the mathematical incompleteness of finite basis sets leads to errors in the molecular calculations. Therefore, extrapolation to the complete basis set (CBS) limit is required for estimation of the binding energies.42,43 However, the formal scaling of CCSD(T) limits its use to small systems with large basis sets or moderate-sized systems with small basis sets. Furthermore, the CCSD(T) level binding energies are found to be extremely sensitive to the basis set used.44,45 It is already known40,44,45 that use of double-ζ basis set yields inaccurate binding energies for both counterpoise-corrected and uncorrected binding energies. In most cases, corrections calculated using double-ζ basis sets are not close to the CBS values, and two-point extrapolation at least using a triple-ζ basis set is needed to get the correction to converge to the right CBS limit.44,45 In view of this, the estimation of benchmark quality binding energies at CCSD(T) level of theory requires calculations at either the aVTZ or aVQZ basis sets. On the contrary, owing to the requirement of large computational resources in terms of time and memory, the conventional CCSD(T) calculations using either of the above basis sets are not routinely possible. As a notable exception to this, Xantheas and co-workers25 have reported CCSD(T)/aVQZ energies up to (H2O)8 cluster isomers and CCSD(T)/aVTZ energies for (H2O)16 and (H2O)17 ones, although for larger clusters, calculations with such large basis-sets are conspicuous by their absence. Góra et al.46 recently reported a protocol for accurate estimation of binding energies for large water clusters through stratified approximation many-body approach (SAMBA). According to them, the efficiency of this method is due to rapid convergence of the many-body expansion. They have applied this method to obtain the interaction energies of the water clusters (H2O)n, for n = 6, 16, and 24. This SAMBA approach seems to be a promising economical method for

■

METHODOLOGY: COMPUTATIONAL DETAILS This section presents the computational details for estimating accurate binding energies (BEs) of water clusters at the gold standard CCSD(T)1−3 level of theory employing MTA,20−23 assisted by the grafting procedure.34−38 The latter is a method for accurate and inexpensive estimation of CCSD(T) level BEs of large-water clusters at the CBS limit.50 In the present work, five (viz. 4444a, 4444b, boat-a, boat-b, and antiboat) low-lying isomers of the (H2O)16 cluster and two (sphere and 552/5) energetically favorable isomers of (H2O)17 cluster optimized23,38 at the MP2 level of theory employing Dunning’s correlation consistent basis sets, viz. aVDZ and aVTZ, are used. The MP2/aVDZ level optimized initial geometries are taken from the earlier work of Gadre and coworkers,38 whereas the MP2/aVTZ ones are taken from the work of Xantheas and co-workers.23 The CCSD(T)/ aVDZ 5707

DOI: 10.1021/acs.jpca.6b04519 J. Phys. Chem. A 2016, 120, 5706−5714

Article

The Journal of Physical Chemistry A

Figure 1. Low-lying isomers of (H2O)16 and (H2O)17 clusters at the MP2/aug-cc-pVTZ level of theory. See text for details.

small basis set (sb) is added to the MTA energies computed with a large basis set (lb). The accurate estimation of MP2/ aVTZ and MP2/aVQZ level energies has been done employing aVDZ and aVTZ as sb, respectively (eq 2). The basis setindependent nature of the errors in the MTA-derived molecular properties have been utilized for this purpose.53 On the contrary, the grafting correction for CCSD(T)/aVDZ and CCSD(T)/aVTZ level of theory is performed in accordance with eq 3. Unlike the grafting correction within MP2 level, the extrapolation is done within one basis set but with two levels of theory, namely, MP2 and CCSD(T).38 The CCSD(T) level grafting procedure presented here hinges upon the well-known fact that the MP2 method captures most of the electron correlation energy for water clusters and can be a suitable starting point for approximating to the highly accurate CCSD(T) energies.29,49 Because the MP2 part is already built inside the CCSD(T) level calculation, the entire procedure for the grafting correction within CCSD(T) level involves only a two-step computation (eq 3). For instance, the CCSD(T)/ aVDZ level grafting correction uses one MTA-based CCSD(T)/aVDZ and one FC at MP2/aVDZ level whereas the correction for CCSD(T)/aVTZ requires one MTA and one FC calculation, respectively, at the CCSD(T)/aVTZ and MP2/ aVTZ levels of theory.

level calculations on MP2/aVDZ optimized (H2O)n (n = 16 and 17) clusters and CCSD(T)/aVTZ level calculations on MP2/aVTZ optimized water clusters are performed employing MTA. Figure 1 displays the MP2/aVTZ optimized geometries of all these isomers. This is also followed by the MTA-based CCSD(T)/aVTZ level single-point energy calculation at the MP2/aVDZ optimized geometries. For a purpose of estimating the CCSD(T)/CBS level binding energies, single-point energy calculations at the MP2/aVQZ level of theory are performed on both the MP2/aVDZ and MP2/aVTZ level optimal geometries. All the CCSD(T) calculations employing either the aVDZ or aVTZ basis are performed within the MTA framework. Due to different scaling behavior, for all the clusters, pentameric fragments are utilized for the CCSD(T)/aVDZ and CCSD(T)/ aVTZ level calculations whereas a mixture of hexameric and heptameric fragments is used for the MP2 level calculations. MTA is a fragmentation-based method to make such highlevel calculations on large systems possible with minimal hardware by scissoring a large system into small overlapping fragments.31−33 Ab initio calculations are performed on all the fragments, employing Gaussian 09 package51 and the results obtained are patched together to obtain the energy, E, of whole cluster by the use of the inclusion−exclusion principle,31 viz. N

E=

N

∑ Ef

i

i=1

−

N

∑ E f ∩f i

i>j

j

+ ... + ( −1)K − 1

∑

E fi ∩ f j ∩ ...fk

E(GMTA/lb) = E(MTA/lb) + [E(FC/sb)

i > j > ... > k

− E(MTA/sb)]

(1)

Here Ef i denotes the value of E for the ith fragment, Ef i∩f j stands for the energy of the binary overlap between the ith and jth fragments and so on. The sign (−1)K−1 is determined by the number K of the fragments included at each summation. The details of the fragmentation and working principle of MTAbased calculations can be found from the earlier research publications of our group.31−39 It has been found from the earlier works of Gadre and coworkers52 that the errors in the molecular properties and hence in the energies derived from the MTA-based calculations vis-àvis the FC values for a given level of theory are almost independent of the basis set used. This fact is utilized for minimizing such errors within MTA methodology. In grafting procedure, a “correction” estimated with MTA energies using

(2)

E(GMTA)CCSD(T) = E(FC)MP2 − E(MTA)MP2 + E(MTA)CCSD(T)

(3)

Here aVDZ is used as the sb for aVTZ as lb whereas aVQZ uses aVTZ as the sb. On the contrary, the CCSD(T) level grafting procedure works within a similar basis set.38 It has already been noticed34−38 that the difference (δE) between the FC MP2 and GMTA-MP2 total energies (δE) is dramatically reduced upon grafting by an order of magnitude to just about