scnc
/
conquest
mirror of https://github.com/OrderN/CONQUEST-release.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #18 from OrderN/manual 

Commit pre-release form of manual
This commit is contained in:
David Bowler 2020-01-24 15:55:35 +00:00 committed by GitHub
parent 3a92702224 699095b488
commit 82d87ccc62
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
25 changed files with 5657 additions and 66 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
docs/CONQUEST-manual/MDPlot.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 24 KiB

515
docs/CONQUEST-manual/basissets.rst Normal file
View File

@ -0,0 +1,515 @@
.. _basissets:
==========
Basis sets
==========
As we have mentioned in :ref:`finding the ground state <groundstate>`,
the density matrix is represented by :ref:`support functions
<gs_suppfunc>`. These, in turn, are made up of basis functions, and
the choice of basis set and how it represents the support functions
affects both the accuracy and performance of a calculation.
There are two kinds of basis functions which are used in CONQUEST to
represent the support functions: pseudo-atomic orbitals (keyword
``PAOs``) which are the default; and b-splines (keyword ``blips``)
which allow systematic convergence at the expense of greater
complexity.
The basis set is selected as follows:
::
Basis.BasisSet PAOs
When the basis set is taken to be PAOs, there are three different ways
to construct the support functions, discussed below:
* Each support function is represented by a single PAO (primitive PAOs)
* Multi-site support functions, built from PAOs on several atoms
* on-site support functions, built from PAOs on one atom
Primitive PAOs are efficient for small systems. When using large PAO
basis sets for systems containing more than several hundred atoms,
multi-site support functions and on-site support functions will be
more efficient than primitive PAOs.
Go to :ref:`top <basissets>`.
.. _basis_paos:
Pseudo-atomic orbitals (PAOs)
-----------------------------
PAOs are solutions of the Schrodinger equation for isolated atoms,
using pseudopotentials, with some confinement applied. They consist
of radial functions multiplied by spherical harmonics. For the
valence orbitals, the radial functions are referred to as *zeta*
(:math:`\zeta`), while for unoccupied orbitals, they are termed
*polarisation*.
A minimal PAO basis set is *single-*:math:`\zeta` *(SZ)*, with one
radial function for each angular momentum quantum number in the
valence electrons. While the cost is significantly lower than for
other basis sets, the accuracy will be rather low.
The accuracy of a calculation can be improved by adding polarisation
functions and multiple radial functions for different angular momentum
values, though systematic improvement is rather difficult to achieve
(this is straightforward with a :ref:`blip function basis<basis_blips>`).
The PAO utility included with CONQUEST generates basis sets with
differing sizes and accuracies; full details of the performance of
these basis sets can be found elsewhere :cite:`b-Bowler:2019fv`.
* minimal (single zeta)
* small (single zeta and polarisation, SZP)
* medium (double zeta, single polarisation, DZP)
* large (triple zeta, double polarisation, TZDP)
Go to :ref:`top <basissets>`.
.. _basis_primitivepaos:
Primitive PAOs as support functions
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The easiest way to prepare support functions is to use primitive PAOs
as the support functions without any modifications. In this case, the
input parameters related to the support functions are automatically
set by obtaining the information from the PAO files (``.ion``
files) so long as they are generated by the CONQUEST ``MakeIonFiles``
utility, version 1.0.3 or later. No further input parameters
need to be set in ``Conquest_input``.
Go to :ref:`top <basissets>`.
.. _basis_mssf:
Multi-site support functions
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Since the computational cost of Conquest scales cubically with the
number of support functions, contracting the PAOs into a smaller set
of support functions is an efficient way to
reduce the computational cost when we use large multiple-:math:`\zeta`
PAO basis sets. Multi-site support functions (MSSFs)
:cite:`b-Nakata2014,b-Nakata2015` are constructed for each
atom by taking linear combinations of the atom's
PAOs and the PAOs from neighbouring atoms within a certain range
(set with the parameter ``Atom.MultisiteRange`` in
the :ref:`atom specification block <input_atomic_spec>`).
Multi-site support functions can be selected by setting the following
parameters:
::
Basis.BasisSet PAOs
Basis.MultisiteSF T
Various other parameters need to be set in the
:ref:`atom specification block <input_atomic_spec>`.
The number of support functions for the atoms must be set, and is
normally equivalent to a minimal (single zeta) basis; it is set with
``Atom.NumberOfSupports``.
(To use a number of support functions larger than this minimal number, the
parameter ``Multisite.nonminimal`` needs to be set to ``T``.)
The range for the multi-site support functions (the PAOs of any atom within this
distance of the atom will be included in the support functions)
is set with ``Atom.MultisiteRange``. The accuracy of the MSSF will
improve as this range is increased, though the computational cost will
also increase; careful tests must be made to find an appropriate
range. For a minimal number of MSSF, the range must be large enough
to include other atoms, though this restriction can be removed (see
:ref:`on-site support functions <basis_ossf>` for more details).
As well as setting the range for the MSSFs, we need to specify an
approach for finding the expansion coefficients. A reasonable set of MSSF
coefficients can be found using the *local filter diagonalization
(LFD)* method. For improved accuracy, this should
be followed by variational *numerical optimisation*.
Go to :ref:`top <basissets>`.
.. _basis_mssf_lfd:
Local filter diagonalization (LFD)
++++++++++++++++++++++++++++++++++
In this method, which is selected by setting ``Multisite.LFD T``, the
MSSF coefficients are found by diagonalising the
Hamiltonian in the primitive PAO basis, for a small cluster of atoms
surrounding the target atom. The MSSF coefficients :math:`C` are determined by
projecting the sub-space molecular orbitals :math:`C_{sub}` around each
atom onto localized *trial* vectors :math:`t`,
:math:`C = C_{sub} f(\varepsilon_{sub}) C_{sub}^T S_{sub} t`
The cluster for diagonalisation must be at least as large as the MSSF
range, but larger clusters tend to give better MSSF coefficients (at the
expense of an increased computational cost).
The LFD sub-space region is determined for each atom by setting
``Atom.LFDRange``.
An example set of parameters for an MSSF calculation for bulk Si would be:
::
Basis.BasisSet PAOs
Basis.MultisiteSF T
Multisite.LFD T
%block ChemicalSpeciesLabel
1 28.07 Si
%endblock
%block Si
Atom.NumberOfSupports 4
Atom.MultisiteRange 8.0
Atom.LFDRange 8.0
%endblock
When calculating binding energy curves or optimising cells, a change of
lattice constant can suddenly bring a new set of atoms within the
range of the support functions. In this case, a smearing can be
applied at the edges of the range, by setting ``Multisite.Smear T``.
Further details are :ref:`given below <basis_mssf_advanced>`.
Some form of self-consistency between the MSSF and the charge density
is required (as the MSSF will determine the Hamiltonian and hence the
output charge density). At present, this is performed as a complete
SCF cycle for each set of MSSF coefficients (though this is likely to
be updated soon for improved efficiency). This is selected by default
(but can be turned off by setting the parameter ``Multisite.LFD.NonSCF T``).
This iterative process is not variational, but is terminated when the
absolute energy change between iterations is less than
``Multisite.LFD.Min.ThreshE``, or the residual (defined in
:ref:`self-consistency <gs_scf>`) is less than
``Multisite.LFD.Min.ThreshD``.
An example input block for this process would be as follows:
::
Multisite.LFD T
Multisite.LFD.Min.ThreshE 1.0e-6
Multisite.LFD.Min.ThreshD 1.0e-6
Go to :ref:`top <basissets>`.
.. _basis_sets_numopt:
Numerical optimisation
++++++++++++++++++++++
The MSSF coefficients can also be optimised by minimizing the
DFT energy with respect to the coefficients, in a variational
process. The threshold and the
maximum iteration number of the numerical optimisation are specified
by ``minE.EnergyTolerance`` and ``minE.SupportVariations``. The
optimisation is based on the conjugate gradient (CG) method, and the
initial CG step size can be specified by ``minE.InitStep_paomin``
(default is 5.0).
::
minE.VaryBasis T
minE.EnergyTolerance 1.0e-6
minE.SupportVariations 30
The numerical optimisation provides more accurate coefficients than
the LFD method but is usually more time consuming. Therefore, it is
generally better to start from good initial values, for example, the
coefficients calculated by LFD. When both ``Multisite.LFD``
and ``minE.VaryBasis`` are selected,
the initial coefficients will be calculated by LFD
and the coefficients will then be optimised.
::
Basis.MultisiteSF T
Multisite.LFD T
minE.VaryBasis T
If good initial coefficient values have been found in a previous
calculation, reading these from files (the base name of these files is
``SFcoeffmatrix2``) and performing only the
numerical optimisation is also a good choice.
::
Basis.LoadCoeffs T
Basis.MultisiteSF T
Multisite.LFD F
minE.VaryBasis T
Go to :ref:`top <basissets>`.
.. _basis_mssf_advanced:
Advanced MSSF concepts
++++++++++++++++++++++
*Smearing the edge of the support functions*
Here, we are concerned with changes of lattice constant which may
bring new atoms inside the support function range.
We can set the smearing-function type
``Multisite.Smear.FunctionType`` (default=1:Fermi-Dirac, 2=Error
function), the center position of the function
``Multisite.Smear.Center`` (default is equal to the range of the
support functions), offset of the center position
``Multisite.Smear.Shift`` and the width of the Fermi-Dirac function
``Multisite.Smear.Width`` (default=0.1).
*Selecting states from the sub-space*
Here, we consider how to create the MSSF themselves from the results
of the sub-space diagonalisation.
The Fermi function :math:`f` with :math:`\varepsilon_{sub}`
``Multisite.LFD.ChemP`` and :math:`kT` ``Multisite.LFD.kT`` in the
equation removes the effects of the subspace molecular orbitals in
higher energy region.
In default, :math:`\varepsilon_{sub}` is automatically set to the mean
value of the subspace HOMO and LUMO energies for each subspace. If
users want to modify this, set ``Multisite.LFD.UseChemPsub F`` and the
:math:`\varepsilon_{sub}` value with ``Multisite.LFD.ChemP``.
For the LFD trial functions :math:`t`, when ``Atom.NumberOfSupports``
is equal to the number of SZ or single-zeta plus polarization (SZP),
the PAOs which have the widest radial functions for each spherical
harmonic function are chosen as the trial vectors automatically in
default.
When ``Atom.NumberOfSupports`` is equal to the number of SZP and
``Multisite.nonminimal.offset`` is set, the other PAOs will have the
weight in the trial vectors with the value of
``Multisite.nonminimal.offset``.
The users can also provide the trial vectors from the input file using the ``LFDTrialVector`` block
::
# Trial vectors of Au (element 1) and O (element 2) atoms.
# Au: 15 PAOs (DZP) -> 6 support functions, O: 13 PAOs (DZP) -> 4 support functions.
%block LFDTrialVector
# species sf npao s s x y z d1 d2 d3 d4 d5 d1 d2 d3 d4 d5 for Au
1 1 15 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
1 2 15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
1 3 15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
1 4 15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
1 5 15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
1 6 15 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0
2 1 13 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2 2 13 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
2 3 13 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
2 4 13 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
# species sf npao s s x y z x y z d1 d2 d3 d4 d5 for O
%endblock LFDTrialVector
The first, second and third columns correspond to the indices of
species, support functions for each species, and the number of PAOs
for each species. The other columns provide the initial values of the
trial vectors. For example, in the first line in the above example,
the second *s* PAO is chosen as the trial vector for the first support
function of Au.
*Self-consistent LFD*
Two further conditions are applied to end the LFD self-consistency
process. The maximum number of iterations is set with
``Multisite.LFD.Min.MaxIteration``. It is also possible, as the
process is not variational, that the energy can increase as well as
decrease between iterations. If the energy *increase* is less than
``Multisite.LFD.Min.ThreshEnergyRise`` (which defaults to ten times
``Multisite.LFD.Min.ThreshE``) then convergence is deemed to have been
reached.
Go to :ref:`top <basissets>`.
.. _basis_ossf:
On-site support functions
^^^^^^^^^^^^^^^^^^^^^^^^^
On-site support functions (OSSF) are similar to multi-site support
functions, but are linear combinations of PAOs only on the target atom.
In this case, ``Atom.MultisiteRange`` should be small enough not to
include any neighboring atoms (suggested values between 0.1 to
0.5). The number of support functions must be equivalent to the number
of functions in an SZP basis (if polarisation functions are in the
basis set) or an SZ basis (if there are no polarisation functions).
The parameter ``Multisite.nonminimal`` should be set to true if
polarisation functions are included.
The coefficients can be determined in the same was as for MSSF (with
the LFD method and/or the numerical optimisation described above). It
is likely that significant improvement in accuracy will be found with
numerical optimisation. It is also important to test the effect of
the parameter ``Atom.LFDRange`` which should be large enough to
include several shells of neighbouring atoms.
The OSSF approach is most likely to be useful when linear scaling
calculations with large basis sets are required. An example set of
parmeters is found below.
::
Basis.BasisSet PAOs
Basis.MultisiteSF T
Multisite.LFD T
Multisite.nonminimal T
minE.VaryBasis T
# example of Si
%block Si
Atom.NumberOfSupports 9
Atom.MultisiteRange 0.1
Atom.LFDRange 8.0
%endblock
Go to :ref:`top <basissets>`.
.. _basis_blips:
Blips
-----
Blips (which are a type of piecewise continuous polynomial called a
B-spline) :cite:`b-Hernandez1997` are useful for very accurate calculations, since the basis set
can be systematically improved, in the same way as a planewave basis set. However, the
calculations can be expensive depending on the parameters, and
the code for blip optimisation is under development. The following
description, and possible keywords, may change during development.
The blips are defined on a blip grid, which is a regular cubic grid
centred on the atoms, which also moves with the atoms. The basis set
can be systematically improved, by increasing the support function radius
and/or reducing the spacing of the blip grids. (The support grid
spacing, which defines the grid for the blips, is equivalent to a
plane wave cutoff; for a given support grid spacing the energy
decreases variationally with support function radius.) For each species,
we need to provide these two parameters, as well as the
number of support functions, which should have a minimal basis size.
(At present, the smallest blip-grid spacing is used for all species.)
For a given atom, we would set:
::
%block atom
Atom.NumberOfSupports 4
Atom.SupportFunctionRange 6.0
Atom.SupportGridSpacing 0.3
%endblock
For each atomic species, an ion file with a minimal (SZ) basis set is
required for the charge density and to initialise the blips.
The blip-grid spacing is directly related to the cutoff energy of the
wavefunctions in planewave calculations. For a given cutoff
energy :math:`E_{\rm cutoff}` in Hartree, the blip-grid spacing should
be :math:`\frac{2\pi}{\sqrt{2 E_{\rm cutoff}}}` in bohr. Note that
the grid spacing of integration grids (or FFT grids for the charge
density) should be half the spacing of the blip grid, or smaller.
It is essential to optimise the support functions (blip coefficients)
in the case of blips. The tolerance and maximum number of iterations
can be set with the following keywords:
::
minE.VaryBasis T
minE.EnergyTolerance 0.10E-07
minE.SupportVariations 30
It is not recommended, but if memory problems are encountered for
very accurate blip calculations, you may need to switch off the
preconditioning procedure for length-scale ill conditioning by setting
the parameter ``minE.PreconditionBlips F``
Go to :ref:`top <basissets>`.
.. _basis_readcoeffs:
Reading coefficients from files
-------------------------------
The calculated linear-combination coefficients of the support
functions are stored in ``SFcoeffmatrix2`` files for PAOs or
``blip_coeffs`` files for blips. Those files can be read by setting
``Basis.LoadCoeffs T`` in the subsequent calculations.
Go to :ref:`top <basissets>`.
.. _basis_bsse:
Basis Set Superposition Error
-----------------------------
Basis set superposition error (BSSE) arises when the two monomer
units come closer and the basis set localized on one unit can act as
diffuse functions for the electrons from the other unit, and therefore
could be responsible for the overestimation of the binding energy for
the interacting systems. It is unlikely to affect blip basis
calculations :cite:`b-Haynes:2006qe`.
To correct this BSSE, the Counterpoise (CP) correction method
:cite:`b-Boys:1970aa` is used, where the artificial stabilization is
controlled by enabling the atoms in monomer calculations to improve
their basis sets by including the basis sets from other monomers
(using so-called ghost atoms).
When systems A and B approach and make a new system AB, the typical
interaction energy between A and B is calculated as:
:math:`E_{AB}^{int} = E_{AB}(AB) - E_A(A) - E_B(B).`
where :math:`E_{AB}(AB)` is the energy of system AB and
:math:`E_{A}(A)` and :math:`E_{B}(B)` are the energies of isolated A
and B. The lowerscript and parentheses correspond to the system and
its structure, respectively.
Now, the estimate for the amount of artificial stabilization of A
coming from the extra basis functions from B is:
:math:`E_{A}^{BSSE} = E_{A\bar{B}}(AB) - E_A(A\text{ in }AB),`
where :math:`\bar{A}` and :math:`\bar{B}` are the ghost atoms, which
have basis functions, but no potential or charge density.
:math:`E_{A\bar{B}}(AB)` is the energy of system A with
the basis sets from ghost-atom system B in the AB structure. :math:`E_A(A\text{ in }AB)`
is the energy of system A in the AB structure but without system B
(neither basis functions nor atoms). Therefore, the subtraction
corresponds to how much system A is stabilized by the basis function
of B.
Similarly, for monomer B,
:math:`E_{B}^{BSSE} = E_{\bar{A}B}(AB) - E_B(B\text{ in }AB),`
Subtracting the BSSE part of A and B units from the typical
interaction energy mentioned above, the counterpoise corrected
interaction energy without BSSE :math:`(E_{AB}^{int,CP})` will be:
:math:`E_{AB}^{int,CP} = E_{AB}^{int} - E_{A}^{BSSE} - E_{B}^{BSSE}.`
Practically, to calculate :math:`E_{A\bar{B}}(AB)`, the basis
functions of B should be placed on atomic centers of B, however with
zero nuclear charge and mass. This can be performed in CONQUEST by
specifying negative masses for the ghost atoms in B in the ``block
ChemicalSpeciesLabel`` of the input file:
::
%block ChemicalSpeciesLabel
1 1.01 A
2 -1.01 B
%endblock
Go to :ref:`top <basissets>`.
.. bibliography:: references.bib
:cited:
:labelprefix: B
:keyprefix: b-
:style: unsrt
Go to :ref:`top <basissets>`.

BIN
docs/CONQUEST-manual/compare.jpg Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 98 KiB

6
docs/CONQUEST-manual/conf.py
View File

@ -50,8 +50,8 @@ master_doc = 'index'
# General information about the project.
project = u'CONQUEST'
copyright = u'2018, CONQUEST Developers'
author = u'J. Kane Shenton'
copyright = u'2018-2019, CONQUEST Developers'
author = u'CONQUEST Developers'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the
@ -146,7 +146,7 @@ latex_elements = {
# author, documentclass [howto, manual, or own class]).
latex_documents = [
(master_doc, 'CONQUEST.tex', u'CONQUEST Documentation',
u'J. Kane Shenton', 'manual'),
u'CONQUEST Developers', 'manual'),
]

228
docs/CONQUEST-manual/convergence.rst Normal file
View File

