LCAO Mode

GPAW supports an alternative mode of calculation, LCAO mode [LCAO-article], which will use a basis set of atomic orbital-like functions rather than grid-based wave functions. This makes calculations considerably cheaper, although the accuracy will be limited by the quality of the chosen basis.

The sections below explain briefly how LCAO mode works, how to generate basis sets and how to use them in calculations. LCAO mode is available for TD-DFT via the LCAOTDDFT module.

Introduction

In the LCAO mode the Kohn-Sham pseudo wave functions \(\tilde{\psi}_n(\mathbf r)\) are expanded onto a set of atomic-like orbitals \(\Phi_{nlm}(\mathbf r)\), in the same spirit as the SIESTA method [Siesta] :

\[\tilde{\psi}_n(\mathbf r) = \sum_\mu c_{\mu n} \Phi_{\mu}(\mathbf r)\]

The basis functions are constructed as products of numerical radial functions and spherical harmonics

\[\Phi_{nlm}(\mathbf{r}) = \Phi_{nlm}(\mathbf r^a + \mathbf R^a) = \varphi_{nl}(r^a) Y_{lm}(\hat{\mathbf{r}}^a)\]

where \(\mathbf R^a\) is the position of nucleus \(a\), and \(\mathbf r^a = \mathbf r - \mathbf R^a\).

In this approximation the variational parameters are the coefficients \(c_{\mu n}\) rather than the real space wave function. The eigenvalue problem then becomes

\[\sum_\nu H_{\mu\nu} c_{\nu n} = \sum_{\nu} S_{\mu\nu} c_{\nu n} \epsilon_n\]

which can be solved by directly diagonalizating the Hamiltonian in the basis of the atomic orbitals.

Some detailed information can be found in the master theses 1 and 2.

Basis-set generation

In order to perform an LCAO calculation, a basis-set must be generated for every element in your system. This can be done by using the gpaw-basis tool, located in your gpaw-directory/tools directory. For example, typing:

$ gpaw-basis H Cl

will generate the basis-set files H.dzp.basis and Cl.dzp.basis for hydrogen and chlorine with sensible default parameters. Note that dzp stands for double zeta polarized which is the default basis-set type. The basis-set should be placed in the same directory as the GPAW setups (see Installation of PAW datasets for details). For a complete list of the parameters do:

$ gpaw-basis --help

For technical reasons, the basis set generator always generates the corresponding PAW, even if the latter exists on the user’s system. Use the --save-setup option to save the calculated setup along with the basis set.

Running a calculation

In order to run an LCAO calculation, the lcao mode and a basis-set should be set in the calculator:

>>> calc = GPAW(mode='lcao',
>>>             basis='dzp',
>>>             ...)

The calculator can then be used in the usual way. The basis keyword accepts the same types of values as the setups keyword, such as basis={'H' : 'dzp', 'O' : 'mine', 'C' : 'szp'}.

For larger systems, to get good performance, be sure to enable ScaLAPACK to parallelize the cubic-scaling diagonalization step and distribute many matrices. If possible, install and enable Elpa [Elpa] to further save time. See the parallel keyword.

Example

The following example will relax a water molecule using the LCAO calculator. The QuasiNewton minimizer will use the forces calculated using the localized basis set.

from ase import Atoms
from ase.optimize import QuasiNewton
from gpaw import GPAW

a = 6
b = a / 2

mol = Atoms('H2O',
            [(b, 0.7633 + b, -0.4876 + b),
             (b, -0.7633 + b, -0.4876 + b),
             (b, b, 0.1219 + b)],
            cell=[a, a, a])

calc = GPAW(nbands=4,
            h=0.2,
            mode='lcao',
            basis='dzp')

mol.calc = calc
dyn = QuasiNewton(mol, trajectory='lcao_h2o.traj')
dyn.run(fmax=0.05)

It is possible to switch to the Finite Difference mode and further relax the structure simply by doing:

>>> calc.set(mode='fd')
>>> dyn.run(fmax=0.05)

More on basis sets

A minimal basis set consists of one atomic orbital-like function for each valence state of the atom. Extra radial functions can be added to improve the span of the basis; basis sets are called single-zeta (sz), double-zeta (dz) and so on, depending on the number of such radial functions per valence state.

