Fast and accurate predictions of covalent bonds in chemical space

15 Pages·2017·2.19 MB·English

Fast and accurate predictions of covalent bonds in chemical space
,
K. Y. Samuel Chang, Stijn Fias, Raghunathan Ramakrishnan, and O. Anatole von Lilienfeld
Citation: The Journal of Chemical Physics 144, 174110 (2016); doi: 10.1063/1.4947217
View online: http://dx.doi.org/10.1063/1.4947217
View Table of Contents: http://aip.scitation.org/toc/jcp/144/17
Published by the American Institute of Physics
Articles you may be interested in
Communication: Understanding molecular representations in machine learning: The role of uniqueness and
target similarity
The Journal of Chemical Physics 145, 161102161102 (2016); 10.1063/1.4964627
Guiding ab initio calculations by alchemical derivatives
The Journal of Chemical Physics 144, 104103104103 (2016); 10.1063/1.4943372
Electronic spectra from TDDFT and machine learning in chemical space
The Journal of Chemical Physics 143, 084111084111 (2015); 10.1063/1.4928757
Accurate ab initio energy gradients in chemical compound space
The Journal of Chemical Physics 131, 164102164102 (2009); 10.1063/1.3249969
Tuning dissociation using isoelectronically doped graphene and hexagonal boron nitride: Water and other
small molecules
The Journal of Chemical Physics 144, 154706154706 (2016); 10.1063/1.4945783
THEJOURNALOFCHEMICALPHYSICS144,174110(2016)
Fast and accurate predictions of covalent bonds in chemical space
K.Y.SamuelChang,1StijnFias,2RaghunathanRamakrishnan,1
andO.AnatolevonLilienfeld1,2,a)
1DepartmentofChemistry,InstituteofPhysicalChemistryandNationalCenterforComputationalDesign
andDiscoveryofNovelMaterials(MARVEL),UniversityofBasel,4056Basel,Switzerland
2GeneralChemistry(ALGC),FreeUniversityBrussels(VUB),Pleinlaan2,1050Brussels,Belgium
(Received12January2016;accepted8April2016;publishedonline4May2016)
We assess the predictive accuracy of perturbation theory based estimates of changes in covalent
bondingduetolinear alchemicalinterpolationsamongmolecules. Wehaveinvestigated σ bonding
to hydrogen, as well as σ and π bonding between main-group elements, occurring in small sets of
iso-valence-electronicmoleculeswithelementsdrawnfromsecondtofourthrowsinthe p-blockof
theperiodictable.NumericalevidencesuggeststhatfirstorderTaylorexpansionsofcovalentbonding
potentialscanachievehighaccuracyif(i)thealchemicalinterpolationisvertical(fixedgeometry),(ii)
itinvolveselementsfromthethirdandfourthrowsoftheperiodictable,and(iii)anoptimalreference
geometry is used. This leads to near linear changes in the bonding potential, resulting in analytical
predictionswithchemicalaccuracy(∼1kcal/mol).Secondorderestimatesdeterioratetheprediction.
Ifinitialandfinalmoleculesdiffernotonlyincompositionbutalsoingeometry,allestimatesbecome
substantiallyworse,withsecondorderbeingslightlymoreaccuratethanfirstorder.Theindependent
particle approximation based second order perturbation theory performs poorly when compared to
the coupled perturbed or finite difference approach. Taylor series expansions up to fourth order
of the potential energy curve of highly symmetric systems indicate a finite radius of convergence,
as illustrated for the alchemical stretching of H+. Results are presented for (i) covalent bonds to
2
hydrogenin12moleculeswith8valenceelectrons(CH ,NH ,H O,HF,SiH ,PH ,H S,HCl,GeH ,
4 3 2 4 3 2 4
AsH , H Se, HBr); (ii) main-group single bonds in 9 molecules with 14 valence electrons (CH F,
3 2 3
CH Cl,CH Br,SiH F,SiH Cl,SiH Br,GeH F,GeH Cl,GeH Br);(iii)main-groupdoublebondsin
3 3 3 3 3 3 3 3
9moleculeswith12valenceelectrons(CH O,CH S,CH Se,SiH O,SiH S,SiH Se,GeH O,GeH S,
2 2 2 2 2 2 2 2
GeH Se);(iv)main-grouptriplebondsin9moleculeswith10valenceelectrons(HCN,HCP,HCAs,
2
HSiN, HSiP, HSiAs, HGeN, HGeP, HGeAs); and (v) H+ single bond with 1 electron.Published by
2
AIPPublishing.[http://dx.doi.org/10.1063/1.4947217]
I. INTRODUCTION coupling to the problem of efficiently estimating the PES of
newmoleculesusingTaylorseriesexpansionsinCCS,rather
Solving Schrödinger’s time independent equation for thanempiricism.
the unperturbed electronic ground-state within the Born- Thealchemicalcouplingapproachcanberelatedtogrand-
Oppenheimer approximation yields the potential energy canonical ensemble theory (Widom insertion)23–25 and has
surface(PES)ofanymoleculeasafunctionofnuclearcharges beenwellestablishedforempiricalforce-fieldbasedmolecular
{Z }(stoichiometry),nuclearpositions{R }(geometry),and dynamicsstudies.26–30UsingQM,alchemicalchangesareless
I I
number of electrons N (molecular charge).1,2 The PES plays commondespiteE.B.Wilson’searlyproposalofvariable Z,
afundamentalroleinchemistryandelsewherebecausemany backin1962.31WithinQM,anytwoiso-electronicmolecules
propertiescanbederivedfromit.Whileonecanstudyefficient in CCS can be coupled “alchemically” through interpolation
waysofpredictingthePESofsinglecompounds,3–5efficient of their external potentials. Here, we have investigated if
estimatesofPESofensemblesofmoleculesaremoreuseful alchemical predictions can be used to model the PES of
(and challenging) in the context of virtual compound design covalent bonds occurring in small closed-shell molecules
efforts.6–10Theseeffortstypicallyattempttosearchchemical madeupfrommaingroupelements.Wehavelimitedourselves
compound space (CCS) spanned by {{Z },{R },N}11,12 tocovalentbondstohydrogen,aswellassingle,double,and
I I
for novel materials with desirable properties. As such, triple bonds in molecules with no more than 14 valence
accurate yet efficient quantum mechanics (QM) based PES electrons.Wepresentanddiscussnumericalevidenceforthe
estimatesholdthekeyforsuccessfulrationalcompounddesign following set of observations: First order Taylor-expansions
applications.6,7,13–15 While many inexpensive semi-empirical of covalent bonding potentials can reach chemical accuracy
QMmethodsareavailable,forthisstudywerestrictourselves (∼1kcal/mol)ifthreeconditionsaremet.First,thealchemical
to first principles in the spirit of Refs. 12 and 16–22. More change has to be “vertical,” meaning that initial reference
specifically, we investigate the application of “alchemical” moleculeaswellasfinaltargetmoleculehavetopossessthe
same number of atoms located at the exact same positions.
a)[email protected] Second,allelementsinvolvedinthealchemicalchange,i.e.,all
0021-9606/2016/144(17)/174110/14/$30.00 144,174110-1 PublishedbyAIPPublishing.
174110-2 Changetal. J.Chem.Phys.144,174110(2016)
{Z } destined to vary, have to occur late in the periodic Given a pair of iso-electronic reference/target systems,
I
table. Third, the off-set of the bond potential, determined describedby {{ZR},{RR},N} and {{ZT},{RT},N},respec-
I I I I
by reference distance, is optimized. Second order Taylor- tively,onecancouplethetwosystemssuchthatcertainZRand
I
expansionbasedpredictionsarelessaccuratethanfirstorder ZT are paired. Note that ZR or ZT can be scaled down to/up
I I I
predictionsiftheseconditionsaremet.Ifreferenceandtarget fromzeroifthenumberofatomsinonemoleculeissmaller.
molecules have different geometries, the predictive power Underiso-electronicconditions,theλ-dependenttermsinthe
of the first order Taylor expansion substantially deteriorates, coupling Hamiltonian (Eq. (2)) are the electron-nucleus and
whilesecondorderestimatesbasedoncoupledperturbed(CP) nucleus-nucleusinteractionoperators,
Kohn-Sham equations offer some improvement, however,
without reaching chemical accuracy. Second order estimates NI ( (1−λ)ZR λZT )
v (r)= − I − I ,
basedontheindependentparticleapproximation(IPA)result λ |r−RR| |r−RT|
in Taylor expansion estimates that are even worse than I I I (3)
NI ((1−λ)ZRZR λZTZT )
first order estimates. For highly alchemical changes with V = I J + I J .
symmetric endpoints, such as the dissociation of H+, a finite λ |RR−RR| |RT−RT|
2 I<J I J I J
radiusofconvergenceisfound.
Since different pairing schemes result in different v (r) and
In Sec. II we briefly summarize the framework of λ
V ,itisobviousthatthealchemicalperturbationisalignment
alchemical derivatives within Hartree-Fock and density λ
dependent. To investigate the behaviour of higher order
functional theory (DFT) as well as our notations. Numerical
correctionsandtheeffectsofvaryinggeometry/stoichiometry,
estimates of covalent bond stretching energies of small
we neglect all relaxation effects for vertical iso-valence-
moleculesarepresentedanddiscussedinSec.III:Extending
electronicchanges(seeSec.IIG).
previous work on alchemical perturbation,19,20,32 we discuss
alchemical energy derivatives with respect to vertical
transmutation, interpolating only the identity of the atoms
B. Firstorderderivative
whilekeepingthegeometryfixed.Estimatesforsingle,double,
andtriplebondsareincludedasanapplication.Wealsoreport The first order derivative of the energy with respect
numericalresultsforalchemicalstretchingofchemicalbonds to an alchemical interpolation parameter connecting any
using non-vertical transmutations. Conclusions are drawn in two iso-electronic molecules can be computed according to
Sec.IV. the Hellmann-Feynman theorem,33 as has been shown for
molecularHOMOeigenvalues,20

