# Basis-set convergence

This lesson was created for an earlier version of BigDFT,
from before the change to the yaml input format.
A new version of the tutorial, in the form of a Jupyter Notebook, is
available at [1].
The corresponding notebook is also available in the source of BigDFT,
at `BIGDFT_ROOT/PyBigDFT/source/tutorials`.

## Contents |

## Running a wavelet computation on a CH_{4} molecule

The purpose of this lesson is to get familiar with the basic variables needed to run a wavelet computation in isolated boundary conditions. At the end of the lesson, one can run a wavelet run, check the amount of needed memory and understand the important part of the output.

### Introduction: running the code

This lesson is based on this skeleton
`input.yaml` file.

Beside this input file, BigDFT requires the atomic positions
for the studied system and optionaly the pseudo-potential files. For the
following tutorial, a methane molecule will be used. The position file
is a simple XYZ file named `posinp.xyz`:

5 angstroemd0 # methane molecule free C 0 0 0 H -0.63169789 -0.63169789 -0.63169789 H +0.63169789 +0.63169789 -0.63169789 H +0.63169789 -0.63169789 +0.63169789 H -0.63169789 +0.63169789 +0.63169789

The pseudo-potential files are following the ABINIT structure and are of GTH or HGH types (see the pseudo-potential file] page on the ABINIT website for several LDA and GGA files and the page of M. Krack on the CP2K server for HGH pseudo for several functionals). The following files may be used for this tutorial: psppar.C and psppar.H.

Running BigDFT is done using the `bigdft`
executable in a standard Unix way, the output being by default the
standard output, it must be redirected to a file or to a pipe, like with the unix command
`tee`:

user@garulfo:~/CH4/$ ls bigdft psppar.C psppar.H input.yaml posinp.xyz user@garulfo:~/CH4/$ ./bigdft | tee screenOutput ...

Warning, to run properly, the pseudo-potential files must be
`psppar.XX` where XX is the symbol used in the position
file. The other files can have user-defined names, as explained in "First runs with BigDFT" lesson.

If the code has been compiled with MPI capabilities (which is enabled by default), running BigDFT on several cores is as easy as run it as a serial job. There is no need to change anything in the input files. The following example shows how to run it on a Debian system with installed OpenMPI on a 4 core machine:

user@garulfo:~/CH4/$ ls bigdft psppar.C psppar.H input.yaml posinp.xyz user@garulfo:~/CH4/$ mpirun -np 4 ./bigdft | tee screenOutput ...

### The wavelet basis set, a convergence study

The wavelet is a systematic basis set (as plane waves are), which means than one can increase arbitrarily the accuracy of the results by varying some parameters.

#### hgrid

The second and third lines of `input.yaml` are used to set
up the basis set. In free boundary conditions, the basis set is characterised
by a spatial expansion and a grid step, as shown in the side
figure.

There is *one float value* on the second line describing the
*grid steps* in the three space directions (*i.e.* x, y and z). These
values are in bohr unit and typically range from 0.3 to 0.65. The
harder the pseudo-potential, the lower value should be set
up. These values are called `hgrid`.

#### crmult, frmult

The third line contains two float values that are two
multiplying factors. They multiply quantities that are chemical
species dependent. The first factor is the most important since it
describes *the spatial expansion* of the basis set (in
yellow on the figure beside). Indeed the basis set
is defined as a set of real space points with non-zero values. These points are
on a global regular mesh and located inside spheres centered on
atoms. The first multiplying factor is called `crmult` for
Coarse grid Radius MULTiplier. Increasing it means that further
spatial expansion is possible for the wavefunctions. Typical values
are 5 to 7.

**Exercise**: run BigDFT for the following
values of `hgrid` and `crmult` and plot the
total energy convergence versus `hgrid`. The final total energy (keyword EKS)
can be retrieved at the end of the screen output, or using this command
``grep FINAL screenOutput``, the value is in Hartree. A
comprehensive explanation of the screen output
will be given later in this tutorial.