It is normally desirable to add a basis function corresponding to the lowest unoccupied angular momentum quantum number. This is called a polarization function.

Double-zeta polarized basis sets are normally required and sufficient to obtain results of reasonable accuracy; they are also the basis type generated by default.

A dzp basis set for nitrogen will have 2s and 2p valence states, each with two radial functions, plus a polarization function of type d, for a total of 5 distinct radial functions. Each will be degenerate by \(2l+1\), meaning that GPAW will use a total of 13 basis functions to represent the atom during a calculation. Transition metals, having s- and d-type valence states, get a p-type polarization function and thus a total of 15 basis functions.

To plot already generated basis functions, use the gpaw-analyse-basis command like:

$ gpaw-analyse-basis -f H.dzp.basis O.dzp.basis

This will plot the basis functions in the specified files. If the -f option is not included, the script will look for the first matching file in the GPAW setups paths, rather than the precise specified files. Run gpaw-analyse-basis --help for more options.

Ghost atoms and basis set superposition errors

In the vicinity of a surface with many basis functions, an adsorbate can “benefit” from the degrees of freedom from the surface basis functions, resulting in a lower energy compared to a calculation on the isolated adsorbate and thus too strong binding energy. This error referred to as the basis set superposition error. It can be eliminated by adding ghost atoms to the calculation on the isolated adsorbate. Ghost atoms possess basis functions as normal but do not otherwise affect the calculation (no projectors, compensation charges and so on), thereby ensuring that the same degrees of freedom are available to the wave functions in any calculation.

A calculation with ghost atoms is performed precisely like a normal calculation, in the sense that the ASE atoms object should contain all the involved atoms including those which are ghosts, with the only difference being that ghost atoms have their setup type set to ghost. It is stressed that the only difference between an ordinary atom and the corresponding ghost atom is the setup type.

Perform a calculation using ghost copper atoms and ordinary oxygen and hydrogen atoms:

>>> GPAW(setups={'Cu' : 'ghost', 'O' : 'paw', 'H' : 'paw'},
         basis='dzp',
         mode='lcao',
         ...)

Perform a calculation where atom 17 and atom 42 (designated by their indices in the Atoms object) use ordinary setups, while all other atoms are ghosts:

>>> GPAW(setups={'default': 'ghost', 17: 'paw', 42:'paw'},
         basis='dzp',
         mode='lcao',
         ...)

Notes on performance

For larger LCAO calculations, it is crucial to use ScaLAPACK and recommended to also use Elpa. See the dedicated section on ScaLAPACK for more information. Below are some hints on how to obtain good performance for operations not related to ScaLAPACK.

The only difference between the FD (grid-based finite-difference) and LCAO modes is the way in which pseudo wave functions are represented. The usual real-space grid methods are still used for the density and potential. The associated computations will therefore take a larger percentage of the CPU time compared to FD mode, where operations on the wave functions usually dominate. Thus it makes sense to pay some attention to the performance of these operations.

This example shows the most important parameters to achieve good performance with LCAO. The example is actually much too small to make much use of parallelism, but these parameters will provide good performance for large-scale systems.

from ase.build import molecule
from ase.optimize import QuasiNewton
from gpaw import GPAW

atoms = molecule('CH3CH2OH', vacuum=4.0)
atoms.rattle(stdev=0.1)  # displace positions randomly a bit

calc = GPAW(mode='lcao',
            basis='dzp',
            nbands='110%',
            parallel=dict(band=2,  # band parallelization
                          augment_grids=True,  # use all cores for XC/Poisson
                          sl_auto=True,  # enable ScaLAPACK parallelization
                          use_elpa=True))  # enable Elpa eigensolver
atoms.calc = calc

opt = QuasiNewton(atoms, trajectory='opt.traj')
opt.run(fmax=0.05)

Note

The following paragraph refers to the old FDPoissonSolver. This has since been replaced by FastPoissonSolver which always performs well, and for which the paragraph does not apply.

