.. _lcao:

=========
LCAO Mode
=========

.. highlight:: bash

GPAW supports an alternative mode of calculation,
:dfn:`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
:ref:`LCAOTDDFT <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]_ :

.. math::

 \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

.. math::

  \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

.. math::

 \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`_.

.. _1: ../../_static/askhl_master.pdf
.. _2: ../../_static/marco_master.pdf

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
:command:`gpaw-basis` tool, located in your
:file:`{gpaw-directory}/tools` directory. For example, typing::

  $ gpaw-basis H Cl

will generate the basis-set files :file:`H.dzp.basis` and
:file:`Cl.dzp.basis` for hydrogen and chlorine with sensible default
parameters. Note that :file:`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 :ref:`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.

.. literalinclude:: lcao_h2o.py

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
:command:`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:

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',
           ...)

.. _poisson_performance:

Notes on performance
--------------------

For larger LCAO calculations, it is crucial to use ScaLAPACK
and recommended to also use Elpa.
See the dedicated section on :ref:`manual_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
:ref:`most important parameters <manual_parallel>` 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.

.. literalinclude:: lcao_opt.py

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

.. autoclass:: gpaw.atom.basis.BasisMaker

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.

.. literalinclude:: basisgeneration.py

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)