LCAO Mode

GPAW supports an alternative mode of calculation, LCAO mode, 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'}.

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

The multigrid method used in the 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 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%).

The convergence criterion of the Poisson solver in FD mode, eps=2e-10, is very strict. A value of around eps=1e-7 can reduce the required Poisson iteration count considerably without increasing the required number of SCF steps. Larger values like eps=1e-5 tend to increase the number of SCF steps, possibly making the calculation take longer.

Example:

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

a = 6
b = a / 2

mol = Atoms('H2O',
             positions=[(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,
            mode='lcao',
            basis='dzp',
            gpts=(32, 32, 32),
            poissonsolver=PoissonSolver(relax='GS', eps=1e-7))

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

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

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