The multigrid method used in the FD Poisson solver relies on alternating interpolations and restrictions of the density on grids of different sizes. Make sure that the grid point count along each axis is divisible by 8, by specifying e.g. gpts=(96, 96, 96) when creating the calculator – this will dramatically reduce the number of required Poisson iterations in large or very oblong systems in those cases where the code would otherwise have chosen a grid point count not divisible by 8. By default, the FD Poisson solver uses the Jacobi method. To increase performance further use the Gauss-Seidel method instead, which usually reduces the Poisson iteration count by around 40% (ideally 50%). Again, please note that none of the above applies to the FastPoissonSolver which is now default.

Advanced basis generation

The class gpaw.atom.basis.BasisMaker is the backend of the basis generation programme. Use this to create basis sets with specialized parameters that cannot be set using the command line interface. In particular, the basis generator relies on the setup generator to define the basis functions; therefore, any parameters that apply to setup generation will equally apply to basis set generation.

class gpaw.atom.basis.BasisMaker(generator, name=None, run=True, gtxt='-', non_relativistic_guess=False, xc='PBE', save_setup=False)[source]

Class for creating atomic basis functions.

This example shows how to generate an RPBE double-zeta basis set for gold, in which the otherwise empty p-state is considered a valence state, and using a non-standard size of the augmentation sphere.

from gpaw.atom.generator import Generator
from gpaw.atom.configurations import parameters
from gpaw.atom.basis import BasisMaker

symbol = 'Au'
args = parameters[symbol]  # Dictionary of default setup parameters
args['rcut'] = 2.6  # Set cutoff of augmentation sphere

generator = Generator(symbol, 'RPBE')
generator.N *= 2  # Increase grid resolution
generator.run(write_xml=False, **args)

bm = BasisMaker(generator, name='special', run=False)

# Create double-zeta RPBE basis where p orbital is considered a valence state
# (ordinary dzp basis would use a smaller p-type Gaussian for polarization)
# The list jvalues indicates which states should be included, and the
# ordering corresponds to the valence states in the setup.
basis = bm.generate(zetacount=2, polarizationcount=0,
                    energysplit=0.1, jvalues=[0, 1, 2],
                    rcutmax=12.0)

basis.write_xml()  # Dump to file 'Au.special.dz.basis'

Miscellaneous remarks

In FD or PW mode, a single LCAO iteration is used to initialize the wave functions and density. Specifying a basis to the calculator in FD or PW mode can be used to increase the quality of the initial guess, but does not in any other way affect the subsequent iterations:

>>> calc = GPAW(mode='fd', basis='dzp', ...)

In either mode, if a basis is not specified to the calculator, the calculator will use the pseudo partial waves \(\tilde \phi_i^a(\mathbf r)\), smoothly truncated to 8 Bohr radii, as a basis. This corresponds roughly to a single-zeta basis in most cases. Depending on the unoccupied states defined on the PAW setups, it may be roughly equivalent to a single-zeta polarized basis set for certain elements.

Local Orbitals

In LCAO mode, it is possible to obtain a reduced basis set of localised orbitals that can be used to define effective tight-binding Hamiltonians. Contrary to Wannier functions (WFs), the local orbital (LO) construction is not based on a projection of the Kohn-Sham states and does not require any physical input such as the initial guesses for the WFs. In fact, the LOs are obtained directly from a sub-diagonalization of the LCAO Hamiltonian. The LOs are constructed for any atom in the system through a sub-diagonalization of the Hamiltonian block of its AOs. This procedure yields a set of LOs whichare atomic-like functions and are by construction atom-centred and orthogonal within the same atom (but not among different atoms). Furthermore, the LO representation can coexist with the original AO one ne, in the sense that one can sub-diagonalize only a subset of atoms in the system. This is useful if one is particularly interested in a limited part of a system, such as a molecular bridge in a quantum junction, or an adsorbate on a substrate. More details and examples can be found in Local Orbitals tutorial.

[LCAO-article]

A. H. Larsen, M. Vanin, J. J. Mortensen, K. S. Thygesen, and K. W. Jacobsen, Phys. Rev. B 80, 195112 (2009)

[Siesta]

J.M. Soler et al., J. Phys. Cond. Matter 14, 2745-2779 (2002)

[Elpa]

A Marek et al., J. Phys.: Condens. Matter 26 213201 (2014)