@ -0,0 +1,228 @@
.. _convergence:
=====================
Converging Parameters
=====================
There are various important parameters in CONQUEST that affect the
convergence of the total energy, and need to be tested. Integrals are
calculated on a grid; the density matrix is found approximately; a
self-consistent charge density is calculated; and support functions
are, in some modes of operations, optimised. These parameters are
described here.
.. _conv_grid:
Integration Grid
----------------
While many integrals are calculated analytically or on fine grids that
move with the atoms, there are still some integrals that must be found
numerically, and CONQUEST uses an orthorhombic, uniform grid to
evaluate these integrals (this grid is also used for the Fourier
transforms involved in finding the Hartree potential). The
spacing of the grid will affect the accuracy of the calculation, and
it is important to test the convergence of the total energy with the
grid spacing.
The grid spacing can be set intuitively using an energy (which
corresponds to the kinetic energy of the shortest wavelength wave that
can be represented on the grid). In atomic units, :math:`E = k^2/2`
with :math:`k = \pi/\delta` for grid spacing :math:`\delta`. The
cutoff is set with the parameter:
::
Grid.GridCutoff E
where ``E`` is an energy in Hartrees. The grid spacing can also be
set manually, by specifying the number of grid points in each
direction:
::
Grid.PointsAlongX N
Grid.PointsAlongY N
Grid.PointsAlongZ N
If setting the grid in this manner, it is important to understand a
little more about the internal workings of CONQUEST. The grid is divided up into
*blocks* (the default size is 4 by 4 by 4), and the number of grid
points in any direction must correspond to an integer multiple of the
block size in that direction. The block size can be set by the user:
::
Grid.InBlockX N
Grid.InBlockY N
Grid.InBlockZ N
Note that the blocks play a role in parallelisation and memory use, so
that large blocks may require larger memory per process; we recommend
block sizes no larger than 8 grid points in each direction.
There is also, at present, a restriction on the total number of grid
points in anuy direction, that it must have prime factors of only 2, 3 and 5. This will be
removed in a future release.
Go to :ref:`top <convergence>`.
.. _conv_dm:
Finding the density matrix
--------------------------
As discussed in the section on :ref:`finding the ground
state<groundstate>`,
the density matrix is found either
with exact diagonaliation, or the linear scaling approach. These
two methods require different convergence tests, and are described separately.
.. _conv_dm_bz:
Diagonalisation: Brillouin Zone Sampling
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The sampling of the Brillouin zone must be tested for convergence, and
the parameters are described :ref:`here <gs_diag_bz>`. The
convergence of charge density will be faster than detailed electronic
structure such as density of states (DOS), and it will be more
accurate for these types of calculations to generate a converged charge
density, and then run non self-consistently (see the section on
:ref:`self consistency <gs_scf>`) with appropriate k-point sampling.
Go to :ref:`top <convergence>`.
.. _conv_on:
Linear Scaling
~~~~~~~~~~~~~~
The range applied to the density matrix (``DM.L_range``) determines
the accuracy of the calculation, as well as the computational time
required (as the number of non-zero elements will increase based on a
sphere with the radius of the range, the time will increase roughly
proportional to the cube of the range). In almost all circumstances,
it is best to operate with a range which converges energy
*differences* and forces, rather than the absolute energy. Testing
for this convergence is an essential part of the preparation for
production calculations.
The tolerance applied to the density matrix optimisation
(``minE.Ltolerance``) must be
chosen to give adequate convergence of the energy and forces. The
tolerance is applied to the residual in the calculation, defined as:
.. math::
R = \sqrt{\sum_{i\alpha j\beta} \partial E/\partial L_{i\alpha j\beta}
\cdot \partial E/\partial L_{i\alpha j\beta} }
The dot product uses the inverse of the overlap matrix as the metric.
The approximate, sparse inversion of the overlap matrix is performed
before the optimisation of the density matrix. The method used,
Hotelling's method (a version of a Newton-Raphson approach) is
iterative and terminates when the characteristic quantity
:math:`\Omega` increases. On termination, if :math:`\Omega` is below
the tolerance ``DM.InvSTolerance`` then the inverse is accepted;
otherwise it is set to the identity (the density matrix optimisation
will proceed in this case, but is likely to be inefficient). We
define:
.. math::
\Omega = (Tr[I - TS])^2
where :math:`T` is the approximate inverse. The range for the inverse
must be chosen (``Atom.InvSRange`` in the species block); by default
it is same as the support function range
(which is then doubled to give the matrix range) but can be
increased. The behaviour of the inversion with range is not simple,
and must be carefully characterised if necessary.
Go to :ref:`top <convergence>`.
.. _conv_scf:
Self-consistency
----------------
The standard self-consistency approach uses the Pulay RMM method, and
should be robust in most cases. It can be monitored via the residual,
which is currently defined as the standard RMS difference in charge
density:
.. math::
R = \sqrt{\int \mathrm{d}\mathbf{r}\mid \rho^{out}(\mathbf{r}) -
\rho^{in}(\mathbf{r})\mid^2}
where :math:`\rho^{in}` is the input charge density for an iteration,
and :math:`\rho^{out}` is the resulting output charge density. The
SCF cycle is terminated when this residual is less than the parameter
``minE.SCTolerance``. The maximum number of iterations is set with
``SC.MaxIters`` (defaults to 50).
There are various further approaches and parameters which can be used
if the SCF cycle is proving hard to converge. As is standard, the
input for a given iteration is made by combining the charge density
from a certain number of previous steps (``SC.MaxPulay``, default 5).
The balance between input and output charge densities from these
previous steps is set with ``SC.LinearMixingFactor`` (default 0.5;
N.B. for spin polarised calculations,
``SC.LinearMixingFactor_SpinDown`` can be set separately). Reducing
this quantity may well improve stability, but slow down the rate of
convergence.
Kerker-style preconditioning (damping long wavelength charge
variations) can be selected using ``SC.KerkerPreCondition T`` (this is
most useful in metallic and small gap systems). The preconditioning
is a weighting applied in reciprocal space:
.. math::
K = \frac{1}{1+q^2_0/q^2}
where :math:`q_0` is set with ``SC.KerkerFactor`` (default 0.1).
This is often very helpful with slow convergence or instability.
Go to :ref:`top <convergence>`.
.. _conv_suppfunc:
Support Functions
-----------------
The parameters relevant to support functions depend on the basis set
that is used. In the case of pseudo-atomic orbitals (PAOs), when
support functions are primitive PAOs, the only relevant parameter is
the basis set size, which is set when the ion files are generated. It
is important to test the accuracy of a given basis set carefully for
the problem that is to be modelled.
When using multi-site support functions (MSSF), the key parameter is
the radius of the MSSF (``Atom.MultisiteRange`` in
the :ref:`atomic specification <input_atomic_spec>` block).
As this is increased, the accuracy of the
calculation will also increase, but with increased computational
effort. Full details of the MSSF (and related OSSF) approach are
given in the section on :ref:`multi-site support functions
<basis_mssf>`.
For the blip basis functions, the spacing of the grid where the blips
are defined is key (``Atom.SupportGridSpacing`` in
the :ref:`atomic specification <input_atomic_spec>` block),
and is directly related to an equivalent plane
wave cutoff (via :math:`k_{bg} = \pi/\delta` and :math:`E_{PW} =
k_{bg}^2/2`, where :math:`\delta` is the grid spacing in Bohr radii
and :math:`E_{PW}` is in Hartrees). For a particular grid spacing,
the energy will converge monotonically with support function radius
(``Atom.SupportFunctionRange`` in
the :ref:`atomic specification <input_atomic_spec>` block).
A small support function radius will introduce some approximation to
the result, but improve computational performance. It is vital to
characterise both blip grid spacing and support function radius in any
calculation. A full discussion of the blip function basis is found
:ref:`here <basis_blips>`.
Go to :ref:`top <convergence>`.

5
docs/CONQUEST-manual/errors.rst
View File

@ -4,4 +4,7 @@
Error codes
===========
These are to be implemented, but when done, this file will list them all...
Error codes generated by different CONQUEST routines will be collated
here, along with explanations of the root cause and suggested fixes.
This will form part of the full release, but is not implemented in the
pre-release.

294
docs/CONQUEST-manual/examples.rst Normal file
View File

@ -0,0 +1,294 @@
.. _examples:
====================
Example calculations
====================
All example calculations here use diagonalisation and PAO basis sets
(with a simple one-to-one mapping between PAOs and support functions).
.. _ex_static:
Static calculation
------------------
We will perform a self-consistent electronic structure calculation on
bulk silicon. The coordinate file that is needed is:
::
10.36 0.00 0.00
0.00 10.36 0.00
0.00 0.00 10.36
8
0.000 0.000 0.000 1 T T T
0.500 0.500 0.000 1 T T T
0.500 0.000 0.500 1 T T T
0.000 0.500 0.500 1 T T T
0.250 0.250 0.250 1 T T T
0.750 0.750 0.250 1 T T T
0.250 0.750 0.750 1 T T T
0.750 0.250 0.750 1 T T T
You should save this in an appropriate file (e.g. ``coords.dat``).
The inputs for the ion file can be found in ``pseudo-and-pao/PBE/Si``
(for the PBE functional). Changing to that directory and running the
``MakeIonFiles`` utility (in ``tools``) will generate the file
``SiCQ.ion``, which should be copied to the run directory, and renamed
to ``Si.ion``. The ``Conquest_input`` file requires only a few simple
lines at its most basic:
::
AtomMove.TypeOfRun static
IO.Coordinates coords.dat
Grid.GridCutoff 50
Diag.MPMesh T
Diag.GammaCentred T
Diag.MPMeshX 2
Diag.MPMeshY 2
Diag.MPMeshZ 2
General.NumberOfSpecies 1
%block ChemicalSpeciesLabel
1 28.086 Si
%endblock
The parameters above should be relatively self-explanatory; the grid
cutoff (in Hartrees) sets the integration grid spacing, and can be
compared to the *charge density* grid cutoff in a plane wave code
(typically four times larger than the plane wave cutoff). The
Monkhorst-Pack k-point mesh (``Diag.MPMeshX/Y/Z``) is a standard
feature of solid state codes; note that the grid can be forced to be
centred on the Gamma point.
The most important parameters set the number of species and give
details of what the species are (``ChemicalSpeciesLabel``). For each
species label (in this case ``Si``) there should be a corresponding
file with the extension ``.ion`` (again, in this case ``Si.ion``).
CONQUEST will read the necessary information from this file for
default operation, so no further parameters are required. This block
also allows the mass of the elements to be set (particularly important
for molecular dynamics runs).
The output file starts with a summary of the calculation requested,
including parameters set, and gives details of papers that are
relevant to the particular calculation. After brief details of the
self-consistency, the total energy, forces and stresses are printed,
followed by an estimate of the memory and time required. For this
calculation, these should be close to the following:
::
Harris-Foulkes Energy : -33.792210321858057 Ha
Atom X Y Z
1 -0.0000000000 0.0000000000 0.0000000000
2 -0.0000000000 0.0000000000 0.0000000000
3 -0.0000000000 0.0000000000 -0.0000000000
4 0.0000000000 0.0000000000 0.0000000000
5 -0.0000000000 0.0000000000 -0.0000000000
6 0.0000000000 0.0000000000 0.0000000000
7 -0.0000000000 0.0000000000 0.0000000000
8 -0.0000000000 -0.0000000000 0.0000000000
Maximum force : 0.00000000(Ha/a0) on atom, component 2 3
X Y Z
Total stress: -0.01848219 -0.01848219 -0.01848219 Ha
Total pressure: 0.48902573 0.48902573 0.48902573 GPa
The output file ends with an estimate of the total memory and time
used.
You might like to experiment with the grid cutoff to see how the
energy converges (note that the
number of grid points is proportional to the square root of the energy,
while the spacing is proportional to one over this, and
that the computational effort will scale with the *cube* of the number
of grid points); as with all DFT
calculations, you should ensure that you test the convergence with
respect to all parameters.
Go to :ref:`top <examples>`.
.. _ex_relax:
Relaxation
----------
.. _ex_relax_atoms:
Atomic Positions
~~~~~~~~~~~~~~~~
We will explore structural optimisation of the methane molecule (a
very simple example). The coordinates required are:
::
20.000 0.000 0.000
0.000 20.000 0.000
0.000 0.000 20.000
5
0.500 0.500 0.500 1 F F F
0.386 0.500 0.500 2 T F F
0.539 0.607 0.500 2 T T F
0.537 0.446 0.593 2 T T T
0.537 0.446 0.407 2 T T T
The size of the simulation cell should, of course, be tested carefully
to ensure that there are no interactions between images. We have
fixed the central (carbon) atom, and restricted other atoms to prevent
rotations or translations during optimisation.
The ``Conquest_input`` file changes only a little from before, as
there is no need to specify a reciprocal space mesh (it defaults to
gamma point only, which is appropriate for an isolated molecule). We
have set the force tolerance (``AtomMove.MaxForceTol``) to a
reasonable level (approximately 0.026 eV/A). Note that the ion files
can be generated in the same way :ref:`as before <ex_static>`, and
that we assume that the ion files are renamed to ``C.ion`` and ``H.ion``.
::
IO.Coordinates CH4.in
Grid.GridCutoff 50
AtomMove.TypeOfRun lbfgs
AtomMove.MaxForceTol 0.0005
General.NumberOfSpecies 2
%block ChemicalSpeciesLabel
1 12.00 C
2 1.00 H
%endblock
The progress of the optimisation can be followed by searching for the
string ``Geom`` (using ``grep`` or something similar). In this case,
we find:
::
GeomOpt - Iter: 0 MaxF: 0.04828504 E: -0.83676760E+01 dE: 0.00000000
GeomOpt - Iter: 1 MaxF: 0.03755566 E: -0.83755762E+01 dE: 0.00790024
GeomOpt - Iter: 2 MaxF: 0.02691764 E: -0.83804002E+01 dE: 0.00482404
GeomOpt - Iter: 3 MaxF: 0.00613271 E: -0.83860469E+01 dE: 0.00564664
GeomOpt - Iter: 4 MaxF: 0.00126136 E: -0.83862165E+01 dE: 0.00016958
GeomOpt - Iter: 5 MaxF: 0.00091560 E: -0.83862228E+01 dE: 0.00000629
GeomOpt - Iter: 6 MaxF: 0.00081523 E: -0.83862243E+01 dE: 0.00000154
GeomOpt - Iter: 7 MaxF: 0.00073403 E: -0.83862303E+01 dE: 0.00000603
GeomOpt - Iter: 8 MaxF: 0.00084949 E: -0.83862335E+01 dE: 0.00000316
GeomOpt - Iter: 9 MaxF: 0.00053666 E: -0.83862353E+01 dE: 0.00000177
GeomOpt - Iter: 10 MaxF: 0.00033802 E: -0.83862359E+01 dE: 0.00000177
The maximum force reduces smoothly, and the structure converges well.
By adjusting the output level (using ``IO.Iprint`` for overall output,
or ``IO.Iprint_MD`` for atomic movement) more information about the
structural relaxation can be produced (for instance, the force
residual and some details of the line minimisation will be printed for
``IO.Iprint_MD 2``).
Go to :ref:`top <examples>`.
.. _ex_relax_cell:
Cell Parameters
~~~~~~~~~~~~~~~
We will optimise the lattice constant of the bulk silicon cell that we
studied for the static calculation. Here we need to change the type
of run, and add one more line:
::
AtomMove.TypeOfRun cg
AtomMove.OptCell T
Adjust the simulation cell size to 10.26 Bohr radii in all three
directions (to make it a little more challenging). If you run this
calculation, you should find a final lattice constant
of 10.372 after 3 iterations. The progress of the optimization can be
followed in the same way as for structural relaxation, and gives:
::
GeomOpt - Iter: 0 MaxStr: 0.00011072 H: -0.33790200E+02 dH: 0.00000000
GeomOpt - Iter: 1 MaxStr: 0.00000195 H: -0.33792244E+02 dH: 0.00204424
GeomOpt - Iter: 2 MaxStr: 0.00000035 H: -0.33792244E+02 dH: -0.00000017
Go to :ref:`top <examples>`.
.. _ex_md:
Simple Molecular Dynamics
-------------------------
We will perform NVE molecular dynamics for methane, CH4, as a simple
example of how to do this kind of calculation. You should use the
same coordinate file and ion files as you did for the structural
relaxation, but change the atomic movement flags in the coordinate
file to allow all atoms to move (the centre of mass is fixed during MD
by default). Your coordinate file should look like this:
::
20.00000000000000 0.00000000000000 0.00000000000000
0.00000000000000 20.00000000000000 0.00000000000000
0.00000000000000 0.00000000000000 20.00000000000000
5
0.500 0.500 0.500 1 T T T
0.386 0.500 0.500 2 T T T
0.539 0.607 0.500 2 T T T
0.537 0.446 0.593 2 T T T
0.537 0.446 0.407 2 T T T
The input file should be:
::
IO.Coordinates CH4.in
AtomMove.TypeOfRun md
AtomMove.IonTemperature 300
AtomMove.NumSteps 100
General.NumberOfSpecies 2
%block ChemicalSpeciesLabel
1 12.00 C
2 1.00 H
%endblock
where the default timestep (0.5fs) is necessary for simulations
involving light atoms like hydrogen. The file ``md.stats`` contains
details of the simulation, while the trajectory is output to
``trajectory.xsf`` which can be read by VMD among other programs. In
this simulation, the conserved quantity is the total energy (the sum
of ionic kinetic energy and potential energy of the system) which is
maintained to better than 0.1mHa in this instance. More importantly,
the variation in this quantity is much smaller than the variation in
the potential energy. This can be seen in the plot below.
.. image:: MDPlot.png
Go to :ref:`top <examples>`.
.. _ex_tut:
Tutorials
---------
We recommend that you work through, in order, the tutorials included
in the distribution in the ``tutorials/`` directory
to become familiar with the modes of operation of the code.
**NOTE** In the initial pre-release of CONQUEST (January 2020) we have
not included the tutorials; they will be added over the coming months.
Go to :ref:`top <examples>`.
.. _ex_next:
Where next?
-----------
While the tutorials have covered the basic operations of Conquest,
there are many more subtle questions and issues, which are given in
the User Guide.
Go to :ref:`top <examples>`.

426
docs/CONQUEST-manual/ext-tools.rst Normal file
View File