hgrid = 0.55bohr / crmult = 3.5 hgrid = 0.50bohr / crmult = 4.0 hgrid = 0.45bohr / crmult = 4.5 hgrid = 0.40bohr / crmult = 5.0 hgrid = 0.35bohr / crmult = 5.5 hgrid = 0.30bohr / crmult = 6.0 hgrid = 0.20bohr / crmult = 7.0

This precision plot shows the systematicity of the wavelet basis set: by improving the basis set, we improve the value of the total energy.

hgrid = 0.55bohr / crmult = 3.5 --> -8.025214Ht hgrid = 0.50bohr / crmult = 4.0 --> -8.031315Ht hgrid = 0.45bohr / crmult = 4.5 --> -8.032501Ht hgrid = 0.40bohr / crmult = 5.0 --> -8.033107Ht hgrid = 0.35bohr / crmult = 5.5 --> -8.033239Ht hgrid = 0.30bohr / crmult = 6.0 --> -8.033300Ht hgrid = 0.20bohr / crmult = 7.0 --> -8.033319Ht

To go further, one can vary `hgrid` and
`crmult` independently. This is shown in the previous
figure with the grey line. The shape of the convergence curve shows that both these
parameters should be modified simoultaneously in order to increase accuracy.
Indeed, there are two kind of errors arising from the
basis set. The first one is due to the fact the basis set can't
account for quickly varying wavefunctions (value of `hgrid` should
be decreased). The second error is the fact that the wavefunctions are
constrained to stay inside the defined basis set (output values are
zero). In the last case `crmult` should be raised.

#### Fine tuning of the basis set

The multi-scale property of the wavelets is used in BigDFT and
a two level grid is used for the calculation. We've seen previously
the coarse grid definition using the the multiplying factor
`crmult`. The second multiplying value on this line of the
input file is used for the fine grid and is called
`frmult`. Like `crmult`, it defines a factor for
the radii used to define the fine grid region where the number of degrees of freedom
is indeed eight times the one of the coarse grid. It allows to define region
near the atoms where the wavefunctions are allowed to vary more
quickly. Typical values for this factor are 8 to 10. It's worth to
note that even if the value of the multiplier is greater than
`crmult` it defines a smaller region due to the fact that
the units which are associated to these radii are significantly different.

The physical quantities used by `crmult` and
`frmult` can be changed in the pseudo-potential by adding
an additional line with two values in bohr. The two values that the
code is using (either computed or read from the pseudo-potential
files) are output in the following way in the screen output:

------------------------------------------------------------------ System Properties Atom N.Electr. PSP Code Radii: Coarse Fine CoarsePSP Calculated File Si 4 10 1.80603 0.43563 0.93364 X H 1 10 1.46342 0.20000 0.00000 X

## Analysing the output

The output of BigDFT is divided into four parts:

- Input values are printed out, including a summary of the different input files (DFT calculation parameters, atom positions, pseudo-potential values...);
- Input wavefunction creation, usually called "input guess";
- The SCF (Self-Consistent Field) loop itself;
- The post SCF calculations including the forces calculation and other possible treatment like a finite size effect estimation or a virtual states determination.

### The system parameters output

All the read values from the different input files are printed out at the program startup. Some additional values are provided there also, like the memory consumption. Values are given for one process, which corresponds to one core in an MPI environment.

#-------------------------------------------------------- Estimation of Memory Consumption Memory requirements for principal quantities (MiB.KiB): Subspace Matrix : 0.1 # (Number of Orbitals: 4) Single orbital : 0.302 # (Number of Components: 38568) All (distributed) orbitals : 2.363 # (Number of Orbitals per MPI task: 4) Wavefunction storage size : 30.617 # (DIIS/SD workspaces included) Nonlocal Pseudopotential Arrays : 0.59 Full Uncompressed (ISF) grid : 8.852 Workspaces storage size : 0.710 Accumulated memory requirements during principal run stages (MiB.KiB): Kernel calculation : 102.802 Density Construction : 81.7 Poisson Solver : 112.694 Hamiltonian application : 81.650 Estimated Memory Peak (MB) : 112 Ion-Ion interaction energy : 9.51544109841576E+00 The overall memory requirement needed for this calculation is thus: 112 MB

