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