Optimization of Gaussian basis sets within the Hartree-Fock Method

Question

I am revisiting some exercises in Thijssen's Computational Physics book, particularly chapter 4 on the Hartree-Fock method. I am interested in the method of nonlinear optimisation for its own purposes, but if there is a better way to approach the problem in this context, I'd be interested in that as well. I know this is discussed in Numerical Recipes (Thijssen cites it), but either I do not understand or am looking at the wrong section of the book.

The problem, in brief, is to calculate the ground state wave function of the helium atom. He does this by expanding the wave function in 4 gaussian basis functions, $$ \phi({\bf r}) =\sum_{p=1}^4 C_p e^{-\alpha_p r^2}, $$ but he provides all $\alpha_p$s so that we have a linear variational problem, which I've already successfully solved (in Fortran).

However, I am wondering how he obtains these $\alpha_p$s and how I would approach such problems more generally. It seems to me that the problem gets much more messy, as we have 2 variables per basis function and of course the vector ${\bf C}$ is not a sufficient representation. The $Q$ tensor that he defines, $$ Q_{prqs} = \int d^3r_1d^3r_2\chi_p(r_1)\chi_r(r_2)\frac 1{|r_1-r_2|}\chi_q(r_1)\chi(r_2) = \frac{2\pi^{5/2}}{(\alpha_p+\alpha_q)(\alpha_r+\alpha_s)\sqrt{\alpha_p+\alpha_q+\alpha_r+\alpha_s}} $$ where the $\chi$s are gaussians, also becomes much more complicated since it contains non-constant $\alpha$s now. More generally, it seems like the problem is $$ \min_{{\bf C}, {\bf \alpha}} 2{\bf C}^Th{\bf C}+\sum_{pqrs}Q_{prqs}C_pC_qC_rC_s $$ where $h$ is the independent-particle Hamiltonian, all subject to a highly nonlinear constraint $H{\bf C}=E'S{\bf C}$ where $H$ contains these messy $\alpha$s. There doesn't seem to be any easy way to solve this, so I was wondering where to go from here (do we simplify further, use a dual problem, or is there some kind of really powerful convex solver that I am unaware of?).

Answer 1

There are a lot of different methods, but a good starting point to understand how this might be done is the method of steepest descent. First we take the generalised eigenvalue expression you cite and rearrange to express the energy expectation value as

$$E=\frac{\phi^\dagger \mathrm{H}\phi}{\phi^\dagger \mathrm{S}\phi},$$

where I've generalised to allow complex wavefunctions for completeness, so the Hermitian conjugate (denoted by $\dagger$) has replaced the transpose in your expression.

The variational principle tells us that the ground-state energy we seek is the minimum possible value of this expectation value. If we start with some trial $\phi$ and compute the energy expectation value, we may also ask how this expectation value would change if we made a small change to $\phi$. Differentiating the expectation value with respect to $\phi$ using the quotient rule (where the matrix-derivative here may be interpreted as being element-wise) we obtain:

$$ \begin{array}{cl} \frac{\delta E}{\delta \phi^\dagger} &= \frac{\left(\phi^\dagger \mathrm{S}\phi\right)\mathrm{H}\phi - \left(\phi^\dagger \mathrm{H}\phi\right)\mathrm{S}\phi}{\left(\phi^\dagger \mathrm{S}\phi\right)^2}\\ &= \frac{\mathrm{H}\phi - E\mathrm{S}\phi}{\phi^\dagger \mathrm{S}\phi}. \end{array} $$ If we assume that $\phi$ is S-orthonormalised, such that $\phi^\dagger \mathrm{S}\phi=I$ (where $I$ is the identity) then we can write this more compactly as $$ \begin{array}{cl} \frac{\delta E}{\delta \phi^\dagger} &= \mathrm{H}\phi - E\mathrm{S}\phi. \end{array} $$

This expression tells us the direction in which to change $\phi^\dagger$ in order to increase $E$ most rapidly; since we wish to decrease $E$, we should change $\phi^\dagger$ in the opposite direction.

Now at this point I should note that we don't actually want to change $\phi^\dagger$, we want to change $\phi$ and there's a complication here which arises from the fact that this is a generalised eigenvalue problem. The gradient I've derived is covariant (it transforms like $S\phi$ under a change of basis), and we need a contravariant quantity (transforms like $\phi$ under a change of basis) if we're going to use it to update $\phi$. There are several avenues for getting this, but probably the simplest and most computationally efficient is to pre-multiply by $S^{-1}$.

We may now define a direction,

$$D = -S^{-1}\frac{\delta E}{\delta \phi^\dagger} = -S^{-1}\mathrm{H}\phi - E\phi,$$

in your example above, there are two ways in which we may change $\phi$: we can change $C$; and we can change $\alpha$. Our search direction $D$ contains changes to both,

$$ \begin{array}{rl} D = \left(\begin{array}{c} D_C\\ D_\alpha \end{array}\right) &= -\left(\begin{array}{c} \frac{\delta E}{\delta C^\dagger}\\ \frac{\delta E}{\delta \alpha} \end{array}\right)\\ &= -\frac{\delta E}{\delta \phi^\dagger}\left(\begin{array}{c} \frac{\delta \phi^\dagger}{\delta C^\dagger}\\ \frac{\delta \phi^\dagger}{\delta \alpha}\\ \end{array}\right)\\ &= -(H-ES)\phi\left(\begin{array}{c} \frac{\delta \phi^\dagger}{\delta C^\dagger}\\ \frac{\delta \phi^\dagger}{\delta \alpha} \end{array}\right) \end{array} $$