In this example, the memory requirement is given for one process run and the peak of memory will be in the initialisation during the Poisson solver kernel creation, while the SCF loop will reach 112MB during the Poisson solver calculation. For bigger systems, with more orbitals, the peak of memory is usually reached during the Hamiltonian application.

**Exercise**: run a small utility program provided with
BigDFT called `bigdft-tool` to estimate the memory requirement
of a run before submitting it to the queue system of a
super-computer. It reads the same input file than the
`bigdft` executable, and is thus convenient to validate inputs.

The executable take one mandatory argument that is the number of cores to run BigDFT on. Try several values from 1 to 6 and discuss the memory distribution.

user@garulfo:~/CH4/$ ls bigdft-tool psppar.C psppar.H input.yaml posinp.xyz user@garulfo:~/CH4/$ ./bigdft-tool -n 2 ...

BigDFT distributes the orbitals over the available processes (the value W does not decrease anymore after 4 processes since there are only 4 bands in our example). This means that running a parallel job with more processors than orbitals will result in a bad speedup. The number of cores involved in the calculation might be however increased via OMP parallelisation, as it is indicated in "Scalability with MPI and OpenMP" lesson.

### The input guess

The initial wavefunctions in BigDFT are calculated using the atomic orbitals for all the electrons of the s, p, d shells, obtained from the solution of the PSP self-consistent equation for the isolated atom.

#----------------------------------- Wavefunctions from PSP Atomic Orbitals Initialization Input Hamiltonian: Total No. of Atomic Input Orbitals : 8 Atomic Input Orbital Generation: - {Atom Type: C, Electronic configuration: {s: [ 2.00], p: [ 2/3, 2/3, 2/3]}} - {Atom Type: H, Electronic configuration: {s: [ 1.00]}}

The corresponding hamiltonian is then diagonalised and the
n_band (`norb` in the code notations) lower eigenfunctions are used to start the SCF loop. BigDFT outputs the
eigenvalues, in the following example, 8 electrons were used in the
input guess and the resulting first fourth eigenfunctions will be used
for a four band calculation.

Input Guess Overlap Matrices: {Calculated: Yes, Diagonalized: Yes} #Eigenvalues and New Occupation Numbers Orbitals: [ {e: -6.501330042644E-01, f: 2.0000}, # 00001 {e: -3.636213579778E-01, f: 2.0000}, # 00002 {e: -3.636197175998E-01, f: 2.0000}, # 00003 {e: -3.636197175998E-01, f: 2.0000}, # 00004 <- Last InputGuess eval, H-L IG gap: 20.6959 eV {e: 3.947505813172E-01, f: 0.0000}, # 00005 <- First virtual eval {e: 3.947523909004E-01, f: 0.0000}, # 00006 {e: 3.947523909004E-01, f: 0.0000}, # 00007 {e: 5.960254618828E-01, f: 0.0000}] # 00008 IG wavefunctions defined : Yes

### The SCF loop

The SCF loop follows a direct minimisation scheme and is made of the following steps:

- Calculate the charge density from the previous wavefunctions.
- Apply the Poisson solver to obtain the Hartree potential from the charges and calculate the exchange-correlation energy and the energy of the XC potential thanks to the chosen functional.
- Apply the resulting hamiltonian on the current wavefunctions.
- Precondition the result and apply a steepest descent or a DIIS
history method. This depends on
`idsx`, not specified in the present`input.yaml`file. It is therefore set to the default value, which is 6 (see the bigdft output). To perform a SD minimisation one should add "`idsx: 0`" to the`input.yaml`file. - Orthogonalise the new wavefunctions.

Then, BigDFT outputs a summary of the parts of the energy:

Energies: {Ekin: 6.66763910179E+00, Epot: -1.04233602201E+01, Enl: 4.35158411655E-01, EH: 1.51771192882E+01, EXC: -3.08539890002E+00, EvXC: -4.03428148190E+00},

