# Optimization of Gaussian basis sets within the Hartree-Fock Method

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?).

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.

• Thanks for the detailed response. One thing I am confused about, how does any of this update $\alpha$? It looks like this only works by modifying $C$. Also, as a side note, it looks like your derivative is a gradient, is that correct? I have not seen such notation using $\delta$ before.
– tmph
Nov 12 '20 at 7:10
• Sorry that wasn't clear, nothing about this procedure is unique to C -- you simply compute dE/dalpha as well and change alpha in the same way as C. In general you need a functional derivative, for which the conventional symbol is delta. Nov 12 '20 at 15:48
• I've rewritten it in terms of phi, to try to clarify it. Is that better? Would you like d/dC and d/dalpha explicitly? Nov 12 '20 at 16:02
• Also updated now to include explicitly how you would treat the C and alpha changes. Is that clearer? Nov 12 '20 at 17:33
• @tmph OK, I've reworked it to be more explicit. The short answer is that in steepest descent you do a single 1D optimisation which changes both $C$ and $\alpha$ simultaneously - there are better methods! I've also tried to briefly explain what co- and contravariance mean here, but in essence $(H-ES)\phi$ transforms like $S\phi$ if we change the basis, and if we're going to use this as a direction to add to $\phi$ it needs to transform in the same way as $\phi$. Nov 13 '20 at 2:01

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).

• Interesting answer; I shall definitely visit some of these papers. I am (re-)studying Hartree-Fock in the interest of developing post-Hartree-Fock methods so I find your fifth paragraph interesting. I'm not familiar with tempering, but by what you are describing it sounds almost like to avoid unreliable nonlinear optimisation, you simply "overfit" with a redundant amount of exponentials. Is that right? And to make sure I understand, when we use, say, GTO basis sets, we do not incorporate the $\alpha$s, but just take them from atomic calculations, which I suppose is where Thijssen got them.
– tmph
Nov 13 '20 at 17:07
• @tmph not really; it's just that to get a reasonable approximation you have to have many exponents, and even if you optimize them fully, the exponents will form a quasi-linear sequence in $\log \alpha$. If you just want to do method development, then you should just use established basis sets. GTO basis set development is an extremely well-established field, and published sequences like the Dunning-style {,aug-}cc-p{,C,wC}VXZ series are sufficient for all applications, especially if you stick to the light elements. If you want to go heavy or very accurate, then they might not suffice. Nov 13 '20 at 19:33