@ -0,0 +1,426 @@
.. _ext-tools:
==============
External tools
==============
.. _et_md_scripts:
Molecular dynamics analysis
---------------------------
Several scripts that may be helpful with postprocessing molecular dynamics are
included with CONQUEST. The can be found in the ``tools`` directory, and the
executables are ``plot_stats.py``, ``md_analysis.py`` and ``heat_flux.py``. They
have the following dependencies:
* Python 3
* Scipy/Numpy
* Matplotlib
If Python 3 is installed the modules can be added easily using ``pip3 install
scipy`` etc.
These scripts should be run in the calculation directory, and will automatically
parse the necessary files, namely ``Conquest_input``, ``input.log``,
``md.stats`` and ``md.frames`` assuming they have the default names. They will
also read the CONQUEST input flags to determine, for example, what ensemble is
used, and process the results accordingly.
Go to :ref:`top <ext-tools>`.
.. _et_plot_stat:
Plotting statistics
+++++++++++++++++++
::
usage: plot_stats.py [-h] [-c] [-d DIRS [DIRS ...]]
[--description DESC [DESC ...]] [--skip NSKIP]
[--stop NSTOP] [--equil NEQUIL] [--landscape]
[--mser MSER_VAR]
Plot statistics for a CONQUEST MD trajectory
optional arguments:
-h, --help show this help message and exit
-c, --compare Compare statistics of trajectories in directories
specified by -d (default: False)
-d DIRS [DIRS ...], --dirs DIRS [DIRS ...]
Directories to compare (default: .)
--description DESC [DESC ...]
Description of graph for legend (only if using
--compare) (default: )
--skip NSKIP Number of equilibration steps to skip (default: 0)
--stop NSTOP Number of last frame in analysis (default: -1)
--equil NEQUIL Number of equilibration steps (default: 0)
--landscape Generate plot with landscape orientation (default:
False)
--mser MSER_VAR Compute MSER for the given property (default: None)
Running ``plot_stats.py --skip 200`` in your calculation will generate a plot
which should resemble the example below, skipping the first 200 steps. This
example is a molecular dynamics simulation of 1000 atoms of bulk silicon in the
NPT ensemble, at 300 K and 0.1 GPa.
.. image:: stats.jpg
The four plots are respectively the breakdown of energy contributions, the
conserved quantity, the temperature and the pressure, the last of which is only
included for NPT molecular dynamics. Several calculations in different
directories can be compared using ``plot_stats.py --compare -d dir1
dir2 --description "dir1 description" "dir2 description"``. The following
example compares the effect of changing the L tolerance in the above simulation.
Note that the contents of the description field will be in the legend of the
plot.
.. image:: compare.jpg
Go to :ref:`top <ext-tools>`.
.. _et_md_ana:
MD analysis
+++++++++++
::
usage: md_analysis.py [-h] [-d DIRS [DIRS ...]] [--skip NSKIP]
[--stride STRIDE] [--snap SNAP] [--stop NSTOP]
[--equil NEQUIL] [--vacf] [--msd] [--rdf] [--stress]
[--nbins NBINS] [--rdfwidth RDFWIDTH] [--rdfcut RDFCUT]
[--window WINDOW] [--fitstart FITSTART] [--dump]
Analyse a CONQUEST MD trajectory
optional arguments:
-h, --help show this help message and exit
-d DIRS [DIRS ...], --dirs DIRS [DIRS ...]
Directories to compare (default: .)
--skip NSKIP Number of equilibration steps to skip (default: 0)
--stride STRIDE Only analyse every nth step of frames file (default:
1)
--snap SNAP Analyse Frame of a single snapshot (default: -1)
--stop NSTOP Number of last frame in analysis (default: -1)
--equil NEQUIL Number of equilibration steps (default: 0)
--vacf Plot velocity autocorrelation function (default:
False)
--msd Plot mean squared deviation (default: False)
--rdf Plot radial distribution function (default: False)
--stress Plot stress (default: False)
--nbins NBINS Number of histogram bins (default: 100)
--rdfwidth RDFWIDTH RDF histogram bin width (A) (default: 0.05)
--rdfcut RDFCUT Distance cutoff for RDF in Angstrom (default: 8.0)
--window WINDOW Window for autocorrelation functions in fs (default:
1000.0)
--fitstart FITSTART Start time for curve fit (default: -1.0)
--dump Dump secondary data used to generate plots (default:
False)
The script ``md_analysis.py`` script performs various analyses of the trajectory
by parsing the `md.frames`` file. So far, these include the radial distribution
function, the velocity autocorrelation function, the mean squared deviation, and
plotting the stress. For example, the command,
``md_analysis.py --rdf --stride 20 --rdfcut 8.0 --nbins 100 --dump --skip 200 --stop 400``
computes the radial distribution function of the simulation in the first example
from every 20th time step (every 10 fs in this case), stopping after 400 steps,
with a cutoff of 8.0 A, and the histogram is divided into 100 bins.
.. image:: rdf.jpg
Go to :ref:`top <ext-tools>`.
.. _et_cq_struc:
CONQUEST structure file analysis
++++++++++++++++++++++++++++++++
::
usage: structure.py [-h] [-i INFILE] [--bonds] [--density] [--nbins NBINS]
[-c CUTOFF [CUTOFF ...]] [--printall]
Analyse a CONQUEST-formatted structure
optional arguments:
-h, --help show this help message and exit
-i INFILE, --infile INFILE
CONQUEST format structure file (default:
coord_next.dat)
--bonds Compute average and minimum bond lengths (default:
False)
--density Compute density (default: False)
--nbins NBINS Number of histogram bins (default: 100)
-c CUTOFF [CUTOFF ...], --cutoff CUTOFF [CUTOFF ...]
Bond length cutoff matrix (upper triangular part, in
rows (default: None)
--printall Print all bond lengths (default: False)
The script ``structure.py`` can be used to analyse a CONQUEST-formatted
structure file. This is useful to sanity-check the bond lengths or density,
since an unphysical structure is so often the cause of a crash. For example, the
bond lengths can be computed with
``structure.py --bonds -c 2.0 3.0 3.0``
where the ``-c`` flag specifies the bond cutoffs for the bonds 1-1, 1-2 and 2-2,
where 1 is species 1 as specified in ``Conquest_input`` and 2 is species 2. The
output will look something like this:
::
Mean bond lengths:
O-Si: 1.6535 +/- 0.0041 (24 bonds)
Minimum bond lengths:
O-Si: 1.6493
Go to :ref:`top <ext-tools>`.
.. _et_ase:
Atomic Simulation Environment (ASE)
-----------------------------------
ASE is a set of Python tools for setting up, manipulating, running, visualizing
and analyzing atomistic simulations. ASE contains a CONQUEST interface, so that
it can be used to calculate energies, forces and stresses for calculations that
CONQUEST can't do (yet). Detailed instructions on how to install and invoke it
can be found on its `website <https://wiki.fysik.dtu.dk/ase/>`_, but we provide
some details and examples for the CONQUEST interface here.
Note that the script will need to set environmental variables specifying the
locations of the CONQUEST executable ``Conquest``, and if required, the basis
set generation executable ``MakeIonFiles`` and pseudopotential database.
::
import os
# The command to run CONQUEST in parallel
os.environ["ASE_CONQUEST_COMMAND"] = "mpirun -np 4 /path/to/Conquest_master"
# Path to a database of pseudopotentials (for basis generation tool)
os.environ["CQ_PP_PATH"] = "~/Conquest/PPDB/"
# Path to the basis generation tool executable
os.environ["CQ_GEN_BASIS_CMD"] = "/path/to/MakeIonFiles"
Go to :ref:`top <ext-tools>`.
.. _et_ase_input:
Keywords for generating the Conquest_input file
+++++++++++++++++++++++++++++++++++++++++++++++
The calculator object contains a dictionray containing a small number of
mandatory keywords, listed below:
::
default_parameters = {
'grid_cutoff' : 100, # DFT defaults
'kpts' : None,
'xc' : 'PBE',
'scf_tolerance' : 1.0e-6,
'nspin' : 1,
'general.pseudopotentialtype' : 'Hamann', # CONQUEST defaults
'basis.basisset' : 'PAOs',
'io.iprint' : 2,
'io.fractionalatomiccoords' : True,
'mine.selfconsistent' : True,
'sc.maxiters' : 50,
'atommove.typeofrun' : 'static',
'dm.solutionmethod' : 'diagon'}
The first five key/value pairs are special DFT parameters, the grid cutoff, the
k-point mesh, the exchange-correlation functional, the SCF tolerance and the
number of spins respectively. The rest are CONQUEST-specific input flags.
The atomic species blocks are handled slightly differently, with a dictionary of
their own. If the ``.ion`` files are present in the calculation directory, they
can be specified as follows:
::
basis = {"H": {"valence_charge": 1.0,
"number_of_supports": 1,
"support_fn_range": 6.9},
"O": {"valence_charge": 6.0,
"number_of_supports": 4,
"support_fn_range": 6.9}}
If the basis set ``.ion`` files are present in the directory containing the ASE
script are pressent and are named ``element.ion``, then the relevant parameters
will be parsed from the ``.ion`` files and included when the input file is
written and this dictionary can be omitted. It is more important when, for
example, setting up a multisite calculation, when the number of contracted
support functions is different from the number in the ``.ion`` file.
ASE can also invoke the CONQUEST basis set generation tool, although care should
be taken when generating basis sets:
::
basis = {"H": {"basis_size": "minimal",
"pseudopotential_type": hamann",
"gen_basis": True},
"O": {"basis_size": "minimal",
"pseudopotential_type": hamann",
"gen_basis": True}}
Finally, non-mandatory input flags can be defined in a new dictionary, and
passed as an expanded set of keyword arguments.
::
conquest_flags = {'IO.Iprint' : 1, # CONQUEST keywords
'DM.SolutionMethod' : 'ordern',
'DM.L_range' : 8.0,
'minE.LTolerance' : 1.0e-6}
Here is an example, combining the above. We set up a cubic diamond cell
containing 8 atoms, and perform a single point energy calculation using the
order(N) method (the default is diagonalisation, so we must specify all of the
order(N) flags). We don't define a basis set, instead providing keywords that
specify that a minimal basis set should be constructed using the MakeIonFiles
basis generation tool.
::
from ase.build import bulk
from ase.calculators.conquest import Conquest
os.environ["ASE_CONQUEST_COMMAND"] = "mpirun -np 4 Conquest_master"
os.environ["CQ_PP_PATH"] = "/Users/zamaan/Conquest/PPDB/"
os.environ["CQ_GEN_BASIS_CMD"] = "MakeIonFiles"
diamond = bulk('C', 'diamond', a=3.6, cubic=True) # The atoms object
conquest_flags = {'IO.Iprint' : 1, # Conquest keywords
'DM.SolutionMethod' : 'ordern',
'DM.L_range' : 8.0,
'minE.LTolerance' : 1.0e-6}
basis = {'C': {"basis_size" : 'minimal', # Generate a minimal basis
"gen_basis" : True,
"pseudopotential_type" : "hamann"}}
calc = Conquest(grid_cutoff = 80, # Set the calculator keywords
xc="LDA",
self_consistent=True,
basis=basis,
nspin=1,
**conquest_flags)
diamond.set_calculator(calc) # attach the calculator to the atoms object
energy = diamond.get_potential_energy() # calculate the potential energy
Go to :ref:`top <ext-tools>`.
.. _et_ase_mssf:
Multisite support functions
+++++++++++++++++++++++++++
Multisite support functions require a few additional keywords in the atomic
species block, which can be specified as follows:
::
basis = {'C': {"basis_size": 'medium',
"gen_basis": True,
"pseudopotential_type": "hamann",
"Atom.NumberofSupports": 4,
"Atom.MultisiteRange": 7.0,
"Atom.LFDRange": 7.0}}
Note that we are constructing a DZP basis set (size medium) with 13 primitive
support functions using ``MakeIonFiles``, and contracting it to multisite basis
of 4 support functions. The calculation requires a few more input flags, which
are specified in the ``other_keywords`` dictionary:
::
other_keywords = {"Basis.MultisiteSF": True,
"Multisite.LFD": True,
"Multisite.LFD.Min.ThreshE": 1.0e-7,
"Multisite.LFD.Min.ThreshD": 1.0e-7,
"Multisite.LFD.Min.MaxIteration": 150,
}
Go to :ref:`top <ext-tools>`.
.. _et_ase_load_dm:
Loading the K/L matrix
++++++++++++++++++++++
Most calculation that involve incrementally moving atoms (molecular dynamics,
geometry optimisation, equations of state, nudged elastic band etc.) can be made
faster by using the K or L matrix from a previous calculation as the initial
guess for a subsequent calculation in which that atoms have been moved slightly.
This can be achieved by first performing a single point calculation to generate
the first K/L matrix, then adding the following keywords to the calculator:
::
other_keywords = {"General.LoadL": True,
"SC.MakeInitialChargeFromK": True}
These keywords respectively cause the K or L matrix to be loaded from file(s)
``Kmatrix.i**.p*****``, and the initial charge density to be constructed from
this matrix. In all subsequent calculations, the K or L matrix will be written
at the end of the calculation and used as the initial guess for the subsequent
ionic step.
Go to :ref:`top <ext-tools>`.
.. _et_eos:
Equation of state
+++++++++++++++++
The following code computes the equation of state of diamond by doing single
point calculations on a uniform grid of the ``a`` lattice parameter. It then
interpolates the equation of state and uses ``matplotlib`` to generate a plot.
::
import scipy as sp
from ase.build import bulk
from ase.io.trajectory import Trajectory
from ase.calculators.conquest import Conquest
# Construct a unit cell
diamond = bulk('C', 'diamond', a=3.6, cubic=True)
basis = {'C': {"basis_size": 'minimal',
"gen_basis": True,
"pseudopotential_type": "hamann"}}
calc = Conquest(grid_cutoff = 50,
xc = "LDA",
basis = basis,
kpts = [4,4,4]}
diamond.set_calculator(calc)
cell = diamond.get_cell()
traj = Trajectory('diamond.traj', 'w') # save all results to trajectory
for x in sp.linspace(0.95, 1.05, 5): # grid for equation of state
diamond.set_cell(cell*x, scale_atoms=True)
diamond.get_potential_energy()
traj.write(diamond)
from ase.io import read
from ase.eos import EquationOfState
configs = read('diamond.traj@0:5')
volumes = [diamond.get_volume() for diamond in configs]
energies = [diamond.get_potential_energy() for diamond in configs]
eos = EquationOfState(volumes, energies)
v0, e0, B = eos.fit()
import matplotlib
eos.plot('diamond-eos.pdf') # Plot the equation of state
Go to :ref:`top <ext-tools>`.

73
docs/CONQUEST-manual/faq.rst Normal file
View File

@ -0,0 +1,73 @@
.. _faq:
==========================
Frequently Asked Questions
==========================
.. _faq_when_cq:
When should I use CONQUEST?
---------------------------
You can use CONQUEST for any DFT simulations that you need to
perform. It is efficient for small problems, though may not be as
efficient as other codes (e.g. plane wave codes) because it has been
designed for massively parallel operation, which brings some
overhead. If you need to perform DFT calculations on large systems
(several hundred atoms or beyond) or want to perform highly parallel
calculations, you should definitely consider CONQUEST.
CONQUEST uses `Hamann`_ optimised norm-conserving Vanderbilt (ONCV)
pseudopotentials, which can also be used by `PWSCF`_
and `Abinit`_ which allows direct comparisons between the codes.
.. _Hamann: http://www.mat-simresearch.com
.. _PWSCF: https://www.quantum-espresso.org
.. _Abinit: https://www.abinit.org
.. _faq_when_on:
When should I use linear scaling?
---------------------------------
You should use linear scaling if you need to model systems with more
than about 5,000 atoms, though gains are often found for smaller
systems (from 1,000 atoms upwards).
Linear scaling calculations offer the prospect of scaling to
significantly larger systems than traditional DFT calculations;
however, they make approximations and require some care and
characterisation. In particular, instead of solving for eigenvalues
and eigenstates, linear scaling methods solve for the density matrix,
so that energy-resolved information (e.g. DOS and band energies) are
not available. To enable linear scaling, a range is also imposed on the
density matrix and it is important to test the effect of this range.
.. _faq_implement:
Will you implement a specific feature for me?
---------------------------------------------
We cannot guarantee to implement specific features, though we are
always happy to take suggestions. We also welcome new developers: if
there is something that you would like to see in the code, please do
talk to us about joining the development effort.
.. _faq_bug:
How do I report a bug?
----------------------
Please use the `GitHub issues`_ page. Include details of the compiler
and libraries used, the version of CONQUEST, and the input and output
files (if possible). We will do our best to check the bug and fix it,
but cannot guarantee to help on any timescale.
.. _GitHub issues: http://github.com/OrderN/CONQUEST-release/issues
.. _faq_help:
How do I get help?
------------------
The Conquest mailing list (**details**) is the best place to get
help. However, the developers cannot guarantee to answer any
questions, though they will try. Bug reports should be made through
the `GitHub issues`_ page.
.. _GitHub issues: http://github.com/OrderN/CONQUEST-release/issues

439
docs/CONQUEST-manual/groundstate.rst Normal file
View File