II. METHOD ∂λEλ =⟨∂λHλ⟩λ = dr ρλ(r)∂λvλ(r)+∂λVλ, (4)
A. TaylorexpansioninCCS where ρ (r) denotes the electron density, dependent on λ.
λ
ATaylorexpansionofthepotentialenergyinCCScanbe At λ =0 we have ρλ(r)= ρR(r), which is independent of
constructedwiththeexclusiveknowledgeacquiredbysolving the target system. As such, the first order derivative can
Schrödinger’s equation for some reference molecule, with be calculated with a single reference density and without
Hamiltonian H , additionalself-consistentfield(SCF)calculationforanytarget
R
system.Inseveralcircumstances,Taylorexpansionestimates
E(∆λ)= ER+∆λ∂λEλ(cid:12)(cid:12)(cid:12)λ=0+ ∆2λ2∂λ2Eλ(cid:12)(cid:12)(cid:12)λ=0+···. (1) uacscinugracfiyrstfoorrdtheer raalpchidempirceadlicdtieornivaotfivpesrophearvteiesshtohwronugghoooudt
Derivativesofthetotalpotentialenergycanbeobtainedby CCS.12,21,32,34,35 In general, however, first order derivatives
couplingareferenceHamiltoniantosometargetHamiltonian, might not be sufficient. Taking higher order derivatives
H ,suchthat H transforms H into H , into account might offer higher accuracy, assuming Eq. (1)
T λ R T
convergesrapidly.
H =(1−λ)H +λH , (2)
λ R T
as the coupling parameter λ goes from 0 to 1. And
consequently, ∂mE =∂m⟨H ⟩, with ∂ H = H −H = H′ C. Secondorderderivative
λ λ λ λ λ λ T R
being the alchemical perturbing Hamiltonian. If these Differentiation of Eq. (4), based on linear interpolated
derivatives can be computed, E can be estimated according HamiltonianinEq.(2),yields
T
to Eq. (1) by setting ∆λ =1. Note that we couple reference 
(cid:0) (cid:1)(cid:0) (cid:1)
and target systems in a linear and global fashion. This is an ∂2E = dr ∂ ρ(r) ∂ v (r) , (5)
λ λ λ λ λ
arbitrary choice; non-linear and local interpolation functions
could have been chosen just as well. In fact, in Ref. 20, requiring the density response due to the alchemical
an empirical quadratic interpolation function is found to perturbation. Again, at λ =0 this amounts to the density
yield superior results for first order predictions of highest responseofthereferencesystem.EvaluationofEq.(5)implies
occupiedmolecularorbital(HOMO)eigenvalues.Inthisstudy a differing density response for each target system. We have
of alchemical changes of covalent bonding, we begin with considered three approximations to ∂ ρ including second
λ
linearandglobalinterpolations;futureworkmightdealwith order perturbation theory with IPA,36 CP approaches,37,38 as
alternativefunctions. wellasfinitedifference(FD)approximation.NotethatEq.(5)
174110-3 Changetal. J.Chem.Phys.144,174110(2016)

