Atomic effective potentials for starting molecular electronic structure calculations

TL;DR

This paper introduces a set of atomic effective one-electron potentials, expressed as Gaussian charge distributions, optimized for all elements from Hydrogen to Nobelium, to improve initial guesses in molecular electronic structure calculations.

Contribution

It presents a compact, one-parameter-per-atom model of atomic effective potentials optimized across isolated atoms and molecules, enhancing starting points for electronic structure computations.

Findings

Effective potentials are accurate for a wide range of elements.

The model simplifies initial molecular calculations.

Potential improves computational efficiency and accuracy.

Abstract

Atomic effective one-electron potentials in a compact analytic form in terms of a few Gaussian charge distributions are developed, for Hydrogen through Nobelium, for starting molecular electronic structure calculations by a simple diagonalization. For each element, all terms but one are optimized in an isolated-atom Hartree--Fock calculation, and the last one is parametrized on a set of molecules. This one-parameter-per-atom model gives a good starting guess for typical molecules and may be of interest even on its own.

Tables

Table 1. Table 1: Atomic parameters.

Atoms	$n$	$a_{n+1}$	$ℰ$a	$\Delta E$a,b
H	0	$1/3$
He	1	$1/3$	0.000024	0.0004
Li	1	$1/16$	0.000154	0.0050
Be	2	$1/16$	0.000075	0.0040
B–Ne	2	$1/3$	0.000320	0.0094
Na	2	$1/32$	0.000379	0.0135
Mg	3	$1/32$	0.000099	0.0069
Al–Ar	3	$1/8$	0.000143	0.0103
K–Ca	3	$1/32$	0.000275	0.0138
Sc–Zn	4	$1/6$	0.000029	0.0191
Ga–Kr	4	$1/12$	0.000180	0.0210
Rb–Sr	4	$1/32$	0.000085	0.0159
Y–Cd	4	$1/8$	0.000391	0.0326
In–Xe	5	$1/12$	0.000035	0.0225
Cs–Yb	5	$1/32$	0.000205	0.0519
Lu–Hg	5	$1/12$	0.000239	0.0537
Tl–Rn	6	$1/12$	0.000069	0.0450
Fr–No	6	$1/32$	0.000297	0.0835

Equations

$$v_{\mathrm{a}}({\rm\bf r})=\sum\limits_{i=1}^{n}c_{i}\,\frac{\operatorname{\mathrm{erf}}\big{(}\sqrt{a_{i}}\,|{\rm\bf r}|\big{)}}{|{\rm\bf r}|}$$

$$\sum _{i=1}^n{c}_i=Q-1,$$

$$\mathcal{E}=\tfrac{1}{2}\iint\frac{\big{(}\rho_{v}({\rm\bf r}_{1})-\rho_{0}({\rm\bf r}_{1})\big{)}\big{(}\rho_{v}({\rm\bf r}_{2})-\rho_{0}({\rm\bf r}_{2})\big{)}}{\left|{\rm\bf r}_{1}-{\rm\bf r}_{2}\right|}\,\mathrm{d}^{3}{\rm\bf r}_{1}\,\mathrm{d}^{3}{\rm\bf r}_{2}$$

$$f_{\mathrm{E}}=\sum _i\left(⟨{\varphi }_i⟩{\hat{F}}^0{\varphi }_i-⟨{\varphi }_i^0⟩{\hat{F}}^0{\varphi }_i^0\right),$$

$$f_{\mathrm{S}}=\sum _j\left(1-\sum _i{\left|⟨{\varphi }_i⟩{\varphi }_j^0\right|}^2\right),$$