@ -0,0 +1,439 @@
.. _groundstate:
========================
Finding the ground state
========================
Finding the electronic ground state is the heart of any DFT code. In
CONQUEST, we need to
consider several linked stages: the density
matrix (found using :ref:`diagonalisation <gs_diag>` or :ref:`linear scaling <gs_on>`);
:ref:`self-consistency between charge and potential <gs_scf>`; and the
:ref:`support functions <gs_suppfunc>` (though these are not always optimised).
The basis functions in CONQUEST are :ref:`support functions <gs_suppfunc>` (localised
functions centred on the atoms), written as
:math:`\phi_{i\alpha}(\textbf{r})` where :math:`i` indexes an atom and
:math:`\alpha` a support function on the atom. The support functions
are used as basis functions for the density matrix and the Kohn-Sham
eigenstates:
.. math::
\psi_{n\mathbf{k}}(\mathbf{r}) = \sum_{i\alpha} c^{n\mathbf{k}}_{i\alpha}
\phi_{i\alpha}(\mathbf{r})\\
\rho(\mathbf{r}, \mathbf{r}^\prime) = \sum_{i\alpha j\beta}
\phi_{i\alpha}(\mathbf{r}) K_{i\alpha, j\beta} \phi_{j\beta}(\mathbf{r}^\prime)
where :math:`n` is an eigenstate index and :math:`\mathbf{k}` is a
point in the Brillouin zone (see :ref:`here <gs_diag_bz>` for more on
this). The total energy can be written in terms of the density
matrix, as:
.. math::
E_{KS} = \mathrm{Tr}[HK] + \Delta E_{Har} + \Delta E_{XC}
for the Hamiltonian matrix :math:`H` in the basis of support
functions, with the last two terms the standard Harris-Foulkes
:cite:`g-Harris1985,g-Foulkes1989` correction terms.
For diagonalisation, the density matrix is made from the coefficients
of the Kohn-Sham eigenstates, :math:`c^{n\mathbf{k}}_{i\alpha}`, while
for :ref:`linear scaling <gs_on>` it is found directly during the variational
optimisation of the energy.
The question of whether to find the density matrix via diagonalisation
or linear scaling is a complex one, depending on the system size,
the accuracy required and the computational resources available. The
simplest approach is to test diagonalisation before linear scaling.
.. _gs_diag:
Diagonalisation
---------------
Exact diagonalisation in CONQUEST uses the ScaLAPACK library which
scales reasonably well in parallel, but becomes less efficient with
large numbers of processes. The computational time
will scale as :math:`N^3` with the number of atoms :math:`N`, but will
probably be more efficient than linear scaling for systems up to a few
thousand atoms. (Going beyond a thousand atoms with diagonalisation
is likely to require the :ref:`multi-site support function
<basis_mssf>` technique.)
To choose diagonalisation, the following flag should be set:
::
DM.SolutionMethod diagon
It is also essential to test relevant parameters, as described below:
the k-point grid in reciprocal space (to sample the Brillouin zone
efficiently); the occupation smearing approach; and the
parallelisation of k-points.
Go to :ref:`top <groundstate>`
.. _gs_diag_bz:
Brillouin zone sampling
~~~~~~~~~~~~~~~~~~~~~~~
We need to specify a set of discrete points in reciprocal space to
approximate integrals over the Brillouin zone. The simplest approach
is to use the Monkhorst-Pack approach :cite:`g-Monkhorst:1976kf`,
where a grid of points is specified in all directions:
::
Diag.MPMesh T
Diag.MPMeshX 2
Diag.MPMeshY 2
Diag.MPMeshZ 2
This grid can be forced to be centred on the gamma point (often an
important point) using the parameter ``Diag.GammaCentred T``.
The origin of the Monkhorst-Pack grid may also be offset by an
arbitrary vector from the origin of the Brillouin zone, by specifying:
::
Diag.MPShiftX 0.0
Diag.MPShiftY 0.0
Diag.MPShiftZ 0.0
Alternatively, the points in reciprocal space can be specified
explicitly by giving a number of points and their locations and weights:
::
Diag.NumKpts 1
%block Diag.Kpoints
0.00 0.00 0.00 1.00
%endblock Diag.Kpoints
where there must be as many lines in the block as there are k-points.
It is important to note that CONQUEST does not consider space group
symmetry when integrating over the Brillouin zone.
Go to :ref:`top <groundstate>`.
.. _gs_diag_para:
K-point parallelization
~~~~~~~~~~~~~~~~~~~~~~~~
It is possible to parallelise over k-points: to split the processes
into sub-groups, each of which is responsible for a sub-set of the
k-points. This can be very efficient, and is specified by the
parameter ``Diag.KProcGroups N``, where it is important that the number
of processes is an integer multiple of the number of groups ``N``. It
will be most efficient when the number of k-points is an integer
multiple of the number of groups.
Go to :ref:`top <groundstate>`.
.. _gs_diag_smear:
Electronic occupation smearing
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The occupation numbers of the eigenstates are slightly smeared near
the Fermi level, following common practice. The default smearing type
is Fermi-Dirac smearing with a temperature (in Hartrees) set with the
flag ``Diag.kT`` which defaults to 0.001Ha.
The Methfessel-Paxton approach :cite:`g-Methfessel:1989ny` to occupations allows much higher
smearing temperatures with minimal effect on the free energy (and
hence accuracy) of the energy. This generally gives a similar accuracy
with fewer k-points, and is selected as:
::
Diag.SmearingType 1
Diag.MPOrder 0
where ``Diag.MPOrder`` specifies the order of the Methfessel-Paxton
expansion. It is recommended to start with the lowest order and
increase gradually, testing the effects.
Go to :ref:`top <groundstate>`.
.. _gs_on:
Linear Scaling
--------------
A linear scaling calculation is selected by setting
``DM.SolutionMethod ordern``. There are two essential parameters that must be
set: the range of the density matrix, and the tolerance on the
optimisation.
::
DM.L_range 16.0
minE.Ltolerance 1.0e-6
The tolerance is applied to the residual (the RMS value of the
gradient of the energy with respect to the density matrix). The
maximum number of iterations in the density matrix optimisation can
be set with ``DM.LVariations`` (default 50).
At present, CONQUEST can only operate efficiently in linear scaling
mode with a restricted number of support functions (though this is an
area of active development). PAO basis sets of SZ and SZP size
(minimal and small in the ion file generator) will run without
restrictions. For larger PAO basis sets, the :ref:`OSSF <basis_ossf>`
approach must be used, and is effective. With a blip basis there are
no restrictions, though efficient optimisation is still under active
development.
It is
almost always more efficient to update the charge density while
optimising the density matrix, avoiding the need for a separate
self-consistency loop. This is set by choosing
``minE.MixedLSelfConsistent T``.
An essential part of a linear scaling calculation is finding the
approximate, sparse inverse of the overlap matrix. Normally this will
happen automatically, but it may require some tests. The key
parameters are the range for the inverse (see the
:ref:`input_atomic_spec` block, and specifically the
:ref:`advanced_atomic_spec_tags` block) and the tolerance applied
to the inversion.
::
Atom.InvSRange R
DM.InvSTolerance R
A tolerance of up to 0.2 can give convergence without significantly
affecting the accuracy. The range should be similar to the radius of
the support functions, though increasing it by one or two bohr can
improve the inversion in most cases.
The input tags are mainly found in the :ref:`input_dm` section of the
:ref:`input_tags` page.
Go to :ref:`top <groundstate>`.
.. _gs_scf:
Self-consistency
----------------
The normal mode of operation for CONQUEST involves an iterative search
for self-consistency between the potential and the charge density.
However, it is also possible to run in a non-self-consistent manner,
either with a converged charge density for electronic structure
analysis, or for dynamics, which will be considerably more efficient
than a self-consistent calculation, but less accurate.
Self consistency is set via the following parameters:
::
minE.SelfConsistent T
minE.SCTolerance 1E-7
SC.MaxIters 50
The tolerance is applied to the RMS value of the residual,
:math:`R(\mathbf{r}) = \rho^{out}(\mathbf{r}) - \rho^{in}(\mathbf{r})`,
integrated over all space:
.. math::
R_{RMS} = \sqrt{\Omega \sum_l \left(R(\mathbf{r}_l)\right)^2 }
where :math:`\mathbf{r}_l` is a grid point and :math:`\Omega` is the
grid point volume (integrals are performed
on a grid explained in :ref:`conv_grid`). The maximum number
of self-consistency cycles is set with ``SC.MaxIters``, defaulting
to 50.
For non-self-consistent calculations, the main flag should be set as
``minE.SelfConsistent F``. The charge density at each step will
either be read from a file (if the flag ``General.LoadRho T`` is set),
or constructed from a superposition of
atomic densities. The Harris-Foulkes functional will be used to
find the energy.
Go to :ref:`top <groundstate>`.
.. _ gs_scf_adv:
Advanced options
~~~~~~~~~~~~~~~~
Instabilities during self-consistency are a well-known issue in
electronic structure calculations. CONQUEST performs charge mixing
using the Pulay approach, where the new charge density is prepared by
combining the charge densities from a number of previous iterations.
In general, we write:
.. math::
rho_{n+1}^{in} = \sum_{i} \alpha_i \left[ \rho_{i}^{in} + A R_{i}
\right]
where :math:`R_{i}` is the residual at iteration :math:`i`, defined above. The
fraction of the output charge density that is included is governed by
the variable :math:`A`, which is set by the parameter
``SC.LinearMixingFactor`` (default 0.5). If there is instability
during the self consistency, reducing :math:`A` can help (though will likely
make convergence a little slower).
It is also advisable to apply Kerker preconditioning to the residual
when the system is large in any dimension. This removes long
wavelength components of the residual, reducing charge sloshing. This
is controlled with the following parameters:
::
SC.KerkerPreCondition T
SC.KerkerFactor 0.1
where the Kerker factor gives the wavevector at which preconditioning
starts to reduce. The Kerker preconditioning is applied to the
Fourier transform of the residual, :math:`\tilde{R}` as:
.. math::
\tilde{R} \frac{q^2}{q^2 - q^2_0}
where :math:`q^2_0` is the square of the Kerker factor and :math:`q` is a
wavevector. You should test values of :math:`q_0` around
:math:`\pi/a` where :math:`a` is the longest dimension of the simulation
cell (or some important length scale in your system).
Go to :ref:`top <groundstate>`.
.. _gs_suppfunc:
Support functions
-----------------
Support functions in CONQUEST represent the density matrix, and can be
simple (pseudo-atomic orbitals, or PAOs) or compound, made from simple
functions (either PAOs or blips). If they are compound, made from other
functions, then the search for the ground state involves the
construction of this representation. Full details of how the support
functions are built and represented can be found in the manual section on
:ref:`basis sets <basissets>`.
Go to :ref:`top <groundstate>`
.. _gs_charged:
Charged systems
-----------------
CONQUEST uses periodic boundary conditions, which require overall
charge neutrality. However, charged systems can be modelled:
if an excess of electrons is specified by the user, a uniform
positive background charge is added automatically to restore overall
neutrality. At present, there are no correction schemes implemented,
so it is important to test the convergence of the energy with unit
cell size and shape. Electrons are added by setting the parameter
``General.NetCharge``.
::
General.NetCharge 1.0
This gives the number of extra electrons to be added to the unit cell,
beyond the valence electrons.
Go to :ref:`top <groundstate>`.
.. _gs_spin:
Spin polarisation
-----------------
CONQUEST performs collinear spin calculations only. A spin-polarised
calculation is performed by setting the parameter
``Spin.SpinPolarised`` to T.
Users need to specify *either* the total initial number of spin-up and spin-down electrons in
the simulation cell (using the parameters ``Spin.NeUP`` and
``Spin.NeDN``), *or* the difference between the number of spin-up and
spin-down electrons (using the parameter ``Spin.Magn``).
The number of electrons for each spin channel can be fixed during SCF
calculations by setting the parameter ``Spin.FixSpin`` to T (default is F).
It is possible to specify the spin occupation in the atomic charge
densities (i.e. the number of spin-up and spin-down electrons used to
build the density). This is done in the :ref:`input_atomic_spec`
part of the ``Conquest_input`` file. Within the atom block for
each species, the numbers of electrons should be set with
``Atom.SpinNeUp`` and ``Atom.SpinNeDn``. Note that these numbers
*must* sum to the number of valence electrons for the atom.
Go to :ref:`top <groundstate>`.
.. _gs_spin_example:
Examples: FM and AFM iron
~~~~~~~~~~~~~~~~~~~~~~~~~
A two atom ferromagnetic iron simulation might be set up using the
parameters below. Note that the net spin here is S=1 :math:`\mu_B`
(i.e. two more electrons in the up channel than in the down), and
that the net spin is not constrained.
::
# example of ferro bcc Fe
Spin.SpinPolarised T
Spin.FixSpin F
Spin.NeUP 9.0 # initial numbers of up- and down-spin electrons,
Spin.NeDN 7.0 # which will be optimised by a SCF calculation when Spin.FixSpin=F
%block ChemicalSpeciesLabel
1 55.845 Fe
%endblock ChemicalSpeciesLabel
An equivalent anti-ferromagnetic calculation could be set up as
follows (though note that the initial specification of spin for the
atoms does *not* guarantee convergence to an AFM ground state). By
defining two species we can create spin-up and spin-down atoms (note
that both species will require their own, appropriately labelled, ion
file).
::
# example of anti-ferro bcc Fe
Spin.SpinPolarised T
Spin.FixSpin F
Spin.NeUP 8.0 # initial numbers of up- and down-spin electrons in an unit cell
Spin.NeDN 8.0 # are set to be the same
%block ChemicalSpeciesLabel
1 55.845 Fe1
2 55.845 Fe2
%endblock ChemicalSpeciesLabel
%block Fe1 # up-spin Fe
Atom.SpinNeUp 5.00
Atom.SpinNeDn 3.00
%endblock Fe1
%block Fe2 # down-spin Fe
Atom.SpinNeUp 3.00
Atom.SpinNeDn 5.00
%endblock Fe2
When using multi-site or on-site support functions in spin-polarised
calculations, the support functions can be made spin-dependent
(different coefficients for each spin channel) or not by setting
``Basis.SpinDependentSF`` (T/F, default is T).
Go to :ref:`top <groundstate>`.
.. bibliography:: references.bib
:cited:
:labelprefix: G
:keyprefix: g-
:style: unsrt
Go to :ref:`top <groundstate>`.

32
docs/CONQUEST-manual/important.rst
View File

@ -1,32 +0,0 @@
.. _important:
================
Important Topics
================
Diagonalisation
---------------
Self-consistency
----------------
Basis functions
---------------
Support functions
-----------------
Grids and integration
---------------------
Tolerances
----------
Functionals
-----------
Restarting
----------
Load balancing and partitions
-----------------------------

106
docs/CONQUEST-manual/index.rst
View File

@ -3,27 +3,101 @@
CONQUEST: a local orbital, large-scale DFT code
===============================================
CONQUEST is a massively parallel DFT code, using a local orbital basis
to represent the Kohn-Sham eigenstates. The code can operate in full
diagonalisation mode or linear scaling, with a variety of basis sets.
CONQUEST is a local orbital density functional theory (DFT) code,
capable of massively parallel operation with excellent scaling. It
uses a local orbital basis
to represent the Kohn-Sham eigenstates or the density matrix.
CONQUEST can be applied to atoms, molecules, liquids and solids, but
is particularly efficient for large systems. The code can find the ground
state using exact diagonalisation of the Hamiltonian or via a linear
scaling approach. The code has demonstrated scaling to over 2,000,000
atoms and 200,000 cores when using linear scaling, and over 3,400
atoms and 850 cores with exact diagonalisation.
CONQUEST can perform structural relaxation (including unit cell
optimisation) and molecular dynamics (in NVE, NVT and NPT ensembles
with a variety of thermostats).
Getting Started
---------------
* :doc:`why-conquest`
* :doc:`faq`
* :doc:`quick-overview`
* :doc:`installing`
* :doc:`examples`
.. toctree::
:maxdepth: 1
:hidden:
:caption: Getting Started
why-conquest
faq
quick-overview
installing
examples
User Guide
----------
* :doc:`input-output`
* :doc:`groundstate`
* :doc:`convergence`
* :doc:`basissets`
* :doc:`strucrelax`
* :doc:`moldyn`
* :doc:`ext-tools`
* :doc:`errors`
* :doc:`input_tags`
.. toctree::
:maxdepth: 1
:hidden:
:caption: User Guide
input-output
groundstate
convergence
basissets
strucrelax
moldyn
ext-tools
errors
input_tags
Theory
------
* :doc:`theory-strucrelax`
* :doc:`theory-md`
.. toctree::
:maxdepth: 2
:caption: Contents:
:maxdepth: 1
:hidden:
:caption: Theory
starting
input-output
important
errors
theory-strucrelax
theory-md
Get in touch
------------
Indices and tables
==================
- If you have suggestions for developing the code, please
use `GitHub issues`_. The developers cannot guarantee to offer
support, though we will try to help.
- Report bugs, suggest features or view the source code on `GitHub
issues`_.
- You can ask for help and discuss any problems that you may have on
the Conquest mailing list (**to be completed**)
.. only:: builder_html
.. _GitHub issues: http://github.com/OrderN/CONQUEST-release/issues
* :ref:`genindex`
* :ref:`search`
Licence
-------
.. only:: not builder_html
CONQUEST is available freely under the open source `MIT Licence`__.
We ask that you acknowledge use of the code by citing appropriate
papers, which will be given in the output file (a BiBTeX file
containing these references is also output).
* :ref:`genindex`
__ https://choosealicense.com/licenses/mit/

183
docs/CONQUEST-manual/input-output.rst
View File

