Using the Pseudoatomic Solver to Generate NGWFs

Author:: Nicholas D.M. Hine, University of Warwick
Date:: September 2011
Date:: Updated February 2016

What is being calculated?

When the atomic solver is used [Ruiz-Serrano2012], a Kohn-Sham DFT calculation is performed for a “pseudoatom”. This means that the pseudopotential of a single isolated ion is used as the external potential, and the single-electron Kohn-Sham states are solved self-consistently, for a given XC functional. The resulting states can be expected to form an ‘ideal’ atomic orbital basis for a given calculation [Soler2002], [Artacho1999], [Blum2009], [Tarralba2008], [Chen2010]. In practice they are a good starting point for initial NGWFs, or a passable fixed basis for calculations without NGWF optimisation, in which context they should be at least comparable to the basis sets generated in SIESTA, for example (though note that a high cutoff energy is required for accurate representation on the grid in such cases).

The pseudoatomic orbitals we are looking for solve the Kohn-Sham equation:

\[\hat{H}_{\mathrm{KS}}\psi_{n}(\mathbf{r})=\epsilon_{n}\psi_{n}(\mathbf{r})\]

where

\[\hat{H}_{\mathrm{KS}}=-\frac{1}{2}\nabla^{2}+V_{\mathrm{loc}}(r)+\sum_{i}|p_{i}\rangle D_{ij}\langle p_{j}|\]

is the Hamiltonian in terms of kinetic, local potential and nonlocal potential contributions. For a norm conserving pseudopotential the matrix \(D_{ij}\) is diagonal and there is one term per angular momentum channel, so there is generally only one contributing nonlocal projector for each wavefunction. In PAW and USPs, there may be more, and nonzero cross-terms \(D_{ij}\) for \(i\neq j\). Because it is a spherical potential, solutions of this form of the Kohn-Sham equation are easily separable into an angular part, solved by the spherical harmonics (\(Y_{lm}(\hat{\mathbf{r}})\), or equivalently the real spherical harmonics \(S_{lm}(\hat{\mathbf{r}})\)) multiplied by radial part, which we will be solving explicitly in terms of a basis.

The atomic orbitals are solved in a sphere of radius \(R_{c}\). The most appropriate basis to use therefore consists of normalised spherical Bessel functions of given angular momentum \(l\). The radial parts of the basis functions, \(B_{l,\nu}(r)\) are

\[B_{l,\nu}(r)=j_{l}(q_{l,\nu}r)\:/\,\Big[\int_{0}^{R_{c}}|j_{l}(q_{l,\nu}r)|^{2}r^{2}dr\Big]^{\frac{1}{2}}\;,\]

with \(q_{l,\nu}\) such that \(q_{l,\nu}R_{c}\) are the zeros of the spherical Bessel functions. Thence \(B_{l,\nu}(R_{c})=0\) and \(\int_{0}^{R_{c}}|B_{l,\nu}(r)|^{2}r^{2}dr=1\) for all \(\nu\). The basis would be complete if \(\nu\) were infinite: in practice it must be truncated, and the number of functions included is determined by a kinetic energy cutoff \(E_{\mathrm{cut}}\). The criterion \(\frac{1}{2}q_{l,\nu}^{2}<E_{\mathrm{cut}}\) determines the largest \(\nu\) for each \(l\).

In the pseudoatom calculation we are therefore calculating Kohn-Sham states of the form

\[\psi_{n}(\mathbf{r})=\sum_{\nu}c_{n,\nu}\, B_{l_{n},\nu}(r)\, S_{l_{n}m_{n}}(\hat{\mathbf{r}})\;,\]