$$\bar{c}_{i}=\left\{\begin{array}[]{cl}0,&i<\bar{n},\\ \sum_{j=1}^{\bar{n}}c_{j}-Q_{0},&i=\bar{n},\\ c_{i},&i>\bar{n},\end{array}\right.$$

Full text

††journal: Theoretical Chemistry Accounts

1 Introduction

Every molecular electronic structure calculation needs a good starting guess. Most likely the earliest and still widely used extended Hückel H1963 theory gives a minimal-basis solution that can be then projected KSKWP1967 onto the working basis set. The underlying minimal basis should be given in some way — some workers SBBEGJKMNSWDM1993 used short three-Gaussians-per-shell sets HAKRST1984 for atoms up to Radon, with known shortcomings. We thought it would be handy to have a tool to build a better minimal set straight from the working set — clearly, atomic Hartree–Fock calculations could be done for that, but even then a starting guess is needed. We have found that the optimized effective potentials SH1953 ; TS1976 for atoms can be given, accurately enough, in a very compact form in terms of Gaussian charge distributions — they can be pre-tabulated and used to get the starting guess for an atom with one diagonalization of an easily-computable matrix. One may wonder whether this has already been done for all elements, but we see only works of this kind MSBG2011 with rather complicated, and even piecewise, analytic form of the potential.

Another well-known starting guess for molecules uses an effective one-electron Hamiltonian of density-functional theory KS1965 within the approximation H1985 of overlapped but frozen atomic densities. The nonlinearity of the exchange-correlation functionals makes the matrix elements of the potential tractable only by numerical integration, for which good schemes B1988 ; LKO2018 are known but need to be implemented. A linearized version thereof is simply the sum of fixed neutral-atom potentials taken from a density-functional calculation, which, remarkably, has been introduced and benchmarked for the first time only a few months ago Lehtola2019 . There is also a Hartree–Fock analog LZDG2006 that needs up to three-center two-electron integrals, or its spherical approximation AC2001 with less costly models of the direct and exchange terms.

Of course, the idea to approximate the molecular effective one-electron potential by a sum of atom-centered functions is not new and dates back to the early density-functional works SF1975 , where the Coulomb part was calculated using density fitting and the exchange-correlation potential was also fitted, either on the grid DCS1979 , or, in the special case of Hartree–Fock–Slater S1951 theory, by an elegant grid-free variational method D1986 with analytical evaluation of all three-center integrals. Though the approximation of both Coulomb and exchange-correlation terms by an electrostatic potential arising from a sum of atom-centered (Gaussian-type) effective charge distributions was put forward by one of us L1997 more than 20 years ago, there seems to be a new interest in using this functional form to derive transferable atomic potentials for molecules in educational work NW2017 and beyond W2019 ; W2019b .

Our atom-optimized compact representations can also be adopted and adapted for molecules by “capping” each atom with one more adjustable potential term to cancel the $-1/r$ asymptotic behavior and make it short-ranged, this way we get a new tunable model worth being studied. Our tests show it to be at least as good as an extended Hückel guess, we also compare it with the sum of tailless Hartree–Fock–Slater-based neutral-atom potentials and we see that our tunable model often works much better, so we report it here in the hope that it may be helpful for others.

2 Theory

For an atom, our effective potential $v({\rm\bf r})$ is the sum of its bare-nuclear potential $v_{\mathrm{n}}({\rm\bf r})$ and an approximate (direct and exchange) screening term

[TABLE]

in the form of Coulomb potential of Gaussian B1950 charge distributions. We optimize the exponents $a_{i}$ and coefficients $c_{i}$, with the constraint

[TABLE]

$Q$ being the nuclear charge, to minimize the measure

[TABLE]

of the difference between the Hartree–Fock $\rho_{0}$ and the effective-potential-based $\rho_{v}$ densities. Spherically-averaged configurations L2019b are taken for atoms with open-shell ground states. Eq. (2) ensures the right asymptotic behavior of $v({\rm\bf r})$ as $|{\rm\bf r}|\to\infty$. We use the four-component scalar-relativistic D1994 Hamiltonian and a finite Gaussian VD1997 nucleus model for all atoms Hydrogen through Nobelium.

Table 1 shows the number of terms $n$ and the value of $\mathcal{E}$ for the groups of atoms, the overall accuracy is quite good, $\mathcal{E}<0.0004$; $n-1$ terms would raise $\mathcal{E}$ by about 10 times, while $n+1$ can hardly do any better and would often lead to run-away solutions. We also check the atomic Hartree–Fock energy errors $\Delta E$ within a two-component scalar-relativistic approximation L2019a thanks to its boundedness from below.

For a molecule with $N$ atoms, the atomic potential tails of $-1/r$ would have summed up to an unphysical $-N/r$; as a lesser evil, we decided to “cap” each atom with one more term having $c_{n+1}=1$ and a to-be-optimized $a_{n+1}$. (Setting $c_{n+1}=(N-1)/N$ may be kept as an option.) This also gives us the freedom to adjust the potential well depth in the valence region to mimic a typical molecular environment.

For the optimization, we tried two objective functions: (a) the energy-based

[TABLE]

or (b) the overlap-based

[TABLE]

where the wavefunctions $\phi_{j}^{0}$ and the Fock operator $\hat{F}^{0}$ are taken from the reference self-consistent field calculation, whereas $\phi_{i}$ are our approximate solutions, the (spin and space) labels $i$ and $j$ run over the occupied valence set. The functions, either $f_{\mathrm{E}}$ or $f_{\mathrm{S}}$, are to be summed up over a training set of molecules and minimized with respect to the atomic parameters.

Our thus derived potentials can work equally well with any all-electron formalism, both two- and four-component. They can also be easily adapted for use with effective core potentials KBT76 by subtracting the core charge $Q_{0}$ from $c_{i}$ referring to the largest $a_{i}$,

[TABLE]

with $\bar{n}$ such that $0\leq\bar{c}_{i}\leq c_{i}$, where $a_{i}>a_{i+1}$ for all $i$.

3 Calculations

The problem of Eq. (3) was solved by a code written by one of us (DNL) within the atomic structure program L2019b , while the further optimization of parameters on molecules was done by an original program111To be found at https://github.com/briling/aepm. written by the other (KRB), and then the starting guess for molecules was double-checked by each of us independently by one’s own molecular electronic structure code.

The Hartree–Fock-optimized atomic parameters for Helium through Nobelium are tabulated in the Supplementary material. The reference molecular calculations were done with the PBE PBE1996 density functional, the scalar-relativistic approximation L2019a , and the L1 basis L2005 ; L2019b .

First, we optimized the $a_{n+1}$ parameters for atoms H through Cl (without He and Ne) on the neutral closed-shell subset of molecules L2011 . For H and B–F we got the values strikingly close to $1/3$ — this, together with the flatness near the minimum, led us to have all parameter values within groups of atoms constrained to be equal during the optimization and rounded afterwards to simple fractions. For the heavier elements, skipping Pr–Tm and Np–Md with their open-shell states, we have built up a smaller set of prototypical molecules and gotten the rest of values shown in Table 1. On average, we get the energy error per atom $f_{E}/N\approx 0.005$ au for the lighter, and $\approx 0.009$ au for the heavier, subsets.

Taking atomic Hartree–Fock–Slater densities as $\rho_{0}$ in Eq. (3) and setting the constraint with $Q$ on the right hand side of Eq. (2), we have optimized another set of potentials of Eq. (1) now with $n$ simply equal to the highest occupied principal quantum number in the atom, they are also tabulated in the Supplementary material for all 102 elements. Symptomatically, we get $a_{1}\approx 0.08$ for Hydrogen, rather diffuse compared to our optimized value of $1/3$. On our set of molecules, we get systematically higher error measures $f_{E}$ in most cases when the parametrized potentials are replaced with their Hartree–Fock–Slater-based analogs, typically by $\sim 10$ times for the organic subset (a full listing can be found in the Supplementary material), thus our parametrization is indeed helpful.

4 Conclusions

It would be straightforward to implement our starting guess, as an option, into any electronic structure code — for Gaussian basis sets, all the integrals can be easily computed analytically B1950 ; MD1978 . On the other hand, this elegant one-parameter-per-atom model may be of interest even on its own. Followed by an energy correction in the spirit of second-order perturbation theory HK1983 ; WP2003 ; DGG2009 ; DGG2010 , it may try to replace the iterative self-consistent solution in lower-accuracy applications.

An unnamed reader has warned us of the dangers of black boxes MBP2009 , and though we find our potentials to work well in the raised case of ferrocene, one should always keep in mind that a guess is a guess, however good it may be on average.

5 Supplementary material

File format description (readme.txt), atomic parameters based on Hartree–Fock (ac.txt) and Hartree–Fock–Slater (ac0.txt) theory, and error measures for all molecules (fm.txt).