@ -4,30 +4,213 @@
Input and output
================
.. _io_files:
Input files
-----------
.. _io_cq_in:
Conquest_input
++++++++++++++
All necessary input parameters should be specified in the ``Conquest_input`` file,
including the names of the coordinate file and the ion files. This
file controls the run; there are many sensible default values for
input parameters, but you should ensure that you understand what
they mean. After a run, the full set of relevant input parameters
(whether specified by the user, or default, are available in the file ``input.log``).
The most common input tags are listed briefly here.
Full documentation can be found in :ref:`input_tags`.
* ``AtomMove.TypeOfRun`` takes ``static``, ``md``, ``lbfgs``, ``cg``
* ``IO.Coordinates`` File name
* ``DM.SolutionMethod`` diagon
* ``Diag.MPMesh`` T/F
* ``Diag.MPMeshX`` (and ``Y`` and ``Z``) N
* ``Diag.GammaCentred`` T/F
* ``General.NumberOfSpecies`` N
* ``%block ChemicalSpeciesLabel`` Specifies element number, mass and
ion file name
* ``General.LoadDM`` T/F (to reload the density matrix)
* ``IO.FractionalAtomicCoords`` T/F
* ``Spin.SpinPolarised`` T/F
* ``Spin.FixSpin`` T/F
* ``Spin.Magn`` Difference between spin channel occupations
* ``Grid.GridCutoff`` Energy in Ha
* ``minE.SCTolerance`` Fractional tolerance
* ``AtomMove.NumSteps`` N
* ``AtomMove.MaxForceTol`` in Ha/bohr
* ``AtomMove.OptCell`` T/F
* ``SC.KerkerPreCondition`` T/F (for Kerker preconditioning of SCF)
* ``SC.MaxIters`` N (maximum number of SCF iterations
Go to :ref:`top <input-output>`
.. _io_ion:
Ion files
+++++++++
The ion files contain data on the different species being modelled:
valence charge, pseudopotentials, pseudo-atomic orbitals (PAOs) etc. Full
details on how the PAOs are used as basis functions for CONQUEST can
be found in the manual section on :ref:`basis sets <basissets>`. A
utility for generating these files is provided with CONQUEST, but Siesta ion
files can also be read. The CONQUEST utility uses the
pseudopotentials generated by the `ONCVPSP`_ code (though note that to
generate new files for CONQUEST, you will need a small patch).
A set of input files for all elements in the `PseudoDojo`_ library for
both LDA and PBE exchange-correlation functionals is provided in the
directory ``pseudo-and-pao``. This will allow you to generate ion
files for these elements easily.
The utility for generating ion files is called `MakeIonFiles`, and its
source code is found in the ``tools/BasisGeneration`` directory. It
uses the same ``system.make`` file as CONQUEST, and following
compilation the executable will be moved to the ``bin`` directory.
The key parameters to be set in the ``Conquest_input`` file for the
ion file generation are:
* ``General.NumberOfSpecies`` to specify number of species
* ``%block SpeciesLabels`` to specify what the species are
* In the species block (set with ``%block XX`` for species XX):
* ``Atom.PseudopotentialFile`` to specify the input file for the ONCVPSP code
* ``Atom.VKBFile`` to specify the file that CONQUEST needs to read
(included in the library of inputs)
* ``Atom.BasisSize`` to specify the size of the basis; at present
this can take the values: ``minimal``; ``small``; ``medium``; and
``large``.
Further fine-grained control can be applied to the basis functions;
this will be documented after the pre-release of CONQUEST.
.. _ONCVPSP: http://http://www.mat-simresearch.com
.. _PseudoDojo: https://www.pseudo-dojo.org/
Go to :ref:`top <input-output>`
.. _io_coords:
Coordinates
+++++++++++
The coordinates are specified in a separate file with relatively
simple format. The coordinates can be specified in fractional form
(default) or cartesian (set the input tag ``IO.FractionalAtomicCoords T``).
Distance units can be Bohr radii (default) or Angstroms (set the input tag
``General.DistanceUnits`` to ``Ang``). At present,
CONQUEST only handles *orthorhombic* unit cells.
The coordinate file is formatted as follows:
::
a 0.0 0.0
0.0 b 0.0
0.0 0.0 c
NAtoms
x y z species MoveX MoveY MoveZ
.
.
.
Note that the flags ``MoveX`` etc take values T/F and indicate whether
atoms are free to move in x, y and z, respectively. The flag
``species`` is an integer, and selects based on species defined in the
:ref:`atomic specification <input_atomic_spec>` section of the
``Conquest_input`` file.
Go to :ref:`top <input-output>`
.. _io_output:
Output files
------------
.. _io_output_main:
Main output
+++++++++++
By default, CONQUEST writes output to the ``Conquest_out`` file
(though the filename can be set with the parameter ``IO.OutputFile``,
and the flag ``IO.WriteOutToFile`` (T/F) selects output to file
or ``stdout``). This file contains all details of the calculation,
including energies, forces and information on the different stages of
the calculation. The output verbosity is controlled by the
``IO.Iprint`` family of parameters, which allows different levels of
output detail to be set for different areas of the code.
Note that, during the pre-release period of CONQUEST, we will be
tidying the output and making it more consistent.
Go to :ref:`top <input-output>`
.. _io_output_elec:
Electronic structure
++++++++++++++++++++
Different electronic structure outputs are available; in each case,
the key output flag is given. Further output flags are described in :ref:`input_tags`.
* Charge density
* Band-resolved charge density (``IO.outputWF``)
* Density of states (``IO.writeDOS``)
* Atom-projected density of states (``IO.write_proj_DOS``)
* Atomic charges, using the Mulliken approach (``IO.AtomChargeOutput``)
The Kohn-Sham eigenvalues are output in the ``Conquest_out`` file.
The charge densities need post-processing to convert from the
standard output format to a file compatible with visualisation
(current supported formats include Gaussian CUBE file and OpenDX
files).
Note that Becke charges can be calculated if the following parameters
are set:
::
SC.BeckeWeights T
SC.BeckeAtomicRadii T
IO.Iprint_SC 3
This method of output will be refined soon.
Go to :ref:`top <input-output>`
.. _io_output_atoms:
Atomic structure
++++++++++++++++
During structural relaxation and molecular dynamics, the atomic
structure at the end of each step is saved in the output file
``coord_next.dat``. This is in the same format as the input.
Go to :ref:`top <input-output>`
.. _io_md:
Molecular dynamics
++++++++++++++++++
A molecular dynamics run will generate a number of additional plain text output
files:
* ``md.stats`` --- summarises thermodynamic quantities at each steps
* ``md.frames`` --- contains the complete physical state of the system (lattice
parameters, atomic positions, velocities, forces, stress).
* ``md.checkpoint`` --- data required for MD restart, namely atomic velocities
and extended system variables.
* ``md.positions`` --- Atomic coordinates saved at the moment of checkpointing
* ``trajectory.xsf`` --- atomic coordinates save in .xsf format, which can be
visualised using (for example) VMD, if ``AtomMove.WriteXSF`` is true..
Full details are available in :ref:`moldyn`.
Go to :ref:`top <input-output>`

1651
docs/CONQUEST-manual/input_tags.rst Normal file
View File

File diff suppressed because it is too large Load Diff

68
docs/CONQUEST-manual/installing.rst Normal file
View File

@ -0,0 +1,68 @@
.. _install:
============
Installation
============
You will need to download and compile the code before you can use it;
we do not supply binaries.
.. _install_down:
Downloading
-----------
CONQUEST is accessed from `the GitHub repository
<https://github.com/OrderN/CONQUEST-release/>`_;
it can be cloned:
``git clone https://github.com/OrderN/CONQUEST-release destination-directory``
where ``destination-directory`` should be set by the user.
Alternatively, it can be downloaded from GitHub as a zip file and
unpacked:
`<https://github.com/OrderN/CONQUEST-release/archive/master.zip>`_
.. _install_compile:
Compiling
---------
Once you have the distribution, you will need to compile the main
Conquest code (found in the ``src/`` directory), along with the ion file
generation code (found in the ``tools/`` directory). Conquest requires
a working MPI installation including a Fortran90 compiler (often
``mpif90`` but this can vary), along with a few standard libraries:
* BLAS and LAPACK (normally provided by the system vendor)
* FFTW 3.x (more detail can be found at `http://www.fftw.org/ <http://www.fftw.org/>`_)
* ScaLAPACK (often provided as part of an HPC system; the source code
can be obtained from `the netlib repository <http://www.netlib.org/scalapack/>`_ if
you need to compile it)
Additionally, Conquest can use LibXC if it is available (v2.x or
later).
The library locations are set in the ``system.make`` file in the ``src/``
directory, along with other parameters needed for compilation.
* ``FC`` (typically ``FC=mpif90`` will be all that is required)
* ``COMPFLAGS`` (set these to specify compiler options such as
optimisation)
* ``BLAS`` (specify the BLAS and LAPACK libraries)
* ``SCALAPACK`` (specify the ScaLAPACK library)
* ``FFT_LIB`` (must be left as FFTW)
* ``XC_LIBRARY`` (choose ``XC_LIBRARY=CQ`` for the internal Conquest
library, otherwise ``XC_LIBRARY=LibXC_v2`` LibXC v2. or ``XC_LIBRARY=LibXC``
for LibXC v3.x or higher)
* Two further options need to be set for LibXC:
+ ``XC_LIB`` (specify the XC libraries)
+ ``XC_COMPFLAGS`` (specify the location of the LibXC include and
module files, e.g. ``-I/usr/local/include``)
Once these are set, you should make the executable using ``make``.
The ion file generation code is compiled using the same options
required for the main code.

245
docs/CONQUEST-manual/moldyn.rst Normal file
View File

@ -0,0 +1,245 @@
.. _moldyn:
==================
Molecular Dynamics
==================
CONQUEST can perform molecular dynamics both when the density matrix is computed
using diagonalisation and O(N), the latter allowing dynamical simulations of
(but not limited to) tens of thousands of atoms. The equations of motion are
integrated using the velocity Verlet method in the case of the microcanonical
ensemble (NVE), and modifications thereof for the canonical (NVT) and
isobaric-isothermal (NPT) ensembles, the details of which can be found in
:ref:`theory-md`. In addition to converging the parameters for the electronic
structure calculations, the following points must also be considered.
Go to :ref:`top <moldyn>`.
.. _md_scf:
Self-consistency tolerance and XL-BOMD
--------------------------------------
The convergence of the electronic structure is important in MD, as
insufficient convergence can be responsible for "drift" in the
conserved quantity of the dynamics. Although the molecular dynamics
integrators used in CONQUEST are time reversible, *the SCF procedure
is not*. Therefore tight convergence (``minE.SCTolerance`` for
diagonalisation, ``minE.LTolerance`` for linear scaling) is
necessary. In the case of diagonalisation, SCF tolerance of ``1E-6`` is
typically enough to negate the drift. However, extended-Lagrangian
Born-Oppenheimer MD (XL-BOMD) :cite:`md-Niklasson2008`, currently only
implemented for O(N), essentially makes the SCF component of the MD
time-reversible by adding the electronic degrees of freedom to the
Lagrangian, relaxing the constraint on ``minE.LTolerance`` ---
although it is still somewhat dependent on the ensemble. In the NVE
and NVT ensembles, a L-tolerance of ``1E-5`` has been found to be
sufficient to give good energy conservations, decreasing to ``1E-6``
in the NPT ensemble. The following flags are required for XL-BOMD:
::
DM.SolutionMethod ordern
AtomMove.ExtendedLagrangian T
Go to :ref:`top <moldyn>`.
.. _md_restart:
Restarting
----------
Assuming the calculation ended gracefully, it can easily be restarted by
setting,
::
AtomMove.RestartRun T
This will do several things: it will read the atomic coordinates from
``md.position`` and read the ``md.checkpoint`` file, which contains the
velocities and extended system (Nose-Hoover chain and cell) variables. Depending
on the value of ``DM.SolutionMethod``, it will read the K-matrix files
(``diagon``) or the L-matrix files (``ordern``), and if XL-BOMD is being used,
the X-matrix files. Finally, it will *append* new data to the ``md.stats`` and
``md.frames`` files, but it will overwrite all other files, including
``Conquest_out``. Note that this flag is equivalent to setting the following:
::
General.LoadL T
SC.MakeInitialChargeFromK T
XL.LoadL T
In addition to the files mentioned above, CONQUEST will try to read the K-matrix
from ``Kmatrix2.i00.*`` when using diagonalisation or the L-matrix from
``Lmatrix2.i00.*`` when using O(N), and ``Xmatrix2.i0*.*`` if the
extended-Lagrangian formalism is used. Note that metadata for these files is
stored in ``InfoGlobal.i00.dat`` which is also required when restarting. If the
calculation ended by hitting the walltime limit, the writing of these matrix
files may have been interrupted, rendering them unusable. In this case, the
calculation can be restarted by setting the above flags to ``F`` *after* setting
``AtomMove.RestartRun T``. Setting the flag ``General.MaxTime`` to some number
of seconds less (say 30 minutes) than the calculation wall time limit will force
the calculation to stop gracefully, preventing the aforementioned situation.
Go to :ref:`top <moldyn>`.
.. _md_vis:
Visualising the trajectory
--------------------------
Setting the flag ``AtomMove.WriteXSF T`` dumps the coordinates to the file
``trajectory.xsf`` every ``AtomMove.OutputFreq`` steps. The .xsf file can be
read using `VMD <https://www.ks.uiuc.edu/Research/vmd/>`_. A small VMD script,
``view.vmd`` is included with the code, and can be invoked using,
``vmd -e view.vmd``
assuming the vmd executable is in your path.
Go to :ref:`top <moldyn>`.
.. _md_tdep:
TDEP output
-----------
CONQUEST molecular dynamics data can be used to perform lattice dyanmical
calculations using the `Temperature Dependent Effective Potential (TDEP)
<https://ollehellman.github.io/index.html>`_ code. Setting the flag ``MD.TDEP
T`` will make conquest dump configurations, forces and metadata in a format
readable by TDEP.
Go to :ref:`top <moldyn>`.
.. _md_nonh:
Non-Hamiltonian dynamics
------------------------
.. _md_nvt:
Canonical (NVT) ensemble
++++++++++++++++++++++++
The thermostat is set using the ``MD.Thermostat`` flag, and can take the values
``svr`` (stochastic velocity rescaling) and ``nhc`` (Nose-Hoover
chain). These thermostats generate the correct canonical ensemble
phase space distribution, and both give a conserved quantity that
allows the quality of the dynamics to be monitored.
1. Stochastic velocity rescaling
::
AtomMove.IonTemperature 300.0
MD.Ensemble nvt
MD.Thermostat svr
MD.tauT 10
While the NHC uses chaotic sensitivity to initial conditions to achieve better
ergodicity, the SVR thermostat :cite:`md-Bussi2007` uses a judiciously chosen stochastic force
coupled to a weak scaling thermostat to correctly generate the
canonical phase space distribution. The ``MD.tauT`` parameter gives
the coupling timescale; the velocity scaling factor is modified by a
factor :math:`\Delta t/\tau`, so a larger :math:`\tau` results in a
more slowly varying temperature. While some characterisation of the
system is recommended, values of :math:`\tau` around 20--200fs are
reasonable. To reproduce a simulation, the random number
generator seed can be set with the ``General.RNGSeed <integer>`` flag.
2. Nose-Hoover chain
::
AtomMove.IonTemperature 300.0
MD.Ensemble nvt
MD.Thermostat nhc
MD.nNHC 5
MD.nYoshida 5
MD.tauT 30
When thermostatting using a Nose-Hoover chain :cite:`md-Nose1984,md-Hoover1985,md-Martyna1992`, it may be necessary to set a
couple more flags. ``MD.nNHC`` sets the number of thermostats in the chain (the
default of 5 is generally sensible), and ``MD.nYoshida`` determines the order of
Yoshida-Suzuki integration. This is essentially a higher level integration
scheme that *can* improve energy conservation in cases when rapid changes in the
Nose-Hoover thermostat velocity is causing integration errors. **Note that**
``MD.tauT`` **means something different to the SVR case**. A good guess is
the time period of the highest frequency motion of the system in fs; however, in
the NVT ensemble, the energy conservation is not very sensitive to this value.
The NHC masses can also be set manually using the following block.
::
MD.CalculateXLMass F
MD.nNHC 5
%block MD.NHCmass
5 1 1 1 1
%endblock
Go to :ref:`top <moldyn>`.
.. _md_npt:
Isobaric-Isothermal (NPT) ensemble
++++++++++++++++++++++++++++++++++
There is one implemented barostat at present, the extended
system, Parrinello-Rahman :cite:`md-Parrinello1981`. At present the
barostat should be treated as a beta-version implementation, which
will be fully characterised and made robust for the full release of
the code.
1. Parrinello-Rahman
::
AtomMove.IonTemperature 300.0
AtomMove.TargetPressure 10.0
MD.Ensemble npt
MD.Thermostat nhc
MD.Barostat pr
MD.nNHC 5
MD.nYoshida 5
MD.tauT 100
MD.tauP 200
MD.PDrag 10.0
The Parrinello-Rahman barostat generates the correct ensemble, but can
be subject to low frequency "ringing" fluctuations in the
temperature and pressure that can destabilise the system or slow equilibration.
Unlike in the NVT ensemble, this combination of barostat and thermostat is
*very* sensitive to the choice of both ``MD.tauT`` and ``MD.tauP``; note that
their values are somewhat higher in this case, since integration errors in the
NHC tend to be more severe due to coupling of the cell and atomic motions. They
are dependent on the system, so it is advised that you find a combination of
these parameters that gives the best energy conservation. The cell is
thermostatted using a separate Nose-Hoover chain to the atoms by default, but
they can be controlled with the same chain by setting ``MD.CellNHC F``. An *ad
hoc* drag factor specified by ``MD.PDrag`` reduces the thermostat and cell
velocities at every timestep to damp out the ringing fluctuations. In this case,
they are reduced by :math:`10/200 \simeq 5\%`, which strictly speaking breaks the NPT
dynamics, but not significantly, and the stability is significantly improved.
Note that the NPT ensemble can also be generated correctly by thermostatting
using the SVR thermostat, although the meaning of the parameter ``MD.tauT`` is
different in this case, as in NVT dynamics.
Postprocessing tools
--------------------
Details of Python post-processing tools for CONQUEST can be found in :ref:`et_md_scripts`.
Go to :ref:`top <moldyn>`.
.. bibliography:: references.bib
:cited:
:labelprefix: MD
:keyprefix: md-
:style: unsrt
Go to :ref:`top <moldyn>`.

82
docs/CONQUEST-manual/quick-overview.rst Normal file
View File

@ -0,0 +1,82 @@
.. _quick_over:
==============
Quick Overview
==============
.. _setting_up:
Setting up a calculation
------------------------
CONQUEST requires three types of file for a calculation:
* A coordinate file
* Ion files (pseudopotentials)
* The input file (``Conquest_input``)
.. _qo_coords:
Coordinates
===========
CONQUEST works with orthorhombic unit cells (i.e. with angles between
lattice vectors at ninety degrees). The coordinate file is laid out
simply: lattice vectors, number of atoms, atom coordinates (along with
species and movement flags). Either fractional or Cartesian
coordinates can be read (the default is fractional; Cartesian
coordinates require a flag to be set in the
input file). CONQUEST also reads and writes PDB format
coordinate files for biomolecular simulations. More information can
be found in :ref:`io_coords`.
.. _qo_ions:
Ion files
=========
The ion files contain the pseudopotentials and pseudo-atomic orbitals
for the elements, and follow a format similar to the ion files from Siesta
(CONQUEST can read Siesta ion files). A set of default inputs to
generate ion files is available in the directory ``pseudo-and-pao``.
These contain pseudopotentials based on the `PseudoDojo`_ library, and
allow ion files to be produced with the basis set generation code that
is included with CONQUEST in the ``tools/BasisGeneration`` directory.
Full details are found :ref:`here <io_ion>`.
.. _PseudoDojo: https://www.pseudo-dojo.org/
.. _qo_cqinput:
Conquest\_input
===============
The ``Conquest_input`` file contains all of the input flags to
control a CONQUEST run. At a minimum, the file must specify: the run
type (e.g. ``static`` or ``md``); the coordinate file name; and the number
of species and the ion file names. For a well characterised
calculation, further options must be given (for instance setting
details for the calculation of the density matrix). Simple examples
are given in :ref:`examples` and full documentation of all options can
be found in :ref:`input_tags`.
Go to :ref:`top <quick_over>`
.. _qo_output:
Output from a calculation
-------------------------
The main output from CONQUEST is in a single file, named
``Conquest_out`` by default (this can be changed, and output can be
written to ``stdout`` rather than a file). This file
contains details of the calculation, energies, forces and stresses and
the various electronic structure and atomic movement calculations
performed. The most important files that are produced during a run are:
* ``Conquest_out`` The output file
* ``Conquest_warnings`` A list of any warnings issued by the code
(also in ``Conquest_out``)
* ``coord_next.dat`` The updated set of atomic positions
* ``conquest.bib`` References suggested for the calculation performed
* ``input.log`` A log of input options (both set by user and
defaults)
Other files are produced by different run types, and are discussed
elsewhere.

BIN
docs/CONQUEST-manual/rdf.jpg Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 54 KiB

619
docs/CONQUEST-manual/references.bib Normal file
View File

@ -0,0 +1,619 @@
%% This BibTeX bibliography file was created using BibDesk.
%% http://bibdesk.sourceforge.net/
%% Created for David Bowler at 2019-12-20 09:46:12 +0000
%% Saved with string encoding Unicode (UTF-8)
@article{Methfessel:1989ny,
Author = {Methfessel, M. and Paxton, A. T.},
Date-Added = {2019-12-20 09:45:53 +0000},
Date-Modified = {2019-12-20 09:46:11 +0000},
Doi = {10.1103/PhysRevB.40.3616},
Issue = {6},
Journal = {Phys. Rev. B},
Month = {Aug},
Numpages = {0},
Pages = {3616--3621},
Publisher = {American Physical Society},
Title = {High-precision sampling for Brillouin-zone integration in metals},
Url = {https://link.aps.org/doi/10.1103/PhysRevB.40.3616},
Volume = {40},
Year = {1989},
Bdsk-Url-1 = {https://link.aps.org/doi/10.1103/PhysRevB.40.3616},
Bdsk-Url-2 = {https://doi.org/10.1103/PhysRevB.40.3616}}
@article{Monkhorst:1976kf,
Author = {H. J. Monkhorst and J. D. Pack},
Date-Added = {2019-12-19 07:55:12 +0000},
Date-Modified = {2019-12-19 08:11:08 +0000},
Doi = {10.1103/PhysRevB.13.5188},
Journal = {Phys. Rev. B},
Pages = {5188},
Title = {Special points for Brillouin-zone integrations},
Volume = {13},
Year = {1976}}
@article{Wu:2006cu,
Author = {Wu, Z. and Cohen, R. E.},
Date = {2006/06/20/},
Date-Added = {2019-05-23 14:18:32 +0100},
Date-Modified = {2019-05-23 14:20:59 +0100},
Day = {20},
Doi = {10.1103/PhysRevB.73.235116},
Id = {10.1103/PhysRevB.73.235116},
J1 = {PRB},
Journal = {Phys. Rev. B},
Journal1 = {Phys. Rev. B},
Month = {06},
Number = {23},
Pages = {235116},
Publisher = {American Physical Society},
Title = {More accurate generalized gradient approximation for solids},
Ty = {JOUR},
Volume = {73},
Year = {2006},
Bdsk-Url-1 = {https://link.aps.org/doi/10.1103/PhysRevB.73.235116},
Bdsk-Url-2 = {https://doi.org/10.1103/PhysRevB.73.235116}}
@article{Zhang:1998oq,
Author = {Zhang, Y. and Yang, W.},
Date = {1998/01/26/},
Date-Added = {2019-05-23 14:16:34 +0100},
Date-Modified = {2019-05-23 14:16:34 +0100},
Day = {26},
Doi = {10.1103/PhysRevLett.80.890},
Id = {10.1103/PhysRevLett.80.890},
J1 = {PRL},
Journal = {Phys. Rev. Lett.},
Journal1 = {Phys. Rev. Lett.},
Month = {01},
Number = {4},
Pages = {890--890},
Publisher = {American Physical Society},
Title = {Comment on ``Generalized Gradient Approximation Made Simple''},
Ty = {JOUR},
Url = {https://link.aps.org/doi/10.1103/PhysRevLett.80.890},
Volume = {80},
Year = {1998},
Bdsk-Url-1 = {https://link.aps.org/doi/10.1103/PhysRevLett.80.890},
Bdsk-Url-2 = {https://doi.org/10.1103/PhysRevLett.80.890}}
@article{Troullier1991a,
Author = {Troullier, N. and Martins, J. L.},
Doi = {10.1103/PhysRevB.43.1993},
Journal = {Phys. Rev. B},
Pages = {1993},
Title = {Efficient pseudopotentials for plane-wave calculations},
Volume = {43},
Year = {1991},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.43.1993}}
@article{Troullier1991b,
Author = {Troullier, N. and Martins, J. L.},
Doi = {10.1103/PhysRevB.43.8861},
Journal = {Phys. Rev. B},
Pages = {8861},
Title = {Efficient pseudopotentials for plane-wave calculations},
Volume = {43},
Year = {1991},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.43.8861}}
@article{Harris1985,
Author = {Harris, J.},
Doi = {10.1103/PhysRevB.31.1770},
Journal = {Phys. Rev. B},
Pages = {1770},
Title = {Simplified method for calculating the energy of weakly interacting fragments},
Volume = {31},
Year = {1985},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.31.1770}}
@article{Foulkes1989,
Author = {W. Foulkes and R. Haydock},
Doi = {10.1103/PhysRevB.39.12520},
Journal = {Phys. Rev. B},
Pages = {12520},
Title = {Tight-binding models and density-functional theory},
Volume = {39},
Year = {1989},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.39.12520}}
@article{Bowler2002,
Author = {Bowler, D. R. and Miyazaki, T. and Gillan, M. J.},
Doi = {10.1088/0953-8984/14/11/303},
Journal = {J. Phys.: Condens. Matter},
Pages = {2781},
Title = {Recent progress in linear scaling {{\em ab initio}} electronic structure techniques},
Volume = {14},
Year = {2002},
Bdsk-Url-1 = {https://doi.org/10.1088/0953-8984/14/11/303}}
@article{Bowler2006,
Author = {Bowler, D. R. and Choudhury, R. and Gillan, M. J. and Miyazaki, T.},
Doi = {10.1002/pssb.200541386},
Journal = {Phys. Stat. Solidi B},
Pages = {989},
Title = {Recent progress with large-scale ab initio calculations: the CONQUEST code},
Volume = {243},
Year = {2006},
Bdsk-Url-1 = {https://doi.org/10.1002/pssb.200541386}}
@article{Miyazaki2004,
Author = {Miyazaki, T. and Bowler, D. R. and Choudhury, R. and Gillan, M. J.},
Doi = {10.1063/1.1787832},
Journal = {J. Chem. Phys.},
Pages = {6186},
Title = {Atomic force algorithms in density functional theory electronic-structure techniques based on local orbitals},
Volume = {121},
Year = {2004},
Bdsk-Url-1 = {https://doi.org/10.1063/1.1787832}}
@article{Bowler1998a,
Author = {Bowler, D. R. and Gillan, M. J.},
Doi = {10.1016/S0010-4655(98)00061-7},
Journal = {Comp. Phys. Commun.},
Pages = {103},
Title = {Length-scale ill conditioning in linear-scaling DFT},
Volume = {112},
Year = {1998},
Bdsk-Url-1 = {https://doi.org/10.1016/S0010-4655(98)00061-7}}
@article{McWeeny1960,
Author = {McWeeny, R.},
Doi = {doi.org/10.1103/RevModPhys.32.335},
Journal = {Rev. Mod. Phys.},
Pages = {335},
Title = {Some recent advances in density matrix theory},
Volume = {32},
Year = {1960},
Bdsk-Url-1 = {https://doi.org/10.1103/RevModPhys.32.335}}
@article{Goringe1997,
Author = {Goringe, C. M. and Hern{\'{a}}ndez, E. and Gillan, M. J. and Bush, I. J.},
Doi = {doi.org/10.1016/S0010-4655(97)00029-5},
Journal = {Comput. Phys. Commun.},
Pages = {1},
Title = {Linear-scaling DFT-pseudopotential calculations on parallel computers},
Volume = {102},
Year = {1997},
Bdsk-Url-1 = {https://doi.org/10.1016/S0010-4655(97)00029-5}}
@article{Bowler2001,
Author = {Bowler, D. R. and Miyazaki, T. and Gillan, M. J.},
Doi = {10.1016/S0010-4655(01)00164-3},
Journal = {Comput. Phys. Commun.},
Pages = {255-273},
Title = {Parallel sparse matrix multiplication for linear scaling electronic structure calculations},
Volume = {137},
Year = {2001},
Bdsk-Url-1 = {https://doi.org/10.1016/S0010-4655(01)00164-3}}
@article{Bowler1999,
Author = {Bowler, D. R. and Gillan, M. J.},
Doi = {10.1016/S0010-4655(99)00221-0},
Journal = {Comp. Phys. Commun.},
Pages = {95},
Title = {Density matrices in O(N) electronic structure calculations: theory and applications},
Volume = {120},
Year = {1999},
Bdsk-Url-1 = {https://doi.org/10.1016/S0010-4655(99)00221-0}}
@article{Hernandez1997,
Author = {Hern\'andez, E. and Gillan, M. J. and Goringe, C. M.},
Doi = {10.1103/PhysRevB.55.13485},
Journal = {Phys. Rev. B},
Pages = {13485--13493},
Title = {Basis functions for linear-scaling first-principles calculations},
Volume = {55},
Year = {1997},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.55.13485}}
@article{Sankey1989,
Author = {Sankey, O. F. and Niklewski D. J.},
Doi = {10.1103/PhysRevB.40.3979},
Journal = {Phys. Rev. B},
Pages = {3979},
Title = {Ab initio multicenter tight-binding for molecular-dynamics simulations and other applications in covalent systems},
Volume = {40},
Year = {1989},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.40.3979}}
@article{Soler2002,
Author = {Soler, J. M. and Artacho, E. and Gale, J. D. and García, A. and Junquera, J. and Ordejón, P. and Sánchez-Portal, D.},
Doi = {10.1088/0953-8984/14/11/302},
Journal = {J. Phys.:Condens. Matter},
Pages = {2745-2770},
Title = {The SIESTA method for ab initio order-N materials simulation},
Volume = {14},
Year = {2002},
Bdsk-Url-1 = {https://doi.org/10.1088/0953-8984/14/11/302}}
@article{Kenny2000,
Author = {Kenny, S. D. and Horsfield, A. P. and Fujitani, H.},
Doi = {10.1103/PhysRevB.62.4899},
Journal = {Phys. Rev. B},
Pages = {4899},
Title = {Transferable atomic-type orbital basis sets for solids},
Volume = {62},
Year = {2000},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.62.4899}}
@article{Junquera2001,
Author = {Junquera, J. and Paz, O. and Sanchez-Portal, D. and Artacho, E.},
Doi = {10.1103/PhysRevB.64.235111},
Journal = {Phys. Rev. B},
Pages = {235111},
Title = {Numerical atomic orbitals for linear-scaling calculations},
Volume = {64},
Year = {2001},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.64.235111}}
@article{Ozaki2003,
Author = {Ozaki, T.},
Doi = {10.1103/PhysRevB.67.155108},
Journal = {Phys. Rev. B.},
Pages = {155108},
Title = {Variationally optimized atomic orbitals for large-scale electronic structures},
Volume = {67},
Year = {2003},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.67.155108}}
@article{Ozaki2004,
Author = {Ozaki, T. and Kino, H.},
Doi = {10.1103/PhysRevB.69.195113},
Journal = {Phys. Rev. B.},
Pages = {195113},
Title = {Numerical atomic basis orbitals from H to Kr},
Volume = {69},
Year = {2004},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.69.195113}}
@article{Goedecker1999,
Author = {Goedecker, S.},
Doi = {10.1103/RevModPhys.71.1085},
Journal = {Rev. Mod. Phys.},
Number = {4},
Pages = {1085},
Title = {Linear scaling electronic structure methods},
Volume = {71},
Year = {1999},
Bdsk-Url-1 = {https://doi.org/10.1103/RevModPhys.71.1085}}
@article{Goedecker1994,
Author = {Goedecker, S. and Colombo, L.},
Doi = {10.1103/PhysRevLett.73.122},
Journal = {Phys. Rev. Lett.},
Pages = {122},
Title = {Efficient linear scaling algorithm for tight-binding molecular dyanmics},
Volume = {73},
Year = {1994},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevLett.73.122}}
@article{Li1993,
Author = {Li, X.-P. and Nunes, R. W. and Vanderbilt, D.},
Doi = {10.1103/physrevb.47.10891},
Journal = {Phys. Rev. B},
Pages = {10891},
Title = {Density-matrix electronic-structure method with linear system-size scaling},
Volume = {47},
Year = {1993},
Bdsk-Url-1 = {https://doi.org/10.1103/physrevb.47.10891}}
@article{Palser1998,
Author = {Palser, A. H. R. and Manolopoulos, D. E.},
Doi = {10.1103/PhysRevB.58.12704},
Journal = {Phys. Rev. B},
Number = {19},
Pages = {12704},
Title = {Canonical purification of the density matrix in electronic-structure theory},
Volume = {58},
Year = {1998},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.58.12704}}
@article{Perdew1981,
Author = {Perdew, J. P. and Zunger, A.},
Doi = {10.1103/PhysRevB.23.5048},
Journal = {Phys. Rev. B},
Pages = {5048},
Title = {Self-interaction correction to density-functional approximations for many-electron systems},
Volume = {23},
Year = {1981},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.23.5048}}
@article{Goedecker1996,
Author = {Goedecker, S. and Teter, M. and Hutter, J.},
Doi = {10.1103/PhysRevB.54.1703},
Journal = {Phys. Rev. B},
Pages = {1703},
Title = {Separable dual-space Gaussian pseudopotentials},
Volume = {54},
Year = {1996},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.54.1703}}
@article{Perdew1996,
Author = {Perdew, J. P. and Burke, K. and Ernzerhof, M.},
Doi = {10.1103/PhysRevLett.77.3865},
Journal = {Phys. Rev. Lett.},
Number = {18},
Pages = {3865},
Title = {Generalized Gradient Approximation Made Simple},
Volume = {77},
Year = {1996},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevLett.77.3865}}
@article{Perdew1998,
Author = {Perdew, J. P. and Burke, K. and Wang, Y.},
Doi = {10.1103/PhysRevB.57.14999},
Journal = {Phys. Rev. B},
Pages = {14999},
Title = {Erratum: Generalized gradient approximation for the exchange-correlation hole of a many-electron system [Phys. Rev. B 54, 16 533 (1996)]},
Volume = {57},
Year = {1998},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.57.14999}}
@article{Hammer1999,
Author = {Hammer, B. and Hansen, L. B. and Nørskov, J. K.},
Doi = {10.1103/PhysRevB.59.7413},
Journal = {Phys. Rev. B},
Pages = {7413},
Title = {Improved adsorption energetics within density-functional theory using revised Perdew-Burke-Ernzerhof functionals},
Volume = {59},
Year = {1999},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.59.7413}}
@article{Perdew1992,
Author = {Perdew, John P. and Wang, Yue},
Doi = {10.1103/PhysRevB.45.13244},
Journal = {Phys. Rev. B},
Pages = {13244},
Title = {Accurate and simple analytic representation of the electron-gas correlation energy},
Volume = {45},
Year = {1992},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.45.13244}}
@article{Hirakawa2017,
Author = {Hirakawa, T. and Suzuki, T. and Bowler, D. R. and Miyazaki, T.},
Doi = {10.1088/1361-648X/aa810d},
Journal = {J. Phys.: Condens. Matter},
Pages = {405901},
Title = {Canonical-ensemble extended Lagrangian Born-Oppenheimer molecular dynamics for the linear scaling density functional theory},
Volume = {29},
Year = {2017},
Bdsk-Url-1 = {https://doi.org/10.1088/1361-648X/aa810d}}
@book{Tuckerman2010,
Author = {Tuckerman, M. E.},
Publisher = {Oxford Graduate Texts},
Title = {Statistical mechanics: theory and molecular simulations},
Year = {2010}}
@book{Frenkel2002,
Author = {Frenkel, D. and Smit, B.},
Publisher = {Academic Press},
Title = {Understanding molecular simulation: from algorithms to application},
Year = {2002}}
@article{Nose1984,
Author = {Nosé, S.},
Doi = {10.1063/1.447334},
Journal = {J. Chem. Phys.},
Pages = {511},
Title = {A unified formulation of the constant temperature molecular dynamics methods},
Volume = {81},
Year = {1984},
Bdsk-Url-1 = {https://doi.org/10.1063/1.447334}}
@article{Martyna1992,
Author = {Martyna, G. J. and Klein, M. L. and Tuckerman, M.},
Doi = {10.1063/1.463940},
Journal = {J. Chem. Phys.},
Pages = {2635},
Title = {Nosé{\textendash}Hoover chains: The canonical ensemble via continuous dynamics},
Volume = {97},
Year = {1992},
Bdsk-Url-1 = {https://doi.org/10.1063/1.463940}}
@article{Martyna1996,
Author = {Martyna, G. J. and Tuckerman, M. E. and Tobias, D. J. and Klein, M. L.},
Doi = {10.1080/002689799163235},
Journal = {Mol. Phys.},
Pages = {1117},
Title = {{Explicit reversible integrators for extended systems dynamics}},
Volume = {87},
Year = {1996},
Bdsk-Url-1 = {https://doi.org/10.1080/002689799163235}}
@article{Berendsen1984,
Author = {Berendsen, H. J. C. and Postma, J. P. M. and van Gunsteren, W. F. and Dinola, A. and Haak, J. R.},
Doi = {10.1063/1.448118},
Journal = {J. Chem. Phys.},
Pages = {3684},
Title = {{Molecular dynamics with coupling to an external bath}},
Volume = {81},
Year = {1984},
Bdsk-Url-1 = {https://doi.org/10.1063/1.448118}}
@article{Hoover1985,
Author = {Hoover, W. G.},
Doi = {10.1103/PhysRevA.31.1695},
Journal = {Phys. Rev. A},
Pages = {1695},
Title = {{Canonical dynamics: Equilibrium phase-space distributions}},
Volume = {31},
Year = {1985},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevA.31.1695}}
@article{Niklasson2014,
Author = {Niklasson, A. M. N. and Cawkwell, M. J.},
Doi = {10.1063/1.4898803},
Journal = {J. Chem. Phys.},
Pages = {164123},
Title = {{Generalized extended Lagrangian Born-Oppenheimer molecular dynamics}},
Volume = {141},
Year = {2014},
Bdsk-Url-1 = {https://doi.org/10.1063/1.4898803}}
@article{Niklasson2008,
Author = {Niklasson, A. M. N.},
Doi = {10.1103/PhysRevLett.100.123004},
Journal = {Phys. Rev. Lett.},
Pages = {123004},
Title = {{Extended Born-Oppenheimer Molecular Dynamics}},
Volume = {100},
Year = {2008},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevLett.100.123004}}
@article{Arita2014,
Author = {Arita, M. and Bowler, D. R. and Miyazaki, T.},
Doi = {10.1021/ct500847y},
Journal = {J. Chem. Theor. Comput.},
Pages = {5419},
Title = {{Stable and Efficient Linear Scaling First-Principles Molecular Dynamics for 10000+ Atoms}},
Volume = {10},
Year = {2014},
Bdsk-Url-1 = {https://doi.org/10.1021/ct500847y}}
@article{Andersen1980,
Author = {Andersen, H. C.},
Doi = {10.1063/1.439486},
Journal = {J. Chem. Phys.},
Pages = {2384},
Title = {{Molecular dynamics simulations at constant pressure and/or temperature}},
Volume = {72},
Year = {1980},
Bdsk-Url-1 = {https://doi.org/10.1063/1.439486}}
@article{Ray1983,
Author = {Ray, J. R.},
Doi = {10.1063/1.445636},
Journal = {J. Chem. Phys.},
Pages = {5128},
Title = {{Molecular dynamics equations of motion for systems varying in shape and size}},
Volume = {79},
Year = {1983},
Bdsk-Url-1 = {https://doi.org/10.1063/1.445636}}
@article{Wentzcovitch1991,
Author = {Wentzcovitch, R. M.},
Doi = {10.1103/physrevb.44.2358},
Journal = {Phys. Rev. B},
Pages = {2358},
Title = {{Invariant molecular-dynamics approach to structural phase transitions}},
Volume = {44},
Year = {1991},
Bdsk-Url-1 = {https://doi.org/10.1103/physrevb.44.2358}}
@article{Artacho1999,
Author = {Artacho, E. and Sanchez-Portal, D. and Ordejon, P. and Garcia, A. and Soler, J. M.},
Doi = {10.1002/(SICI)1521-3951(199909)215:1<809::AID-PSSB809>3.0.CO;2-0},
Journal = {Phys. Stat. Solidi B},
Pages = {809},
Title = {Linear-scaling ab-initio calculations for large and complex systems},
Volume = {215},
Year = {1999},
Bdsk-Url-1 = {https://doi.org/10.1002/(SICI)1521-3951(199909)215:1%3C809::AID-PSSB809%3E3.0.CO;2-0}}
@article{Shinoda2004,
Author = {Shinoda, W. and Shiga, M. and Mikami, M.},
Doi = {10.1103/PhysRevB.69.134103},
Journal = {Phys. Rev. B},
Pages = {4396},
Title = {{Rapid estimation of elastic constants by molecular dynamics simulation under constant stress}},
Volume = {69},
Year = {2004},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevB.69.134103}}
@article{Parrinello1981,
Author = {Parrinello, M. and Rahman, A.},
Doi = {10.1063/1.328693},
Journal = {J. Appl. Phys.},
Month = dec,
Pages = {7182--7190},
Title = {{Polymorphic transitions in single crystals: A new molecular dynamics method}},
Volume = {52},
Year = {1981},
Bdsk-Url-1 = {https://doi.org/10.1063/1.328693}}
@article{Bussi2007,
Author = {Bussi, G. and Donadio, D. and Parrinello, M.},
Doi = {10.1063/1.2408420},
Journal = {J. Chem. Phys.},
Pages = {014101},
Title = {{Canonical sampling through velocity rescaling}},
Volume = {126},
Year = {2007},
Bdsk-Url-1 = {https://doi.org/10.1063/1.2408420}}
@article{Bussi2009,
Author = {Bussi, G. and Zykova-Timan, T. and Parrinello, M.},
Doi = {10.1063/1.3073889},
Journal = {J. Chem. Phys.},
Pages = {074101},
Title = {{Isothermal-isobaric molecular dynamics using stochastic velocity rescaling}},
Volume = {130},
Year = {2009},
Bdsk-Url-1 = {https://doi.org/10.1063/1.3073889}}
@article{Bitzek2006,
Author = {Bitzek, E. and Koskinen, P. and Gähler, F. and Moseler, M. and Gumbsch, P.},
Doi = {10.1103/PhysRevLett.97.170201},
Journal = {Phys. Rev. Lett.},
Pages = {2897},
Title = {{Structural Relaxation Made Simple}},
Volume = {97},
Year = {2006},
Bdsk-Url-1 = {https://doi.org/10.1103/PhysRevLett.97.170201}}
@article{Pfrommer1997,
author = {Pfrommer, B. G. and C\^ot\'e, M. and Louie, S. and Cohen, M. L.},
title = {{Relaxation of Crystals with the Quasi-Newton Method}},
journal = {J. Comput. Phys.},
year = {1997},
volume = {131},
pages = {233},
doi = {10.1006/jcph.1996.5612}}
@article{Bowler:2019fv,
Author = {David R. Bowler and Jack S. Baker and Jack T. L. Poulton and Shereif Y. Mujahed and Jianbo Lin and Sushma Yadav and Zamaan Raza and Tsuyoshi Miyazaki},
Doi = {10.7567/1347-4065/ab45af},
Journal = {Jap. J. Appl. Phys.},
Pages = {100503},
Title = {Highly accurate local basis sets for large-scale {DFT} calculations in CONQUEST},
Volume = {58},
Year = 2019}
@article{Nakata2015,
Author = {Nakata, Ayako and Bowler, David and Miyazaki, Tsuyoshi},
Journal = {Phys. Chem. Chem. Phys.},
doi = {10.1021/ct5004934},
Pages = {31427},
Title = {{Optimized multi-site local orbitals in the large-scale DFT program CONQUEST}},
Volume = {17},
Year = {2015}}
@article{Nakata2014,
Author = {Nakata, Ayako and Bowler, David R. and Miyazaki, Tsuyoshi},
Journal = {J. Chem. Theory Comput.},
doi = {10.1039/C5CP00934K},
Pages = {4813},
Title = {{Efficient Calculations with Multisite Local Orbitals in a Large-Scale DFT Code CONQUEST}},
Volume = {10},
Year = {2014}}
@article{Haynes:2006qe,
Author = {P.D. Haynes and C.-K. Skylaris and A.A. Mostofi and M.C. Payne},
Doi = {10.1016/j.cplett.2006.02.086},
Journal = {Chem. Phys. Lett.},
Pages = {345 - 349},
Title = {Elimination of basis set superposition error in linear-scaling density-functional calculations with local orbitals optimised in situ},
Volume = {422},
Year = {2006}}
@article{Boys:1970aa,
author = {S.F. Boys and F. Bernardi},
title = {The calculation of small molecular interactions by the differences of separate total energies. Some procedures with reduced errors},
journal = {Mol. Phys.},
volume = {19},
year = {1970},
doi = {10.1080/00268977000101561},
pages = {553}}

14
docs/CONQUEST-manual/starting.rst
View File

@ -1,14 +0,0 @@
.. _starting:
===============
Getting started
===============
Installation
------------
Tutorials
---------
Where next?
-----------

BIN
docs/CONQUEST-manual/stats.jpg Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 160 KiB

171
docs/CONQUEST-manual/strucrelax.rst Normal file
View File

@ -0,0 +1,171 @@
.. _strucrelax:
=====================
Structural relaxation
=====================
This section describes how to find the zero-Kelvin equilibrium structure, given
a starting structure with non-zero forces and/or stresses. CONQUEST
can employ a variety of algorithms
algorithm to minimise energy with respect to
atomic positions, including: L-BFGS; conjugate gradients; and damped
molecular dynamics (both MDMin and FIRE approaches). The minimisation
of energy or enthalpy with respect to cell vectors is restricted to
conjugate gradients at present, though L-BFGS will be implemented.
Setting ``AtomMove.WriteXSF T`` for all flavours of optimisation will dump the
trajectory to the file ``trajectory.xsf``, which can be visualised using `VMD
<https://www.ks.uiuc.edu/Research/vmd/>`_. Setting ``AtomMove.AppendCoords T``
will append the structure at each step to ``UpdatedAtoms.dat`` in the format of a
CONQUEST structure input.
For the L-BFGS and conjugate gradients relaxations, the progress of the calculation can be
monitored by searching for the word ``GeomOpt``; grepping will print the
following:
::
$ grep GeomOpt Conquest_out
GeomOpt - Iter: 0 MaxF: 0.00329282 H: -0.14168571E+03 dH: 0.00000000
GeomOpt - Iter: 1 MaxF: 0.00331536 H: -0.14168995E+03 dH: 0.00424155
GeomOpt - Iter: 2 MaxF: 0.00350781 H: -0.14168997E+03 dH: 0.00001651
GeomOpt - Iter: 3 MaxF: 0.00504075 H: -0.14169161E+03 dH: 0.00164389
GeomOpt - Iter: 4 MaxF: 0.00725611 H: -0.14169172E+03 dH: 0.00010500
GeomOpt - Iter: 5 MaxF: 0.01134145 H: -0.14169329E+03 dH: 0.00157361
GeomOpt - Iter: 6 MaxF: 0.01417229 H: -0.14169385E+03 dH: 0.00056077
GeomOpt - Iter: 7 MaxF: 0.01434628 H: -0.14169575E+03 dH: 0.00190304
GeomOpt - Iter: 8 MaxF: 0.01711197 H: -0.14170001E+03 dH: 0.00425400
GeomOpt - Iter: 9 MaxF: 0.02040556 H: -0.14170382E+03 dH: 0.00381110
GeomOpt - Iter: 10 MaxF: 0.01095167 H: -0.14170752E+03 dH: 0.00370442
In this example, MaxF is the maximum single force component, H is the enthalpy and dH is the
change in enthalpy.
Go to :ref:`top <strucrelax>`.
.. _sr_ions:
Ionic relaxation
----------------
To optimise the ionic positions with respect to the DFT total energy, the
following flags are essential:
::
AtomMove.TypeOfRun lbfgs
AtomMove.MaxForceTol 5e-4
AtomMove.ReuseDM T
The parameter ``AtomMove.TypeOfRun`` can take the value ``lbfgs`` or
``cg`` for iterative optimisation. Both algorithms are robust and
relatively efficient in most instances; L-BFGS is preferred. The
parameter ``AtomMove.MaxForceTol`` specifies the force
convergence criterion in Ha/bohr, i.e. the calculation will terminate
when the largest force component on any atom is below this value.
The parameter
``AtomMove.ReuseDM`` specifies that the density matrix (the K-matrix for
diagonalisation or L-matrix for O(N) calculations) from the
previous step will be used as an initial guess for the SCF cycle after
propagating the atoms; this should generally decrease the number of SCF cycles
per ionic step.
If the self-consistency tolerance is too low, the optimisation may fail to
converge with respect to the force tolerance; this may necessitate a tighter
``minE.SCTolerance`` for diagonalisation (also possibly
``minE.LTolerance`` for O(N) calculations). A grid which is too
coarse can also cause problems with structural relaxation to high tolerances.
For problematic cases where the conjugate gradients algorithm fails to find a
downhill search direction, it may be worth trying quenched molecular dyanamics,
which propagates the equations of motion following a simple NVE
approach, but resets the velocities to zero when the dot product of
force and velocity is zero.
::
AtomMove.TypeOfRun md
AtomMove.QuenchedMD T
AtomMove.MaxForceTol 5e-4
AtomMove.ReuseDM T
The FIRE algorithm :cite:`sr-Bitzek2006` is a variant of quenched MD
that has been shown to outperform conjugate gradients in some
circumstances.
::
AtomMove.TypeOfRun md
AtomMove.FIRE T
AtomMove.MaxForceTol 5e-4
AtomMove.ReuseDM T
Go to :ref:`top <strucrelax>`.
.. _sr_cell:
Cell optimisation
-----------------
The unit cell can be optimised with respect to enthalpy *with fixed fractional
coordinates* (``AtomMove.OptCellMethod 1``) using the following input:
::
AtomMove.TypeOfRun cg
AtomMove.OptCell T
AtomMove.OptCellMethod 1
AtomMove.TargetPressure 1.0
AtomMove.ReuseL T
AtomMove.EnthalpyTolerance 1E-6
AtomMove.StressTolerance 0.01
Here, we specify the target pressure in GPa and two new tolerances, the enthalpy
tolerance in Ha and the stress tolerance in GPa.
Go to :ref:`top <strucrelax>`.
.. _sr_both:
Combined optimisation
---------------------
For simple crystals, the fractional ionic coordinates vary trivially with
changes in the lattice vectors; however for more complicated systems such as
molecular crystals and amorphous materials, it is necessary simultaneously relax
the ionic positions and lattice vectors. This can be done by setting
``AtomMove.OptCellMethod 3``
::
AtomMove.TypeOfRun cg
AtomMove.OptCell T
AtomMove.OptCellMethod 3
AtomMove.TargetPressure 1.0
AtomMove.ReuseL T
AtomMove.MaxForceTol 5e-4
AtomMove.EnthalpyTolerance 1E-6
AtomMove.StressTolerance 0.01
Note that the enthalpy will generally converge much more rapidly than the force
and stress, and that it may be necessary to tighten ``minE.SCTolerance``
(diagonalisation) or ``minE.LTolerance`` (order(N)) to reach the force
tolerance, if it is even possible.
Due to the nature of the complex partitioning system, large and sudden changes in volume
may cause the calculation to crash, particlularly in the case of combined
optimisation. In such cases, it may help to try ``AtomMove.OptCellMethod 2``,
which uses a simple but robust double-loop minimisation: a full ionic conjugate
gradients relaxation for the inner loop and a single cell steepest descent
relaxation for the outer loop. This is considerably less efficient, but
may help in particularly problematic cases.
Go to :ref:`top <strucrelax>`.
.. bibliography:: references.bib
:cited:
:labelprefix: SR
:keyprefix: sr-
:style: unsrt
Go to :ref:`top <strucrelax>`.

364
docs/CONQUEST-manual/theory-md.rst Normal file
View File

@ -0,0 +1,364 @@
.. _theory-md:
==========================
Molecular Dynamics: Theory
==========================
.. _th_md_nve:
Microcanonical (NVE) ensemble
-----------------------------
The Hamiltonian for the microcanonical ensemble is,
.. math::
\mathcal{H} = \sum_{i=1}^N \frac{\mathbf{p}_i^2}{2m_i} + U(\mathbf{r}_i)
where :math:`\mathbf{p}_i` and :math:`\mathbf{r}_i` are the position and
momentum of particle :math:`i` and :math:`U` is the DFT total (potential)
energy. Hamilton's equations can be solved to give the following equations of
motion:
.. math::
\mathbf{\dot{r}}_i &= \frac{\mathbf{p}_i}{m_i} \\
\mathbf{\dot{p}}_i &= \frac{\partial U(\mathbf{r}_i)}{\partial\mathbf{r}_i} = \mathbf{F_i}
In order to construct a time-reversible algorithm from these equations, the
Liouvillian formulation is employed :cite:`t-Frenkel2002` (trivially, in this
case). The Liouville operator :math:`L` can be defined in terms of position and
momentum components:
.. math::
iL = \mathbf{\dot{r}}\frac{\partial}{\partial\mathbf{r}} + \mathbf{\dot{p}}\frac{\partial}{\partial\mathbf{p}} = i(L_r + L_p).
The Liouvillian can be used to construct the classical propagator, which relates
the state :math:`f` of the system at time 0 to its state at time :math:`t`:
.. math::
f[\mathbf{p}^N(t),\mathbf{r}^N(t)] = e^{iLt}f[\mathbf{p}^N(0),\mathbf{r}^N(0)]
Taking the individual position and momentum parts of the Liouvillian :math:`L_r`
and :math:`L_p`, it can be shown that applying it to the state :math:`f` result
in a simple linear shift in coordinates and a simple linear shift in momentum
respectively:
.. math::
iL_rf(t) &= f[\mathbf{p}^N(0),\mathbf{r}^N(0) + \mathbf{\dot{r}}^N(0)t] \\
iL_pf(t) &= f[\mathbf{p}^N(0) + \mathbf{F}^N(0)t,\mathbf{r}^N(0)]
However, we cannot simply replace :math:`e^{iLt}` with :math:`e^{iL_rt}` because
:math:`iL_r` and :math:`iL_p` are non-commuting operators, so we must employ the
Trotter-Suzuki identity:
.. math::
e^{A+B} = \lim_{P\rightarrow\infty}\left(e^{A/2P}e^{B/P}e^{A/2P}\right)^P
Thus for a small enough time step :math:`\Delta t = t/P` and to first order, a
discrete time step corresponds to the application of the discrete time
propagator :math:`G`,
.. math::
G(\Delta t) = U_1\left(\frac{\Delta t}{2}\right)U_2\left(\Delta t\right)U_1\left(\frac{\Delta t}{2}\right) = e^{iL_1\frac{\Delta t}{2}}e^{iL_2\Delta t}e^{iL_1\frac{\Delta t}{2}},
which can be shown to be unitary and therefore time-reversible. Applying the
operators :math:`U` in the sequence determined by the Trotter decomposition
generates the velocity Verlet algorithm, which is used to integrate
microcanonical molecular dynamics in CONQUEST. For a detailed derivation of the
algorithm, refer to Frenkel & Smit :cite:`t-Frenkel2002`.
Go to :ref:`top <theory-md>`.
.. _th_md_xlbomd:
Extended Lagrangian Born-Oppenheimer MD (XL-BOMD)
-------------------------------------------------
If the electronic density from the previous ionic step is used as an initila
guess for the next SCF cycle, a problem arises because this process breaks the
time-reversibility of the dynamics. This is manifested as a gradual drift in the
total energy in the case of a NVE simulation, or the conserved quantity in the
case of non-Hamiltonian dynamics. The solution proposed by Niklasson
:cite:`t-Niklasson2008,t-Niklasson2014` is to introduce auxilliary electronic
degrees of freedom into the Lagrangian, which can be propagated via
time-reversible integrators.
The extended Lagrangian used in CONQUEST is :cite:`t-Arita2014`,
.. math::
\mathcal{L}^\mathrm{XBO}\left(\mathbf{X}, \mathbf{\dot{X}}, \mathbf{R}, \mathbf{\dot{R}}\right) = \mathcal{L}^\mathrm{BO}\left(\mathbf{R}, \mathbf{\dot{R}}\right) + \frac{1}{2}\mu\mathrm{Tr}\left[\mathbf{\dot{X}}^2\right] - \frac{1}{2}\mu\omega^2\mathrm{Tr}\left[(\mathbf{LS} - \mathbf{X})^2\right],
where :math:`\mathbf{X}` is a sparse matrix with the same range as
:math:`\mathbf{LS}`, :math:`\mu` is the fictitious electron mass and
:math:`\omega` is the curvature of the auxiliary harmonic potential. The
Euler-Lagrange equations of motion are then,
.. math::
m_i\mathbf{\ddot{r}_i} &= -\frac{\partial U[{{\mathbf{R;LS}}}]}{\partial\mathbf{r}_i} = \mathbf{F_i} \\
\mathbf{\ddot{X}} &= \omega^2(\mathbf{LS} - \mathbf{X}),
The first of these is simply Newton's second law, and the velocity update
equation of motion in the microcanonical ensemble. The second can be integrated
using a time-reversible algorithm, the velocity Verlet scheme in the case of
CONQUEST :cite:`t-Arita2014`:
.. math::
\mathbf{X}(t+\delta t) &= 2\mathbf{X}(t) -\mathbf{X}(t-\delta t) + \delta t^2\omega^2\left[\mathbf{L}(t)\mathbf{S}(t)-\mathbf{X}(t)\right] \\
&+ a\sum_{m=0}^M c_m\mathbf{X}(t-m\delta t)
i.e. the trajectory of :math:`\mathbf{X}(t)` is time-reversible, and evolves in
a harmonic potential centred on the ground state density
:math:`\mathbf{L}(t)\mathbf{S}(t)`. The matrix :math:`\mathbf{XS}^{-1}` is a
good guess for the :math:`\mathbf{L}` matrix in the Order(N) scheme.
Despite the time-reversitibility, the :math:`\mathbf{X}` matrix tends in
practice to gradually drift from the harmonic centre over time, increasing the
number of SCF iterations required to reach the minimum over the course of the
simulation. To remove such numerical errors, the final dissipative term is
included, and is found to have a minimal effect on the time-reversibility. We
note that since the auxiliary variable :math:`X` is used to generate an intial
guess for the SCF process, it does not appear in the conserved
(pseudo-Hamiltonian) quantity for the dynamics.
Go to :ref:`top <theory-md>`.
.. _th_md_nonH:
Non-Hamiltonian dynamics
------------------------
.. _th_md_ext:
Extended system method
~~~~~~~~~~~~~~~~~~~~~~
Hamiltonian dynamics generally describes systems that are isolated from their
surroundings, but in the canonical and isobaric-isothermal ensembles, we need to
couple the system to an external heat bath and/or stress. It is possible to
model such systems by positing a set of equations of *non-Hamiltonian* equations
of motion, and proving that they generate the correct statistical ensemble
:cite:`t-Tuckerman2010`. This is the extended system approach: we modify the
Hamiltonian to include the thermostat and/or barostat degrees of freedom, derive
the (pseudo-) Hamiltonian equations of motion, and demostrate that the correct
phase space distribution for the ensemble is recovered.
Go to :ref:`top <theory-md>`.
.. _th_md_nvt:
Canonical (NVT) ensemble
~~~~~~~~~~~~~~~~~~~~~~~~
In the Nose-Hoover formulation :cite:`t-Nose1984,t-Hoover1985`, the Hamiltonian
for a system in the canonical ensemble can be written,
.. math::
\mathcal{H} = \sum_i \frac{1}{2}m_i s^2\mathbf{\dot{r}}_i^2 + U(\mathbf{r}_i) + \frac{1}{2}Q\dot{s}^2 - (n_f + 1)k_B T \ln s,
where :math:`\mathbf{r}_i` and :math:`\mathbf{\dot{r}_i}` are respectively the
position and velocity of particle :math:`i`, :math:`U` is the potential energy,
in this case the DFT total energy, :math:`s` is a dimensionless quantity that
can be interpreted post-hoc as a time step scaling factor, :math:`Q` is the
fictitious mass of the heat bath and :math:`n_f` is the number of ionic degrees
of freedom. Hamilton's equations can be solved to generate the Nose-Hoover
equations of motion. However Martyna *et al*. demonstrate that this method does
not generate an ergodic trajectory, and proposed an alternative formulation
:cite:`t-Martyna1992` in which the temperature is controlled by a chain of
:math:`M` coupled thermostats of mass :math:`Q_k`, notional position
:math:`\eta_k` and conjugate momentum :math:`p_{\eta_k}`:
.. math::
\mathbf{\dot{r}_i} &= \frac{\mathbf{p}_i}{m_i} \\
\mathbf{\dot{p}_i} &= -\frac{\partial U(\mathbf{r})}{\partial \mathbf{r}_i} - \frac{p_{\eta_1}}{Q_1}\mathbf{p}_i \\
\dot{\eta}_k &= \frac{p_{\eta_k}}{Q_k} \\
\dot{p}_{\eta_1} &= \left(\sum_{i=1}^N\frac{\mathbf{p}_i}{m_i} - n_fk_BT\right) - \frac{p_{\eta_{2}}}{Q_{\eta_{2}}}p_{\eta_1} \\
\dot{p}_{\eta_k} &= \left(\frac{p^2_{\eta_{k-1}}}{Q_{k-1}} - k_BT\right) - \frac{p_{\eta_{k+1}}}{Q_{k+1}}p_{\eta_k} \\
\dot{p}_{\eta_M} &= \left(\frac{p^2_{\eta_{M-1}}}{Q_{M-1}} - k_BT\right)
The Liouvillian for these equations of motion can be non-uniquely decomposed
into components of ionic position (:math:`iL_r`) and momentum (:math:`iL_p`) as
in the microcanonical case, the extended Lagrangian (:math:`iL_\mathrm{XL}`, and
a Nose-Hoover chain component (:math:`iL_\mathrm{NHC}`)
.. math::
iL = iL_\mathrm{NHC} + iL_p + iL_{\mathrm{XL}} + iL_r,
which is directly translated into an algorithm with the Trotter-Suzuki
expansion,
.. math::
\exp(iL\Delta t) = &\exp\left(iL_\mathrm{NHC}\frac{\Delta t}{2}\right)\exp\left(iL_p\frac{\Delta t}{2}\right) \times \\
&\exp\left(iL_\mathrm{XL}\frac{\Delta t}{2}\right)\exp\left(iL_r\Delta t\right)\exp\left(iL_\mathrm{XL}\frac{\Delta t}{2}\right) \times \\
&\exp\left(iL_p\frac{\Delta t}{2}\right)\exp\left(iL_\mathrm{NHC}\frac{\Delta t}{2}\right)
This is recognisable as the velocity Verlet algorithm with extended Lagrangian
integration which can be reduced to a single step, as described in
:ref:`Extended Lagrangian Born-Oppenheimer MD (XL-BOMD)`, with a half time step
integration of the Nose-Hoover chain equations of motion before and after. For
full details of the integration scheme, see Hirakawa *et al*.
:cite:`t-Hirakawa2017`.
Go to :ref:`top <theory-md>`.
.. _th_md_npt:
Isobaric-Isothermal (NPT) ensemble
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The Parinello-Rahman equations of motion :cite:`t-Parrinello1981` extend the
fixed cell equations of motion to include the cell degrees of freedom in the
extended system approach. We use the Martyna-Tobias-Tuckerman-Klein modification
:cite:`t-Martyna1996`, which couples the variable cell equations of motion to a
Nose-Hoover chain the themrostat the system, recovering the isobaric-isothermal
(NPT) ensemble. For an unconstrained cell (i.e. the lattice vectors can change
freely), the equations of motion are,
.. math::
\mathbf{\dot{r}}_i &= \frac{\mathbf{p}_i}{m_i} + \frac{\mathbf{p}_g}{W_g}\mathbf{r}_i \\
\mathbf{\dot{p}}_i &= \mathbf{F}_i - \frac{\mathbf{p}_g}{W_g}\mathbf{p}_i - \left(\frac{1}{N_f}\right)\frac{\mathrm{Tr}[\mathbf{p}_g]}{W_g}\mathbf{p}_i - \frac{p_\xi}{Q}\mathbf{p}_i \\
\mathbf{\dot{h}} &= \frac{\mathbf{p}_g\mathbf{h}}{W_g} \\
\mathbf{\dot{p}_g} &= V(\mathbf{P}_\mathrm{int}-\mathbf{I}P_\mathrm{ext}) + \left[\frac{1}{N_f}\sum_{i=1}^N\frac{\mathbf{p}_i^2}{m_i}\right]\mathbf{I} - \frac{p_\xi}{Q}\mathbf{p}_g \\
\dot{\xi} &= \frac{p_\xi}{Q} \\
\mathbf{\dot{p}}_g &= \sum_{i=1}^N\frac{\mathbf{p}_i^2}{m_i} + \frac{1}{W_g}\mathrm{Tr}[\mathbf{p}_g^T\mathbf{p}_g] - (N_f + d^2)kT
Here, :math:`\mathbf{r}_i`, :math:`\mathbf{p}_i` and :math:`m_i` are the
position, momentum and mass of particle :math:`i` respectively, :math:`\xi`,
:math:`p_\xi` and :math:`Q` are the position, momentum and mass of the
thermostat, and :math:`\mathbf{h}`, :math:`\mathbf{p}_g` and :math:`W_g` are the
matrix of lattice vectors, matrix of cell velocities and cell mass respectively.
Note that these equations only include one Nose-Hoover thermostat for
simplicity. Conquest uses the Shinoda-Shiga-Mikami splitting of the Liouvillian
:cite:`t-Shinoda2004` to propagate the system. The Liouvillian is decomposed as,
.. math::
iL = iL_r + iL_h + iL_v + iL_\mathrm{bath},
which can be further split,
.. math::
iL_\mathrm{bath} &= iL_\mathrm{box} + iL_\mathrm{particles} \\
iL_\mathrm{box} &= iL_\mathrm{vbox} + iL_\xi + iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}} \\
iL_\mathrm{particles} &= iL_\mathrm{vpart} + iL_\xi + iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}}
Using Liouville's theorem, we have,
.. math::
iL_r &= \sum_{i=1}^N[\mathbf{v}_i + \mathbf{v}_g\mathbf{r}_i]\cdot\nabla_{\mathbf{r}_i} \\
iL_h &= \sum_{\alpha,\beta}\mathbf{v}_{g,\alpha\beta}\mathbf{h}_{\alpha\beta}\frac{\partial}{\partial\mathbf{h}_{\alpha\beta}} \\
iL_v &= \sum_{i=1}^N\left(\frac{\mathbf{F}_i}{m_i}\right)\cdot\nabla_{\mathbf{v}_i} \\
iL_\mathrm{bath} &= iL_\mathrm{vpart} + iL_\mathrm{vbox} + iL_\xi + iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}} \\
&= \sum_{i=1}^N\left[-\left\{\mathbf{v}_g + \frac{1}{N_f}\mathrm{Tr}(\mathbf{v}_g) + v_{\xi_1}\right\}\mathbf{v}_i\right]\nabla_{\mathbf{v}_i} \\
&+ \sum_{\alpha,\beta}\left[\frac{F_\mathrm{box}}{W} - v_{\xi_1}\mathbf{v}_{g,\alpha\beta}\right]\frac{\partial}{\partial\mathbf{v}_{g,\alpha\beta}} \\
&+ \sum_{k=1}^M v_{\xi_k}\frac{\partial}{\partial\xi_k} \\
&+ \left[\frac{F_{\mathrm{NHC}_1}}{Q_1} - v_{\xi_1}v_{\xi_2}\right]\frac{\partial}{\partial v_{\xi_1}} \\
&+ \sum_{k=2}^M\left[\frac{1}{Q_k}(Q_{k-1}v_{\xi_{k-1}}^2 - kT_\mathrm{ext}) - v_{\xi_k}v_{\xi_{k+1}}\right]\frac{\partial}{\partial v_{\xi_k}} \\
&+ \left[\frac{1}{Q_M}(Q_{M-1}v_{\xi_{M-1}}^2 - kT_\mathrm{ext})\right]\frac{\partial}{\partial v_{\xi_M}}
Here we use :math:`M` heat baths in a Nose-Hoover chain. The Trotter-Suzuki
expansion is,
.. math::
e^{iL\Delta t} = e^{iL_\mathrm{bath}\frac{\Delta t}{2}}e^{iL_v\frac{\Delta t}{2}}e^{iL_h\frac{\Delta t}{2}}e^{iL_r\Delta t}e^{iL_h\frac{\Delta t}{2}}e^{iL_v\frac{\Delta t}{2}}e^{iL_\mathrm{bath}\frac{\Delta t}{2}}.
The Liouvillian for the heat baths can be further expanded:
.. math::
e^{iL_\mathrm{particles}\frac{\Delta t}{2}} = e^{\left(iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}}\right)\frac{\Delta t}{4}}e^{\left(iL_\xi + iL_\mathrm{vpart}\right)\frac{\Delta t}{2}}e^{\left(iL_\xi + iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}}\right)\frac{\Delta t}{4}}
Finally, expanding the first propagator in the previous expression, we have,
.. math::
e^{\left(iL_{v_{\xi_1}} + iL_{v_{\xi_k}} + iL_{v_{\xi_M}}\right)\frac{\Delta t}{4}} &= e^{-i\left(-v_{\xi_1}v_{\xi_2}\frac{\partial}{\partial \xi_1} - \sum_{k=2}^Mv_{\xi_k}v_{\xi_{k+1}}\frac{\partial}{\partial \xi_k} - v_{\xi_{M-1}}v_{\xi_M}\frac{\partial}{\partial \xi_M}\right)\frac{\Delta t}{8}} \\
&\times e^{i\left(F_{\mathrm{NHC}_1}\frac{\partial}{\partial v_{\xi_1}} + F_{\mathrm{NHC}_k}\frac{\partial}{\partial v_{\xi_k}} + F_{\mathrm{NHC}_M}\frac{\partial}{\partial v_{\xi_M}}\right)\frac{\Delta t}{4}} \\
&\times e^{-i\left(-v_{\xi_1}v_{\xi_2}\frac{\partial}{\partial \xi_1} - \sum_{k=2}^Mv_{\xi_k}v_{\xi_{k+1}}\frac{\partial}{\partial \xi_k} - v_{\xi_{M-1}}v_{\xi_M}\frac{\partial}{\partial \xi_M}\right)\frac{\Delta t}{8}}
These expressions are directly translated into the integration algorithm.
Go to :ref:`top <theory-md>`.
.. _th_md_weak:
Weak coupling thermostat/barostat
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Instead of modifying the Hamiltonian, the Berendsen-type weak coupling method
:cite:`t-Berendsen1984` involves coupling the ionic degrees of freedom to a an
external temperature and/or pressure bath via "the principle of least local
perturbation consistent with the required global coupling." Thermostatting is
acheived via a Langevin-type equation of motion, in which the system is globally
coupled to a heat bath and subjected to random noise:
.. math::
m_i\ddot{\mathbf{r}}_i = \mathbf{F}_i + m_i \gamma\left(\frac{T_0}{T}-1\right)\dot{\mathbf{r}}_i,
where :math:`\gamma` is a global friction constant chosen to be the same for all
particles. This can be acheived in practice by rescaling the velocities
:math:`\mathbf{v}_i \rightarrow \lambda\mathbf{v}_i`, where :math:`\lambda` is,
.. math::
\lambda = \left[ 1 + \frac{\Delta t}{\tau_T}\left(\frac{T_0}{T}-1\right)\right]^{\frac{1}{2}}
A similar argument can be applied for weak coupling to an external pressure
bath. In the isobaric-isoenthalpic ensemble, the velocity of the particles can
be expressed,
.. math::
\dot{\mathbf{r}} = \mathbf{v} - \frac{\beta(P_0 - P)}{3\tau_P}\mathbf{r},
i.e. the fractional coordinates are scaled by a factor determined by the
difference between the internal and external pressures, the isothermal
compressibility :math:`\beta` and a pressure coupling time constant $\tau_P$.
In the isotropic case, the cell scaling factor :math:`\mu` can be expressed,
.. math::
\mu = \left[ 1 - \frac{\Delta t}{\tau_P}(P_0 - P)\right]^{\frac{1}{3}},
where the compressibility is absorbed into the time time constant
:math:`\tau_P`. Allowing for fluctuations of all cell degrees of freedom, the
scaling factor becomes,
.. math::
\mathbf{\mu} = \mathbf{I} - \frac{\beta\Delta t}{3\tau_P}(\mathbf{P}_0 - \mathbf{P})
While trivial to implement and in general stable, the weak-coupling method does
not recover the correct phase space distribution for the canonical or
isobaric-isothermal ensembles, for which the extended system method is
required. It is no longer supported in CONQUEST, but the concepts are useful.
Go to :ref:`top <theory-md>`.
.. _th_md_svr:
Stochastic velocity rescaling
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Stochastic velocity rescaling (SVR) :cite:`t-Bussi2007` is a modification of the
weak coupling method, in which a correctly constructed random force is added to
enforce the correct NVT (or NPT) phase space distribution. The kinetic energy is
rescaled such that the change in kinetic energy between thermostatting steps is,
.. math::
dK = (\bar{K} - K)\frac{dt}{\tau} + 2\sqrt{\frac{K\bar{K}}{N_f}}\frac{dW}{\sqrt{\tau}}
where :math:`\bar{K}` is the target kinetic energy (external temperature),
:math:`dt` is the time step, :math:`\tau` is the time scale of the thermostat,
:math:`N_f` is the number of degrees of freedom and :math:`dW` is a Wiener
process. Practically, the particle velocities are rescaled by a factor of
:math:`\alpha`, defined via,
.. math::
\alpha^2 = e^{-\Delta t/\tau} + \frac{\bar{K}}{N_fK}\left(1-e^{-\Delta t/\tau}\right)\left(R_1^2 + \sum_{i=2}^{N_f}R_i^2\right) + 2e^{-\Delta t/2\tau}\sqrt{\frac{\bar{K}}{N_fK}\left(1-e^{-\Delta t/\tau}\right)R_1}
Where :math:`R_i` is a set of :math:`N_f` normally distributed random numbers
with unitary variance. This method can be applied to thermostat the NPT ensemble
by barostatting the system with the Parinello-Rahman method, and using the above
expressions, but with additional :math:`R_i`'s for the cell degrees of freedom,
and thermostatting the cell velocities as well as the particle velocities
:cite:`t-Bussi2009`.
.. bibliography:: references.bib
:cited:
:labelprefix: Ta
:keyprefix: t-
:style: unsrt
Go to :ref:`top <theory-md>`.

134
docs/CONQUEST-manual/theory-strucrelax.rst Normal file
View File

@ -0,0 +1,134 @@
.. _theory-strucrelax:
=============================
Structural relaxation: Theory
=============================
Structural relaxation involves optimisation of the ionic coordinates,
optimisation of the simulation cell, or both, with respect to the DFT total
energy or the enthalpy if the cell is not fixed.
.. _th_sr_ions:
Ionic relaxation
----------------
.. _th_sr_lbfgs:
L-BFGS:
~~~~~~~
To be written...
.. _th_sr_cg:
Conjugate gradients
~~~~~~~~~~~~~~~~~~~
The most naive geometry optimisation algorithm is steepest descent: we calculate
the gradient of the DFT total energy (i.e. the force) and propagate the system
in the direction of the steepest gradient (the direction of the force vector)
until the energy stops decreasing. We choose the direction (largest gradient in
this case) and perform a line search. This will be sufficient if the potential
energy surface is well-behaved, but in most cases convergence will require many
iterations. Conjugate gradients is a well-established method the improves upon
steepest descent in the choice of search direction. Without going into too much
detail, we choose a new search direction that is orthogonal to all previous
search directions using the *conjugacy ratio* :math:`\beta`. At iteration
:math:`n`, it is given by,
.. math::
\beta_n = \beta_{n-1} + \frac{\mathbf{f}_n^T\mathbf{f}_n}{\mathbf{f}_{n-1}^T\mathbf{f}_{n-1}}
This is the Fletcher-Reeves formulation; note that :math:`\beta_0 = 0`. We can
then construct the search direction at step :math:`n`, :math:`D_n`,
.. math::
D_n = \beta_n D_{n-1} + \mathbf{f}_n,
and peform the line minimisation in this direction. This process is repeated
until the maximum force component is below some threshold.
Go to :ref:`top <theory-strucrelax>`.
.. _th_sr_qmd:
Quenched MD
~~~~~~~~~~~
The system is propagated in the direction of steepest descent as determined by
the DFT forces, and the velocity is scaled down as the system approaches its
zero-temperature equilibrium configuration.
Go to :ref:`top <theory-strucrelax>`.
.. _th_sr_fire:
FIRE Quenched MD
~~~~~~~~~~~~~~~~
The system is propagated using the modified equation of motion :cite:`t2-Bitzek2006`,
.. math::
\mathbf{\dot{v}}(t) = \mathbf{F}(t)/m -
\gamma(t)|\mathbf{v}(t)|[\mathbf{\hat{v}}(t) - \mathbf{\hat{F}}(t)]
which has the effect of introducing an acceleration in a direction that is
steeper than the current direction of motion. If the power :math:`P(t) =
\mathbf{F}(t)\cdot\mathbf{v}(t)` is positive then the system is moving
"downhill" on the potential energy surface, and the stopping criterion is when
it becomes negative (moving "uphill").
Go to :ref:`top <theory-strucrelax>`.
.. _th_sr_cell_opt:
Cell optimisation
-----------------
When optimising the cell with *fixed fractional ionic coordinates*, the same
conjugate gradients method is used as above, but minimising the enthalpy with
respect to the cell vectors.
Go to :ref:`top <theory-strucrelax>`.
.. _th_sr_comb_opt:
Combined optimisation
---------------------
The ionic and cell degrees of freedom can be relaxed simultaneously by combining
all of their coordinates into a single vector and optimising them with respect
to the enthalpy of the system. However, this atomic forces and total stresses
having numerical values of the same order of magnitude, and changes in ionic
coordinates and cell vectors being of the same order of magnitude. Using the
method of Pfrommer *et al*. :cite:`t2-Pfrommer1997`, the latter can be enforced
by using fractional coordinates for the ionic positions, and fractional lattice
vectors of the form :math:`h = (1 + \epsilon)h_0` where h is the matrix of
lattice vectors, :math:`h_0` is the matrix for some reference configuration and
epsilon is the strain matrix. The "fractional" force on the *i* th atom is then
:math:`\mathbf{F}_i = g\mathbf{f}_i` where :math:`\mathbf{f}_i` is the
DFT-calculated force multiplied by the metric tensor :math:`g = h^Th`. The
"fractional" stress is,
.. math::
f^{(\epsilon)} = -(\sigma + p\Omega)(1 + \epsilon^T)
where :math:`\sigma` is the DFT-calculated stress, :math:`p` is the target
pressure and :math:`\Omega` is the volume. The resulting vector is optimised
using the same conjugate gradients algorithm as before, minimising the enthalpy.
.. bibliography:: references.bib
:cited:
:labelprefix: Tb
:keyprefix: t2-
:style: unsrt
Go to :ref:`top <theory-strucrelax>`.

68
docs/CONQUEST-manual/why-conquest.rst Normal file
View File

@ -0,0 +1,68 @@
.. _overview:
=======================
Overview: Why CONQUEST?
=======================
There are already many DFT codes which are available under open-source
licences. Here we give reasons why you might choose to use CONQUEST.
.. _over_large:
Large-scale simulations
-----------------------
CONQUEST is designed to scale to large systems, either using exact
diagonalisation (with the multisite support function approach, we have
demonstrated calculations on over 3,000 atoms) or with linear scaling
(where calculations on over 2,000,000 atoms have been demonstrated).
Moreover, the same code and basis sets can be used to model systems
from 1 atom to more than 1,000,000 atoms.
.. _over_para:
Efficient parallelisation
-------------------------
CONQUEST is an inherently parallel code, with scaling to more than 800
cores demonstrated for exact diagonalisation, and nearly 200,000 cores
with linear scaling. This scaling enables efficient use of HPC
facilities. CONQUEST (in linear scaling mode, as well as to a certain
extent for exact diagonalisation) scales best with weak scaling:
fixing the number of atoms per core (or thread) and choosing a number
of cores based on the number of atoms.
CONQUEST also offers some OpenMP parallelisation in linear scaling
mode, with relatively low numbers of MPI threads per node, and further
parallelisation performed with OpenMP.
.. _over_on:
Linear scaling
--------------
The ideas of linear scaling have been current for more than twenty
years, but it has proven challenging to make efficient, accurate codes
to implement these ideas. CONQUEST has demonstrated effective linear
scaling (with excellent parallel scaling), though is still somewhat
restricted in the basis sets that can be used. For calculations
beyond 5,000-10,000 atoms with DFT, linear scaling is the only option.
.. _over_basis:
Basis sets
----------
CONQUEST expresses the Kohn-Sham eigenstates or the density matrix
(which are equivalent) in terms of local orbitals called *support
functions*. These support functions are made from one of two basis
sets: pseudo-atomic orbitals (PAOs) or blip functions (B-splines);
the main basis functions in use in CONQUEST are the PAOs. A PAO
generation code is included with the CONQUEST distribution, with
well-defined and reliable default basis sets for most elements.
The simplest choice is to use one PAO for each support function (typically
this allows calculations up to 1,000 atoms). For diagonalisation
beyond this system size, a composite basis is used,
where PAOs from several are combined into a smaller set of support functions
(multi-site support functions, or MSSF). With MSSF, calculations on
3,000+ atoms are possible on HPC platforms. For linear scaling, more
care is required with basis sets (more details can be found :ref:`here
<gs_on>`).