At what chain length do unbranched alkanes prefer folded conformations?
aa r X i v : . [ phy s i c s . c h e m - ph ] F e b At what chain length do unbranched alkanes prefer folded conformations?
Jason N. Byrd, a) Rodney J. Bartlett, and John A. Montgomery, Jr. Quantum Theory Project, University of Florida, Gainesville, FL 32611 Department of Physics, University of Connecticut, Storrs, CT 06269
Short unbranched alkanes are known to prefer linear conformations, while long unbranched alkanes are folded.It is not known with certainty at what chain length the linear conformation is no longer the global minimum.To clarify this point, we use ab initio and density functional methods to compute the relative energies of thelinear and hairpin alkane conformers for increasing chain lengths. Extensive electronic structure calculationsare performed to obtain optimized geometries, harmonic frequencies and accurate single point energies for theselected alkane conformers from octane through octadecane. Benchmark CCSD(T)/cc-pVTZ single point cal-culations are performed for chains through tetradecane, while approximate methods are required for the longerchains up to octadecane. Using frozen natural orbitals to unambiguously truncate the virtual orbital space,we are able to compute composite CCSD FNO(T) single point energies for all the chain lengths. This approx-imate composite method has significant computational savings compared to full CCSD(T) while retaining ∼ .
15 kcal/mol accuracy compared to the benchmark results. More approximate dual-basis resolution-of-the-identity double-hybrid DFT calculations are also performed and shown to have reasonable 0 . − . Angew. Chem. Int. Ed.
I. INTRODUCTION
Unbranched alkane chains (C n H n +2 ) are of funda-mental importance in organic chemistry. They are con-stituents of fossil fuels and polymers, as well as importantstructural motifs in lipids and other biomolecules. It isclearly important to understand their conformational andthermochemical properties.The conformer potential energy surface of an un-branched alkane is characterized by torsional twistswhich lead from linear chains to highly deformed struc-tures dominated by intramolecular dispersion forces. Attemperatures less than 300 Kelvin, short alkanes ( n =4 −
8) in the gas phase are well known to prefer the lin-ear all-trans ( T = 180 ◦ , X = 90 ◦ and G = 60 ◦ for trans,cross and gauche dihedral angles respectively) conforma-tion. However, as the length of the alkane grows theremust be a point where the attractive intramolecular in-teractions will cause the chain to self-solvate into a foldedconformer. A cross-gauche-cross rotation combination( T . . . XGX . . . T ) is sufficient to fold the chain, but thiscreates an energetically unfavorable syn -pentane like con-formation. In addition, the chain ends are not parallelin this conformation, reducing the possible stabilizationdue to van der Waals attraction. A hairpin conformationwith four gauche rotations ( T . . . GGT GG . . . T ) mini-mizes the number of strained bonds and allows an en-ergetically favorable parallel arrangement of the chainends, leading it to be the suggested global minimum forlonger alkanes.
These three conformational structuresare illustrated in Figure 1. a) Electronic mail: [email protected]fl.edu
It is well known that the computation of relative en-ergies involving weak interactions presents a significantchallenge for computational studies. In the case of shortalkanes ( n = 4 − while ab initio methods past secondorder perturbation theory are necessary for a full de-scription of the more configurationally complicated tran-sition state conformer structures. As the length of thealkane chain increases so does the importance of a proper
FIG. 1. Illustrative optimized alkane structures. In the case of octane ( n = 8)second order perturbation theory calculations underesti-mate the “bowl” conformer ( GGT GG ) energy differencewhile coupled cluster theory (CCSD) will overestimatethe energy difference by ∼ . us-ing force field and semi-empirical calculations suggesteda turning point anywhere between n = 12 to n = 26,with subsequent force field calculations also pointingtowards n = 18 as the critical chain length. This prob-lem was experimentally addressed in the recent work ofL¨uttschwager et al. Their experiment used Raman spec-troscopy with a supersonic jet expansion apparatus andconcluded that the critical chain length is between n = 16and n = 18 at temperatures of 100K. The accompanyingtheoretical work uses a local coupled cluster approach toalso suggest a critical chain length of n = 18.In this paper we will first obtain benchmark ab initio structures, electronic and harmonic vibrational energiesfor the linear (all-trans) and hairpin alkane conformers ofincreasing length starting with octane ( n = 8) throughtetradecane ( n = 14). These benchmark values will beused to characterize various approximate methods thatare extendable to longer chains ( n >
14) with which wepredict the critical alkane chain length. Although theyare not the focus of this study, entropic effects becomeimportant as the temperature increases, and this is dis-cussed briefly in the conclusions.
II. ELECTRONIC STRUCTURECALCULATIONS
We perform ab initio and DFT electronic structure cal-culations on the n-alkane ( n = 8 , −
18) linear andhairpin conformers using the GAMESS and ACES III quantum chemistry packages running on the Universityof Florida HPC, HiPerGator and University of Connecti-cut BECAT clusters. All calculations in this work useDunning’s correlation consistent family of basis sets (cc-pVnZ, n=T,Q). For resolution-of-the-identity calcula-tions the triple-zeta fitting basis set of Weigend et al. (cc-pVTZ-RI) are used in conjunction with the standardcc-pVTZ basis set. Unless explicitly stated all correlationcalculations in this work assume the frozen-core approx-imation where all 1 s carbon orbitals are dropped fromthe correlation space while the corresponding last virtualorbital is retained.For shorter ( n = 2 −
7) alkane chains it has beenshown that the quality of the final optimized ge-ometry is more strongly dependent on the level of thecorrelation theory than on the basis set. Additional testsshow that including correlation beyond that of second or- der Møller-Plesset perturbation theory (MP2) is unneces-sary. For the longer alkane species we find that the use ofa small basis set (such as the Pople split valence 6-311G*basis set) that has insufficient polarization functions willlead to erroneous hairpin structures. Therefore geome-try optimizations of the alkane conformers used in thiswork are obtained using the MP2 level of theory (withanalytic gradients) and the cc-pVTZ basis set. Harmoniczero point energy (ZPE) shifts are computed numericallyusing analytic first derivatives at the MP2/cc-pVTZ all-electron level of theory. Due to the cost of numericalMP2 hessians, the massively parallel ACES III programis used to perform the necessary first derivatives with thecaveat that only all-electron MP2 gradients are available.This change of theory between geometries and hessian isfound to introduce a negligible error of 0 .
03 kcal/molwhen considering relative conformer ZPE shifts.High level single point energy calculations are com-puted using coupled cluster theory with singles, dou-bles, and perturbative triples (CCSD(T)) and thecc-pVTZ basis set. Higher order effects such as contribu-tions from the core-valence correlation energy or higherorder excitations (full triples, quadruples etc.) are notincluded as their effects are small and cancel nearly iden-tically in conformational energy differences. It is conve-nient to analyze the calculated energies by orders of per-turbation theory. In this way the MP2 correlation energyis given by ∆MP2 = E (MP2) − E (SCF), while higher or-der contributions can be conveniently given as ∆CCSD = E (CCSD) − E (MP2), ∆CCSD(T) = E (CCSD(T)) − E (MP2), and ∆(T) = E (CCSD(T)) − E (CCSD). Be-cause the basis set convergence of the post-MP2 correla-tion is much faster than the second order contribution ,we can estimate the effects of going to the complete basisset (CBS) limit by combining large basis MP2 energiesand small basis coupled cluster correlation energies. Us-ing the cc-pVTZ and cc-pVQZ basis sets, we separatelyextrapolate the SCF and MP2 correlation energy usingthe linear extrapolation formulas of Schwenke E ∞ (SCF) = E n − (SCF)+ F n − ,n ( E n (SCF) − E n − (SCF))(1)(with F , = 1 . et al. ∆ ∞ MP2 = n ∆ n MP2 − ( n − ∆ n − MP2 n − ( n − (2)respectively. Here E n (SCF) and ∆ n MP2 refers to theSCF and MP2 correlation energy computed with the cc-pVnZ basis set. Adding in the coupled cluster correlationenergy to form a composite CBS energy E (CCSD(T)) / CBS = E ∞ (SCF)+∆ ∞ MP2+∆CCSD(T)(3)we obtain a 0 . n = 8 − n =14 starts to become prohibitively expensive (over 18002 ABLE I. Benchmark ab initio conformer energy differ-ences (hairpin − linear, in kcal/mol) for the n = 8 − H C H C H C H cc-pVTZSCF 5.19 5.50 5.88 9.76∆MP2 -4.10 -4.76 -5.59 -10.47∆CCSD 0.83 1.00 1.18 2.11∆(T) -0.49 -0.58 -0.70 -1.41ZPE 0.58 0.60 0.58 0.69cc-pVQZSCF 5.24 5.56 5.95 9.89∆MP2 -4.07 -4.70 -5.52 -10.50CCSD(T) 1.44 1.18 0.79 0.03CCSD(T)/CBS 1.57 1.36 0.99 0.14 basis functions and 130 electrons) we instead omit theCBS correction and include a 0 . O ( o v ) (for o oc-cupied and v virtual orbitals) scaling of the perturbativetriples quickly becomes problematic. With our availablecomputational resources the perturbative triples contri-bution could be computed only for alkane chains up to n = 14. In order to alleviate the cost of includingtriples it is possible to truncate the virtual space bysome amount p , providing a p prefactor which can en-able larger calculations to be done. To do so system-atically and unambiguously we use the frozen naturalorbital (FNO) method which uses the MP2 densitymatrix to make new virtual orbitals. The Hartree-Fockvirtual orbitals can then be replaced with the appropri-ately transformed MP2 virtual natural orbitals, resultingin a set of virtual orbitals sorted by their contributionto the correlation energy. The virtual space can thenbe truncated by examining the MP2 virtual occupationnumbers (eigenvalues of the MP2 density matrix) anddropping orbitals with an occupation smaller than somepredetermined threshold. In this work we take a thresh-old of 1 × − , which results in 40% of the virtual or-bitals being dropped (a prefactor p of 0 . E (FNO CCSD(T)) = E (MP2) + ∆ FNO
CCSD(T) (4)and E (CCSD FNO(T)) = E (CCSD) + ∆ FNO (T) . (5)both of which have comparable accuracies for smaller sys-tems (see Table II). We choose to use the latter composite method (Equation 5), which makes approximations onlyin the triples calculation. Our particular implementationin GAMESS and ACES III uses the converged T and T amplitudes from an ∆ FNO
CCSD calculation to computethe perturbative triples contribution. The error associ-ated with using these amplitudes compared to completevirtual space CCSD amplitudes which are then truncatedby the FNO prescription is small. Similar compositemethods have been used by DePrince and Sherrill withcomparable accuracy obtained. An approximate theoretical method that scales bet-ter than the usual CCSD O ( o v ) calculation with veryreasonable accuracy for short alkanes is double-hybriddensity functional theory (DH-DFT). This methodmixes SCF and DFT exchange with DFT and MP2correlation energy then corrects the dispersion energyempirically using Grimme’s D3 correction. To facili-tate calculations of even larger molecules, we have re-cently implemented in GAMESS the dual-basis SCF method. Here the SCF energy is approximated bya converged small (truncated) basis energy calculation.The large (with polarization functions) basis contribu-tion is then approximated by a single new Fock matrixconstructed from the projected small basis density ma-trix. This approximate SCF method is much faster thana full SCF calculation with errors comparable to standarddensity-fitted SCF methods. We evaluate the DH-DFTmethod in GAMESS using dual-basis DFT married withthe resolution-of-the-identity MP2 method (referred to-gether with the dual-basis SCF to as DB-RI), a furtherapproximation that adds trivial errors while shifting allof the leading computational cost to the single large basisFock matrix build. III. COMPUTATIONAL RESULTS ANDDISCUSSION
When computing relative conformer energies for short( n ≤
6) alkane chains, values taken from MP2 level calcu-lations are sufficient to give 0 .
15 kcal/mol RMS accuracy (taking CCSD(T) values as the correlation benchmark).Even for more demanding structures such as transitionstates (pentane for example) an RMS of 0 . . − . With our available computing resources we are ableto compute accurate ab initio coupled cluster conformerenergies (CCSD(T)/cc-pVTZ) for the alkane chains n =8 −
14, requiring 25,000 CPU hours for each n = 12 con-former and 65,000 CPU hours for each n = 14 conformer3 ABLE II. Computed ab initio conformer energy differences (hairpin − linear, in kcal/mol) for the n = 8 −
18 alkane chainsusing the cc-pVTZ basis set. Both FNO composite energies are shown, as well as the temperature dependent enthalpy changes.Method C H C H C H C H C H C H SCF 5.19 5.50 5.88 9.76 10.11 11.06∆MP2 -4.10 -4.76 -5.59 -10.47 -11.70 -13.77∆CCSD 0.83 1.00 1.18 2.11 2.34 2.75∆
FNO
MP2 -3.80 -4.32 -5.09 -9.62 -10.71 -12.61∆
FNO
CCSD -0.13 0.94 1.16 2.05 2.27 2.66∆
FNO (T) -0.59 -0.49 -0.59 -1.23 -1.39 -1.64CCSD FNO(T) 1.51 1.25 0.89 0.18 -0.63 -1.60FNO CCSD(T) 1.50 1.20 0.86 0.12 -0.70 -1.69∆ H K H K H K (see Figure 2. Obtaining accurate values for this widerange of chain lengths allows us to benchmark any furtherapproximate methods that we choose to use in order toextend our analysis to longer alkane chains. The break-down of the post SCF correlation contributions throughCCSD(T) for these medium length alkanes can be foundin Table I. Two systematic trends can be noticed here:the MP2 conformer energy consistently underestimatesby 0 . . This illustrates the size of bothinfinite order singles and doubles, providing a much morecomplete dispersion contribution, and the significant roleof connected triples in these extended systems.As mentioned earlier, while we are able to performCCSD energy calculations for all the chain lengths underconsideration, the O ( o v ) cost of the perturbative triplesbecomes untenable for chains longer than n = 14 on avail-able computational resources. Therefore we opted for anapproximate treatment using the CCSD FNO(T) virtualspace truncation scheme. Numerical tests show thatan occupation number threshold of 1 × − is sufficient todrop approximately 40% of the virtual space while retain-ing ∼ . . n = 14. A beneficial cance-lation of error can be observed here where the FNO com-posite method predicts conformer energies between thatof the CCSD(T) value and the CCSD(T)/CBS extrapo-lated CBS estimate. Because of computational cost con-siderations we include this ± . n = 8 −
14) at the all-electron MP2 level of theory,with the results presented in Table I. Thermodynamicconsiderations are taken into account by computing thechange in internal energy (vibration, rotation and trans-lation) as a function of temperature. The increase in number of degrees of freedom as the chains lengthen sig-nificantly increases the computational cost, necessitatingan approximate ZPE shift for the longer alkanes. Ourfinal CCSD FNO(T) relative conformer energies withZPE and temperature dependent shifts (simply referredto as ∆ H ) are given in Table II. Noting the nearly con-stant ZPE and temperature shift in ∆ H for each of theshorter alkanes, we take as an approximation that the∆ H shift for the longer ( n >
14) alkanes is the same asfor the n = 14 alkane. Because of this added approxi-mation to the final relative conformer energy we increasethe estimated error bars to 0 . FIG. 2. Total CPU hours used to compute the correlationenergy for a single alkane conformer. Also listed above eachtiming bar is the number of processors used. Calculations forthe n = 8 −
12 and n = 16 chains were performed on theHPC cluster while the n = 14 and n = 18 chain calculationswere performed on the new HiPerGator cluster. The differ-ence between clusters being the increased number of avaliableprocessors and inter-node communication bandwidth. IG. 3. Calculated enthalpy differences, ∆ H = ∆ H hairpin − ∆ H linear , using CCSD FNO(T)/cc-pVTZ single point energiesand MP2/cc-pVTZ harmonic vibrational frequencies includ-ing temperature dependent shifts. alkane chains. Our final ∆ H values (plotted in Figure3) show that the n = 18 alkane chain the hairpin con-former is definitely preferred over the linear structure bya full kcal/mol. For low temperatures (0 ∼
100 K) the n = 16 hairpin conformer is possibly preferred with theestimated error bars extending on either side of the 0 line,however as the temperature increases our ∆ H K valuestrongly suggests that the n = 16 hairpin is preferred.These results are in complete agreement with the gas jetexperimental (performed at 100K) work of L¨uttschwager et al. where the hairpin preference is found to be be-tween n = 16 and 18.With the FNO coupled cluster calculations costing12% of the corresponding full virtual space calculations,we were able to perform CCSD FNO(T) calculations forall the alkane chains through n = 18. Even so, the O ( o v ) scaling remains such that the n = 18 chain re-quired 60,000 CPU hours for each conformer to obtainthe relative energy. Clearly the attractiveness of an ap-proximate theory with a reduced computational scaling isgreat. Qualitatively some force field and semi-empiricalcalculations perform well, with OPLS-AA and MM2 both predicting that the n = 18 conformer energeticallyprefers the hairpin with PM3 predicting n = 12 (this isexcluding ZPE and other thermodynamic shifts). How-ever other commonly used force field methods do not per-form as well such as MM3 and AMBER, which predict the hairpin turning point at n = 25 and n = 26 respec-tively. Taking advantage of the polymer like repetitivestructure of the long alkanes L¨uttschwager et al. haveperformed a composite local coupled cluster (local-CC)correlation calculation for the n = 14 −
22 series, includ- ing ZPE and thermal shifts, and compute that n = 18is the first lowest energy hairpin length. Qualitativelythis is in good agreement with our own results, thoughquantitatively we find the local-CC results to be too highcompared to our CCSD FNO(T) values, but very consis-tent with the DH-DFT results (see Table III and Figure4).Density functional theory is not typically a good choicefor non-covalently bonded and weakly interacting con-former studies (or any system where dispersion is animportant consideration). However with the inclusionof London type dispersion through an empirically de-rived additive correction (in this case using Grimme’s-D3 function ) DFT methods can be quantitative within the limits of the training set. In our previouswork the accuracy of the B2-PLYP , B2GP-PLYP and DSD-BLYP double-hybrid DFT composite func-tionals using the dual-basis resolution-of-the-identity approximation (see the previous methods section) weretested for a variety of non-covalently bonded dimers andsmall alkane conformers (so called S22 and ACONF test set). The results of which were promising, with anaverage RMS error of 0 . . n = 4 − n = 8 − ab initio results very well, the best curve comingfrom the DSD-BLYP functional with a nearly constant0 . IV. SUMMARY AND CONCLUSIONS
As the length of unbranched alkane chains reachessome critical length, intramolecular dispersion forcescause a self-solvation effect in which the chains assume afolded conformation. To accurately determine this crit-ical chain length, linear and hairpin alkane conformerstructures were optimized using the MP2/cc-pVTZ levelof theory for chains of length up through n = 18. Bench-mark CCSD(T)/cc-pVTZ single point energy calcula-tions were then performed for octane through tetrade-cane using the ACES III massively parallel quantumchemistry package. Harmonic zero point energies andtemperature shifts were computed using the MP2/cc-pVTZ level of theory.For chains longer than n = 14 it was necessary touse more approximate methods to obtain conformer en-ergy differences. It was found that our CCSD FNO(T)5 ABLE III. Relative conformer energies (hairpin − linear, in kcal/mol) for the n = 8 −
18 alkane chains computed using variouscomposite methods with the cc-pVTZ basis set. Also shown are the local-CC values from L¨uttschwager et al. BenchmarkCCSD(T) energies are included for avaliable chain lengths. Here * denotes the reported hairpin/linear crossing.Method C H C H C H C H C H C H CCSD(T) 1.44 1.18 0.79 0.03CCSD FNO(T) 1.51 1.25 0.89 0.18 -0.63 -1.60B2-PLYP-D3/DB-RI 1.79 1.55 1.17 0.47 -0.39 -1.38B2GP-PLYP-D3/DB-RI 1.81 1.59 1.25 0.71 -0.10 -1.03DSD-BLYP-D3/DB-RI 1.65 1.40 1.03 0.35 -0.49 -1.49local-CC *PM3 * method which takes the full CCSD correlation energy andadds the perturbative triples correlation energy takenfrom a frozen natural orbital calculation where retain-ing only 60% of the virtual space is required to ob-tain results comparable to full CCSD(T) results (Ta-ble II). We have also explored the effectiveness of thedual-basis resolution-of-the-identity double-hybrid den-sity functional theory approach to this problem. Witha computational cost several orders of magnitude lessthan the more rigorous ab initio methods considered herewe find that these approximate DFT methods performedwell with errors ∼ . n = 18 show that the temperaturedependent hairpin preference takes place at n ≥ ∼ . FIG. 4. Composite relative single point energies, E hairpin − E linear , compared against benchmark CCSD(T)/cc-pVTZ val-ues. L¨uttschwagger et al. As the temperature increases, entropic effects will be-come important. The Gibbs free energy difference be-tween linear and folded chains will tend to decrease dueto the increased entropy of the folded conformation. Anaccurate assessment of this effect would require confor-mational sampling of the CCSD(T) energy surface, acalculation that greatly exceeds the computational re-sources available to us. However, the close agreement ofour calculated results with the low temperature resultsof L¨uttschwagger et al. suggests that our conclusionswould be significantly altered by entropic effects only atmuch higher temperatures. V. ACKNOWLEDGEMENTS
JNB and RJB would like to acknowledge funding sup-port from the United States Army Research Office grantnumber W911NF-12-1-0143 and ARO DURUP grantnumber W911-12-1-0365.
REFERENCES S. Grimme, J. Antony, T. Schwabe, and C. M¨uck-Lichtenfeld, Organic & Biomolecular Chemistry , 741(2007). J. Goodman, J. Chem. Inf. Comp. Sci. , 876 (1997). L. L. Thomas, T. J. Christakis, and W. L. Jor-gensen, The Journal of Physical Chemistry B ,21198 (2006). N. L¨uttschwager, T. Wassermann, R. Mata, andM. Suhm, Angew. Chem. Int. Ed. , 463 (2013). D. Gruzman, A. Karton, and J. M. L. Martin, J. Phys.Chem. A , 11974 (2009). J. M. L. Martin, J. Phys. Chem. A , 3118 (2013). M. W. Schmidt and K. K. Baldridge and J. A. Boatzand S. T. Elbert and M. S. Gordon and J. J. Jensenand S. Koseki and N. Matsunaga and K. A. Nguyenand S. Su and T. L. Windus and M. Dupuis and J. A.Montgomery, Jr. , J. Comput. Chem. , 1347 (1993).6 M. S. Gordon and M. W. Schmidt, in
Theory and Ap-plications of Computational Chemistry, the first fortyyears , edited by C. E. Dykstra, G. Frenking, K. S. Kim,and G. E. Scuseria (Elsevier, Amsterdam, 2005) pp.1167–1189. V. Lotrich, N. Flocke, M. Ponton, A. Yau, A. Perera,E. Deumens, and R. Bartlett, J. Chem. Phys. ,2722 (2008). T. Dunning Jr, J. Chem. Phys. , 1007 (1989). F. Weigend, A. K¨ohn, and C. H¨attig, J. Chem. Phys. , 3175 (2002). J. B. Klauda, B. R. Brooks, A. D. MacKerell, R. M.Venable, and R. W. Pastor, The Journal of PhysicalChemistry B , 5300 (2005). G. D. Purvis III and R. J. Bartlett, J. Chem. Phys. ,1910 (1982). K. Raghavachari, G. W. Trucks, J. A. Pople, andM. Head-Gordon, Chem. Phys. Lett. , 479 (1989). R. Bartlett and M. Musia l, Rev. Mod. Phys. , 291(2007). J. A. Platts, J. G. Hill, K. E. Riley, J. ˇRez´aˇc, andP. Hobza, J. Chem. Theory Comput. , 330 (2013). D. Schwenke, J. Chem. Phys. , 014107 (2005). T. Helgaker, W. Klopper, H. Koch, and J. Noga, J.Chem. Phys. , 9639 (1997). A. G. Taube and R. J. Bartlett, Coll. Czech. Chem.Comm. , 837 (2005). A. G. Taube and R. J. Bartlett, J. Chem. Phys. ,164101 (2008). A. E. DePrince III and C. D. Sherrill, J. Chem. TheoryComput. , 293 (2012). S. Grimme, J. Chem. Phys. , 034108 (2006). T. Schwabe and S. Grimme, Acc. Chem. Res. , 569(2008). S. Grimme, J. Antony, S. Ehrlich, and H. Krieg, J.Chem. Phys. , 154104 (2010). J. N. Byrd and J. A. Montgomery, Jr., To be submitted. K. Wolinski and P. Pulay, J. Chem. Phys. , 9497 (2003). W. Liang and M. Head-Gordon, J. Phys. Chem. A ,3206 (2004). T. Nakajima and K. Hirao, J. Chem. Phys. , 184108(2006). R. Steele, R. DiStasio, Y. Shao, J. Kong, and M. Head-Gordon, J. Chem. Phys. , 074108 (2006). R. Steele, R. DiStasio, and M. Head-Gordon, J. Chem.Theory Comput. , 1560 (2009). M. Katouda and S. Nagase, Int. J. Quantum Chem. , 2121 (2009). R. Steele, Y. Shao, R. DiStasio, and M. Head-Gordon,J. Phys. Chem. A , 13915 (2006). T. Takatani, E. Hohenstein, M. Malagoli, M. Marshall,and C. Sherrill, J. Chem. Phys. , 144104 (2010). L. ˇSimov´a, J. ˇRez´aˇc, and P. Hobza, J. Chem. TheoryComput. , 3420 (2013). M. Pitoˇn´ak, P. Neogr´ady, J. ˇCern´y, S. Grimme, andP. Hobza, ChemPhysChem , 282 (2009). J. O. Hirschfelder, C. F. Curtiss, and R. B. Bird,
Molecular Theory of Gases and Liquids (Wriley, NewYork, 1954). L. Goerigk and S. Grimme, J. Chem. Theory Comput. , 107 (2010). L. Goerigk and S. Grimme, Chem. Theory Comput. ,291 (2011). R. Podeszwa and K. Szalewicz, J. Chem. Phys. ,161102 (2012). A. Karton, A. Tarnopolsky, J.-F. Lame´ere, G. Schatz,and J. M. L. Martin, J. Phys. Chem. A , 12868(2008). S. Kozuch, D. Gruzman, and J. M. L. Martin, J. Phys.Chem. C , 20801 (2010). P. Jureˇcka, J. ˇSponer, J. ˇCern´y, and P. Hobza, Phys.Chem. Chem. Phys.8