水團簇構象穩定性起源和本質的密度泛函理論與量子分子動力學研究

王友娟 趙東波 榮春英,* 劉述斌,2,*

(1湖南師范大學化學化工學院,資源精細化與先進材料湖南省高校重點實驗室,化學生物學及中藥分析教育部重點實驗室,長沙 410081; 2 Research Computing Center,University of North Carolina,Chapel Hill,North Carolina 27599-3420,U.S.A.）

1 Introduction

For a given polyatomic molecule,there often exist a few experimentally accessible conformations.As the number of atoms in a molecule increases,the total number of local minima skyrockets exponentially,thus impossible to enumerate them exhaustively.Natural questions to ask are which one is most stable,why,and what factor or factors dictate the relative stability of these local conformational minima.A convincing answer to these questions is not easy,even for the simplest molecules like ethane and hydrogen peroxide.1With the tools available it is now straightforward to identify which conformation has a lower energy,but to find out what factor or factors contribute to or dominate in its stability is controversial.2-6This is often where disagreements arise.From the physicochemical viewpoint,however,to have a definite and well-accepted answer is essential for our understanding.

In this work,using the octamer water cluster as an example,we investigate the conformational stability of water clusters,trying to understand the nature and origin of their stability.To that end,we employ quantum molecular dynamics to generate a large number of conformations for octamer water clusters and then employ two energy partition schemes recently established from density functional theory(DFT）to pinpoint the principles governing the stability of these species.The key question we want to answer is which interaction or interactions determine the molecular stability.Water clusters are bound together through hydrogen bonds.It is generally believed that the nature of hydrogen bonding is predominantly electrostatic,even though quantum contributions through covalent bonding could also be important.Is it really true that the electrostatic interaction is the dominant factor in a water cluster?Do other effects such as steric and exchange-correlation contributions play a role as well?We will provide our answer to these questions in this study.

2 Methodology and computational details

From the theoretical point of view,using virial theorem,the energy difference ΔE between two stable isomers should satisfy7

where T and V denote the kinetic and potential energies,respectively,of the system in concern.Eq.(1）suggests that the stability difference between two conformers is equal to either the entire kinetic energy difference or half of the total potential energy difference.These energy components are,however,not chemically meaningful.We often wish to obtain insights from such effects as steric,electrostatic,or quantum,which are missing in Eq.(1）.More importantly,Eq.(1）does not work for density functional theory,7-9because a portion of the kinetic energy,Tc[ρ],has already been incorporated in the exchange-correlation energy Exc,making the DFT version of the viral theorem much more complicated.10-13

In DFT,the total energy of a system comes from five different contributions:

where Ts,Vne,J,Exc,and Vnnrepresent the non-interacting kinetic,nuclear-electron attraction,classical electron-electron repulsion,exchange-correlation,and nuclear-nuclear repulsion interactions,respectively.Since Vne,J,and Vnnare electrostatic in nature,these three components can be bundled together,yielding Ee[ρ]=Vne[ρ]+J[ρ]+Vnn.Therefore,the conventional approach to perform the decomposition for the total energy difference in DFT is the following,7-9

where Eestands for the electrostatic energy components.

Recently,we proposed an alternative scheme to perform energy difference partition in the framework of DFT,14

where the total energy difference ΔE comes from the contribution of three independent effects,steric ΔEs,electrostatic ΔEe,and fermionic quantum ΔEq,which results from the exchange and correlation effects among electrons.It has been shown that the energy contribution from the steric effect can simply be expressed by the Weizs?cker kinetic energy,Es≡Tw,with

where ρ(r）and ?ρ(r）are the total electron density and its gradient,respectively.Also,the fermionic quantum energy contribution due to the exchange-correlation effect(because electrons are fermions）,Eq[ρ],is the sum of the conventional exchangecorrelation energy Exc[ρ],which includes a kinetic counterpart of the dynamic electron correlation,and the Pauli energy,15,16EPauli[ρ],which is the contribution to the kinetic energy from the antisymmetric requirement of the many-body wave function required by the Pauli Exclusion Principle.10,11,13It is known in the literature that the Pauli energy is the difference between the non-interacting kinetic energy Tsand the Weizs?cker kinetic energy Tw.17Therefore,

The physical meaning of this new quantification of the steric contribution,Eq.(5）,is based on the introduction of a new reference state,where electrons in atoms and molecules are assumed to behave like bosons.If the density of the hypothetical boson state is the same as that of the fermionic state,ρ(r）,the total wave function of the hypothetical state will be(ρ(r）/N）1/2,where N is the number of electrons.The total kinetic energy of the hypothetical state,from which Weisskopf′s“kinetic energy pressure”14,18for the steric effect is calculated,is simply Eq.(5）.

A few prominent features and properties of this novel quantification of steric contribution have been revealed.For instance,the integrand of Eq.(5）is non-negative everywhere,and thus repulsive in nature.It vanishes for the case of a homogeneous electron gas.It is extensive because it is homogeneous of degree one in density scaling,11,19,20so the larger the system,the larger the steric repulsion.If Bader′s atoms-in-molecules approach is adopted,the steric energy can be partitioned at the atomic and functional group levels as well.Its corresponding steric potential,steric charge,and steric force have been defined and evaluated.14,21Its relationship with information theory has been investigated.22,23This approach has been applied to a number of systems,such as conformational changes of small molecules,1,24,25SN2 reactions,26chained and branched alkanes,27and other systems.28-30Reasonably good trends and linear relationships between theoretical and experimental scales (by Taft）of the steric effect have recently been observed at both group and entire molecular levels.31

In this work,we take the octamer water cluster,(H2O）8,as an example and investigate the conformational stability of water clusters,trying to understand the nature and origin of their stability and to address which interaction or interactions determine the molecular stability of these species.To that end,a large number of conformations for octamer water clusters are needed to perform the two energy partition schemes discussed above.To generate as many different local minima as possible,we employ the quantum molecular dynamics(QMD）approach,which has successfully been used for other purposes elsewhere.32In QMD simulations,which were performed with the NWChem package,33atomic nuclei are treated as Newtonian particles whose forces are obtained from the fully converged electronic structure calculation in the Born-Oppenheimer approximation.The simulation protocol is the following.We start the QMD simulations with a few structures from the literature.After an initial structure optimization of 120 steps,each structure is undergone QMD simulations under 300 K for 100 ps with a step size of 0.5 fs and the leapfrog integration algorithm.We employ a constant temperature ensemble using Berendsen′s thermostat with the temperature relaxation time set to be 2 fs.The cutoff radius for short range interactions is 2.8 nm.SHAKE is disabled and thus all bonding interactions are treated according to the force calculated from quantum mechanics.Trajectories are saved in every 50 femtoseconds.A shell script has been written to extract distinct conformations with a total of at least 0.8 nm derivations from the last selected local minimum.The initial structure was from the literature.34A total of 185 distinct isomers have been obtained from these processes.A few selected low-energy local minima from QMD simulations are shown in Scheme 1.For each structure extracted by the post-processing script,a full geometrical optimization is performed at the level of M062X/aug-cc-PVTZ theory35using the Gaussian 09 package36with tight SCF convergence and ultrafine integration grids.Energy partition analyses are ensued after the optimized structure is obtained at the same level of theory.Different approximate exchange-correlation functionals and basis sets were tested and no substantially different results were obtained(results not shown）.

3 Results and discussion

Fig.1 shows three strong correlations obtained for these systems between the total energy difference ΔE of different octamer water clusters and their energy components.Our first observation is that all these two energy components,ΔExc,ΔEe,and ΔEs,are negative in sign,indicating that they are contributing positively to the molecular stability,because ΔE＜0.In Figs.1(a）and 1(c）,we find that the relative stability ΔE of the water cluster is proportional to the exchange-correlation energy difference ΔExcand to the steric effect difference ΔEs,with the correlation coefficient equal to 0.954 and 0.987,respectively.A less significant correlation between the relative stability and the total electrostatic interaction difference ΔEewith the correlation coefficient R2=0.767 was also observed.The positive slope in these three relationships suggests that these energy components all contribute positively to the relative stability.The less than unit slope in these correlations indicates that these energy components are larger in magnitude than the total energy difference itself.

Scheme 1 A few low-energy structures obtained from QMD simulations in this study

Fig.1 Three strong linear correlations between the total energy difference,ΔE,and three energy components,

A few working principles about the relative stability of this molecular system can be obtained from Fig.1.Notice that in magnitude Eq＞0,Es＞0,Exc＜0,Ee＜0,and E＜0.First,a strong linear correlation between ΔE and ΔEsin Fig.1(c）shows that the more stable a water octamer cluster,the smaller its steric repulsion,suggesting that for a lower energy cluster structure,its steric repulsion should be smaller.Since a smaller steric repulsion also implies smaller size,29this result suggests that more stable clusters are often compact and possess smaller sizes.This working principle of molecular stability can be called the minimum steric repulsion principle.Another principle is from Fig.1(a）,where the exchange-correlation difference ΔExcis proportional to the relative stability ΔE,meaning that the more stable an isomer,the larger its exchange-correlation interaction.This result indicates that stable water clusters prefer to have strong exchange-correlation interactions.This can be called the maximum exchange-correlation interaction principle.For the relationship in Fig.1(b）,the correlation is not as strong as the other two energy components,yet it still appears that the electrostatic interaction in a lower energy structure possesses a stronger electrostatic interaction.This latter point answers the question where or not the electrostatic interactions in water cluster is predominant.What we observe in this study is that the electrostatic interaction is indeed a strong,positive contribution to the stability of water clusters,but its correlation with the relative stability,ΔE,is not as strong as the steric repulsion ΔEsand the exchange-correlation interaction ΔExc.These results also provide inputs for other questions.For example,is the quantum effect(exchange-correlation interactions）important?The answer is certainly yes,as illustrated in Fig.1(a）.In addition,Fig.1(c）adds another factor into the picture of our consideration,that is,the steric effect.This effect has not been previously taken into consideration,but our present results clearly showcased its relevance.Put together,our results in Fig.1 suggest that more stable structures of water clusters prefer to have smaller size and smaller steric repulsion,and at the same time,strong exchange-correlation and electrostatic interactions.

Shown in Fig.2 are two strong correlations between energy components.The first one is between the electrostatic interaction energy difference ΔEeand the total noninteracting kinetic energy difference ΔTswith the correlation coefficient equal to 0.976,and the other is between the Fermionic quantum energy difference ΔEqand the steric energy difference ΔEswith R2=0.999.The second correlation has already been discovered elsewhere,24-28whereas the first one is peculiar only to this system.The two relationships are converse correlations,each with a negative slope,meaning that(i）the non-interacting kinetic energy difference ΔTsand the Fermionic quantum energy difference ΔEqare both positive quantities,contributing negatively to the molecular stability,and(ii）the two energy components involved are canceling one another because ΔEeand ΔEsare negative values.

With the fitted formulas from Fig.2,we have

Together with Eqs.(3）and(4）,there result

Eq.(8）shows that ΔEeis the dominant contributor to ΔE＜0 because the second quantity in this equation is positive(since ΔEs＜0）,whereas in Eq.(7）the contribution comes from both terms,with the governing contributor from ΔExcbut the remnant of ΔEe/ΔTsalso contributing positively to ΔE.The correlation coefficients for Eqs.(7）and(8）are found to be 0.96 and 0.77,respectively.These equations provide us with two different approaches to find out which energy component is the dictating factor in governing the relative molecular stability.

Fig.2 Strong linear correlations(a）between total electrostatic energy difference ΔE eand the kinetic energy difference ΔT s,and(b）between Fermionic quantum energy difference ΔE qand the steric repulsion ΔE s

Fig.3 Correlations between the calculated relative stability of water clusters and the two fitted models using two-variable least-square fitting from the energy decomposition schemes in density functional theory

Given the strong correlations in Fig.2,another way to simplify Eqs.(3）and(4）is to use two of the three quantities in Fig.1,which are found to be positively proportional to the relative stability ΔE,to perform least-square fittings.Fig.3 shows the twovariable fitting results in this manner.Using ΔEeplus ΔExcor ΔEs,much better fits can be obtained,with all quantities contributing positively to ΔE and R2equal to or better than 0.99.In Fig.3(a）,the fitted formula is

where we find that the dominant contribution is from the exchange-correlation interaction with the latter possessing a larger coefficient,whereas in Fig.3(b）,

where we see that the electrostatic terms possesses a larger coefficient than the steric repulsion term and thus ΔEeis the dominant contributor.These results are consistent with what we found in Eqs.(7）and(8）,where ΔExcand ΔEewere shown to play dominant roles in the two energy partition schemes,respectively.

Put together,our present results unambiguously show that there exist clear working principles governing the relative stability for such molecular systems as water clusters.Three energy components,electrostatic,steric,and exchange-correlation,are found to all contribute positively to the molecular stability,with the correlation coefficient of the last two correlations better than 0.95.These relationships demonstrate that a more stable structure possesses less steric repulsion,and stronger Fermionic and exchange-correlation interactions.We also found that there exist strong correlations between energy components,such as ΔEevs ΔTs,and ΔEsvs ΔEq.These relationships enable us to simplify the two energy partition schemes in Eqs.(3）and(4）and to obtain either Eqs.(7）and(8）or Eqs.(9）and(10）,where ΔEeand ΔExcare found to be the dominant contributor,respectively.

Our current results also shed new light on how to account for the origin of molecular stability for systems like water clusters.Same as other systems,1the relative stability of an isomer comes from the net contribution from all energetic effects involved.These effects,including electrostatic,steric,kinetic,exchange-correlation,and Fermionic quantum interactions,have different values for different isomers and they follow different trends in the conformation space.Some effects contribute positively to the molecular stability,while others do so negatively,canceling contributions from other interactions.One of the main results in this work is the finding that exchange-correlation interaction and steric repulsion are strongly correlated to the relative stability of water clusters,whereas for the electrostatic interaction,a less strong correlation has been observed.Even though in Eqs.(8）and(10）,the electrostatic interaction is dominant,Fig.1(b）shows that its correlation with molecular stability is weaker than the exchange-correlation interaction or steric repulsion.Using Eq.(7）or(9）,where ΔExcis dominant,much stronger correlation with relative molecular stability can be obtained.

4 Conclusions

To summarize,in this work,we employ quantum molecular dynamics to obtain a large number of distinct structures for the octamer water cluster and then perform energy partition studies using two approaches from density functional theory to identify working principles governing the relative molecular stability for these water clusters.We find that the exchange-correlation interaction and steric repulsion are two strong indicators of their relative conformation stability.We also identify strong correlations between energy components.It appears that a more stable structure possesses a smaller size and less steric repulsion,and at the meantime it has stronger electrostatic and exchange-correlation interactions.Two strong linear correlations using two different quantities are subsequently proposed to account for their relative stability,each with the correlation coefficient larger than 0.99.This work should shed new light to our fundamental understanding about the origin and nature of molecular stability for systems like water clusters as well as other similar molecular complexes formed through intermolecular interactions.

Finally,we mention in passing that our present approach is different from others scheme in performing energy decomposition analysis,such as the one by Morokuma,37where its focus is on the total interaction energy.In our case,we consider the total energy of the system instead.Also,what we have obtained in this work is only for the octamer.Are our conclusions applicable to other sizes of the water cluster as well?How sensitive are they to the choice of basis sets or density functionals?More interestingly,even though our approach is different from the other energy partition scheme(by Morokuma and others）in the literature,is there any correlation from the terms obtained these different approaches applied to the same systems?More systematic studies are in progress.These and other questions will be addressed elsewhere.

(2）Pophristic,V.;Goodman,L.Nature 2001,411,565.doi:10.1038/35079036

(3）Bickelhaupt,F.M.;Baerends,E.J.An gew.Chem.Int.Edit.2003,42,4183.

(4）Weinhold,F.Angew.Ch em.Int.E dit.2003,42,4188.

(5）Mo,Y.R.Nat.Chem.2010,2,666.doi:10.1038/nchem.721

(6）Mo,Y.;Gao,J.Accounts Chem.R es.2007,40,113.doi:10.1021/ar068073w

(7）Parr,R.G.;Yang,W.Density Functional Theory of Atoms Molecules;Oxford University Press:New York,1989.

(8）Geerlings,P.;De Proft,F.;Langenaeker,W.Chem.Rev.2003,103,1793.doi:10.1021/cr990029p

(9）Liu,S.B.Acta P hys.-Chim.Sin.2009,25,590.[劉述斌.物理化學學報,2009,25,590.]doi:10.3866/PKU.WHXB20090332

(10）Levy,M.;Perdew,J.P.Ph ys.Rev.A 1985,32,2010.doi:10.1103/PhysRevA.32.2010

(11）Liu,S.B.;Parr,R.G.P hys.R ev.A 1996,53,2211.doi:10.1103/PhysRevA.53.2211

(12）Liu,S.B.;Nagy,A.;Parr,R.G.Phys.Rev.A 1999,59,1131.doi:10.1103/PhysRevA.59.1131

(13）Liu,S.B.;Morrison,R.C.;Parr,R.G.J.Ch em.Phys.2006,125,174109.doi:10.1063/1.2378769

(14）Liu,S.B.J.Chem.Ph ys.2007,126,244103.doi:10.1063/1.2747247

(15）March,N.H.Phys.L ett.A 1986,113,476.doi:10.1016/0375-9601(86）90123-4

(16）Holas,A.;March,N.H.Phys.R ev.A 1991,44,5521.doi:10.1103/PhysRevA.44.5521

(17）von Weizs?cker,C.F.Z.Phys.1935,96,431.doi:10.1007/BF01337700

(18）Weisskopf,V.F.Science 1975,187,605.doi:10.1126/science.187.4177.605

(19）Liu,S.B.Phys.R ev.A 1996,54,4863.doi:10.1103/PhysRevA.54.4863

(20）Liu,S.B.;Parr,R.G.P hys.Rev.A 1997,55,1792.doi:10.1103/PhysRevA.55.1792

(21）Tsirelson,V.G.;Stash,A.I.;Liu,S.B.J.Chem.P hys.2010,133,114110.doi:10.1063/1.3492377

(22）Liu,S.B.J.Chem.P hys.2007,126,191107.doi:10.1063/1.2741244

(23）Esquivel,R.O.;Liu,S.B.;Angulo,J.C.;Dehesa,J.S.;Antolín,J.;Molina-Espíritu,M.J.Phys.Chem.A 2011,115,4406.doi:10.1021/jp1095272

(24）Liu,S.B.;Govind,N.J.Phys.Chem.A 2008,112,6690.doi:10.1021/jp800376a

(25）Liu,S.B.;Govind,N.;Pedersen,L.G.J.Chem.Phys.2008,129,094104.doi:10.1063/1.2976767

(26）Liu,S.B.;Hu,H.;Pedersen,L.G.J.Phys.Chem.A 2010,114,5913.doi:10.1021/jp101329f

(27）Ess,D.H.;Liu,S.B.;De Proft,F.J.Phys.Chem.A 2010,114,12952.doi:10.1021/jp108577g

(28）Huang,Y.;Zhong,A.G.;Yang,Q.;Liu,S.B.J.Chem.P hys.2011,134,084103.doi:10.1063/1.3555760

(29）Zhao,D.B.;Rong,C.Y.;Jenkins,S.;Kirk,S.R.;Yin,D.L.;Liu,S.B.Acta Phys.-Chim.Sin.2013,29,43.[趙東波,榮春英,蘇 曼,蘇 文,尹篤林,劉述斌.物理化學學報,2013,29,43.]doi:10.3866/PKU.WHXB201211121

(30）Tsirelson,V.G.;Stash,A.I.;Karasiev,V.V.;Liu,S.B.Comp.T heor.Chem.2013,1006,92.doi:10.1016/j.comptc.2012.11.015

(31）Torrent-Sucarrat,M.;Liu,S.B.;De Proft,F.J.Ph ys.Ch em.A 2009,113,3698.doi:10.1021/jp8096583

(32）Liu,S.B.J.Chem.Sci.2005,117,477;Zhong,A.G.;Rong,C.Y.;Liu,S.B.J.Phys.Chem.A 2007,111,3132.doi:10.1007/BF02708352

(33）Valiev,M.;Bylaska,E.J.;Govind,N.;Kowalski,K.;Straatsma,T.P.;Van Dam,H.J.J.;Wang,D.;Nieplocha,J.;Apra,E.;Windus,T.L.;de Jong,W.Comput.Phys.Commun.2010,181,1477.

(34）Maeda,S.;Ohno,K.J.P hys.Chem.A 2007,111,4527.doi:10.1021/jp070606a

(35）Zhao,Y.;Truhlar,D.G.T h eor.Ch em.A cc.2008,120,215.doi:10.1007/s00214-007-0310-x

(36）Frisch,M.J.;Trucks,G.W.;Schlegel,H.B.;et al.Gaussian 09,Revision C.01;Gaussian,Inc.:Wallingford,CT,2009.

(37）Kitaura,K.;Morokuma,K.Int.J.Quantum Chem.1976,10,325.