which have eigenvalues \(\epsilon_{n}\) and occupancies \(f_{n}\) which include spin-degeneracy. The occupancies are fixed, and determined before the main calculation, such that they sum to the number of valence electrons. Spherical symmetry of the density is assumed, so the occupancies of all members of a given set of \(m\)-degenerate orbitals are always equal — and in fact in practice the \(m\) states for a given \(l\), \(n\) pair are amalgamated into one state with \(f_{n}\) summed over all the degenerate \(m\)’s. For example, for a nitrogen ion with valence configuration \(2s^{2}\,2p^{3}\), we would have \(f_{2s}=2\), \(f_{2p}=3\). For this, we therefore need to find the lowest-energy self-consistent eigenstate of each of \(l=0\) and \(l=1\). Henceforth we will only consider the radial dependence \(\psi_{n}(r)\). All radial quantities will be considered to have been integrated over solid angle already, so factors of \(4\pi\) are omitted and \(\int|\psi_{n}(r)|^{2}r^{2}dr=1\) for a normalised orbital.

We define the local potential through

\[V_{\mathrm{loc}}(r)=V_{\mathrm{psloc}}(r)+V_{H}[n](r)+V_{XC}[n](r)+V_{\mathrm{conf}}(r)\]

where for a spherical charge distribution \(n(r)=\sum_{n}f_{n}|\psi_{n}(r)|^{2}\), the Hartree potential is given by

\[V_{H}(r)=\frac{1}{r}\int_{0}^{r}n(r')r'^{2}dr'+\int_{r}^{\infty}n(r')r'\, dr'\;.\]

and the XC potential is \(V_{XC}[n](r)=\frac{\partial E_{XC}[n]}{\partial n(r)}\). \(V_{\mathrm{conf}}(r)\) is an optional confining potential whose specific form will be discussed later [Blum2009].

For each \(l\) we can define the Hamiltonian matrix

\[H_{\nu,\nu'}^{l}=\int_{0}^{R_{c}}B_{l,\nu}(r)\left[\hat{H}B_{l,\nu'}(r)\right]r^{2}dr\]

and the overlap matrix

\[S_{\nu,\nu'}^{l}=\int_{0}^{R_{c}}B_{l,\nu}(r)B_{l,\nu'}(r)r^{2}dr\]

We then solve the secular equation

\[\mathbf{H}^{l}.\mathbf{c}_{n}=\epsilon_{n}\mathbf{S}^{l}.\mathbf{c}_{n}\]

to give the coefficients \(c_{n,\nu}\) which describe the orbitals. The orbitals are generated on the real-space grid and density mixing with a variable mixing parameter \(\alpha\) is then used until self-consistency is obtained. The result is deemed to be converged once a) the Harris-Foulkes estimate of the total energy (the bandstructure energy) matches the total energy as determined from the density to within a given tolerance (\(10^{-5}\) Ha) and the energy has stopped changing each iteration, to within a given tolerance (\(10^{-7}\) Ha).

Performing a Calculation with the Pseudoatomic Solver

The atomic solver is the default approach to NGWF initialisation, so if you do not need to change any settings for any species, simply omit the %block species_atomic_set block.

If there are any tweaks to be made to the default, this block is required, and for each element symbol the string “SOLVE” should appear in the entry. If you want to use automatic initialisation of the number of NGWFs, then specify that to be -1 in the species block. The code will attempt to determine how many orbitals to use, which orbitals constitute the valence, and what their default occupancies should be. To illustrate what will happen, we present some simple examples.

Let us imagine setting up a calculation with only nitrogen, for which \(Z_{\mathrm{atom}}=7\) and \(Z_{\mathrm{ion}}=5\). The valence manifold consists of \(4\) NGWFs of radius \(R_{c}=8.0\) per atom, so we would have the following blocks in our input file:

%block species

N N 7 4 8.0

%endblock species

%block species_atomic_set

N “SOLVE”

%endblock species_atomic_set

Note that because we may well want to add extra options to this string later, it’s best to always use the “” quotes around SOLVE. These settings will activate the pseudoatomic solver and it will attempt to guess a default configuration for the atom. Since \(Z_{\mathrm{ion}}=5\), the code will count back five electrons from the end of the default neutral atom occupancy, which is \(1s^{2}\,2s^{2}\,2p^{3}\), and will discover that the valence states are \(2s^{2}\,2p^{3}\). Since we have asked for \(N=4\) NGWFs, the solver will then count forward from the start of the valence states and determine that by including the whole first set of \(s\) and \(p\) states it has enough to span the valence space and create four orbitals (and thus four NGWFs). The solver will therefore solve for one state with \(l=0\), \(f=2\) and one state with \(l=1\), \(f=3\), all with radius \(R_{c}=8.0\), and from these states will produce one \(s\)-like NGWF and the three degenerate \(p_{x}\), \(p_{y}\) and \(p_{z}\) NGWFs.

