Atomic basis functions for molecular electronic structure calculations
aa r X i v : . [ phy s i c s . c h e m - ph ] A ug Atomic basis functions for molecular electronic structure calculations
Dimitri N. Laikov a) Chemistry Department, Moscow State University, 119991 Moscow, Russia (Dated: 17 August 2018)
Electronic structure methods for accurate calculation of molecular properties have a high cost that growssteeply with the problem size, therefore, it is helpful to have the underlying atomic basis functions that areless in number but of higher quality. Following our earlier work [Chem. Phys. Lett. , 116 (2005)] wheregeneral correlation-consistent basis sets are defined, for any atom, as solutions of purely atomic functionalminimization problems, and which are shown to work well for chemical bonding in molecules, we take a furtherstep here and define a new kind of atomic polarization functionals, the minimization of which yields addi-tional sets of diffuse functions that help to calculate better molecular electron affinities, polarizabilities, andintermolecular dispersion interactions. Analytical representations by generally-contracted Gaussian functionsof up to microhartree numerical accuracy grades are developed for atoms Hydrogen through Nobelium withinthe four-component Dirac-Coulomb theory and its scalar-relativistic approximation, and also for Hydrogenthrough Krypton in the two-component nonrelativistic case. The convergence of correlation energy with thebasis set size is studied, and complete-basis-set extrapolation formulas are developed.
I. INTRODUCTION
The idea that the molecular electronic structure prob-lem can be solved in terms of atom-centered basis func-tions is as old as the quantum theory of electronitself , but its computational realization has gone along way from qualitative pictures to accurate quanti-tative predictions of molecular properties and reactiv-ity. Atomic functions of exponential type were thefirst to be used for molecular calculations at the Hartree-Fock (HF) and limited configuration interaction (CI) level; within the density-functional theory , theyare still widely used thanks to the numerical integrationand density-fitting schemes speeding up the calcula-tions, and are optimized for all atoms within the zeroth-order regular approximation . Gaussian-type functionsare unique among all classes of elementary functions inthat all multicenter molecular integrals can be computedanalytically , and they have become the standard prim-itive basis in correlated molecular calculations. Otherkinds of functions can be least-squares fitted by sumsof primitive Gaussians, this was done first for the ex-ponentials and was popular for some time, but itsoon became clear that energy-optimized Gaussian ap-proximations of atomic HF wavefunctions work muchbetter and are a must-have part of a good basis set. Thelow variational flexibility of the minimal basis was under-stood, so radial and angular polarization functionsbegan to be added, often tweaked by hand to get bettercomputed properties of a favorite set of small molecules.Later, the so-called diffuse functions were found to beimportant, as were multiple polarization functions, alsoof higher angular momentum, and it was around thattime that the need for a more general construction wasfelt. a) E-Mail: [email protected]; Homepage:http://rad.chem.msu.ru/˜laikov/
The many-body perturbation theory at second (MP2), third , fourth , and seldom fifth or-der, as well as the coupled-cluster theory with sin-gle and double (CCSD) and perturbative triple (CCSD(T)) substitutions are a family of systematicmethods approaching chemical accuracy in the completebasis limit, their fifth- to seventh-power scaling with thesystem size soon makes the integral evaluation, with itsfourth-power scaling, a small part of the work, so onecan be generous in the choice of the primitive expansionlength. Moreover, there are cubic-scaling integral eval-uation techniques with down to quadratic scaling withthe primitive set size: the pseudospectral decomposi-tion can even reduce the scaling of the correla-tion energy calculation; the density-fitting (resolution-of-identity) approximation speeds up the MP2 cal-culations , greatly lowers the memory usage, and ishighly parallelizable.Atomic natural orbitals (ANO) were introduced forthe general contraction of Gaussian basis sets, and thisis indeed a general method as, ideally, only the well-defined atomic solutions would be enough; for the Hy-drogen atom, however, one had to use the H molecule,and the same would have to be done for all atoms withone valence electron. No general way to add the diffusefunctions within the ANO method seems to be found, theauthors resorted to uncontracting or adding primitiveGaussians whose exponents may be quite arbitrary. Un-fortunately, the natural occupation numbers have no di-rect connection to the correlation energy contribution .The overwhelming breakthrough in the field is the de-velopment of the now-classic correlation-consistent ba-sis sets — systematic sequences of energy-optimizedsets of single-Gaussian polarization functions added togenerally-contracted minimal HF sets, with an optionfor core-core and core-valence correlation; the con-vergence with growing set size allows the extrapola-tion to the (apparent) complete basis set limit;single-Gaussian diffuse functions optimized on atomic an-ions can also be added — even though some atomshave a too small or no electron affinity, the limitedsingle-Gaussian functional form helps to get meaningfulexponents for typical molecular applications, more dif-fuse even-tempered sets are used for electrical responseproperties. It is remarkable, how well the single-Gaussianform works for the lighter atoms (up to Ne), althoughslightly less so already for the second row . But it is alsoclear, that a somewhat greater accuracy can be achievedwith the same number of functions if they are made upfrom longer primitive sets .New approaches are evolving: a completeness pro-file can be used to set up a black-box minimizationprocedure that yields completeness-optimized basissets; bound virtual states of multiply-ionized atoms areused as polarization functions in the numerical grid-basedatomic basis sets for density-functional calculations;polarization-consistent basis sets are optimized, on aset of prototypical molecules, for faster convergence of HFenergy, and also for density-functional calculations;aiming at specific intermolecular potentials, interaction-optimized basis sets can be constructed by direct min-imization of the counterpoise-corrected interaction en-ergy for the dimer.We have also found a new way to define atomic ba-sis functions as solutions of atomic variational problems— the closed-shell MP2 correlation energy expression isused as a general minimization functional to drive the op-timization of virtual space. It is generalized to be appliedto all atoms (even Hydrogen) by the use of an effectiveHamiltonian constructed from the spherical average-of-configurations Fock operator, which can be understoodas a model of an atom in a closed-shell-like molecularenvironment. The use of the simplest MP2 functionalseems to be no serious limitation as can be seen fromthe comparison with the classic basis sets of thesame size in molecular CCSD benchmarks. (The molec-ular MP2-optimized virtual space was shown in teststo be close to optimal also for the high-level correlatedmethods such as CCSD.) Our method fits naturally tothe four-component wavefunction theories, such as Dirac-Coulomb Hamiltonian and its scalar-relativistic approxi-mation , and yields the atomically-balanced contractionof kinitically-balanced primitive sets. All chemically-interesting atoms have been covered and the sets havebeen used and worked well in the studies of compoundswith atoms as heavy as actinides . Still, the limi-tations show up when one tries to study intermolecularpotentials, dipole polarizabilities, and negatively-chargedmolecular systems — the functions optimized for the cor-relation energy of the neutral atom are too localized, anddiffuse functions are missing in the set. We sought a gen-eral solution to this problem. Lowest Rydberg stateswere tried in this role , but are too diffuse for a typicaluse. From the known connection between the ionizationpotential and the long-range behavior of the densityfollows that a minimal set of diffuse functions (one foreach angular symmetry) cannot fit for all imaginable an- ionic (or even neutral) species, and a lower bound for theionization potential should be set anyway.After years of thinking, we found new closed-shell-likepolarizability functionals, derived from a simplified (dou-ble) perturbative treatment of atomic electron affinity atthe MP2 level, whose minimization yields a first set of dif-fuse atomic basis functions that already recover most ofthe dipole polarizability and C dispersion coefficient andhelp to get accurate intermolecular potentials for neutralmolecules, and are also helpful for anionic species withstrong enough electron binding. We find a way to get fur-ther diffuse function sets to be used for molecules with aslow ionization potential as needed. All these variationalprocedures can be used to contract the primitive sets ofGaussian (or other well-behaved) functions, and we havealso worked out a weighted least-squares fitting techniqueto optimize the primitive sets of any size needed to geta given accuracy. Here we present these new ideas (aftera short review of our older underlying work ) and showhow they are used to make a database of atomic basissets for all atoms from Hydrogen through Nobelium. Wealso study the convergence and extrapolation towards thecomplete basis limit on a few atoms and molecules. II. THEORY
Our general atomic basis functions can be eitherfour-component spinors of the Dirac-Coulomb theory(or its scalar-relativistic approximation ), or the non-relativistic two-component wavefunctions. They can bedefined either as solutions of integro-differential equa-tions or as linear combinations of some primitive func-tions whose coefficients are solutions of constrained func-tional optimization problems, the latter algebraic formu-lation is more practical and will be given here.We begin with a set of N occupied wavefunctions { φ i } that make the average-of-configurations Hartree-Fock en-ergy E = N X i =1 w i H ii + N X i,j =1 a ij w i R ii,jj w j (1)stationary (the orthogonality constraints are always im-plied). The one-electron { H µν } and antisymmetrizedtwo-electron { R κλ,µν } integrals are computed over thegiven set of N two-component (or 2 N four-component)primitive functions { φ µ } and transformed with coeffi-cients C µi to those in Eq. (1). The occupancies 0 < w i ≤ a ij = a ji are such as keepthe spherical symmetry of a neutral atom, in the closed-shell case w i = 1, a ij = 1. The stationarity conditionsare H ui + N X j =1 a ij R ui,jj w j = 0 (2)where u > N counts all unoccupied states, and also( w i − w j ) H ij + N X k =1 ( a ik w i − a jk w j ) R ij,kk w k = 0 . (3)The occupied energies ǫ i = F o ii (4)are set to the diagonal elements of the Fock matrix F o ij = H ij + N X k =1 ( a ik + a jk ) R ij,kk w k , (5)and if the energy of Eq. (1) is invariant to rotations forsome ij -pair, then the condition F o ij = 0 for i = j, w i = w j , a ik w i = a jk w j , (6)should also be met. For the virtual subspace, the Fockmatrix is taken to be F v µν = H µν + N X i =1 R µν,ii w i , (7)in the four-component case, it is diagonalized in that sub-space to get N − N electronic and N positronic states,the latter being discarded and all further work is donewithin the electronic part.For some atoms, there are low-lying states not includedin the average of Eq. (1) but important for chemicalbonding, so N functions should be added to the occu-pied set { φ i } to account for it, and this is done by thediagonalization of a Fock matrix F + µν = H µν + N X i =1 R µν,ii w + i (8)in the space of N − N virtual electronic wavefunctions, F + ij = ǫ i δ ij for i, j > N , (9) N such states with energies ǫ i go into the occupied setthat from now on has the size N o = N + N . The occu-pancies w + i are for the spherically-averaged configurationof the singly-ionized atom.At this time, we build the effective Fock operatorˆ F = N o X i =1 | φ i i ǫ i h φ i | + N o X i,j =1 (cid:0) − | φ i i h φ i | (cid:1) ˆ F v (cid:0) − | φ j i h φ j | (cid:1) (10)that is spherically symmetric, and together with thetwo-electron interaction it is all that is needed as in-put to a correlation energy calculation. We forget aboutthe fractional occupations — now we work with the(pseudo)atom as if it had the closed-shell configurationwith N o fully filled levels. A set of N v virtual wavefunctions { φ u } (needed forelectron correlation, atomic polarization upon chemicalbonding, electron affinity, and dispersion interaction inmolecules) can be grown stepwise by minimization of thefunctionals (shown below) defined in terms of the lin-early transformed set { φ a } with coefficients { C ua } thatdiagonalize their block of the Fock matrix N v X v = N o +1 F v uv C va = ǫ a C ua (11)and have energies { ǫ a } . The stepwise growth means thata set of N (1)v functions { φ u } is optimized first for onefunctional; next, keeping these N (1)v frozen, N (2)v func-tions are added and optimized for another functional ofall N (1)v + N (2)v functions, and so on. (Sometimes (seeSection III) the very first set of N (0)v functions is takenas the next lowest energy states of Eq. (9) after the lowest N .)The minimization of a model second-order correlationenergy functional E = X ijab p ij | R ai,bj | ǫ i + ǫ j − ǫ a − ǫ b (12)with respect to some members of the set { φ u } yields func-tions nearly optimal for the most part of electron correla-tion in atoms and molecules. The factors p ij in Eq. (12)are set to 1 or 0 to switch the valence-only, core-valence,and core-core correlation. The stationarity conditionscan be derived by chain-rule differentiation with respectto rotations that mix the external { φ x } and the virtual { φ u } functions, and can be written as G xu = X a G xa + X b F xb D ab ! C ua = 0 (13)with G xa = X ijb p ij R xi,bj R ai,bj ǫ i + ǫ j − ǫ a − ǫ b , (14) D ab = X ijc p ij R ai,cj R bi,cj ( ǫ i + ǫ j − ǫ a − ǫ c )( ǫ i + ǫ j − ǫ b − ǫ c ) . (15)The second derivatives of Eq. (12) can also be derivedand used in an efficient quadratically-convergent Newton-Raphson optimization with a careful steps size control.Until now, we have followed our earlier work , butsince then the experience has shown that such atomic ba-sis sets lack diffuse functions needed for accurate calcula-tion of electron affinities, polarizabilities, and dispersioninteractions in molecules. Now, we have found a modelpolarization functional E a = X ia p i (cid:12)(cid:12) ¯ U ai (cid:12)(cid:12) ǫ i − ǫ a , (16)¯ U µν = N X i =1 R µν,ii ¯ w i , (17)whose minimization yields a set of functions, with angu-lar momenta up to those in the occupied set, that helpto account for the changes upon electron attachment (ordetachment). The occupancies ¯ w i add up to one and arespherically-averaged, typically ¯ w i = w i − w + i , the switch-ing factors p i are set to 1 or 0. Eq. (16) can be under-stood as the second-order perturbative correction for theaveraged electron attachment energy at the Hartree-Focklevel, and we believe it to be a better functional for driv-ing the diffuse function optimization than the directlycomputed energy of an atomic anion — some atoms havea tiny or no electron affinity so the functions may becometoo diffuse, but the perturbative first order changes aremore localized and still in the right direction. The sta-tionarity conditions are as in Eq. (13) but now with G xa = X i p i ¯ U xi ¯ U ai ǫ i − ǫ a , (18) D ab = X i p i ¯ U ai ¯ U bi ( ǫ i − ǫ a )( ǫ i − ǫ b ) , (19)and there is a simple explicit solution.To optimize the diffuse functions with higher-than-occupied angular momenta, we model the changes in theMP2 correlation energy upon electron attachment. Theperturbed occupied wavefunctions¯ φ i = φ i + N X a = N o +1 φ a ¯ U ai ǫ i − ǫ a , (20)with the sum running here over the whole N − N o vir-tual space in diagonal representation of the matrix ofEq. (7), are computed first, and then the set { ¯ φ i } is usedin Eq. (12) instead of { φ i } to get the new functional ¯ E .The difference ∆ ¯ E = ¯ E − E (21)is a measure of how the members of the set { φ u } underoptimization help to lower the correlation energy of theatomic anion more than of the neutral atom, and thisis the functional that is minimized to get the set of dif-fuse functions. With the intermediate normalization ofEq. (20) one term in ¯ E is exactly canceled by E , andwe like it.Here we can stop, as we now have enough tools tobuild systematic sequences of atomic basis set for a broadrange of molecular applications. Still, we should be awareof their limitations and would like to take a look aheadtowards a better sampling of the diffuse tail region. Wehave studied the dipole polarizability functional E d = X ia p i | r ai | ǫ i − ǫ a , (22) and the homoatomic C dispersion coefficient functional E = X ijab p ij | r ai | | r bj | ǫ i − ǫ j − ǫ a − ǫ b , (23)where r ai are the dipole moment integrals, minimizationof either of them yields a set of functions that are some-what more diffuse than those based on Eqs. (16) and (21)— these can be added to the set after the former ones,but our molecular tests show this to be of little help formost molecular systems, even for noble gas dimers thereis only a small lowering of the potential energy curve.Higher multipole analogs of Eqs. (22) and (23) could havebeen studied, but we do not feel it to be the right wayforward.A smooth sampling of the diffuse tails can be done witha one-parameter family of model Fock operatorsˆ F ζ = ˆ F + ζ ˆ¯ U, (24)with ˆ¯ U from Eq (17) and 0 ≤ ζ ≤ ζ max , that can bediagonalized (with a frozen core option) to get the wave-functions of the form¯ φ iζ = φ i + N X a = N o +1 φ a ¯ C aiζ (25)with energies ¯ ǫ iζ , and ζ max can be found such that¯ ǫ N o ζ max = ǫ max . A good ǫ max < −√− ǫ max r )for the new members of the set { φ u } , such that if the ˆ F ζ is diagonalized in the subspace to get the wavefunctions˜ φ iζ = φ i + N v X u = N o +1 φ u ˜ C uiζ (26)with energies ˜ ǫ iζ , then the integral˜ ǫ = ζ max Z X i ˜ ǫ iζ d ζ (27)is minimized. It can be seen that the functional ofEq. (16) is a lowest-order perturbative approximation tothat of Eq. (27), and the functions of Eq. (20) are noth-ing else than first-order perturbative analogs of thoseof Eq. (25). In the same way, to get the higher-than-occupied functions, the integral∆ ¯ E = ζ max Z (cid:0) ¯ E ( ζ ) − E (cid:1) d ζ (28)can be minimized, where ¯ E ( ζ ) is built upon the set { ¯ φ iζ } instead of { φ i } in E of Eq. (12).The general atomic basis functions outlined above canbe computed to any meaningful accuracy by numericallysolving the underlying variational problems. Once the(nearly) exact solutions { φ k } are at hand, practical ap-proximations, such as the traditional contracted Gaus-sian or exponential (Slater-type) functions, can be de-veloped for use in molecular calculations. First, we opti-mize the exponents of the primitive functions by weightedleast-squares fitting, minimizing Q = N o + N v X k =1 β k Z (cid:12)(cid:12)(cid:12) ˜ φ k ( r ) − φ k ( r ) (cid:12)(cid:12)(cid:12) w ( r ) d r (29)with respect to both the linear coefficients { ˜ c nk } and theprimitive exponents { α n } of the approximate functions { ˜ φ k } , with the weights β k = (cid:26) β o , k ≤ N o , k > N o (30)made heavier for the occupied set by β o , and the radialweight functions w ( r ) = 1 / | r | . After much experimen-tation, we set β o = 2 that gives a good balance be-tween the HF and correlation energies, and our choice of w ( r ) leads to the exponents that are close to the energy-optimized values in the HF case. We have tried to putthe orthonormality constraints on { ˜ φ k } in the fitting, butfound it to heavily complicate and slow down the compu-tations without giving better results, so we put it aside.There is a known pitfall in the work with the exponents— some of them may be driven towards the same valuewhereas the linear coefficients of opposite sign go towardsinfinity — and this is the true but numerically unstableand impractical solution. To overcome this, we put abound on the closeness and parametrize the exponentsas α n = exp p + n X m =2 q p + p m ! , (31)and { p n } , n ≥
1, are now taken as the optimization vari-ables instead of { α n } , and thus α n +1 /α n ≥ exp( p ), wesettle on p = (ln 2) / III. CALCULATIONS
We have written a computer code for solving theatomic variational problems of section II with full useof spherical symmetry — the angular degrees of freedomare integrated out analytically and only the radial equa-tions are worked with. Extended precision floating-pointarithmetics with 256- or 128-bit mantissa is implementedusing the
X86 64 of the speed of light c = 137 . with Gaus-sian charge distribution with exponent (in au) α = (cid:18) · √ M (cid:19) , (32)where M is the (integer) mass number of the most abun-dant isotope; point nucleus is used in the non-relativisticcase. We solve the variational problems of Section II to avery high accuracy over a huge even-tempered primitiveGaussian basis with exponents α p = 2 p/ (33)where − ≤ p ≤ − − for the occupied set and somewhat less for the virtual. Inthe non-relativistic case, a range − ≤ p ≤
225 wouldbe needed to meet the nuclear cusp condition, but wefind it more practical to use α p = 2 p/ + exp( ap − b ) (34)with a and b optimized for each atom on its hydrogen-likeion with one electron, and a more narrow range of p .The standardized electronic configurations of atomsare shown in Table I, where L + is the angular momentumfor which one electron is removed to get the occupan-cies w + i in Eq. (8), L is for the N functions of Eq. (9)added to the occupied set, L is for the first N (0)v virtualfunctions also from Eq. (9), and the atoms are markedfor which the diffuse functions may be added. To get thecoupling coefficients in Eq. (1), we use the average level formalism in the Dirac-Coulomb case, and also the usualaverage of the highest-spin configurations in the scalar-relativistic and non-relativistic cases.The stepwise optimization begins with the outermostoccupied shell block (of the same principal quantum num-ber n ), for which the functional of Eq. (12) with valence-only p ij is minimized first to get a set of virtual func-tions with N v ln radial parts for each angular momentum l , then a set of diffuse functions may be optimized usingEq. (16) followed by Eq. (21); the next block of innervalence (for transition metals) or core shells may then betaken and the virtuals optimized using Eq. (12) with p ij set to core-core and core-valence correlation, and so onto the innermost core shells. The correlation-consistentvirtual set sizes at each step are taken to be N v ln = max (cid:0) λ n − max( l − l max n − , , (cid:1) (35)where l max n is the highest occupied angular momentum ofthe n -block, and λ n is the set number (1 , , . . . ) for this n — thus l reaches up to l max n + λ n and there are λ n radialparts for l ≤ l max n + 1. The size of the diffuse set is simply N v ln = 1 for l ≤ l max n + λ n . The same λ n is typically usedfor all n , but our code can work with any other settings.For most metal atoms, the outermost core shells shouldbe unfrozen, Table I shows the number of the n -blocksthat are always correlated. One may also note the atomsCa, Sr, Ba, and Ra to have a shell added ( L ) to theunoccupied set, we found it to be the key to get the accu-rate bonding properties of these “subtransition” metals,to cure the known pathology . TABLE I. Atomic electronic configurations. atoms number of electrons a L + L L diffuse b n c H, He 1+ 0 + 1Li, Be 3+ 0 1 2B – Ne 4 ,
1+ 1 + 1Na, Mg 5+ , ,
7+ 1 + 1K 7 ,
12 0 1 2Ca 8 ,
12 0 1 2 2Sc – Zn 8 , ,
1+ 0 1 2Ga – Kr 8 , ,
10 1 + 1Rb 9 , ,
10 0 1 2Sr 10 , ,
10 0 1 2 2Y – Cd 10 , ,
11+ 0 1 2In – Xe 10 , ,
20 1 + 1Cs 11 , ,
20 0 1 2Ba 12 , ,
20 0 1 2 2La 12 , ,
21 0 1 2Ce – Yb 12 , , ,
1+ 0 1 3Lu – Hg 12 , , ,
14 0 1 2Tl – Rn 12 , , ,
14 1 1Fr 13 , , ,
14 0 1 2Ra 14 , , ,
14 0 1 2 2Ac 14 , , ,
14 0 1 2Th – No 14 , , ,
15+ 0 1 3 a Number of electrons for each angular momentum; the plus + meansthe number grows in the row starting with the given value. b Whether the diffuse functions are added. c At least n of the outermost principal quantum number shell blocksare correlated. With the nearly-exact solutions at hand, the least-squares optimization based on Eq. (29) is run to getthe primitive basis sets of growing size — this can bedone either separately and independently for each angu-lar momentum l , or all at once with the exponents sharedbetween all l . The former is more flexible, economical,and natural, so we do most of our work this way; butthe latter is helpful to speed up some electronic structuremethods if the integral evaluation can make use of sharedexponents, so we also do it on a smaller scale. The non-linear optimization of exponents sometimes finds multi-ple minima, and we often had to feed it with several setsof starting values, made by hand, to get either the low-est error or, more seldom, a more regular variation ofthe exponents with the atomic number. We have alwaysseen the close-to-exponential convergence of the fit er-ror with the primitive set size M l , and the same is alsotrue for both the HF- and MP2-energy errors, which weestimate for a given M l by running atomic calculations on the sets of size { ¯ M , . . . , ¯ M l − , M l , ¯ M l +1 , . . . , ¯ M l max } and { ¯ M , . . . , ¯ M l max } and subtracting the energies, where¯ M l are big enough. Thus we grade the primitive setsby their energy errors E l ( M ) for each l , and now wehave to decide on the standard set sizes { M ( κ ) l } so that E l (cid:16) M ( κ ) l (cid:17) ≈ E ( κ ) for all l and for both HF and MP2, thatis, the errors are of the same order, and for κ = 1 , , . . . we should have the errors close to a geometric series E ( κ ) / E ( κ − ≈ ε . As a guide, we look at the MP2 valuesof E l max ( κ ) for κ = 1 , . . . , ε ≈ , and so we align all our sets for all atoms.For λ n = 1 , , , κ = 1 , , , ,
5, with and without the correlation ofthe outermost or all core electrons, we have optimizedthe series of basis sets for atoms H through No (7682 innumber) for both the Dirac-Coulomb Hamiltonian and itsscalar-relativistic approximation, and also for H throughKr for the non-relativistic case. It took us years of hardwork at which we grew old and sick, so we cannot giveall details here. Instead, the files in the supplementarymaterial hold all our data sets (to 64-bit precision) andeveryone is welcome to use them or to study their prop-erties.Our mnemonics for the optimized atomic basis sets are: • L λ κ for the valence-only correlation, • L λ a κ the same with the diffuse functions, • L λλ κ , L λλλ κ with the outermost core shells in-cluded into the correlation, • L λλ a κ , L λλλ a κ the same with the diffuse func-tions, • Lx λ κ for all-electron correlation, • Lx λ a κ the same with the diffuse functions.It would have been good to test all these on a set ofmolecules, but here we will only study the convergenceand extrapolation towards the complete basis limit on afew simple but characteristic examples. TABLE II. MP2 correlation energy of closed-shell atoms.
He Ne Ni λ L λ L λ a L λ L λ a L λ . . . ∞ -0.0373774 -0.32021 -1.2032 The atoms He, Ne, and Ni in Table II are archetypalfor the closed-shell correlation, and we took pains to goup to λ = 9 to come close to the asymptotic behavior.(For Ni, the sets are not for the standard configurationof Table I but for the closed-shell state with 18 outerelectrons correlated.)A natural functional form E λ ≈ E ∞ + P X p =3 A p ( λ + ν ) p (36)with some small P > E λ for a range of λ and thus to get an estimate of E ∞ .The only nonlinear parameter ν in Eq. (36) can be ad-justed each time, but we find ν = to be a good fixedvalue, often very close to the optimal one, and we set itso everywhere in the following. By fitting through P − λ = λ , . . . , λ + P −
2, each time for a higher λ ,we can also estimate the residual error of the last E ∞ ( λ )as E ∞ ( λ ) − E ∞ ( λ − E ∞ in Table II seemingly converged to all digits given. Withthese at hand, we can try to find a simple and practicaltwo-point extrapolation formula of the kind E ∞ ≈ E λ + ( E λ − E λ − ) c λ (37)that would follow if it would hold that E λ ≈ E ∞ + Ab λ , (38)with the universal b λ , so that c λ = b λ / ( b λ − − b λ ) . (39)Computing˜ c λ = ( E ∞ − E λ ) / ( E λ − E λ − ) ≈ c λ (40)over the data set of Table II, we see that ˜ c λ are weaklysystem dependent, and a conservative approximation toEq. (38) is simply the shortest form of Eq. (36), E λ ≈ E ∞ + A/ (cid:0) λ + (cid:1) , (41)and thus c λ = 1 .(cid:16)(cid:0) / ( λ + ) (cid:1) − (cid:17) . (42)We have tried to find better extrapolation formulas bysplitting the two-electron correlation energy into its spincomponents — the same-spin part is known to have A = 0 in Eq. (36) — but could not get a higher overallaccuracy. The CCSD correlation energy can as well beextrapolated in the same way — we have found it to beof no help to split it into the MP2 and higher-order parts,as the latter does not seem to show a regular behavior,at least for smaller λ . The perturbative triples energyof CCSD(T), however, can be extrapolated well enoughusing the fourth power instead of the third in Eqs. (41)and (42).Molecular tests in Tables III and IV show the conver-gence with respect to both the number of basis functions TABLE III. MP2 calculations on H , N and LiF molecules. H N LiFset r ∆ E r ∆ E r ∆ E L1 1
L1 2
L1 3
L1 4
L1 5
L1 6
L2 1
L2 2
L2 3
L2 4
L2 5 -0.167002 -0.372687 -0.229689
L3 1
L3 2
L3 3
L3 4
L3 5 -0.167683 -0.381103 -0.233453
L4 1
L4 2
L4 3
L4 4
L4 5 -0.167828 -0.383306 -0.233674
L5 1
L5 2
L5 3
L5 4 -0.167878 -0.383787 -0.233665
L6 4
L1a 1
L1a 2
L1a 3
L1a 4
L1a 5
L1a 6
L2a 1
L2a 2
L2a 3
L2a 4
L2a 5 -0.167475 -0.377247 -0.230105
L3a 1
L3a 2
L3a 3
L3a 4
L3a 5 -0.167770 -0.383046 -0.233407
L4a 1
L4a 2
L4a 3
L4a 4
L4a 5 -0.167876 -0.383797 -0.233438
L5a 1
L5a 2
L5a 3
L5a 4 -0.167889 -0.383981 -0.233446
The bond lengths r and the binding energies ∆ E are in au,the values in italics are extrapolated from those on the line above. and the quality of their approximation with the under-lying primitive set size. The three molecules, H , N ,and LiF, are prototypical for covalent and ionic bondingand, because of the very regular and consistent structureof our basis sets across the periodic table, Table III can guide the choice of the right set for a given application.At the other end, the weakest bonding in the noble gasdimers He and Ne studied in Table IV should be under-stood as the worst case performance — we see here that,without the diffuse functions, the bond lengths and ener-gies do converge but too slowly, and our diffuse sets helpto recover the most part of the attractive interaction thatis somewhat overestimated and can be (over)corrected bythe counterpoise method so that the two binding ener-gies seem to bracket the “exact” value. TABLE IV. MP2 calculations on He and Ne dimers. He Ne set r ∆ E r ∆ E L1 6
L2 5
L3 5
L4 5
L5 4
L6 4
L7 4
L1a 1
L1a 2
L1a 3
L1a 4
L1a 5
L1a 6
L2a 1
L2a 2
L2a 3
L2a 4
L2a 5 -19.20 -20.63 -90.09 -71.61
L3a 1
L3a 2
L3a 3
L3a 4
L3a 5 -20.10 -20.88 -60.76 -79.39
L4a 1
L4a 2
L4a 3
L4a 4
L4a 5 -21.35 -21.66 -74.72 -82.83
L5a 1
L5a 2
L5a 3
L5a 4 -22.09 -22.22 -81.62 -85.11
L6a 4 E are in microhartrees and computedeither relative to the isolated atoms (first column)or with the counterpoise correction (second column),the values in italics are extrapolated from those on the line above.The bond lengths r are in bohrs. The extrapolation of the correlation energy at leastdoes not hurt these weakest bonds, leading to some un-derbinding, but is a great improvement for the strongchemical bonds as seen in Table III, so it should be help-ful for all molecular systems.The dipole polarizabilities in Table V clearly witnessthe need for the diffuse functions, using Eq. (22) we addone more “ b ”-set in L λ ab which yields the highest accu-racy, but the L λ a are already quite good and should beused to get the accurate intermolecular interactions.The minimal primitive representation of angular polar-ization functions of our L λ and L λ a 1 sets for lighteratoms can be compared one-to-one with that of the clas-sical works , a remarkably close match of the expo-nents for atoms He, B–Ne, and Al-Ar can be seen (less TABLE V. MP2 polarizabilities (au) of atoms and molecules. set He Ne H N L1 6
L2 5
L3 5
L4 5
L5 4
L1a 6
L2a 5
L3a 5
L4a 5
L5a 4
L1ab 6
L2ab 5
L3ab 5
L4ab 5 so for H as we take the atom and not the H molecule),and this must be a very good sign for us all. IV. CONCLUSIONS
The atomic basis sets are now made available thatallow quantitative calculations of the structure and en-ergetics of typical molecular systems. We are somewhatsorry for our belated report, for they have already been,and are being, used by our colleagues in their studies of unusual molecules under the unusual conditions of low-temperature high-energy chemistry .They can also be used to set up a database of ac-curate reference values for the parametrization of elec-tronic structure models or molecular mechanicalforce fields of all kinds. We are looking forward to theirfurther useful applications for the good of chemistry,chemists, and society. Atoms are many, molecules areendless, but we are few, not to say alone. W. Heitler and F. London, Z. Physik , 455 (1927). J. E. Lennard-Jones, Trans. Faraday Soc. , 668 (1929). E. Schr¨odinger, Ann. Phys. , 361 (1926). E. Schr¨odinger, Phys. Rev. , 1049 (1926). P. A. M. Dirac, Proc. R. Soc. Lon. Ser. A , 610 (1928). C. Zener, Phys. Rev. , 51 (1930). J. C. Slater, Phys. Rev. , 57 (1930). D. R. Hartree, P. Camb. Philos. Soc. , 89 (1928). V. Fock, Z. Physik , 126 (1930). C. C. J. Roothaan, Rev. Mod. Phys. , 69 (1951). P.-O. L¨owdin, Phys. Rev. , 1474 (1955). W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965). E. J. Baerends, D. E. Ellis, and P. Ros,Chem. Phys. , 41 (1973). E. V. Lenthe and E. J. Baerends,J. Comp. Chem. , 1142 (2003). E. van Lenthe, E. J. Baerends, and J. G. Snijders,J. Chem. Phys. , 4597 (1993). S. F. Boys, Proc. R. Soc. A , 542 (1950). W. J. Hehre, R. F. Stewart, and J. A. Pople,J. Chem. Phys. , 2657 (1969). R. F. Stewart, J. Chem. Phys. , 431 (1970). R. Ditchfield, W. J. Hehre, and J. A. Pople,J. Chem. Phys. , 5001 (1970). R. Ditchfield, W. J. Hehre, and J. A. Pople,J. Chem. Phys. , 724 (1971). W. J. Hehre, R. Ditchfield, and J. A. Pople,J. Chem. Phys. , 2257 (1972). P. C. Hariharan and J. A. Pople,Theor. Chim. Acta , 213 (1973). M. J. Frisch and J. A. Pople, J. Chem. Phys. , 3265 (1984). C. Møller and M. S. Plesset, Phys. Rev. , 618 (1934). R. J. Bartlett, J. Chem. Phys. , 3258 (1975). J. A. Pople, R. Seeger, and R. Krishnan,Int. J. Quantum Chem. , 149 (1977). R. Krishnan and J. A. Pople,Int. J. Quantum Chem. , 91 (1978). R. Krishnan, M. J. Frisch, and J. A. Pople,J. Chem. Phys. , 4244 (1980). K. Raghavachari, J. A. Pople, E. S. Replogle, and M. Head-Gordon, J. Phys. Chem. , 5579 (1990). F. Coester, Nucl. Phys. , 421 (1958). F. Coester and H. K¨ummel, Nucl. Phys. , 477 (1960). J. ˇCiˇzek, J. Chem. Phys. , 4256 (1966). G. D. Purvis and R. J. Bartlett,J. Chem. Phys. , 1910 (1982). N. C. Handy, J. A. Pople, M. Head-Gordon, K. Raghavachari,and G. W. Trucks, Chem. Phys. Lett. , 185 (1989). K. Ragavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon, Chem. Phys. Lett. , 479 (1989). R. A. Friesner, J. Phys. Chem. , 3091 (1988). M. N. Ringnalda, M. Belhadj, and R. A. Friesner,J. Chem. Phys. , 3397 (1990). B. H. Greeley, T. V. Russo, D. T. Mainz, R. A. Friesner, J. Lan-glois, W. A. Goddard, III, R. E. Donnelly, Jr., and M. N. Ring-nalda, J. Chem. Phys. , 4028 (1994). V. Termath and N. C. Handy,Chem. Phys. Lett. , 17 (1994). R. B. Murphy, W. T. Pollard, and R. A. Friesner,J. Chem. Phys. , 5073 (1997). R. Izs´ak and F. Neese, J. Chem. Phys. , 144105 (2011). T. J. Martinez and E. A. Carter,J. Chem. Phys. , 7081 (1993). D. M. Schrader and S. Prager, J. Chem. Phys. , 1456 (1963). J. L. Whitten, J. Chem. Phys. , 4496 (1973). N. H. F. Beebe and J. Linderberg,Int. J. Quantum Chem. , 683 (1977). C. V. Alsenoy, J. Comp. Chem. , 620 (1988). O. Vahtras, J. Alml¨of, and M. W. Feyereisen,Chem. Phys. Lett. , 514 (1993). M. Feyereisen, G. Fitzgerald, and A. Komornicki,Chem. Phys. Lett. , 359 (1993). F. Weigend and M. H¨aser, Theor. Chem. Acc. , 331 (1997). J. Alml¨of and P. R. Taylor, J. Chem. Phys. , 4070 (1987). R. C. Raffenetti, J. Chem. Phys. , 4452 (1973). J. Alml¨of and P. R. Taylor, J. Chem. Phys. , 551 (1990). N. Suaud and J. P. Malrieu, Mol. Phys. , 2684 (2017). T. H. Dunning, J. Chem. Phys. , 1007 (1989). D. E. Woon and T. H. Dunning,J. Chem. Phys. , 1358 (1993). D. E. Woon and T. H. Dunning,J. Chem. Phys. , 4572 (1995). K. A. Peterson and T. H. Dunning,J. Chem. Phys. , 10548 (2002). W. Lakin, J. Chem. Phys. , 2954 (1965). R. N. Hill, J. Chem. Phys. , 1173 (1985). J. M. L. Martin and P. R. Taylor,J. Chem. Phys. , 8620 (1997). A. Halkier, T. Helgaker, P. Jørgensen, W. Klop-per, H. Koch, J. Olsen, and A. K. Wilson,Chem. Phys. Lett. , 243 (1998). A. J. C. Varandas, J. Chem. Phys. , 244105 (2007). J. G. Hill, K. A. Peterson, G. Knizia, and H.-J. Werner,J. Chem. Phys. , 194105 (2009). R. A. Kendall, T. H. Dunning, and R. J. Harrison, J. Chem. Phys. , 6796 (1992). D. E. Woon and T. H. Dunning,J. Chem. Phys. , 2975 (1994). A. K. W. T. H. Dunning, K. A. Peterson,J. Chem. Phys. , 9244 (2001). T. Hashimoto, K. Hirao, and H. Tatewaki,Chem. Phys. Lett. , 345 (1997). D. P. Chong, Can. J. Chem. , 79 (1995). P. Manninen and J. Vaara, J. Comp. Chem. , 434 (2006). S. Lehtola, P. Manninen, M. Hakala, and K. Hamalainen,J. Chem. Phys. , 044109 (2013). B. Delley, J. Chem. Phys. , 508 (1990). F. Jensen, J. Chem. Phys. , 9113 (2001). F. Jensen, J. Chem. Phys. , 9234 (2002). F. Jensen, J. Chem. Phys. , 7372 (2002). J. VandeVondele and J. Hutter,J. Chem. Phys. , 114105 (2007). J. G. C. M. van Duijneveldt-van de Rijdt and F. B. van Duijn-eveldt, J. Chem. Phys. , 3812 (1999). S. F. Boys and F. Bernardi, Mol. Phys. , 553 (1970). D. N. Laikov, Chem. Phys. Lett. , 116 (2005). L. Adamowicz and R. J. Bartlett,J. Chem. Phys. , 6314 (1987). K. G. Dyall, J. Chem. Phys. , 2118 (1994). G. A. Shamov, G. Schreckenbach, and T. N. Vo,Chem. Eur. J. , 4932 (2007). G. A. Shamov and G. Schreckenbach,J. Phys. Chem. A , 9486 (2006). Y. A. Ustynyuk, I. P. Gloriozov, S. N. Kalmykov, A. A. Mitro-fanov, V. A. Babain, M. Y. Alyapyshev, and N. A. Ustynyuk,Solv. Extr. Ion Exch. , 508 (2014). Y. A. Ustynyuk, S. N. Kalmykov, M. Alyapyshev, V. Babain,H. V. Lavrov, S. S. Zhokhov, P. I. Matveev, I. P. Glorio-zov, L. I. Tkachenko, I. G. Voronaev, and N. A. Ustynyuk,Dalton Trans. , 10926 (2017). E. V. Saenko, D. N. Laikov, I. A. Baranova, and V. I. Feldman,J. Chem. Phys. , 101103 (2011). N. C. Handy, M. T. Marron, and H. J. Silverstone,Phys. Rev. , 45 (1969). M. M. Morrell, R. G. Parr, and M. Levy,J. Chem. Phys. , 549 (1975). T. Aoyama, M. Hayakawa, T. Kinoshita, and M. Nio,Phys. Rev. Lett. , 111807 (2012). D. Hanneke, S. F. Hoogerheide, and G. Gabrielse,Phys. Rev. A , 052122 (2011). L. Visscher and K. G. Dyall,Atom. Data Nucl. Data , 207 (1997). I. P. Grant, D. F. Mayers, and N. C. Pyper,J. Phys. B: At. Mol. Phys. , 2777 (1976). S. S. Wesolowski, E. F. Valeev, R. A. King, V. Baranovski, andH. F. Schaefer, Mol. Phys. , 1227 (2000). Supplementary material. I. S. Sosulin, E. S. Shiryaeva, D. A. Tyurin, and V. I. Feldman,J. Chem. Phys. , 131102 (2017). S. V. Ryazantsev, D. A. Tyurin, and V. I. Feldman,Spectrochim. Acta A , 39 (2017). S. V. Ryazantsev, D. A. Tyurin, V. I. Feldman, and L. Khri-achtchev, J. Chem. Phys. , 184301 (2017). S. V. Kameneva, D. A. Tyurin, and V. I. Feldman,Phys. Chem. Chem. Phys. , 24348 (2017). S. V. Kameneva, D. A. Tyurin, K. B. Nuzhdin, and V. I. Feld-man, J. Chem. Phys. , 214309 (2016). E. S. Shiryaeva, D. A. Tyurin, and V. I. Feldman,J. Phys. Chem. A , 7847 (2016).
I. S. Sosulin, E. S. Shiryaeva, D. A. Tyurin, and V. I. Feldman,J. Phys. Chem. A , 4042 (2018).
D. N. Laikov, J. Chem. Phys. , 134120 (2011).
K. R. Briling, J. Chem. Phys.147