We may now search along the direction $D$ for an improved estimate of the ground state $\phi$. The simplest approach is to use a line-search method, whereby we define

$$\phi^\prime (\lambda) = \phi + \lambda D$$

by which we mean

$$ \left(\begin{array}{c} C^\prime(\lambda)\\ \alpha^\prime(\lambda) \end{array}\right) = \left(\begin{array}{c} C\\ \alpha \end{array}\right) +\lambda \left(\begin{array}{c} D_C\\ D_\alpha \end{array}\right) $$

Now we find the minimum expectation value of the energy with respect to $\lambda$ (a 1D minimisation problem, so fairly straightforward). Note that the choice of function of $\lambda$ above means that $\phi^\prime$ is not orthonormalised in general, even if $\phi$ and $D$ are, which necessitates re-orthonormalisation later on. More sophisticated methods are available which do not have this side-effect.

Once we have the value of $\lambda$ which minimises the energy expectation value, $\lambda_{opt}$, we have an improved estimate of the ground state $\phi^{new}=\phi(\lambda_{opt})$. We may now compute the energy derivative about $\phi^{new}$, and hence a new search direction, and the procedure repeats.

There is a complication in that the sensitivity to changes in $C$ and $\alpha$ will be very different, leading to ill-conditioning; this can be resolved by "preconditioning" the search direction (see later on).

This iterative update procedure is usually terminated when the change in solution falls below some tolerance criterion -- the change usually being measured by the change in the energy expectation value or the magnitude of the gradient.

This method is relatively crude and only performs well when the expectation values associated with the underlying basis states are similar; if this is not satisfied, then the problem becomes ill-conditioned and can take many iterations to converge (it may even diverge).

The conditioning is dependent on the curvature of the energy expectation with respect to $\phi^\dagger$, and if we knew what that was we could compute an improved search direction (what you do is pre-multiply the gradient by the inverse of the curvature tensor). Unfortunately in electronic structure calculations we don't usually know what the curvature is analytically, and it's very expensive to compute. However, two practical improvements are:

• Preconditioning -- premultiply the gradient by an approximation to the inverse curvature; in the case of electronic structure calculations this would often be based on parts of the Hamiltonian which dominate in particular regimes, for example the work of Teter, Payne and Allan (https://doi.org/10.1103/PhysRevB.40.12255)
• Quasi-Newton method -- essentially, use gradient information from previous iterations to estimate the curvature and produce an improved search direction; examples include the popular conjugate gradients and (L-)BFGS methods.

These two approaches may be combined, and preconditioned quasi-Newton methods can be extremely efficient. Quasi-Newton methods are fairly standard across optimisation problems and you'll see the same methods cropping up time and time again; in contrast, finding a good preconditioning approach will require careful thought about the behaviour of your particular problem, so different approaches are employed in different domains.

Answer 2

Short answer: don't do it.

While optimizing the exponents in a basis set for your system is possible in principle, the optimization is very expensive. You may also end up ruining the accuracy of your properties, since the basis set may become biased.

Basis sets of various quality are routinely available in the literature. Instead of reoptimizing the primitive basis, it is much cheaper to just use a bigger basis set i.e. more exponents from the start.

How are the basis sets in the literature formed? Often, the construction of basis sets starts from atomic calculations. At a later stage, the sets are tested in molecular calculations, and higher-angular-momentum basis functions are added to describe polarization effects.

Phil's answer above explains how to do multidimensional optimization in general. However, since solving the self-consistent field problem with respect to the orbital coefficients ${\bf C}$ is often quite tricky, and since you might also be interested in running calculations at post-Hartree-Fock levels of theory, usually the exponent optimization is just taken out of the electronic structure calculation. Namely, for each set of exponents $\boldsymbol{\alpha}$ you can solve for the electronic wave function $\Psi$ and the resulting energy $E$. You can determine the gradient $\nabla_{\boldsymbol{\alpha}}E$ either analytically from the energy expression, or by finite differences. The multidimensional optimization for $\boldsymbol{\alpha}$ is often realized with the L-BFGS method.

Full optimization is, however, often not required. Once you have enough exponents, the ones in the middle start to follow a linear sequence in $\log \alpha$. Even-tempering is often a very good starting point for basis sets; probably the majority of basis sets that have been published use even-tempering.

In the even tempering scheme, you define your sets in terms of just two parameters $\alpha$ and $\beta>1$ as $\alpha_i = \alpha_0 \beta^i$ for $i\in[0,N-1]$. In the limit of a complete basis, $N\to\infty$, one has $\alpha_0 \to -\infty$ and $\beta\to1$. Now instead of optimizing $N$ exponents, it suffices to optimize just the two parameters $\alpha_0$ and $\beta$.

If you want to be more accurate, you can throw in $k\in[1,N]$ more parameters. The exponent optimization is well-behaved if you express them in a basis of Legendre polynomials, see J. Chem. Phys. 118, 1101 (2003). Note, however, that the possibility of landing on a local minimum instead of the global one increases as a function of the number of parameters!

PS. I have recently suggested an even simpler scheme to forming accurate basis sets from a simple physical principle, requiring no electronic structure calculations at all, see J. Chem. Phys. 152, 134108 (2020) (arXiv:2001.04224).