Finally the total energy and the square norm of the residue
(gnrm) are printed out. The gnrm value is the stopping criterion. It
can be chosen using `gnrm_cv` in the `input.yaml` file. The default
value (1e-4) is used here and a good value can reach 1e-5.

iter: 6, EKS: -8.03335831460998051E+00, gnrm: 1.20E-03, D: -3.55E-05,

**Exercise**: run ``grep "EKS" screenOutput`` and
look at the convergence rate for our methane molecule.

The minimisation scheme coupled with DIIS (and thanks to the good preconditioner) is a very efficient way to obtain convergence for systems with a gap, even with a very small one. Usual run should reach the 1e-4 stop criterion within 15 to 25 iterations. Otherwise, there is an issue with the system, either there is no gap, or the input guess is too symmetric due to the LCAO diagonalization, specific spin polarization...

### The post-SCF treatments

At the end of the SCF loop, a diagonalisation of the current hamiltonian is done to obtain Kohn-Sham eigenfunctions. The corresponding eigenvalues are also given.

The forces are then calculated.

Some other post-SCF may be done depending on the
`input.dft` Media:

- One can run an estimation of finite-size effects. This is explained in the manual (which is not yet completely updated to recent BigDFT versions).
- One can run a Davidson treatment on the current hamiltonian to obtain the energies (and virtual wavefunctions) of the first unoccupied levels.

**Exercise**: Before going further, review the
`input.dft` file to identify the meaning of the different
lines as explained previously.

1st line, "0.450 0.450 0.450" hx, hy, hz are the grid spacing in the three directions.

2nd line, "5.0 9.0" crmult, frmult define the basis set real space expansion.

3rd line, "1" defines the exchange correlation functional, following the ABINIT numbering convention.

6th line, "1.e-04" is the stop criterion.

7th line, "50 10" the first value is the maximum number of SCF iteration and the second is the maximum number of restart after a fresh diagonalisation if convergence is not reached.

8th line, "6 6" the second value is the length of the DIIS history and should be put to 0 to use SD instead.

**Exercise**: run `bigdft-tool` when varying the
DIIS history length and discuss the memory consumption.

Reducing the DIIS history is a good way to reduce the memory consumption when one cannot increase the number of processes. Of course this implies more iterations in SCF loops.

## Adding a charge

BigDFT can treat charged system without the requirement to add a compensating background like in plane waves.

The additional charge to add to the system is set in the
`input.dft` file at the fourth line. In the following
example an electron has been added (-1):

-1 0.0 0.0 0.0 ncharge efield

**Exercise**: remove the last hydrogen atom in the previous
methane example and modify `input.dft` to add an
electron. Then run BigDFT for an electronic convergence.

One can notice that the total charge in the
system is indeed -8 thanks to the additional charge. The convergence
rate is still good for this CH_{3}^{-} radical since
it is a closed shell system.

## Running a geometry optimisation

In the previous charged example the geometry of the radical is kept the same than for the methane molecule, while it is likely to change. One can thus optimize the geometry with BigDFT.

To run geometry calculations (molecular dynamics, structure
optimisations...) one should add another input file called `input.geopt`. The first line
of this file contains the method to use. Here, we look for a local
minimum so we can use the keyword `LBFGS`. The third line of
this file contains the stopping criteria. There are two stopping
criteria: the first being the number of loops (force evaluations) and the second is the maximum on
forces. For isolated systems, the first criterion is well adapted
while the second is good for periodic boundary conditions.

**Exercise**: take the CH_{3}^{-} radical
`posinp.xyz` file, add the `input.geopt` and run
a geometry optimisation.

The evolution of the forces during relaxation
can be easily obtained running `grep FORCES screenOutput`. At each
iteration, BigDFT outputs a file `posoutXXX.xyz` in
the directory `data` with the geometry of the iteration
XXX. You can visualize it using v_sim
(select all files in the "Browser" tab).