Basis sets

As we have mentioned in finding the ground state, the density matrix is represented by support functions. These, in turn, are made up of basis functions, and the choice of basis set and how it represents the support functions affects both the accuracy and performance of a calculation.

There are two kinds of basis functions which are used in CONQUEST to represent the support functions: pseudo-atomic orbitals (keyword PAOs) which are the default; and b-splines (keyword blips) which allow systematic convergence at the expense of greater complexity.

The basis set is selected as follows:

Basis.BasisSet PAOs

When the basis set is taken to be PAOs, there are three different ways to construct the support functions, discussed below:

Each support function is represented by a single PAO (primitive PAOs)
Multi-site support functions, built from PAOs on several atoms
on-site support functions, built from PAOs on one atom

Primitive PAOs are efficient for small systems. When using large PAO basis sets for systems containing more than several hundred atoms, multi-site support functions and on-site support functions will be more efficient than primitive PAOs.

Go to top.

Pseudo-atomic orbitals (PAOs)

PAOs are solutions of the Schrodinger equation for isolated atoms, using pseudopotentials, with some confinement applied. They consist of radial functions multiplied by spherical harmonics. For the valence orbitals, the radial functions are referred to as zeta (\(\zeta\)), while for unoccupied orbitals, they are termed polarisation.

A minimal PAO basis set is single-\(\zeta\) (SZ), with one radial function for each angular momentum quantum number in the valence electrons. While the cost is significantly lower than for other basis sets, the accuracy will be rather low.

The accuracy of a calculation can be improved by adding polarisation functions and multiple radial functions for different angular momentum values, though systematic improvement is rather difficult to achieve (this is straightforward with a blip function basis). The PAO utility included with CONQUEST generates basis sets with differing sizes and accuracies; full details of the performance of these basis sets can be found elsewhere [B1].

minimal (single zeta, SZ)
small (single zeta and polarisation, SZP)
medium (double zeta, single polarisation, DZP)
large (triple zeta, double polarisation, TZDP)

Go to top.

Primitive PAOs as support functions

The easiest way to prepare support functions is to use primitive PAOs as the support functions without any modifications. In this case, the input parameters related to the support functions are automatically set by obtaining the information from the PAO files (.ion files) so long as they are generated by the CONQUEST MakeIonFiles utility, version 1.0.3 or later. No further input parameters need to be set in Conquest_input.

Go to top.

Multi-site support functions

Since the computational cost of Conquest scales cubically with the number of support functions, contracting the PAOs into a smaller set of support functions is an efficient way to reduce the computational cost when we use large multiple-\(\zeta\) PAO basis sets. Multi-site support functions (MSSFs) [B2, B3] are constructed for each atom by taking linear combinations of the atom’s PAOs and the PAOs from neighbouring atoms within a certain range (set with the parameter Atom.MultisiteRange in the atom specification block).

Multi-site support functions can be selected by setting the following parameters:

Basis.BasisSet PAOs
Basis.MultisiteSF T

Various other parameters need to be set in the atom specification block. The number of support functions for the atoms must be set, and is normally equivalent to a minimal (single zeta) basis; it is set with Atom.NumberOfSupports. (To use a number of support functions larger than this minimal number, the parameter Multisite.nonminimal needs to be set to T.) The range for the multi-site support functions (the PAOs of any atom within this distance of the atom will be included in the support functions) is set with Atom.MultisiteRange. The accuracy of the MSSF will improve as this range is increased, though the computational cost will also increase; careful tests must be made to find an appropriate range. For a minimal number of MSSF, the range must be large enough to include other atoms, though this restriction can be removed (see on-site support functions for more details).

As well as setting the range for the MSSFs, we need to specify an approach for finding the expansion coefficients. A reasonable set of MSSF coefficients can be found using the local filter diagonalization (LFD) method. For improved accuracy, this should be followed by variational numerical optimisation.

Go to top.

Local filter diagonalization (LFD)

In this method, which is selected by setting Multisite.LFD T, the MSSF coefficients are found by diagonalising the Hamiltonian in the primitive PAO basis, for a small cluster of atoms surrounding the target atom. The MSSF coefficients \(C\) are determined by projecting the sub-space molecular orbitals \(C_{sub}\) around each atom onto localized trial vectors \(t\),

\(C = C_{sub} f(\varepsilon_{sub}) C_{sub}^T S_{sub} t\)