can be rewritten as ∂2E = drdr′(∂ v(r))(∂ v(r′)) δ2E , E. Predictingchangesincovalentbonds
λ λ λ δv(r)δv(r′)
where δv(δr)2δEv(r′) = χ(r,r′)isthestaticlinearresponsefunction Forthestudyofcovalentbonds,wefocusonthechanges
orsusceptibility,wellestablishedwithinconceptualDFT.39–44 inbindingpotentialduetoalchemicalcoupling.Weconsider
Perturbationtheoryprovideswaystoestimate∂λρλ(r).45 thedifferenceintotalpotentialenergybetweentwocovalently
WithinIPA,36,46,47thestaticdensityresponseforaclose-shell bound fragments at two arbitrary interatomic distances
systemisapproximatedby d andd ,
0

∂ ρ (r)≈−4φ(r)φ (r) dr′φi(r′)φa(r′)∂ v (r′), (6) ∆E(d,d0)= E(d)−E(d0). (9)
λ λ ia i a εa−εi λ λ If, for example, d0 is large and d is the geometry minimum,
∆E becomes the bond dissociation energy. We are interested
where{φi,εi}denotetheithoccupiedmolecularorbital(MO) in changes of ∆E(d,d ) as a function of d due to alchemical
anditseigenvalue,while {φ ,ε } denotethe athunoccupied 0
a a changesforafixedd .Morespecifically,wecoupleareference
0
counterparts. IPA neglects the influence of the alchemical
totargetsystemviathecorrespondingHamiltoniansyielding
perturbation on the Hartree and exchange-correlation (xc)
expectationvaluesasafunctionofλ,
potentials.38,42 Note that Eq. (6) becomes numerically exact
for1-electronsystemwithconvergedbasissetwithinHartree- ∆Eλ(d,d0)=Eλ(d)−Eλ(d0)
FockapproximationbecauseoftheabsenceofCoulomband =⟨H (d)+λ(H (d)−H (d))⟩
R T R
xcinteractionbetweenelectrons. −⟨H (d )+λ(H (d )−H (d ))⟩. (10)
R 0 T 0 R 0
Recently, Yang, Cohen, De Proft, and Geerlings derived
an expression of the density response that also includes the As λ goes from 0 to 1, the two components in Eq. (9)
dependenceofCoulombandxcpotential,37theCPapproach,44 changefromreference(ER(d),ER(d0))totarget(ET(d),ET(d0))
compound.ThetruncatedTaylorexpansionbasedestimateof

∂ ρ (r)=−4 φ(r)φ (r) thetargetcompound’spotentialisthenobtainedvia
λ λ i a
ij ab  ∆E (d,d )≈∆E(m)(d,d )
T 0 T 0
×(M−1)ia,jb dr′φj(r′)φb(r′)∂λvλ(r′), (7) m 1
=∆E (d,d )+ ∂k∆E (d,d ), (11)
R 0 k! λ λ 0
wherethematrixelementsofMare k=1
where the superscript m stands for Taylor expansion with m
MXJiiaa,,jjbb === (εadd−rrddεrri)′′φδφii(j(δrr)a)φφba+((|rrr4))−φJφijra((,r′r|j′b′))φφ+b(4(rrX′′))i(a,,jb,δ2Exc ).(8) tcwoλehr,ordamarennkrsgd,s,eutaushn.spelSaetdisonfesucpmenoectn∆=htdieEoe4rTnnwfcioiosysrfetohtbhfenoeodpnt0srdeotw-drple.eietlnIclrngthbytiethnohogfidmsionfsifotttteuHrerddev+yse,f,troa,twrnitchdteaheleuisnapurvlebcetsoshstcteoimrmgifpati=tthceTaid2sl,
 ia,jb i a j b δρ(r)δρ(r′) 2