A slightly more complex example would be if we were generating orbitals for iron (\(Z_{\mathrm{atom}}=26\), \(Z_{\mathrm{ion}}=8)\):

%block species

Fe Fe 26 9 10.0

%endblock species

%block species_atomic_set

Fe “SOLVE”

%endblock species_atomic_set

This time, to find the default configuration, the solver initialisation routines will count back 8 electrons from the neutral atom configuration of \(1s^{2}\,2s^{2}\,2p^{6}\,3s^{2}\,3p^{6}\,3d^{6}\,4s^{2}\) and thus will determine that the valence states are \(3d^{6}\,4s^{2}\). However, this time we have asked for 9 NGWFs, so it will then count forward from \(3d\), include the fivefold-degenerate lowermost \(d\)-like state and the lowest \(s\)-like state. This only makes six, so it then will also have to include the threefold-degenerate \(4p\)-like state. The solver will have to solve for one unoccupied \(p\)-like orbital, which will have \(f=0\) throughout the calculation.

Controlling the configuration

The default neutral-atom configurations for all the elements up to \(Z=92\) are included in the code, and will be used by default to generate the configuration. However, it is also possible to override these default configurations. For example, to generate NGWFs for iron in the 3+ state, we might want to set the occupancies to \(3d^{5}\,4s^{0}\). To do this we use the “conf=” directive after the SOLVE string:

%block species_atomic_set

Fe “SOLVE conf=3d5 4s0”

%endblock species_atomic_set

Any terms in the configuration which are not overridden are left at their default values. Another example might be if we wanted to force the partial occupation of more higher-lying states than would otherwise be occupied for the neutral atom:

%block species_atomic_set

C “SOLVE conf=2s1.5 2p2.5”

%endblock species_atomic_set

Note that the solver counts through the configuration terms strictly in the order \(n,l\), i.e. \(n\) is looped over outermost, then \(l=0\) to \(l=n-1\) for each \(n\) innermost. This means that sometimes a little thought is required to get the terms one actually wants, and not spurious extra ones. For example, if we wanted to run a calculation of oxygen with 9 NGWFs per atom, what we probably wanted would be to run with 1 \(s\)-like NGWF, 3 \(p\)-like NGWFs and 5 \(d\)-like NGWFs. However, this is not by default what one will get if one asks for

%block species

O O 8 9 9.0

%endblock species

%block species_atomic_set

O “SOLVE”

%endblock species_atomic_set

This will identify \(2s^{2}\,2p^{4}\) as the valence orbitals, and counting forward will identify \(2s\), \(2p\), \(3s\), \(3p\) and just 1 of the 5 degenerate \(3d\) states as the NGWFs required. Therefore, we must instruct the atomsolver to ignore the unwanted excited \(3s\) and \(3p\) terms. We do this with an “X”, which instructs the solver to knock out this term:

%block species_atomic_set

O “SOLVE conf=2s2 2p4 3sX 3pX 3d0”

%endblock species_atomic_set

Strictly speaking, the \(2s\), \(2p\) and \(3d\) strings are not needed, as they are the default values anyway, but they are left in for clarity. I find it advisable, so that I can keep track of the terms which will generate the NGWFs, to add explicitly the terms with zero occupancy to the conf string.

Generating larger, non-minimal bases