The cluster for diagonalisation must be at least as large as the MSSF range, but larger clusters tend to give better MSSF coefficients (at the expense of an increased computational cost). The LFD sub-space region is determined for each atom by setting Atom.LFDRange.

An example set of parameters for an MSSF calculation for bulk Si would be:

Basis.BasisSet PAOs
Basis.MultisiteSF T
Multisite.LFD T

%block ChemicalSpeciesLabel
1 28.07 Si
%endblock

%block Si
Atom.NumberOfSupports 4
Atom.MultisiteRange 8.0
Atom.LFDRange 8.0
%endblock

When calculating binding energy curves or optimising cells, a change of lattice constant can suddenly bring a new set of atoms within the range of the support functions. In this case, a smearing can be applied at the edges of the range, by setting Multisite.Smear T. Further details are given below.

Some form of self-consistency between the MSSF and the charge density is required (as the MSSF will determine the Hamiltonian and hence the output charge density). At present, this is performed as a complete SCF cycle for each set of MSSF coefficients (though this is likely to be updated soon for improved efficiency). This is selected by default (but can be turned off by setting the parameter Multisite.LFD.NonSCF T).

This iterative process is not variational, but is terminated when the absolute energy change between iterations is less than Multisite.LFD.Min.ThreshE, or the residual (defined in self-consistency) is less than Multisite.LFD.Min.ThreshD.

An example input block for this process would be as follows:

Multisite.LFD T
Multisite.LFD.Min.ThreshE 1.0e-6
Multisite.LFD.Min.ThreshD 1.0e-6

Go to top.

Numerical optimisation

The MSSF coefficients can also be optimised by minimizing the DFT energy with respect to the coefficients, in a variational process. The threshold and the maximum iteration number of the numerical optimisation are specified by minE.EnergyTolerance and minE.SupportVariations. The optimisation is based on the conjugate gradient (CG) method, and the initial CG step size can be specified by minE.InitStep_paomin (default is 5.0).

minE.VaryBasis T
minE.EnergyTolerance 1.0e-6
minE.SupportVariations 30

The numerical optimisation provides more accurate coefficients than the LFD method but is usually more time consuming. Therefore, it is generally better to start from good initial values, for example, the coefficients calculated by LFD. When both Multisite.LFD and minE.VaryBasis are selected, the initial coefficients will be calculated by LFD and the coefficients will then be optimised.

Basis.MultisiteSF T
Multisite.LFD T
minE.VaryBasis T

If good initial coefficient values have been found in a previous calculation, reading these from files (the base name of these files is SFcoeffmatrix2) and performing only the numerical optimisation is also a good choice.

Basis.LoadCoeffs T
Basis.MultisiteSF T
Multisite.LFD F
minE.VaryBasis T

Go to top.

Advanced MSSF concepts

Smearing the edge of the support functions Here, we are concerned with changes of lattice constant which may bring new atoms inside the support function range.

We can set the smearing-function type Multisite.Smear.FunctionType (default=1:Fermi-Dirac, 2=Error function), the center position of the function Multisite.Smear.Center (default is equal to the range of the support functions), offset of the center position Multisite.Smear.Shift and the width of the Fermi-Dirac function Multisite.Smear.Width (default=0.1).

Selecting states from the sub-space Here, we consider how to create the MSSF themselves from the results of the sub-space diagonalisation.

The Fermi function \(f\) with \(\varepsilon_{sub}\) Multisite.LFD.ChemP and \(kT\) Multisite.LFD.kT in the equation removes the effects of the subspace molecular orbitals in higher energy region. In default, \(\varepsilon_{sub}\) is automatically set to the mean value of the subspace HOMO and LUMO energies for each subspace. If users want to modify this, set Multisite.LFD.UseChemPsub F and the \(\varepsilon_{sub}\) value with Multisite.LFD.ChemP.

For the LFD trial functions \(t\), when Atom.NumberOfSupports is equal to the number of SZ or single-zeta plus polarization (SZP), the PAOs which have the widest radial functions for each spherical harmonic function are chosen as the trial vectors automatically in default. When Atom.NumberOfSupports is equal to the number of SZP and Multisite.nonminimal.offset is set, the other PAOs will have the weight in the trial vectors with the value of Multisite.nonminimal.offset. The users can also provide the trial vectors from the input file using the LFDTrialVector block