forallothermolecules.Forafixedd ,thefirstorderestimate
0
In the limit of J → 0 and X → 0, Eqs. (6) and (7) iscalculatedaccordingto
ia,jb ia,jb (cid:0) (cid:1)
are equivalent. In our implementation, the evaluation of the ∆ET(1)(d)= E(cid:0)R(d)+∂λEλ(d)|λ=0 (cid:1)
CP second order derivative has a computational complexity
similar to a molecular frequency calculation. IPA incurs − ER(d0)+∂λEλ(d0)|λ=0
negligiblecomputationaloverhead. =E (d)+ dr ρ (d,r)[v (d,r)−v (d,r)]
R R T R
Alternatively, one can also introduce an explicit small
perturbation and converge the new density at ∆λ ≪1. The +[VT(d)−VR(d)]−ER(d0)
density response can then be estimated via FD, ∂λρ(r) − dr ρ (d ,r)[v (d ,r)−v (d ,r)]
≈ ρ∆λ(r)−ρR(r). In practice, instead of starting the SCF for R 0 T 0 R 0
∆λ
the perturbed system from atom based initial guesses, we −[V (d )−V (d )], (12)
T 0 R 0
restart with ρ (r) resulting in convergence within few SCF
R andthesecondorderestimatecorrespondingly,
steps.
(cid:0) (cid:1)
∆ET(2)(d)= ER(d)+∂λEλ(d)|λ=0+ 12∂λ2Eλ(d)|λ=0
D. Higherorderderivatives (cid:0) (cid:1)
1
Møller-Plesset (MP) perturbation theory48,49 is used to − ER(d0)+∂λEλ(d0)|λ=0+ 2∂λ2Eλ(d0)|λ=0 , (13)
estimate correlation energy corrections based on converged using above discussed expressions for the CP, IPA, and FD
Hartree-Fock results. The derivation of higher order definitionsofthesecondorderderivative.
corrections in MP theory is equivalent to the mth order Sincedandd inEq.(9)arearbitrary,onecanobtainthe
0
alchemical derivative. Here, instead of the two-particle bindingcurvebyscanning d foranyfixed d .Thepredictive
0
operator for electron-electron interaction as perturbation in power, however, happens to depend on d . For this reason,
0
MPtheory,thealchemicalperturbationoperatorH −H can we optimize d such that the integrated error in dissociation
T R 0
beused.WithinIPA,theMPformulacanbedirectlyapplied regionisminimal.AsshowninFig.3anddiscussedbelow,an
toobtainanymthorderderivative. empiricallinearrelationshipexistsbetweenequilibriumbond
174110-4 Changetal. J.Chem.Phys.144,174110(2016)
lengthoftargetmoleculedT andd , used, resulting in a constant number of valence electrons.12
eq opt
For example, one can couple carbon to silicon using just
dopt≈0.76deTq+0.97Å. (14) four valence electrons. Non-local PPs are widely used to
mimic the presence of core electrons in atoms50 and are
d isdeterminedaccordingtoEq.(14)forallverticalchanges.
0
amenabletothetuningofawiderangeofpropertiesincluding
If dT is not known, it can easily be estimated with semi-
eq dispersion forces, band-gap, or vibrational frequencies.51–53
empiricalquantumchemistrymethods.
Thenon-localexternalpotentialv (r)inEq.(3)thenbecomes
Fornon-verticalchanges,weusetheequilibriumdistance λ
d in the reference molecule as starting reference geometry,
eq
a=n0d.wEeq.a(l1so2)fisximdp0litfioetsheinstahmisecvaasleu,e, resulting in ∆E(m)(deq) vλ(r,r′)=NI ((1−λ)vIR(r,r′)+λvIT(r,r′)), (17)
 I
∆E(1)(d)= dr ρ (d ,r)[v (d,r)−v (d ,r)]
R eq T T eq
where vR and vT are PPs for ZR and ZT, respectively.
I I I I
+[V (d)−V (d )], (15) Note that v (r,r′) in Eq. (17) and v (r) in Eq. (3) result
T T eq λ λ
in different coupling Hamiltonians, and therefore different
andsodoesEq.(13)forsecondorderestimates.
λ-dependenciesoftheenergyanditsderivatives.
All results have been obtained within the Born-
Oppenheimer approximation, where nuclei are clamped;
F. Errormeasures
nuclear repulsion V is decoupled from the electronic
λ
For bond lengths, we quantify the predictive power wavefunction and is added as a geometry- and λ-dependent
of the Taylor expansions by evaluating the deviation of constant to the electronic energy. Nuclear-nuclear repulsion
prediction from the DFT bond length ∆d = d(m)−d , energy is computed automatically by most QM codes.
eq eq eq
where d(m) stands for the predicted equilibrium distance However,itmustberemovedandrecomputedindependently
eq
of ∆E(m). We calculate the deviation of the predicted for V according to Eq. (3) to avoid self-repulsion between
λ
energy at d(m) from the DFT energy at the DFT minimum, transmutating atoms. Throughout the present study, standard
eq
∆E =∆E(m)(d(m))−∆E(d ). The deviation in harmonic atomic and plane-wave basis functions, linearly interpolated
eq eq eq
vibrationfrequency,∆ω =ω(m)−ω,ofthebondstretchingis PPs, as well as the PBE xc potential54 within Kohn-Sham
alsoincludedinordertoquantifytheaccuracyofthestiffness DFTareused.Thescanningof0.5Å ≤ d ≤ 3.0Åiscarried
of the predicted binding potential. The vibration frequency out with increments ∆d =0.1 Å. For each prediction order
is computed from the curvature of cubic spline interpolated m, ∆E(m)(d) are interpolated with cubic splines, from which
binding potential, ω = 1 keq, where k =∂2∆E(d ) and the stiffness ∂d2∆E(m)(de(mq))= keq is computed. All density
2π µ eq d eq volumetricdataareprintedintoGaussianCUBEfiles,from
µ is the reduced mass. Finally, we measure the integrated
whichintegrateddensityslicesarecalculated.
mean absolute error (MAE) for the dissociative tail, defined
as