ONETEP is generally used to create an in-situ-optimised, minimal basis (eg 4 NGWFs/atom for C, N, O etc). However, it is also possible to fix the NGWFs and run with a much larger, unoptimised basis, in a manner akin to other DFT codes designed for large-scale simulations (eg SIESTA). One would then normally want to use multiple NGWFs for each angular momentum channel corresponding to the valence orbitals. This is known as using a “multiple-zeta” basis set, where zeta refers to the radial part of the valence atomic orbitals. For example, a “triple-zeta” basis for carbon would have \(3\) \(s\)-like functions and \(3\) of each of \(p_{x}\), \(p_{y}\), and \(p_{z}\)-like functions. There are two approaches to generating these extra radial functions. This simplest is just generate the higher-lying orbitals of a given angular momentum. For carbon, for example, a double-zeta basis in this scheme would include \(3s\) and \(3p\)-like states. This approach, however, is not very quick to converge with basis size.

It is often better to apply the commonly-used “split-valence” approach. This allows the orbitals that have been generated to be “split” into multiple functions, so as to generate so-called “split-valence multiple-zeta” basis sets. In this formalism, one function \(f(r)\) can be split into two functions \(g_{1}(r)\) and \(g_{2}(r)\) according to the following:

A matching radius \(r_{m}\) is chosen: for \(r>r_{m}\), we set \(g_{2}(r)=f(r)\). For \(r\leq r_{m}\), we set \(g_{2}(r)=r^{l}(a_{l}-b_{l}r^{2})\), where \(a_{l}\) and \(b_{l}\) are chosen such that \(g_{2}(r_{m})=f(r_{m})\) and \(g_{2}'(r_{m})=f'(r_{m})\).
The other function, \(g_{1}(r)\), is set to \(f(r)-g_{2}(r)\), so \(g_{1}(r)=0\) for \(r\geq r_{m}\).
Both functions are renormalised, so \(\int_{0}^{R_{c}}|g_{1}(r)|^{2}r^{2}dr=1\) and \(\int_{0}^{R_{c}}|g_{2}(r)|^{2}r^{2}dr=1\).

Splitting of a term is activated by adding a colon after the term and specifying the “split norm” value. This is the fraction \(p\) of the total norm of the orbital which is beyond the matching radius \(r_{m}\), such that \(\int_{r_{m}}^{R_{c}}|f(r)|^{2}r^{2}dr=p\). If this colon is present, the solver will take into account the total number of orbitals which will result from this term after splitting, when counting forward in the configuration terms to determine which orbitals to solve. For example, if we wished to generate a Double-Zeta Polarisation (DZP) basis for oxygen (\(2\times1\times s\), \(2\times3\times p\),\(1\times5\times d\)), where the last 15% of the norm was matched for the \(s\) and \(p\)-orbitals, we would use the following:

%block species

O O 8 13 9.0

%endblock species

%block species_atomic_set

O “SOLVE conf=2s2:0.15 2p4:0.15 3sX 3pX 3d0”

%endblock species_atomic_set

So that you can tell that it is happening, the code will output a message along the following lines when splitting a given orbital.

Splitting orbital 1, splitnorm= 0.150000000

Splitting orbital 1, splitnorm= 0.060000000

The result of the splitting can be viewed in the “initial_rad_ngwf_xx” files.

Obtaining Polarisation Orbitals through Perturbation

As well as including more flexibility for the valence orbitals, in the form of multiple-zeta basis sets, one frequently also wants to expand the basis by extending it to higher angular momentum channels. This can be done by simply increasing the number of NGWFs requested and ensuring through the ’conf=’ string that the extra functions obtained are of the right angular momentum. However, the resulting high-\(l\) states tend to be unbound in the free atom, and therefore do not necessarily add anything particularly useful to the basis.

There is an alternative means to generate higher-\(l\) states, using perturbation theory. In this, one effectively applies an electric field to the valence states of angular momentum \(l\) and polarises them, resulting in a set of states of angular momentum \(l+1\). Imagine we have an orbital \(\psi_{0}(\mathbf{r})\) of angular momentum \(l\), \(m\) which is an eigenstates of the original Hamiltonian \(\hat{H}_{0}\) with eigenvalue \(\epsilon_{0}\):

\[\psi_{0}(\mathbf{r})=\sum_{\nu}c_{0,\nu}\, B_{l,\nu}(r)\, S_{lm}(\hat{\mathbf{r}})\;.\]

We wish to polarise this orbital by applying an electric field \(\mathcal{E}\) in the \(z\)-direction:

