M. Griebel, J. Hamaekers, F. Heber A bond order dissection ANOVA approach for efficient ...

26 Pages·2015·1.59 MB·English

Wegelerstraße  •  Bonn • Germany
phone +- • fax +-
www.ins.uni-bonn.de
M. Griebel, J. Hamaekers, F. Heber
A bond order dissection ANOVA approach for
efficient electronic structure calculations
INS Preprint No. 1402
March 2014

A bond order dissection ANOVA approach for
efficient electronic structure calculations
MichaelGriebel,JanHamaekers,andFrederikHeber
Abstract In this article, we present a new decomposition approach for the effi-
cientapproximatecalculationoftheelectronicstructureproblemformolecules.Itis
basedonadimension-wisedecompositionofthespacetheunderlyingSchro¨dinger
equationlivesin,i.e.R3(M+N),whereM isthenumberofnucleiandN isthenum-
ber of electrons. This decomposition is similar to the ANOVA-approach (analysis
of variance) which is well-known in statistics. It represents the energy as a finite
sumofcontributionswhichdependonthepositionsofsinglenuclei,ofpairsofnu-
clei, of triples of nuclei, and so on. Under the assumption of locality of electronic
wave functions, the higher order terms in this expansion decay rapidly and may
thereforebeomitted.Furthermore,additionaltermsareeliminatedaccordingtothe
bondingstructureofthemolecule.Thisway,onlythecalculationoftheelectronic
structureoflocalparts,i.e.smallsubsystemsoftheoverallsystem,isnecessaryto
approximate the total ground state energy. To determine the required subsystems,
weemploymoleculargraphtheorycombinedwithmolecularbondingknowledge.
Inprinciple,thelocalelectronicsubproblemsmaybeapproximatelyevaluatedwith
whatevertechniqueisappropriate,e.g.HF,CC,CI,orDFT.Fromtheselocalener-
gies,thetotalenergyoftheoverallsystemisthenapproximatelyputtogetherina
telescopingsumlikefashion.Thus,ifthesizeofthelocalsubproblemsisindepen-
dentofthesizeoftheoverallmolecularsystem,linearscalingisdirectlyobtained.
Wediscussthedetailsofournewapproachandapplyittoboth,varioussmalltest
systemsandinterferonalphaasanexampleofalargebiomolecule.
MichaelGriebel
InstituteforNumericalSimulation,UniversityofBonn,Wegelerstr.6,53115Bonn,Germany,e-
mail:[email protected]
JanHamaekers
FraunhoferInstituteforAlgorithmsandScientificComputingSCAI,SchlossBirlinghoven,53754
SanktAugustin,Germany,e-mail:[email protected]
FrederikHeber
InstituteforNumericalSimulation,UniversityofBonn,Wegelerstr.6,53115Bonn,Germany,e-
mail:[email protected]
1
2 MichaelGriebel,JanHamaekers,andFrederikHeber
1 Introduction
The coupling of the micro- and the mesoscale of chemical reactions is currently
a field of intensive research. Where the microscale is the realm of quantum me-
chanical effects, the mesoscale is described by statistical mechanics and macro-
scopic thermodynamics. Nevertheless, there are additional strong influences onto
the mesoscale by effects from the microscale. Numerically, the microscale is usu-
allytreatedwithHartree-Fock(HF),ConfigurationInteraction(CI),CoupledClus-
ter (CC), or Density Functional Theory (DFT) methods which yield approximate
solutionstotheunderlyingquantum-mechanical(QM)Schro¨dingerequation(SE),
whereasthemesoscaleiscoveredbyclassicalmolecularmechanics(MM)methods
whichuseNewton’smechanicswithempiricallyfittedpotentialfunctions.
The ultimate goal would be a seamless coupling of quantum mechanical com-
putationswhereneededandclassicalmolecularmechanicssimulationswheresuffi-
cient.Suchapproachesaregenerallyreferredtoasmulti-scalemethods,anextensive
overviewisgivenin[28].AnystartingpointmustbethegeneralSchro¨dingerequa-
tionfortheelectronsandnucleiofthesystemunderconsideration.TheSchro¨dinger
equation however lives in 3(M+N) dimensions, where M denotes the number of
nucleiandNdenotesthenumberofelectrons.Thisrendersadirectnumericaltreat-
ment impossible due to the curse of dimension and one has to resort to model ap-
proximations.Asafirststep,intheBorn-Oppenheimermoleculardynamics(MD)
approach,thewavefunctionsofthenucleiandelectronsareseparated,thesubsys-
temofthenucleiistreatedclassicallywithNewton’smechanicsandtheremaining
3N-dimensionalelectronicSchro¨dingerequationisfurtherapproximatedbyoneof
the aforementioned methods. The potential needed for Newton’s mechanics is ob-
tained from the electronic solution by the Hellmann-Feynman theorem. This way,
QMandMMaregloballycoupled.However,aglobalelectronicQMsolutionis,at
leastforlargermolecules,stilltooexpensiveasconventionalmethodsscaleatbest
withO(M3)duetotheunderlyingproblemofmatrixdiagonalization.
Thus,generallinearscalingelectronicstructuremethodsareemployedtoover-
comethedimensionalityproblem.Asafirststep,forthelong-rangeCoulombinter-
action, the use of the fast multipole method [16] has resulted in O(MlogM) com-
plexity.Furthermore, acutoff radiussuch asfor theMP2 theory[5]was usedin a
Divide&Conquerapproachtotakeadvantageoftheexponentialdecaypropertiesof
electronic wave functions. Altogether, this resulted in linear scaling [14]. Another
commonmethodistheDensityMatrixMinimizationtechnique[27,8].There,the
densitymatrixisunconstrainedlyminimizedviaaconjugategradientscheme,using
idempotencyandnormalization.TheFockmatrixistheminimizedoutput,afteroff-
diagonalelementshavealsobeentruncatedatacutoffradius.Theelectroniclocal-
izationinnon-metallicsystemscanalsobeexploitedforplanewavebasissets[34].
Again,acutoffthenallowsforlinearscaling.Notehoweverthatthereisacrossover
pointuptowhichthestandardcubicscalingapproachesstillperformfasterdueto
smallerprefactorsintheircomputationalcomplexitycounts[15].
Inordertoreducetheconstantsandthustoshiftthiscrossoverpoint,onetriesto
somehowfurtherdecomposethefullglobalelectronicHamiltonianintolocalparts
AbondorderdissectionANOVAapproachforefficientelectronicstructurecalculations 3
andemployslocalQMthere.Letusbrieflysummarizethemostcommondecompo-
sitionapproachesinthefollowing.Oneofthefirst,theForce-MatchingMethodby
Ercolessi[13],triestoautomaticallygenerateempiricalpotentialsbyaleast-square
fittingoftheforcesofab-initiocalculationstogeneralmany-bodypotentialforms
suchasthatoftheEmbeddedAtomapproach[11]orthatofAbell-Tersoff[1,36].
Then, there is a range of methods which employ a decomposition1 directly in R3,
like e.g. the SIBFA (Sumof Interactions Between Fragments computed Ab initio)
procedure [17] and its generalization, the so-called Fragmentation Reconstruction
Method(FRM)[2].Furtherschemesaretheso-calledIMOMMansatzproposedby
Morokuma[30],theONIOMapproach[39]andthewell-knownFragmentMolecu-
larOrbital(FMO)method[25].Aschemeformodelingtheelectrostaticimpactofa
passiveMMenvironmentontheactiveQMsystemisdescribedin[26].Moreover,
in[3]and[33]aninterfaceregimebetweenQMandMMwith”link”atomsispro-
posedtoaccountforthecuttingofbonds.Similartechniquesareusedin[37,38].
Thecommonbasicideaofmostapproachesistouseatelescopingsumovertwore-
gionstodescribethetotalenergy,likee.g.Ω andΩ Ω ,wheretheenergyissplit
1 2 1
asEQM/MM =EMM+EQM EMM,wheretoourkno⊆wledgeallproceduresinvolve
Ω1 Ω2 − Ω2
stringent chemical knowledge to choose the regions (or cuts) as best as possible
buttostillkeeptheground-stateelectronicdensityintact.Anotherapproachdivides
time instead of space in order to generate a coupling between QM and MM. One
of these methods is learn-on-the-fly [10], which is similar to the Force-Matching
method,butisrunduringthecomputation:Atintermittenttimestepscertainclus-
tersofthesimulationdomainarelocallycomputedbyQM,andtheobtainedlocal
forcesareusedtocorrecttheMMcalculation.
Whilealloftheabovemethodshavepromisingfeatures,wefeelthattheygener-
allyeitherinvolvetoomanyadditionalparameters,unchemicallycutbondsinsep-
aratingactivefrompassiveregions,orevenworse,addunphysicalpseudo-atomsin
ordertocompensateforthedifferentenergyandtimescalesandtoavoidspill-out
effectsofelectronicdensityorenergy.Moreover,theyareplainlytoosimpleordo
not grasp the problem in its full complexity, since only a matching or interpola-
tionwithrespecttoenergyorforcesbetweentheQMandMMpartsoftheoverall
approachisemployed.
There is one more group of methods that build upon additivity models, well-
knowninchemistry,see[19]andreferencestherein.Thecentralideaistoconstruct
molecularpropertiesofasystembyaddingupthecorrespondingknownproperties
of its fragments. The principal hope is that a high-dimensional system such as a
complex molecule depends strongly only on few input variables. Rabitz et al [19]
describe a High-DimensionalModel Representation (HDMR) thatcan also be un-
derstoodasanANalysisOfVAriance(ANOVA),whichiswell-knownfromstatis-
tics. They address the problem of the estimation of the enthalpy of formation of a
broadrangeoforganicmoleculesbasedonexperimentaldata,buttheydonotassess
thepossibilitiesoftheansatzinthefieldofelectronicstructurecalculations.Deev
andCollins[12,9]usethisadditivitymodelansatzbycalculatingthetotalelectronic
1Ultimately,theaimwouldbeadecompositionofR3(M+N),thespacewherethefullSchro¨dinger
equationlivesin.
4 MichaelGriebel,JanHamaekers,andFrederikHeber
energyoffragmentsofasystemunderconsiderationtoobtainagoodapproximation
oftheenergyofthetotalmolecule.Theydogiveanalgorithmicdescription,how-
everwhichwefeelisnotfullyconsistentwiththemathematicalbasics,governedby
theANOVAorHDMRscheme.
In this article, we propose a more sophisticated algorithm. The additivity mod-
els place their hope on the same grounds as do many-body potential such as Ter-
soff’s [36], where the energy and the forces of an atom are assumed to depend
onitslocalcoordination.Here,foraproof-of-concept,weconcentrateoncovalent
bonding,henceoncharge-neutralmolecularsystemsandsubsystems.2 Wewilluse
thisknowledgeofcoordinationandbondsbetweennucleitodecomposethespace
R3(M+N) oftheunderlyingSchro¨dingerequationinadimension-wisefashion.This
decompositionissimilartotheANOVA-approach.Itrepresentstheenergyasafi-
nitesumofcontributionswhichdependonthepositionsofsinglenuclei,ofpairsof
nuclei,oftriplesofnuclei,andsoon.Undertheassumptionoflocalityofelectronic
wave functions, the higher order terms in this expansion decay rapidly and may
thereforebeomitted.Furthermore,additionaltermsareeliminatedaccordingtothe
bondingstructureofthemolecule.Thisway,onlythecalculationoftheelectronic
structure of local parts, i.e. small overlapping subsystems of the overall molecule
system,isnecessarytoapproximatethetotalgroundstateenergy.Todeterminethe
required subsystems, we employ molecular graph theory combined with molecu-
lar bonding knowledge. Here, modern graph algorithms are used to create proper
local subproblems as overlapping fragments of the overall molecular system. Fur-
thermore,hydrogenizationisusedtocloseshellsandsaturatebondsthathavebeen
cut.Wethusalsoexploitlocality,howevernotbyanexplicitcutoffradiusasmost
conventional methods do, but by implicitly using it in the inherent bond structure
ofthemolecularsystem.Inprinciple,thelocalelectronicsubproblemsmaybeap-
proximatelyevaluatedwithwhateverQMtechniqueisappropriate,e.g.HF,CI,CC,
or DFT. From these local energies, the total energy of the overall system is then
approximatelyputtogetherinatelescopingsumlikefashion.Thus,ifthesizeofthe
localsubproblemsisindependentofthesizeoftheoverallmolecularsystem,linear
scaling is directly obtained. The 3(M+N)-dimensional full global Hamiltonian is
broken down within the Born-Oppenheimer Approximation to O(M) components,
(k) (k)
thei-thofthemwithM degreesoffreedom,withanupperboundmax M con-
i i{ i }
trolledbyasingleparameterkwhichwenamethebondorderoftheapproximation.
Thisansatzspecificallycombinesthesmallerprefactorofthecubicscalingmethods
withagenerallinearscalingbehavior.Asthesizeofeachsubproblemdependson
thebondcoordinationoftheinvolvedatoms,wecoinedthemethodBOSSANOVA
(BondOrderdiSSectionANOVA).
Theremainderofthisarticleisorganizedasfollows:InSect.2webrieflysum-
marize the basics of the underlying Schro¨dinger equation. In Sect. 3 we describe
the ANOVA-like decomposition of the energy of the Schro¨dinger equation in the
contextofmoleculargraphtheory.InSect.4wegivenumericalresultsforabroad
rangeoforganicmolecules.WeendwithsomeconcludingremarksinSect.5.
2Notehoweverthatourapproachshouldworkequallywellalsointhenon-chargeneutralcase.
AbondorderdissectionANOVAapproachforefficientelectronicstructurecalculations 5
2 Schro¨dingerEquationintheBorn-Oppenheimer
Approximation
Let us consider a molecular system consisting of M nuclei and N electrons. Its
time-dependentstatefunctioncanbewritteningeneralas
Ψ =Ψ(R ,...,R ,r ,...,r ,t),
1 M 1 N
where R and r denote positions in three-dimensional space R3 associated to the
i j
ith nucleus and the jth electron, respectively. The variable t denotes the time-
dependency of the state function. The vector space (space of configurations), in
which the coordinates of the particles are given, is therefore of dimension 3(M+
N). In the following we will abbreviate (R ,...,R ) and (r ,...,r ) with the
1 M 1 N
shorter notation R and r, respectively. Also, we assume thatΨ is normalized to
(cid:82)
Ψ (R,r,t)Ψ(R,r,t)dRdr=1.
∗
Nucleiandelectronsarechargedparticles.Theelectrostaticpotential(Coulomb
potential) of a point charge is 1/r in atomic units, where r is the distance from
thepositionofthechargedparticle.Anelectronmovinginthispotentialpossesses
thepotentialenergyV(r)= 1/r.Neglectingspinandrelativisticinteractionsand
−
assumingthatnoexternalforcesactonthesystem,theHamiltonoperatorinposition
representationassociatedtothesystemofnucleiandelectronsisgivenasthesum
overtheoperatorsforthekineticenergyandtheCoulombpotentials,
H(N,M,Z ,m ,...,Z ,m ;R,r):=
1 1 M M
1 N N 1 N M Z M Z Z 1 M 1
∑∆ +∑ ∑ ∑ j +∑ k j ∑ ∆ ,
−2 rk r r − r R R R −2 m Rk
k=1 k<j|| k− j|| k=1j=1|| k− j|| k<j|| k− j|| k=1 k
(cid:124) (cid:123)(cid:122) (cid:125)
H (N,M,Z ,m ,...,Z ,m ;R,r)
e 1 1 M M
(1)
where we use a semicolon to distinguish between parameters (i.e. the number M
of atoms, the number N of electrons, the nuclei mass in atomic units m and the
j
atomic number Z ) and the degrees of freedom (i.e. the positions R and r). Here,
j
r r arethedistancesbetweenelectrons, r R aredistancesbetweenelec-
k j k j
(cid:107) − (cid:107) (cid:107) − (cid:107)
tronsandnucleiand R R aredistancesbetweennuclei.Wewillomitparam-
k j
(cid:107) − (cid:107)
eters from this list if they are clear from the context. This will later especially be
N,M,Z ,m ,...,Z ,m .
1 1 M M
Now,asystemofequationsfortheelectronicandforthenucleidegreesoffree-
domisusuallyderivedwiththeBorn-Oppenheimerapproximation.Tothisend,the
large difference in masses between electrons and atomic nuclei is exploited to de-
couplethemotionoftheelectronsfromthatofthenuclei.3 Then,oneassumesthat
3 TheratioofthevelocityvK ofanucleustothevelocityofanelectronve isingeneralsmaller
than10−2.
6 MichaelGriebel,JanHamaekers,andFrederikHeber
theelectronsadaptinstantaneouslytoachangeinthenuclearconfigurationandare
thusalwaysinthequantummechanicalgroundstatedenotedbyφ (R(t);r),which
0
is associated to the actual position of the nuclei R(t). Note that this allows us to
write H (R(t);r) instead of H (R(t),r) since the movement of the nuclei during
e e
theadaptationoftheelectronpositionsisnegligiblysmallinthesenseofclassical
dynamics. This justifies to setΨ(R,r,t) ΨBO(R,r,t):=∑∞ χ (R,t)φ (R;r),
≈ j=0 j j
which allows to separate the fast from the slow variables. We then obtain the fol-
lowingsetofequations:
(cid:26)(cid:90) (cid:27)
MkR¨k(t)=−∇Rk|φ0(Rm(ti)n;·)|=1 φ0∗(R(t);r)He(R(t);r)φ0(R(t);r)dr (2)
(cid:124) (cid:123)(cid:122) (cid:125)
=:VBO(R(t))
e
H (R(t);r)φ (R(t);r)=E (R(t))φ (R(t);r). (3)
e 0 0 0
In the end, after time discretization, we have to perform in each time step the
followingtasks:First,wehavetocomputeanapproximatesolutionoftheelectronic
Schro¨dinger equation in (3) for fixed positions R of the nuclei, then we have to
computefromitssolutiontheforcesonthenucleiandfinallywehavetocompute
thepositionsofthenucleiatthenexttimestepbye.g.aVerlettimestepforNewton’s
equationsofmotionofthenucleiin(2).Tothisend,weusetheHellmann-Feynman
Theoremtoobtaintheelectronicforces
(cid:90)
Fk(R(t))=−∇Rk φ0∗(R(t))He(R(t))φ0(R(t))dr
actingonthenuclei.VariantsofthisapproacharetheEhrenfestmoleculardynam-
ics and the Car-Parrinello method. For details of the derivation, see [18] and the
referencescitedtherein.
3 ANOVADecompositionScheme
So far, the Born-Oppenheimer molecular dynamics was employed to split the full
Schro¨dingerproblemintotwoparts,i.e.aclassicalNewton’sequationofmotionfor
thenuclei,and,ineachdiscretizedtimestep,theelectroniceigenproblem(3)which
may approximately be solved by e.g. the Hartree Fock, Configuration Interaction,
Coupled Cluster, or Density Functional method, see [35, 31]. However, such an
overall approach is only feasible for small molecules due to the high complexity
of any approximate solution method for the electronic problem. To overcome this
difficulty, the aforementioned coupling techniques and linear scaling methods had
AbondorderdissectionANOVAapproachforefficientelectronicstructurecalculations 7
been developed. They basically all exploit locality of the electronic wave function
inonewayoranothertoreducethecomplexityoftheelectronicproblem.45
Inthefollowing,wealsoresorttoacertainlocalityoftheelectronicwavefunc-
tion.Itisexpressedinthebondstructureofthemolecularsystem.Wedecompose
theoverallelectronicproblemintosmallsubproblemswhichthenmaybehandled
efficiently.Tothisend,weintroduceanANOVAdecompositionschemefortheen-
ergyofamolecularsystemintolocalpartsbymeansofthebondorderofthenuclei
inthesystem.
3.1 ANOVAExpansion
WewillnowdefinetheenergyfunctionforamolecularsystemanditsANOVAse-
ries expansion. To this end, we consider a molecular system which consists of N
electrons and M nuclei, each with coordinate vector R R3 and atomic number
i
Z N,i 1,...,M .Werestrictourselvestocharge-neu∈tralsystems,i.e.thenum-
i
∈ ∈{ }
berofelectronsN isequalto∑MZ forreasonsofsimplicity.Finally,weconsider
i i
the systems only in their electronic ground state in the framework of the Born-
Oppenheimer molecular dynamics. To this end, we separate the time-independent
electronicSchro¨dingerequationasin(3)anddefineatotalgroundstateenergyfunc-
tionEM :(N R3)M R.Itdependsontheparametersthatcompletelyidentifythe
× →
systemunderconsideration,namelythecoordinatesR andtheatomicnumberZ of
i i
eachnucleiwithfixedanduniquelabeli 1,...,M ,i.e.
∈{ }
EM((Z ,R ),...,(Z ,R )):=
1 1 M M
(cid:124) (cid:123)(cid:122) (cid:125) (cid:124) (cid:123)(cid:122) (cid:125)
=:X1 =:XM
(cid:90) M
(cid:0) (cid:1)
min φ0∗(R(t);r)He N=∑Zi,X1,...,XM)φ0(R(t);r dr, (4)
|φ0(R(t);·)|=1 i=1
wherewefurthersimplifythenotationbydefiningX :=(Z,R).I.e.X combines
i i i i
theatomicnumberandthecoordinatesofthenucleii.Notethat,duetothecharge-
neutrality condition N =∑MZ, the parameter N may now be eliminated from the
i i
parameterlistoftheHamiltonianH.
NowwewilldecomposethefunctionEM inamultivariatetelescopingsum,i.e.
inafiniteseriesexpansioninthenucleicparameters,inasimilarwayastheANOVA
decomposition6[21].ThisdecompositioninvolvesasplittingoftheM-dimensional
4Thisexcludesingeneralmetallicsystems,whoseelectronsmaybedelocalizedduetoavanishing
bandgap.
5Furthermore,thenotionofthelocalityofthewavefunctionisimportantasitleadstothegeneral
chemicalunderstandingofmoleculesfromthegeneralbondstructureuptonucleophilicsites.
6wTithhefAuNdeOpVenAdidnegcoomnlpyoosnititohneovfaraiaMbl-edsiminednisciaotnedalifnuun.cTtiohneffu:n[c0t,io1n]Ms→fusRatriesfaydsthfe=rec∑uur⊆re{1n,c..e.,Mre}lafu-
(cid:82)
tion f0/ =L{1,...,M}(f), fu=L{1,...,M}/u(f)−∑v⊂ufvwithLw(f)= [0,1]|w| f(x1,...,xM)dxw.Thus,
8 MichaelGriebel,JanHamaekers,andFrederikHeber
function into contributions which depend on the positions of single nuclei and as-
sociatedcharges,ofpairsofnucleiandassociatedcharges,oftriplesofnucleiand
charges, and so on. To this end, we consider the subset of the nuclei parameters
X describedbyasetoflabelsI withcardinality I =kandcallitthemolecu-
i i I
{ }∈ | |
larfragment associatedtoI withsizek.Notethatweheredonotneedtoconsider
theelectronicdegreesoffreedomr,asthesystemisassumedtobeingroundstate
and, hence, the electronic state functions are all fixed by the minimum condition
in(4).
First, we define the total electronic ground state energy of lower-dimensional
subsystems of the molecular system under consideration, described by the set of
indicesI= i ,...,i ,
1 k
{ }
(cid:90) k
(cid:0) (cid:1)
E{i1,...,ik}(X1,...,Xk):=|φm0|i=n1 φ0∗(r)He j∑=1Zij,Xi1,...,Xik)φ0(r dr. (5)
Notethatthisisinformverysimilarto(4).Inthenotationoftheelectronicground
statewavefunctionsφ ,thedependencyonR(t)wasdroppedasitisclearfromthe
0
context.
Then, the energy function EM is decomposed analogously to the ANOVA ap-
proachas
EM(X ,...,X )=F
1 M 0/
M
+∑F (X )
i1 {i1} {i1}
M
+ ∑ F (X )
i1<i2 {i1,i2} {i1,i2}
M
+ ∑ F (X )
i1<i2<i3 {i1,i2,i3} {i1,i2,i3}
+...
+F (X )
{i1,...,iM} {i1,...,iM}
=: ∑ F (X ),
U U
U 1,...,M
⊆{ }
whereX denotesthesetofvariables X andU 1,...,M .
U i i U
{ }∈ ⊆{ }
Here,eachtermF isdefinedasfollows:
{i1,...,ik}
f isdecomposedintoaconstant,asumofone-dimensionalfunctions,asumoftwo-dimensional
functions,andsoon.Theinvolvedfunctionsaregeneratedbyproperpartialintegrationandtele-
scopiccorrectionsaccordingtotherecurrencerelation.