# Trial vectors of Au (element 1) and O (element 2) atoms.
# Au: 15 PAOs (DZP) -> 6 support functions, O: 13 PAOs (DZP) -> 4 support functions.
%block LFDTrialVector
# species sf npao   s   s   x   y   z  d1  d2  d3  d4  d5  d1  d2  d3  d4  d5 for Au
        1  1   15 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        1  2   15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
        1  3   15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
        1  4   15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
        1  5   15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
        1  6   15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
        2  1   13 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        2  2   13 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        2  3   13 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
        2  4   13 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
# species sf npao   s   s   x   y   z   x   y   z  d1  d2  d3  d4  d5 for O
%endblock LFDTrialVector

The first, second and third columns correspond to the indices of species, support functions for each species, and the number of PAOs for each species. The other columns provide the initial values of the trial vectors. For example, in the first line in the above example, the second s PAO is chosen as the trial vector for the first support function of Au.

Self-consistent LFD Two further conditions are applied to end the LFD self-consistency process. The maximum number of iterations is set with Multisite.LFD.Min.MaxIteration. It is also possible, as the process is not variational, that the energy can increase as well as decrease between iterations. If the energy increase is less than Multisite.LFD.Min.ThreshEnergyRise (which defaults to ten times Multisite.LFD.Min.ThreshE) then convergence is deemed to have been reached.

Go to top.

On-site support functions

On-site support functions (OSSF) are similar to multi-site support functions, but are linear combinations of PAOs only on the target atom. In this case, Atom.MultisiteRange should be small enough not to include any neighboring atoms (suggested values between 0.1 to 0.5). The number of support functions must be equivalent to the number of functions in an SZP basis (if polarisation functions are in the basis set) or an SZ basis (if there are no polarisation functions). The parameter Multisite.nonminimal should be set to true if polarisation functions are included.

The coefficients can be determined in the same was as for MSSF (with the LFD method and/or the numerical optimisation described above). It is likely that significant improvement in accuracy will be found with numerical optimisation. It is also important to test the effect of the parameter Atom.LFDRange which should be large enough to include several shells of neighbouring atoms.

The OSSF approach is most likely to be useful when linear scaling calculations with large basis sets are required. An example set of parmeters is found below.

Basis.BasisSet PAOs
Basis.MultisiteSF T
Multisite.LFD T
Multisite.nonminimal T

minE.VaryBasis T

# example of Si
%block Si
Atom.NumberOfSupports 9
Atom.MultisiteRange 0.1
Atom.LFDRange 8.0
%endblock

Go to top.

Blips

Blips (which are a type of piecewise continuous polynomial called a B-spline) [B4] are useful for very accurate calculations, since the basis set can be systematically improved, in the same way as a planewave basis set. However, the calculations can be expensive depending on the parameters, and the code for blip optimisation is under development. The following description, and possible keywords, may change during development.

The blips are defined on a blip grid, which is a regular cubic grid centred on the atoms, which also moves with the atoms. The basis set can be systematically improved, by increasing the support function radius and/or reducing the spacing of the blip grids. (The support grid spacing, which defines the grid for the blips, is equivalent to a plane wave cutoff; for a given support grid spacing the energy decreases variationally with support function radius.) For each species, we need to provide these two parameters, as well as the number of support functions, which should have a minimal basis size. (At present, the smallest blip-grid spacing is used for all species.)

For a given atom, we would set:

%block atom
Atom.NumberOfSupports                        4
Atom.SupportFunctionRange                  6.0
Atom.SupportGridSpacing                    0.3
%endblock

For each atomic species, an ion file with a minimal (SZ) basis set is required for the charge density and to initialise the blips.

The blip-grid spacing is directly related to the cutoff energy of the wavefunctions in planewave calculations. For a given cutoff energy \(E_{\rm cutoff}\) in Hartree, the blip-grid spacing should be \(\frac{2\pi}{\sqrt{2 E_{\rm cutoff}}}\) in bohr. Note that the grid spacing of integration grids (or FFT grids for the charge density) should be half the spacing of the blip grid, or smaller.

It is essential to optimise the support functions (blip coefficients) in the case of blips. The tolerance and maximum number of iterations can be set with the following keywords:

minE.VaryBasis              T
minE.EnergyTolerance             0.10E-07
minE.SupportVariations             30

It is not recommended, but if memory problems are encountered for very accurate blip calculations, you may need to switch off the preconditioning procedure for length-scale ill conditioning by setting the parameter minE.PreconditionBlips F

Go to top.

Reading coefficients from files