\[\hat{H}_{1}=\mathcal{E}r\, S_{10}(\hat{\mathbf{r}})\;,\]

(since \(S_{10}(\hat{\mathbf{r}})=z/r\)). Perturbing \(\psi_{0}\) with \(\hat{H}_{1}\) gives us no shift in energy to first-order, since the perturbation is an odd multiplicative function of \(z\), meaning \(\epsilon_{1}=0\). What about the change in the wavefunction? This obeys

\[(\hat{H}_{0}-\epsilon_{0})\psi_{1}(\mathbf{r})=-(\hat{H}_{1}-\epsilon_{1})\psi_{0}(\mathbf{r})\label{eq:pert}\]

In principle, \(\psi_{1}(\mathbf{r})\) could have any angular momentum components, but we can see that in practice it only contains \(L=l\pm1\), since the dipole selection rule excludes all other channels. We already have \(l-1\) states in our basis, so we conclude that \(\psi_{1}(\mathbf{r})\) need only include \(l+1\), and we can expand \(\psi_{1}\) in terms of the \(l+1\) basis functions:

\[\psi_{1}(\mathbf{r})=\sum_{\nu}c_{1,\nu}\, B_{l+1,\nu}\, S_{l+1,m}(\hat{\mathbf{r}})\]

Therefore we can generate a shifted Hamiltonian

(1)\[H_{\nu,\nu'}^{l+1}=\int_{0}^{R_{c}}B_{l+1,\nu}(r)\left[(\hat{H}^{l+1}-\epsilon_{0})B_{l+1,\nu'}(r)\right]r^{2}dr\;,\]

and the components of the RHS of Eq. (1)

\[D_{\nu}=-\int_{0}^{R_{c}}B_{l+1,\nu}(r)r\psi_{0}(r)r^{2}dr\;.\]

