.. _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. .. _lcaotddft: 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 ` 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. 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 `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. 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. 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. 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. 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. .. 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. 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. 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. .. _los in lcao: 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 `los tutorial` 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)