The calculated linear-combination coefficients of the support functions are stored in SFcoeffmatrix2 files for PAOs or blip_coeffs files for blips. Those files can be read by setting Basis.LoadCoeffs T in the subsequent calculations.

Go to top.

Basis Set Superposition Error

Basis set superposition error (BSSE) arises when the two monomer units come closer and the basis set localized on one unit can act as diffuse functions for the electrons from the other unit, and therefore could be responsible for the overestimation of the binding energy for the interacting systems. It is unlikely to affect blip basis calculations [B5].

To correct this BSSE, the Counterpoise (CP) correction method [B6] is used, where the artificial stabilization is controlled by enabling the atoms in monomer calculations to improve their basis sets by including the basis sets from other monomers (using so-called ghost atoms).

When systems A and B approach and make a new system AB, the typical interaction energy between A and B is calculated as:

\(E_{AB}^{int} = E_{AB}(AB) - E_A(A) - E_B(B).\)

where \(E_{AB}(AB)\) is the energy of system AB and \(E_{A}(A)\) and \(E_{B}(B)\) are the energies of isolated A and B. The lowerscript and parentheses correspond to the system and its structure, respectively.

Now, the estimate for the amount of artificial stabilization of A coming from the extra basis functions from B is:

\(E_{A}^{BSSE} = E_{A\bar{B}}(AB) - E_A(A\text{ in }AB),\)

where \(\bar{A}\) and \(\bar{B}\) are the ghost atoms, which have basis functions, but no potential or charge density. \(E_{A\bar{B}}(AB)\) is the energy of system A with the basis sets from ghost-atom system B in the AB structure. \(E_A(A\text{ in }AB)\) is the energy of system A in the AB structure but without system B (neither basis functions nor atoms). Therefore, the subtraction corresponds to how much system A is stabilized by the basis function of B.

Similarly, for monomer B,

\(E_{B}^{BSSE} = E_{\bar{A}B}(AB) - E_B(B\text{ in }AB),\)

Subtracting the BSSE part of A and B units from the typical interaction energy mentioned above, the counterpoise corrected interaction energy without BSSE \((E_{AB}^{int,CP})\) will be:

\(E_{AB}^{int,CP} = E_{AB}^{int} - E_{A}^{BSSE} - E_{B}^{BSSE}.\)

Practically, to calculate \(E_{A\bar{B}}(AB)\), the basis functions of B should be placed on atomic centers of B, however with zero nuclear charge and mass. This can be performed in CONQUEST by specifying negative masses for the ghost atoms in B in the block ChemicalSpeciesLabel of the input file:

%block ChemicalSpeciesLabel
  1   1.01  A
  2  -1.01  B
%endblock

Go to top.

[B1]

David R. Bowler, Jack S. Baker, Jack T. L. Poulton, Shereif Y. Mujahed, Jianbo Lin, Sushma Yadav, Zamaan Raza, and Tsuyoshi Miyazaki. Highly accurate local basis sets for large-scale DFT calculations in conquest. Jap. J. Appl. Phys., 58:100503, 2019. doi:10.7567/1347-4065/ab45af.

[B2]

Ayako Nakata, David R. Bowler, and Tsuyoshi Miyazaki. Efficient Calculations with Multisite Local Orbitals in a Large-Scale DFT Code CONQUEST. J. Chem. Theory Comput., 10:4813, 2014. doi:10.1039/C5CP00934K.

[B3]

Ayako Nakata, David Bowler, and Tsuyoshi Miyazaki. Optimized multi-site local orbitals in the large-scale DFT program CONQUEST. Phys. Chem. Chem. Phys., 17:31427, 2015. doi:10.1021/ct5004934.

[B4]

E. Hernández, M. J. Gillan, and C. M. Goringe. Basis functions for linear-scaling first-principles calculations. Phys. Rev. B, 55:13485–13493, 1997. doi:10.1103/PhysRevB.55.13485.

[B5]

P.D. Haynes, C.-K. Skylaris, A.A. Mostofi, and M.C. Payne. Elimination of basis set superposition error in linear-scaling density-functional calculations with local orbitals optimised in situ. Chem. Phys. Lett., 422:345 – 349, 2006. doi:10.1016/j.cplett.2006.02.086.

[B6]

S.F. Boys and F. Bernardi. The calculation of small molecular interactions by the differences of separate total energies. some procedures with reduced errors. Mol. Phys., 19:553, 1970. doi:10.1080/00268977000101561.

Go to top.