1 dmax 1. Detailsforverticaliso-valence-electronicchanges
MAE= dx|∆E(m)(x)−∆E(x)|, (16)
|dmax−de(mq)| de(mq) Numerical results for vertical iso-valence-electronic
alchemical changes (discussed in Secs. III A and III B)
for vertical iso-valence-electronic changes. Note that while havebeenobtainedwithCPMD,55aplanewavebasiswith100
in principle one would like dmax→ ∞, dmax has been set to Ry cutoff, and Goedecker PPs.56–58 The periodic supercell
correspondroughlytotheinflectionpoint,duetotheissuesofa size is 20×15×15 Å3, and one heavy atom is fixed at
singledeterminantmethod,suchasDFT,fordescribingcova-
(7.5 Å, 7.5 Å, 7.5 Å) while the stretching atom shifts along
lentbond-dissociation.Thisshortcomingisalsoevidentfrom +x-axis.Foreachgeometry,heavyatomsaremutatedtoother
comparisonofDFTtocoupled-clustersingledoublespertur-
elements in the same column of the periodic table while
bative triples (CCSD(T)) curves shown in Fig. 3. Note that
all H’s are fixed at the same location as in the reference
this aspect is irrelevant for alchemical predictions: If a more
compound.
reliablereferencemethodhadbeenused,theerrorintegration
Since Eq. (17) is a non-local operator, Eqs. (4) and
couldhaveeasilybeenexpandedtoincludetheentiredissocia-
(5) need to be converted to wavefunction expressions. The
tivetail.Thesefourquantitiesprovideanumericalindication
first order derivative for the Hamiltonian H is evaluated
R→T
of how good a prediction is. For a perfect prediction, one
using RESTART files in which the reference compound’s
would expect (∆Eeq,∆deq,∆ω,MAE)=(0,0,0,0). Note that density and wavefunctions have been stored: ∂ E =⟨∂ H⟩
λ λ R
wecomparethepredictionstoDFT.Thisisanarbitrarychoice, = E [ρ ]−E [ρ ]. And the second order derivative is eval-
T R R R
anyotherQMmethodcouldhavebeenappliedjustaswell. uatedcorrespondinglyrelyingonFD,∂2E ≈ ⟨∂λH⟩∆λ−⟨∂λH⟩R,
λ ∆λ
with ∆λ =0.05. Wavefunctions of reference compound are
used for ∆E(1), while ∆E(2) is evaluated by FD with linearly
G. Computationaldetails FD
interpolatedPPparameters.
Alchemical interpolations of molecules containing CCSD(T) results obtained for HCl and HBr have been
elementsfromdifferentrowsintheperiodictablecanstillbe computed using Gaussian0959 in aug-cc-pVTZ60 basis and
iso-electronic if effective core or pseudopotentials (PPs) are defaultinputparameters.
174110-5 Changetal. J.Chem.Phys.144,174110(2016)
2. Detailsfornon-verticaliso-electronicchanges
Numericalresultsfornon-verticaliso-electronicalchem-
ical changes have been obtained using atom centered basis-
sets. Restricted open-shell Hartree-Fock calculations have
beencarriedoutusingCartesianaug-cc-pVTZbasisset60for
H+ (discussed in Sec. III D 1). Eq. (6) and higher order
2
derivatives are evaluated analytically by Gaussian expansion
ofMOs.ReferencegeometryisfirstrelaxedbyGaussian0959
and the converged MO coefficients are extracted to evaluate
orbital integrals. NWChem61 is used to scan ∆E as a function
of λ in Fig. 5(d) along alchemical path with discretization
∆λ =0.01. It is done by reassigning nuclear charges in the
system.
Non-verticalalchemicalchangesin10-electronmolecules
(discussed in Sec. III D 2) have been calculated using the
uncontracted Cartesian Def2TZVP basis set.62 Uncontracted
neon basis is used for all second row atoms. Additional
hydrogen basis functions are placed along the stretching
pathway, from d =0.5 Å to d =3.0 Å in increments
∆d =0.1 Å. All systems with integer nuclear charges
have been calculated using Gaussian0959 while systems
with fractional nuclear charges have been calculated using
NWChem61withdiscretization∆λ =0.01.Foreach0 ≤ λ ≤ 1,
theatomicdensityforSCFinitialguessiteratesthrough{C,N,
O,F,Ne}toensureconvergence.InallGaussianandNWChem
calculations,weusedCartesian/realsphericalharmonicbasis
functions.
III. RESULTSANDDISCUSSIONS
A. Verticaliso-valence-electronicchangesofX-H
1. Predictedpotentials
Using Taylor expansions binding potentials have been
estimated for covalent bonds involving hydrogen (X-H)
for the following 12 molecules with 8 valence electrons:
CH , NH , H O, HF (second period); SiH , PH , H S,
4 3 2 4 3 2
HCl (third period); and GeH , AsH , H Se, HBr (fourth
4 3 2
period). Numerical results for vertical first (red) and second
(blue)ordertruncatedTaylorseriesestimates,calculatedwith
Eqs.(12)and(13),featureinFig.1.Theymeasurethechange
in X-H binding energy as one goes from reference to target FIG. 1. ∆E is shown as a function of d in Eq. (11). White background
compound. panels: true (black circles), first (red squares), and second (blue triangles)
order predictions of changes in the covalent bond potential of hydrogen
We first note that the entire potential is reproduced
duetoverticalalchemicalinterpolations.Graybackgroundpanels:thetrue
in semi-quantitative fashion for all combinations of refer-
potentialsofthereferencecompoundsemployedforfirstandsecondorder
ence/target molecules. The precise predictive power strongly predictions.
depends on the choice of reference/target molecule pair, on
the choice of d , and on the expansion going up to first or
0
second order. First order estimates among molecules with HBr as a reference, the second order prediction is more
elements from the third or fourth row are very accurate (see accuratethanfirstorder.Fortheinverseprediction(i.e.,HBr
Fig. 1, bottom and mid-row in mid- and bottom panels, from HF), however, first order is more accurate than second
respectively).Bycontrast,predicting,orstartingwith,second order.
row elements consistently yields worse results. Inclusion TheperformanceoftruncatedTaylorseriesdramatically
of second order corrections does not necessarily lead to variesdependingonthechoiceofthed value.Thetoppanel
0
improved performance. Second order truncated Taylor series in Fig. 2 illustrates this for ∆E(2 Å,d ) for HF→ HBr as a
0
estimatesonlyyieldmoreaccuratepredictionsthanfirstorder function of λ, once with d =0.94 Å= d the equilibrium
0 eq
when the reference molecule contains heavier elements than bondlengthofHF—andoncewithd =1.57Å= d ,avalue
0 opt
the target molecule. For example, if we predict HF using ford whichhappenstolinearize∆Einλ.Whilethecoupling
0
174110-6 Changetal. J.Chem.Phys.144,174110(2016)
FIG.2. AlchemicalcouplingofHF(λ=0)toHBr(λ=1).TOPpanel:E (LEFT)and∆E (RIGHT)wheredeq=0.94Å(red)denotestheequilibriumbond
lengthofreferencemoleculeHF,anddopt=1.57Å(blue)linearizes∆E.BOTTOMpanel:integratedvalenceelectrondensitydifferenceslicesbetweenH-X
at d=2 Å and at dopt (RIGHT) and deq (LEFT), respectively, ∆Pλ(x)= dydz[ρλ(r,d)−ρλ(r,d0)]. Dependence on λ is shown for the same vertical
interpolations,leftandrightcorrespondingtothenon-linear(red)andlinearized(blue)∆Ecurvesinright-handTOPpanel.Densitychangesatheavyatomand
hydrogenpositionsarehighlightedasreddashedanddotted/dashed-dottedlines,respectively.
path of total energies is hardly distinguishable for E(2 Å), (λ → 0), a weak convexity is noticeable in ∆E (blue line),
E(1.57 Å), and E(0.94 Å), ∆E is strongly dependent on the despite the overall concavity of the path. The presence
choice of d . By choosing d =1.57 Å, ∆E(2 Å,1.57 Å) in of inflection points will always lead to a deterioration of
0 0
Eq.(10)becomesnearlylinear,whileplotting∆E(2Å,0.94Å) secondorderpredictions,implyingamoreaccuratefirstorder
reveals substantial curving. This is why, when choosing the estimate.Ontheotherside(λ → 1),nosuchinflectionpoint
rightd ,firstorderpredictionsof∆E canbeverypredictive. exists and the second order term results in the expected
0
The top panel in Fig. 2 also explains why second order improvementoftheprediction.Fortheothercompoundpairs
estimatescanbeworsethanfirstorder,andwhythischanges shown in Fig. 1, similar observations can be made for the
for the reverse coupling: On the side of the lighter element attractive part of the bonding potential (see also ∆E(2 Å,d )
0
FIG.3. LEFT:Scatterplotofoptimizedreferencebondlengthdoptversusequilibriumbondlengthoftargetmolecule,denotedbydeTq.Linearregressiongives
dopt=0.76deq+0.97ÅwithMAE=0.11ÅandRMSE=0.15Å.Firstorder(redemptysquares)andsecondorder(blueemptytriangles)doptofcovalentX-H
bondstretching,aswellasfirstorder(redfilledsquares)andsecondorder(bluefilledtriangles)doptofX-Y,X=Y,X#Ystretchingareshown,where-,=,#
standforsingle,double,andtriplebonds.Someofthealchemicalpathsarehighlightedbyblackarrows.AllnumbersaregiveninTablesSIandSII.72RIGHT:
Alchemicalpredictionscanbemoreaccuratethanapproximateddensityfunctionals.CovalentbindingpotentialsobtainedfromalchemicalPBE0estimate(red
squares)andordinaryPBE(greendiamonds)forHBr(fullsymbols)andHCl(emptysymbols).ThealchemicalestimatecorrespondstoEq.(12)usingPBE0
densityofHBr(HCl)inordertopredictHCl(HBr).Forcomparison,correspondingCCSD(T)resultsareshownaswell(dashed).
174110-7 Changetal. J.Chem.Phys.144,174110(2016)
curves in supplementary material72). We have found that the are more accurate than upward. For example, predicting
inflection point near λ =0 occurs only when the reference HBr using HF as a reference, a better estimate is obtained
molecule has the lighter element. Conversely, no inflection (error=+25.7 kcal/mol) than for predicting HF using HBr
point has been observed for atoms transmutating upward the as a reference (error=+61.1 kcal/mol). Correspondingly,
column.Webelievethatthisbehaviourisduetothespecifics predicting HCl using HBr has an error=+5.4 kcal/mol,
oftheemployedPPs.Futurestudieswillshowwhythisisthe while the prediction of HBr using HCl has only an error of
case,andifsimilartrendsholdforotherPPs. +3.6kcal/mol.Similarobservationsholdforbondlengthsand
force constants. The asymmetry is also illustrated in Fig. 2.
∆E(d,d )isnotnecessarilysymmetricwithrespecttoλ =0.5
2. Integratederror 0
for a given choice of (d,d ). Consequently, truncated Taylor
0
Prediction errors for energy minima, equilibrium bond series based predictions from either ends will not be equally
lengths, force constants, and integrated error in dissociation accurate.
region(calculatedasdescribedinSec.IIF)havebeenobtained
for all predictions in Fig. 1 and are reported in Table SI.72
The results lend quantitative support for the observations 4. Chemicalaccuracy
articulated above. In particular, these results suggest that
We have seen that very accurate, yet inexpensive, first
chemical accuracy can be obtained when using first order
orderalchemicalestimatescanbemadeforverticalalchemical
Taylor series based estimates among compounds containing
changes between third and fourth row elements according to
thirdandfourthrowelements.Secondorderbasedpredictions
Eq. (4)—once the density is converged for a given reference
arealwaysworseexceptwhenamoleculewithheavierelement
molecule. Then, an interesting question is if the alchemical
is used as a reference to predict a molecule with lighter one,
accuracy is on the same order of magnitude as common
forexample,HBr→ HF.
approximations made when solving Schrödinger’s equation.
The best prediction performance is found for first order
WehaveinvestigatedthispointforalchemicalcouplingofHBr
based estimates using reference molecules containing third
andHClusinghybridandgeneralizedgradientapproximated
row elements (nR=3) in order to predict target molecules DFT. When using PBE063 as the method for the reference
madeupoffourthrowelements(n =4).Theoverallaverage
T compound,wefindthefirstorderbasedalchemicalpredictions
deviation from reference bonding potential energies and
accordingtoEq.(12)tobeinbetteragreementwiththePBE0
integratederroris∼2.5kcal/mol.Correspondingpredictions
resultsforthetargetcompoundthantruegeneralizedgradient
of equilibrium distances deviate at most by 0.03 Å, and the basedapproximationPBE.64Fig.3illustratesthispointforthe
vibrationfrequenciesdeviatenotmorethan32cm−1.Second
covalent bindingpotentials ofHCl andHBr calculatedusing
orderestimatesforthesamethirdandfourthrowcombinations
PBE0,PBE0basedverticalfirstorderalchemicalpredictions,
give slightly worse results. The worst predictions are found
andPBE.Forallinteratomicdistancesinthedissociativetail,
if the coupled molecules skip a row, i.e., involve elements
thealchemicalprediction(squares)isclosertoPBE0(circles)
from second and fourth row—for first as well as second
thanPBE(diamonds).Fortherepulsivepartofthepotential,
order truncated Taylor expansions. This is not surprising as
thealchemicalpredictionissubstantiallybetterthanPBEfor
the central atom’s electron density must accommodate the
HBr and slightly worse than PBE for HCl. For comparison,
most severe contractions/expansions for such interpolations.
we also included CCSD(T) results. These results amount
Moreover, second order overcorrections can also be found
to numerical evidence that the predictive power of vertical
(Table SI)72 whenever molecules containing fourth row
alchemical predictions can exceed the accuracy of common
elements are used to predict molecules containing third row
DFT approximations for third or fourth row elements—if
elements: both, predicted energy minimum and equilibrium a sufficiently accurate electron density is provided for the
bondlength,shownegativedeviations.
referencecompound.
3. Alchemicalpredictionsdonotcommute
B. Verticaliso-valence-electronicchangesinvolving
We note the asymmetry in the predictive power of single,double,andtriplebonds
first order based predictions which is due to the lack
1. Predictedpotentials
of commutation: In general ∂λE|λ=0(cid:44)−∂λE|λ=1, except if
reference and target Hamiltonian happen to differ only by Having discussed covalent bonds involving hydrogen,
translation, rotation, or parity (enantiomers, i.e., without we now turn to single (XH -Y), double (XH =Y), and
3 2
accounting for parity violation). Within our restricted triple (HX#Y) bonds among p-block elements. Since third
case of linear interpolations of iso-electronic systems, the row elements can either be alchemically compressed to
perturbing potential does differ only by sign. The integral the corresponding second row (n=2) element in the same
over its product with the electron density, however, differs column, or expanded to the fourth row (n=4) element, we
ingeneral,i.e.,⟨H −H ⟩ = dr ρ (v −v )(cid:44) dr ρ (v chose third row (n=3) based reference systems for single,
A B B B A B A A
−v )=−⟨H −H ⟩ . As such, the error in estimating double, and triple bonds, namely, SiH Cl, SiH S, and HSiP.
B B A A 3 2
A based on B will not be the same as the error in The resulting eight alchemical paths are combinations of
estimating B based on A. Results in Table SI72 suggest changing the Si atom (Si→ C, Si→ Ge) or its binding
that predictions downward the columns in the periodic table partner(Cl→ F,Cl→ Br,S→ O,S→ Se,P→ N,P→ As).
174110-8 Changetal. J.Chem.Phys.144,174110(2016)
InFig.4firstandsecondorderalchemicalpredictions,again reference compound for the eight following molecules with
calculatedwithEqs.(12)and(13),areshownforthebonding 14valenceelectrons:CH F,CH Cl,CH Br(n =2);SiH F,
3 3 3 X 3
potentialusingverticaltransmutationsfromthethreereference SiH Br(n =3);andGeH F,GeH Cl,andGeH Br(n =4).
3 X 3 3 3 X
molecules. For double bonds, we have considered predictions for the
More specifically, single bond predictions have been following eight unsaturated molecules 12 valence electrons
investigated for making predictions using SiH Cl as a and using SiH S as a reference compound: CH O, CH S,
3 2 2 2
CH Se (n =2); SiH O, SiH Se (n =3); and GeH O,
2 X 2 2 X 2
GeH S,andGeH Se(n =4).Andfinallyfortriplebonds,we
2 2 X
have studied the following eight molecules with 10 valence
electrons and using HSiP as a reference compound: HCN,
HCP, HCAs (n =2); HSiN, HSiAs (n =3); and HGeN,
X X
HGeP,andHGeAs(n =4).
X
Numerical results in Fig. 4 indicate qualitatively
correct behavior for all predictions. Regarding quantitative
performance, the accuracy of the alchemical prediction of
∆E(d,d )exhibitssimilartrendsastheonediscussedabovein
0
thecaseofverticalchangesinthehydrogencontainingsingle
bond: First order predictions (red) systematically achieve
strong predictive power whenever the change involves the
coupling of the third row element to a fourth row element.
Corresponding second order predictions (blue) deteriorate
the accuracy due to inflection points near λ =0. If the
coupling involves one lighter element from the second row,
the prediction is no longer quantitative. However, in these
cases, second order predictions provide a slightly superior
prediction. If both atoms are simultaneously transmutated to
lighter atoms from the second row, e.g., SiH Cl→ CH F,
3 3
secondorderestimatesover-correct(changeofsign)thefirst
order prediction. In the case of one element transmutating
upward the column, the other downward, the second order
estimateishardlydistinguishablefromthefirstorderestimate.
We believe that the reason for this is that the coupling to
thelighterelementontheonesiteinthemoleculeyieldsthe
concavebehaviorleadingtoanimprovementintheprediction,
whilethecouplingtotheheavierelementontheothersitein
the molecule yields the convex behavior with the inflection
point,leadingtoadeteriorationoftheprediction.Effectively,
these two effects cancel each other and result in the same
predictive accuracy as the one obtained for the first order
estimate.Thisrationalizationrestsontheassumptionthatthe
discussionofFig.2canalsobeappliedtoalinearcombination
ofeffectsatdifferenttransmutatingsites.
2. Integratederrors
Above observations are consistent with the quantitative
integratedpredictionerrormeasures(definitionsinSec.IIF)
summarized in Table SII.72 All first order based predictions
of target molecules implying a transmutation downward the
periodic table (columns 4/3, 3/4, 4/4) exhibit chemical accu-
racywithatmost1.83kcal/moldeviationinminimalenergy
(GeH S),atmost0.04Ådeviationinbondlength(GeH Se),
2 2
at most −12.8 cm−1 deviation in wavenumber (SiH Se),
2
and at most 1.56 kcal/mol in integrated energy (GeHAs).
FIG.4. Alchemicalpredictionsofsingle(top),double(middle),andtriple The best performance is achieved in the case of changing
(bottom)bondpotentials.Curvesareshownforeighttargetsystems(specified SiH Cl→ GeH Brwithenergyerror∆E =0.6kcal/moland
3 3
asinsets),iso-electronicwithreference(graybackground)moleculeSiH3Cl integrated MAE=0.9 kcal/mol. Corresponding predictions
(upperpanel),SiH2S(middlepanel),andHSiP(bottompanel).Potentialscor- of equilibrium distance deviate 0.03 Å with vibration
respondtotrue(blackcircles),first(redsquares),andsecond(bluetriangles)
orderverticalalchemicalpredictionsofheavyatombonddissociationcurves. frequency deviate −1.1 cm−1. First order predictions do not
174110-9 Changetal. J.Chem.Phys.144,174110(2016)
yield quantitative predictive power for changes involving varying geometry and/or atom types and numbers between
lighter elements (columns 2/3, 3/2, 2/4, 2/2, 4/2). The worst referenceandtargetmolecules.
predictions are found for the simultaneous coupling to two
lighter elements (column 2/2) with 76.29 kcal/mol, 0.35 Å,
+
−568.8 cm−1, and 41.45 kcal/mol deviation in minimum 1. AlchemicalstretchingofH2
energy, bond length, harmonic frequency, and integrated We now turn to the case of alchemical stretching of
energy(HCN). H+ in order to understand the effect of varying geometry
2
Whileforallfirstorderpredictionsallmutualdeviations on alchemical predictions. Since Hartree-Fock is formally
exhibit the same sign, second order corrections introduce exactforone-electronsystems,wehaveemployedanatomic
the sign changes in minimum energy and bond length basis set in an “all-electron” (no PPs) calculation within
alluded to before, namely, both third row elements couple the following alignment scheme: One proton is centered at
to lighter elements from the second row (column 2/2). R =(0,0,0),andtheotherisalignedalongthe+x-axis.The
1
Second order predictions for this column are even worse reference system corresponds to H+ at its equilibrium bond
2
than the corresponding first order predictions. Second order length. Stretching is accomplished not by pulling the atoms
predictions only improve first order predictions in the case apart but rather by simultaneous annihilation and creation
of columns 2/3 and 3/2, alas, not to a degree considered of nuclear charges at RR=(d ,0,0) and RT=(d,0,0),
2 eq 2
satisfying. respectively.OncetheSCFisdoneford ,theentirebinding
eq
In summary, for changes corresponding to columns 4/3, potentialcanbeestimatedupto m=4order,usingEq.(11),
3/4, 4/4, first order based estimates yield chemical accuracy. byscanningthroughvariousd’s,i.e.,settingd = d .
0 eq
For changes corresponding to columns 2/2, first order based Results are shown in Fig. 5(a). Due to the variational
estimates are inaccurate but still better than second order principle for linearly coupled alchemical Hamiltonians,12
estimates. For changes corresponding to columns 2/4 and ∆E(1)>∆E forallinteratomicdistances.Inclusionofsecond
4/2, first order based estimates are similar to second order ordertermimprovesuponthefirstorderprediction,yieldinga
estimates,yetbothareinaccurate.Forchangescorresponding reasonable binding potential. However, when including third
to columns 2/3 and 3/2, second order based estimates are and fourth order, the performance deteriorates again with
inaccuratebutstillbetterthanfirstorderestimates. oscillatingbehaviourforvaryingorder(Fig.5(a)andinsetof
Fig.5(b)),as∆E(3)overshootand∆E(4)over-corrects.Overall
∆E(2)givesthebestprediction.
C. Empiricald
opt
ToexplaintheoscillatingbehaviourinTaylorexpansion
Theaboveobservationshavebeenmadeforoptimizedd0. order,weinvestigateinmoredetailhowthesystemresponds
Itshouldbenotedthatthechoiceof d0inEq.(10)iscrucial toalchemicalperturbation.When λ increasesgraduallyfrom
forlinearizingthepropertyofinterestinalchemicalcoupling 0 to 1, the nuclear charge decreases from 1 to 0 at RR,
2
parameter λ and hence essential for the performance of the while increasing from 0 to 1 at RT. Using the alchemical
2
perturbation based predictions. We have found that the error derivatives at λ =0, truncated Taylor series based estimates
minimizing dopt has an approximately linear dependence on are plotted along with the true energy in Fig. 5(b) as a
the target molecule’s equilibrium bond length deq, no matter function of λ at d =3 Å. ∆E(1), ∆E(2), ∆E(3), and ∆E(4) are
if the reference is the hydrogen containing single bond, or linear, quadratic, third order, and fourth order polynomials,
a single, double, or triple bond involving p-block elements respectively. Clearly, the truncated Taylor series will fail to
from second, third, or fourth row. Furthermore, the linear converge to ∆E at λ =1 due to a sharp change of ∆E at
relationshipispreserved,independentofthefactifpredictions λ ≈0.9. This implies a strong nonlinear electronic response
aremadewithfirstorsecondorderestimates.Thisrelationship occurringlateinthealchemicalcouplingregime,resultingin
is shown in Fig. 3. The parameters of a linear regression are the oscillating behaviour of the predicted PES in Figs. 5(a)
specified as well. The outlier in Fig. 3 at target deq≈1.35Å and 5(b). Note that while the sign of error alternates, the
and dopt≈1.0Å is due to the second order prediction of magnitude of error also increases as one increases the order.
SiH3Cl→ CH3F, i.e., for the above discussed worst case Similarbehaviourcanbeobservedforothervaluesofd.
scenario(column2/2)whereastrongovercorrectionhasbeen The energy gain starting at λ ≈0.9 is due to a rapid
found. rearrangement of electron density for λ >0.9. This is
illustrated in Fig. 5(c) where the integrated electron density
P (x) is plotted as a function of both λ and x at d =3 Å.
D. Non-verticaliso-electronicchanges Cλohen and Mori-Sánchez already pointed out for H+ the
2
Wearenotawareofanymathematicallimitationonhow dramatic changes in electronic structure for infinitesimally
to construct alchemical coupling paths under iso-electronic small changes in nuclear charges at infinite distance.65 One
condition. In addition to the investigation of predicted PES wouldexpectthiseffecttointensifyasmorebasisfunctionsare
of iso-electronic compounds with the same geometry, as taken into account. This behaviour can be seen in Fig. 5(c).
discussed in Secs. III A and III B, we have also investigated The locations of the proton at origin R , the annihilated
1
if one can use only one reference calculation in order to proton at RR, and the created proton at RT, are indicated by
2 2
estimatetheentirePESthrough“non-vertical”interpolations. redlines.
In other words, we have also assessed the applicability of Furtheranalysisshowsthatforλ >0.5,bothgroundand
the Taylor expansion of Eq. (9) to non-vertical changes for firstexcitedstateorbitalsarelocalized:Theelectronicground