To solve for \(c_{1,\nu}\) we just need to invert \(H_{\nu,\nu'}^{l+1}\)and apply it to \(D_{\nu}\), and then renormalise the result to have a norm of 1.

In practice, polarisation of a given configuration term of angular momentum \(l\), to form a perturbative polarisation orbital for \(l+1\), is achieved by adding “|P” to the term, for example:

%block species

O O 8 13 9.0

%endblock species

%block species_atomic_set

O “SOLVE conf=2s2:0.15 2p4:0.15|P”

%endblock species_atomic_set

So that you can tell that it is happening, the code will output a message along the following lines when polarising a given orbital:

Polarising orbital 1 to generate l= 2 function (orbital 3)

Again, the result can be viewed by plotting the relevant “initial_rad_ngwf_xx” file.

Overriding radii

By default, the cutoff radius used for all the orbitals of an atom is the same \(R_{c}\) as defined in the %block species entry for that element. However, we can override this, either for all orbitals, or for certain angular momentum channels.

To override the radius for all channels, for example to 7.0\(a_{0}\), would we add the flag “R=7.0” to the SOLVE string:

%block species_atomic_set

O “SOLVE conf=2s2:0.15 2p4:0.15 3sX 3pX 3d0 R=7.0”

%endblock species_atomic_set

Or leave the default values for all other channels, but override the \(d\)-channel only to 5.0\(a_{0}\), we would use

%block species_atomic_set

O “SOLVE conf=2s2:0.15 2p4:0.15 3sX 3pX 3d0 Rd=5.0”

%endblock species_atomic_set

Adjusting confining potentials

By default, a confining potential is applied, of the form:

\[V_{\mathrm{conf}}(r)=S\,\exp[-w_{l}/(r-R_{c}+w_{l})]/(r-R_{c})^2\]

where \(S\) is the maximum height of the confining potential (at \(r=R_{c}\)), and \(w_{l}\) is the width of the region over which it is applied. By default, \(S=100\) Ha, and\(w_{l}=3.0a_{0}\) for all \(l\)-channels used. These can also be overridden, either all at once or for specific \(l\)-values in the case of \(w\).

For example, to set no confining potential on the confined \(d\)-orbitals in Zinc, but to keep the default one on all the other orbitals, we could set \(w_{d}=0\):

%block species_atomic_set

Zn “SOLVE conf=3d10 4s2 Rd=5.0 wd=0.0”

%endblock species_atomic_set

Or to turn off confinement potentials entirely, and generate \(R_{c}=15.0a_{0}\) orbitals to match those generated by CASTEP’s atomsolver (this should allow direct comparison of energies, given suitable tweaking of the energy cutoffs so that they exactly match:

%block species_atomic_set

O “SOLVE R=15.0 S=0.0”

%endblock species_atomic_set

Note that there can be problems with convergence for certain choices of confining potential. In particular, if you apply different confining potentials to different occupied orbitals, there will be problems obtaining agreement between the Harris-Foulkes estimator and the actual total energy - because the band energy will incorporate the different confining potentials, but the total energy cannot. The confining potential on angular momentum channels with no occupied orbitals can therefore be whatever you like, but those of occupied orbitals must all match. The exception to this is if the cutoff radius ment of one channel is less than the onset radius for the others. In this case, there is no need to apply a confinement to the lower-cutoff channel at all (eg in the example above for Zinc).

Initial Guess Density: Setting initial charges and spins

The atomsolver solutions are by default also used to provide an initial guess for the valence electron density. This is used to generate the initial Hamiltonian, which determines the occupation of the orbitals via Palser-Manolopoulos or other kernel optimisation schemes. Therefore it is important that this initial density is a reasonably good guess to the real density.

A superposition of atomic densities is usually fine for neutral systems without large magnetic moments. Sometimes, however, one needs to adjust the charges and spins of this initial density. It appears to be a rather bad idea to actually adjust the orbital occupations of the pseudoatoms self-consistently: it becomes impossible to disentangle the effect of changing the orbital from that of changing the density.

A better approach, therefore, is to retain the same pseudoatomic solutions for the neutral atom, but adjust the orbital occupations only at the point where they are used to generate the density.

This can be done by specifying “INIT CHARGE=X SPIN=Y” in the SOLVE string. Either CHARGE or SPIN can be omitted if they are not required. For example, for a manganese ion with charge +3 and spin 4, we might set

%block species_atomic_set

Mn “SOLVE conf=3d5 4s2 wd=7.0 INIT SPIN=4 CHARGE=+3”

%endblock species_atomic_set

The default occupation for the neutral atom is \(3d^5\) \(4s^2\). However, we ask it to apply a charge of +3, and this will remove occupation number from the orbitals with the highest energy until the right charge is obtained. In this case the resulting occupation will be \(3d^4\) \(4s^0\). Next, the spin is applied, resulting in \(3d_{\uparrow}^4\) \(3d_{\downarrow}^0\). Note that the charge is applied first, followed by the spin.

[Ruiz-Serrano2012] A. Ruiz-Serrano, N.D.M. Hine and C.-K. Skylaris, in press, (2012).

[Soler2002] J.M. Soler, E. Artacho, J.D. Gale, A. Garcia, J. Junquera, P. Ordejon, and D. Sanchez-Portal,The SIESTA method for ab initio order-N materials simulation, J. Phys. Condens. Matter 14, (2002).

[Artacho1999] E. Artacho, D. Sanchez-Portal, P. Ordejon, A. Garca, and J. M. Soler, Linear-scaling ab-initio calculations for large and complex systems, Phys. Status Solidi B 215, 809 (1999).

[Blum2009] V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter, and M. Scheffler: Ab initio molecular simulations with numeric atom-centered orbitals, Comput. Phys. Commun. 180, 2175 (2009).

[Tarralba2008] A. S. Torralba, M. Todorovic, V. Brazdova, R. Choudhury, T. Miyazaki, M. J. Gillan, and D. R. Bowler. Pseudo-atomic orbitals as basis sets for the O(N) DFT code CONQUEST, J. Phys. Condens. Matt. 20(29), (2008).

[Chen2010] M. Chen, G.-C. Guo, and L. He, Systematically improvable optimized atomic basis sets for ab initio calculations, J. Phys. Condens. Matt. 22, 445501, (2010).