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] :
The basis functions are constructed as products of numerical radial functions and spherical harmonics
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
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 :ref:` parallel keyword <manual_parallel>`.
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.
- 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)