# ABINIT, PAW lesson of the tutorial:

## The PAW method tutorial.

This lesson aims at showing how to perform a calculation in the frame of the PAW method.

You will learn what is a PAW calculation and what are the main input variables to govern convergency and numerical efficiency.
It is supposed that you already know how to use ABINIT in the norm conserving pseudopotential case

This lesson should take about 1 hour to be done.

##### Copyright (C) 2000-2006 ABINIT group (MT,FJ)
This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
For the initials of contributors, see ~abinit/doc/developers/contributors .

Goto : ABINIT home Page | Suggested acknowledgments | List of input variables | Bibliography

Help files : New user's guide | Abinis (main) | Abinis (respfn) | Mrgddb | Anaddb | AIM (Bader) | Cut3D | Optic | Mrgscr

### Content of the PAW lesson

- 0. Summary of the paw method
- 1. The convergence in ecut
- 2. The convergence in pawecutdg
- 3. The mixing in the self-consistent loop
- 4. The other input variables wfoptalg, ortalg
- 5. State of the art of PAW in ABINIT

####

### 0. Summary of the PAW method.

The PAW method means Projector Augmented Wave method. It has been introduced by Peter Blöchl in 1994. As he says, "the projector augmented wave method is an extension of augmented wave methods and the pseudopotential approach, which combines their traditions into a unified electronic structure method".
It is based on a transformation that connects the "true" wavefunctions {psi} with "auxiliary" (or "pseudo") soft wavefunctions {tpsi}:
psi= tpsi+sum_i (phi_i-tphi_i) <p_i|tpsi>,
where the index i runs on all the nuclei of the system. This relation is based on the definition of atomic spheres of radius rc, around the atoms of the system in which the phi_i form a base of atomic wavefunctions, tphi_i are "pseudized" wavefunctions of phi_i, and p_i are dual functions of the tphi_i called projectors.
It is therefore possible to write every quantity depending on psi (density, energy, hamiltonian) as a function of tpsi and to find tpsi solving self-consistent equations.
The PAW method has two main advantages:
- from tpsi, it is always possible to obtain the true ("all electron") wavefunction psi
- if tphi is a soft pseudization of phi, the convergency is comparable to an ultrasoft pseudopotential one.

From a practical point of view, a PAW calculation differs from a norm-conserving pseudopotential one in using a special atomic data file that contain the phi_i, tphi_i and p_i and plays the same role as a pseudopotential file in the norm-conserving case.

It is highly recommanded to read the following papers to understand correctly the bases of the PAW method:
"Projector augmented-wave method", by P. Blöchl, Phys. Rev. B, 50, 17953 (1994)
"Comparison of the projector augmented-wave, pseudopotential, and linearized augmented-plane-wave formalisms for density-functional calculations of solids", by N. Holzwarth et al., Phys. Rev. B, 55, 2005 (1997)
"From ultrasoft pseudopotentials to the projector augmented-wave method", by G. Kresse and D. Joubert, Phys. Rev. B, 59, 1758 (1999)

The implementation of the PAW method in ABINIT follows the notations and formulations of the paper from G. Kresse and D. Joubert.

###

### 1. The convergence in ecut.

1.a Computing the convergence in ecut for diamond in the norm-conserving case

Before beginning, you might work out again steps 1.a to 1.d of the lesson_1.

You might create a subdirectory of the ~abinit/tests/tutorial/Input directory, and use it for the tutorial. In what follows, the name of files are mentioned as if you were in this subdirectory.

The input file ../tpaw1_1.in is an example of a file that contains data for computing the convergence in ecut for diamond. You might use the file ../tpaw1_1.files (with a standard norm-conserving pseudopotential) as "files" file, and get the corresponding output file (it is available as ../Refs/tpaw1_1.out).
Copy the files ../tpaw1_1.in and ../tpaw1_1.files in your Work directory, and issue e.g.

    ../../abinis < tpaw1_1.files > ,log

The run should take around 1 minute on a PC 2.8 GHz.
In the meantime, you can read the input file and see that there is no paw input variable. There are 8 datasets, for which ecut increases from 5 Ha to 40 Ha by step of 5 Ha.

You should obtain the values (output file tpaw1_1.out) :

etotal1 -1.1226497581E+01
etotal2 -1.1840984990E+01
etotal3 -1.2002784186E+01
etotal4 -1.2063510936E+01
etotal5 -1.2074691486E+01
etotal6 -1.2076216642E+01
etotal7 -1.2076541167E+01
etotal8 -1.2076816719E+01

You can check that the energy is converged to 1mHa for ecut=30 Ha (6th dataset)``.

1.b Computing the convergence in ecut for diamond in the PAW case

Use the same input file as in section 1.a, you might use the ../tpaw1_x.files file as "files" file. Copy it in your work directory. The only difference with the tpaw1_1.files file is that the pseudopential filename has been replaced by a PAW atomic data filename for carbon. Issue e.g.

    ../../abinis < tpaw1_x.files > ,log

The run should take about 1 minute also, as for the norm-conserving pseudopotential case.

You should obtain the values (output file tpaw1_2.out) :

etotal1 -1.1027285148E+01
etotal2 -1.1490154915E+01
etotal3 -1.1510603329E+01
etotal4 -1.1511008498E+01
etotal5 -1.1511311237E+01
etotal6 -1.1511357181E+01
etotal7 -1.1511371110E+01
etotal8 -1.1511371801E+01
You can check that the energy is converged to 1mHa for ecut=15 Ha (3rd dataset)``. So with the same input, for the same convergence, a PAW calculation for diamond needs a cutoff divided by 2, compared to a norm-conserving pseudopotential calculation.
1.c Overview of the output file in the PAW case
Compared to an output file for norm-conserving pseudopotentials, an output file for PAW contains the following specific topics:
- some specific default PAW input variables (ngfftdg, pawecutdg, and useylm) are mentionned in the section:
-outvars: echo values of preprocessed input variables --------
- a specific description of the atomic data file (you might follow the tutorial PAW2, devoted to the building of the PAW atomic data, for a complete analysis of the file):
Pseudopotential format is: paw2
basis_size (lnmax)= 4 (lmn_size= 8), orbitals= 0 0 1 1
Spheres core radius: rc_sph= 1.11201554
4 radial meshes are used:
- mesh 1: r(i)=AA*[exp(BB*(i-1))-1], size= 467 , AA= 0.41313E-03 BB= 0.16949E-01
- mesh 2: r(i)=AA*[exp(BB*(i-1))-1], size= 532 , AA= 0.41313E-03 BB= 0.16949E-01
- mesh 3: r(i)=AA*[exp(BB*(i-1))-1], size= 520 , AA= 0.41313E-03 BB= 0.16949E-01
- mesh 4: r(i)=AA*[exp(BB*(i-1))-1], size= 596 , AA= 0.41313E-03 BB= 0.16949E-01
Shapefunction is SIN type: shapef(r)=[sin(pi*r/rc)/(pi*r/rc)]**2
Radial grid used for partial waves is grid 1
Radial grid used for projectors is grid 2
Radial grid used for (t)core density is grid 3
Radial grid used for Vloc is grid 4
- the value of the compensation charge: as in the case of calculation with ultrasoft pseudopotentials, PAW calculations are not norm conserving. The value of the compensation charge is given both when calculated on the (spherical) augmentation grid and on the regular FFT grid. A discussion on these two values will be done in the next section (2)
PAW TEST:
==== Compensation charge inside spheres ============
The following values must be equal...
Compensation charge over spherical meshes = -0.039476479366873
Compensation charge over fine fft grid = -0.039479379648493
- Informations concerning the non-local term (pseudopotential strength Dij ) and the spherical density matrix (augmentation waves occupancies Rhoij)
==== Results concerning PAW augmentation regions ====
Total pseudopotential strength Dij (hartree):
Atom # 1
...
Atom # 2
...
NON-SYMMETRIZED (augmentation) waves occupancies Rhoij:
Atom # 1
...
Atom # 2
...
SYMMETRIZED (augmentation) waves occupancies Rhoij:
Atom # 1
...
Atom # 2
...
- the decomposition of total energy both by direct calculation and double counting calculation:
--------------------------------------------------------------------------------
Components of total free energy (in Hartree) :
Kinetic energy = 6.69582841719232E+00
Hartree energy = 7.54658655606795E-01
XC energy = -3.75586908431682E+00
Ewald energy = -1.28639953292939E+01
PspCore energy = 6.26705130230738E-01
Loc. psp. energy= -3.69753073313453E+00
Spherical terms = 1.21291981613464E+00
>>>>> Internal E= -1.10272831275808E+01
-kT*entropy = -4.72488802703871E-16
>>>>>>>>> Etotal= -1.10272831275808E+01
"Double-counting" decomposition of free energy:
Band energy = 9.51664611278093E-01
Ewald energy = -1.28639953292939E+01
PspCore energy = 6.26705130230738E-01
Dble-C XC-energy= 1.20778231776548E-01
Spherical terms = 1.37562207905920E-01
>>>>> Internal E= -1.10272851481026E+01
-kT*entropy = -4.72488802703871E-16
>>>> Etotal (DC)= -1.10272851481026E+01
>Total energy in eV = -3.00067641174963E+02
>Total DC energy in eV = -3.00067696156160E+02

###

### 2. The convergence in pawecutdg.

In a norm-conserving calculation, the (plane wave) density grid is twice bigger than the (plane wave) wavefunction grid, for each dimension. In a PAW calculation, the (plane wave) density grid is tunable thanks to the input variable pawecutdg (paw ecut double grid). This is because one has to map the compensation density in the spheres onto the regular FFT grid. The number of points of the Fourier grid must be high enough to preserve the accuracy of the compensation density starting from its value on the augmentation (spherical grid). In fact, pawecutdg is used to generate in a convenient way the Fourier grid. An alternative is to use directly the input variable ngfftdg. In the calculations made in previous section, pawecutdg was set to its default value (= 2.0 ecut; it must always be greater than the cutoff energy).

Use now the input file ../tpaw1_3.in . You might still use tpaw1_x.files as "files" file (do not forget to modify it). The only difference with the tpaw1_1.in is that ecut is fixed to 15 Ha, while pawecutdg runs from 15 to 60 Ha. The run should take again about 1 minute on a PC 2.8 GHz.

You should obtain the values (file tpaw1_3.out) :

etotal1 -1.1510584473E+01
etotal2 -1.1510545268E+01
etotal3 -1.1510603146E+01
etotal4 -1.1510603329E+01
etotal5 -1.1510578527E+01
etotal6 -1.1510584467E+01
etotal7 -1.1510584478E+01
etotal8 -1.1510580503E+01
etotal9 -1.1510580510E+01
etotal10 -1.1510580466E+01

The energy is still converged to 1mHa . For ecutdg=35 Ha (5th dataset), the energy is converged to 0.1 mHa``. The steps in the convergency are due to the changes in the grid size (see ngfftdg). The convergency of the compensation charge follows this convergency.

We want now to check the convergence in ecut with a fixed value of 35 Ha for pawecutdg. Modify the file tpaw1_1.in, setting pawecutdg to 35 Ha, and select only the datasets 1 to 7 (in dataset 8 pawecutdg is lower than ecut); you can use tpaw1_4.in instead. You should obtain the values :

etotal1 -1.1027282439E+01
etotal2 -1.1490181274E+01
etotal3 -1.1510578527E+01
etotal4 -1.1511001924E+01
etotal5 -1.1511308631E+01
etotal6 -1.1511354016E+01
etotal7 -1.1511367943E+01
You can check again that the energy is converged to 1mHa for ecut=15 Ha (3rd dataset)``, and to 0.1 mHa for ecut=25 Ha (5th dataset).
Associated with the input variable pawecutdg is the input variable ngfftdg: it gives the size of the FFT grid associated with pawecutdg. Note that pawecutdg is only useful to define the FFT grid for the density (called sometimes fine grid) in a convenient way. You can therefore tune directely ngfftdg to define the size of the FFT grid for the density.
Other input variables associated with pawecutdg : mqgriddg and bxctmindg. They play the same role as mqgrid and boxcutmin, but for the FFT grid associated with the density. It is recommended not to change the default values.

###

3. The mixing in the self-consistent loop.

The mixing in the self-consistent loop is a crucial point to minimize the number of steps to achieve convergency. In a norm-conserving pseudopotentials calculation, the mixing is done on the potential or the density, thanks to the input variable iscf (please, read again the different options of this input variable).
For PAW calculations, the mixing in the self consistent loop includes both a mixing of the density or of the potential on the FFT grid AND a mixing of the augmentation waves occupancies rhoij, that govern the spherical part of the calculation.
In this section we want to compare the mixing on the potential with the mixing on the density. Use now the input file ../tpaw1_5.in . You might still use tpaw1_x.files as "files" file (do not forget to modify it). One dataset is made with iscf=13 and the second one with iscf=3 . You should obtain the values :
etotal1 -1.1510578527E+01
etotal2 -1.1510578527E+01
You can check that the converged value are the same. However, the energy is not computed in the same way:
For the mixing on the potential, the energy is the "direct" energy: etotal=ekin+eloc+ehart+exc+enl+...
For the mixing on the density, the energy is the double counting energy: etotal= eband-ehart+exc-evxc+eloc+enl+...
At convergency, these two values are the same. But it has been observed that the converged value was reached more rapidly by the direct energy, when the mixing is on the potential, and by the double counting energy when the mixing is on the density. So the default choice in the output file is to print the direct energy, when the mixing is on the potential, and the double counting energy when the mixing is on the density.
Mixing on the density:
--------------------------------------------------------------------------------
Components of total free energy (in Hartree) :
Kinetic energy = 6.97039533234181E+00
Hartree energy = 9.50329717667828E-01
XC energy = -3.83170967286506E+00
Ewald energy = -1.28639953292939E+01
PspCore energy = 6.26705130230738E-01
Loc. psp. energy= -5.14605722579266E+00
Spherical terms = 1.78378507234504E+00
>>>>> Internal E= -1.15105469753662E+01
-kT*entropy = -2.91924984698783E-16
>>>>>>>>> Etotal= -1.15105469753662E+01
"Double-counting" decomposition of free energy:
Band energy = 6.76162136549466E-01
Ewald energy = -1.28639953292939E+01
PspCore energy = 6.26705130230738E-01
Dble-C XC-energy= -7.62323364353650E-02
Spherical terms = 1.26781871119307E-01
>>>>> Internal E= -1.15105785278298E+01
-kT*entropy = -2.91924984698783E-16
>>>> Etotal (DC)= -1.15105785278298E+01
>Total energy in eV = -3.13217919552002E+02
>Total DC energy in eV = -3.13218778138220E+02
--------------------------------------------------------------------------------
Mixing on the potential:
--------------------------------------------------------------------------------
Components of total free energy (in Hartree) :
Kinetic energy = 6.97035681631770E+00
Hartree energy = 9.50383771131953E-01
XC energy = -3.83172998567850E+00
Ewald energy = -1.28639953292939E+01
PspCore energy = 6.26705130230738E-01
Loc. psp. energy= -5.14625454651019E+00
Spherical terms = 1.78395561691978E+00
>>>>> Internal E= -1.15105785268824E+01
-kT*entropy = -3.20165780347442E-16
>>>>>>>>> Etotal= -1.15105785268824E+01
"Double-counting" decomposition of free energy:
Band energy = 6.76215506603142E-01
Ewald energy = -1.28639953292939E+01
PspCore energy = 6.26705130230738E-01
Dble-C XC-energy= -7.62822662400389E-02
Spherical terms = 1.26778734045191E-01
>>>>> Internal E= -1.15105782246549E+01
-kT*entropy = -3.20165780347442E-16
>>>> Etotal (DC)= -1.15105782246549E+01
>Total energy in eV = -3.13218778112441E+02
>Total DC energy in eV = -3.13218769888411E+02
--------------------------------------------------------------------------------
The mixing on rhoij (spherical part) is done on the same way as the mixing on the density or the potential. In particular, it is possible to adjust the input variable diemix in case the convergency is difficult to obtain. The default value is diemix=1.0. For big unit cells or complex systems, it is always better to decrease this value (in PAW). You can try to decrease diemix in tpaw1_5.in file and rerun abinit.
If you want to disconnect the mixing on rhoij and the mixing on density or potentiel (not recommended!), you can use the input variable pawsphmix, that governs only the mixing on rhoij. The default value is pawsphmix=diemix.
If you want to save CPU time for complex systems with lot of atoms and/or electrons you can use the input variable pawlmix. Only parts of augmentation occupancies (rhoij quantities) associated with l<=pawlmix are mixed. Other parts of augmentation occupancies are not included in the mixing process. This option is usefull to save CPU time but does influence the numerical accuracy of the results.

### 4. The other input variables

. In a PAW calculation, we have to compute the XC potential and energy inside a sphere, on a spherical grid defined by a radial grid and a discretization of the theta and phi angles. The densities in the spheres can therefore be defined as a function of r, theta and phi. The size of the radial grid is given in the atomic data file. The input variable pawxcdev allows to chose how to compute the angular dependence of the densities in the spheres: if it is set to 0, the XC term in the spherical part of energy is totally computed on the angular mesh. This last one is defined by pawntheta, that gives the number of theta directions, and pawnphi, that gives the number of phi directions. This possibility is available only for a LDA calculation. Another possible choice is to develop the XC term in the spherical part of the energy onto lm-moments. This correspond to pawxcdev=1 and is the default value. For a LDA calculation, it is completely equivalent to the former case but the speed of the calculation is better. For a GGA calculation, a development into lm-moment to the first order is done. This is a current approximation: the angular variation of the density is not taken into account for the GGA contribution.
. When we use a development of the density into lm-moments (pawxcdev=1), it happens that due to the symmetry of the system, several lm-moments are equal to zero. The input variable pawnzlm allows to test the moments at the first step of the electronic iterations and to impose to compute only nonzero moments for the following steps. This is the default value (pawnzlm=1). The moments of the densities in the spheres are printed in the log file when pawprtvol >=2 (see below).
. The expansion of the densities in angular momenta is performed up to l=pawlcutd-1. The choice made for this variable DOES have a bearing on the numerical accuracy of the results, and, as such, should be the object of a convergence study. The default value is 10.
. In the PAW method, the spheres are theoretically supposed not to overlap. However, a small overlap only introduces a small error due to the continuity and derivability of all quantities at the sphere radii. Moreover, a greater sphere radius allows a better softening of the augmentation pseudized wavefunctions, and therefore a decreasing of the cutoff. Finally, when atoms are moving, an overlap may arise during the calculation. That is why the input variable pawovlp allows the spheres to overlap (pawovlp is the voluminal percentage of overlap allowed) or not (when pawovlp=0 %). The default value is 5 %, so that the code stops when the overlap percentage is greater than 5%. There is therefore a compromise to do between accuracy (no overlap and "hard" atomic functions) and speed (significant overlap and "soft" atomic functions).
. The input variable pawprtvol governs the verbosity of the log file.

### 5. State of the art of PAW in ABINIT

From version 4.6, PAW in ABINIT is available for:

- the calculation of energy, forces and stresses (ground state calculation)
- LDA (PZ) or GGA(PBE) calculation
- spin polarized calculation (no spin-orbit coupling for the moment)
- geometry optimisation

As an example, we are going to optimize the geometry of a diamond cell starting from a displaced atom. Use now the input file ../tpaw1_6.in . You might still use tpaw1_x.files as "files" file (do not forget to modify it). The second carbon atom has been put at (0.2; 0.2;0.2).
You can check (file tpaw1_6.out) that the energy at the end of the Broyden steps is:
etotal -1.1510578167E+01

which is in very good agreement with the result etotal1 of tpaw1_5.out (etotal1 -1.1510578527E+01; the tolerance on energy was 1.10-7 Ha). For which concerns the atomic positions, we obtain:

xred -2.5069641107E-02 -2.5069641107E-02 -2.5069641107E-02
2.2506964111E-01 2.2506964111E-01 2.2506964111E-01
that is to say that the position of the second C atom is (0.25014,0.25014,0.25014), with a tolerance on forces tolmxf=5.10-5 ha/bohr.

Goto : ABINIT home Page | Suggested acknowledgments | List of input variables | Bibliography

Help files : New user's guide | Abinis (main) | Abinis (respfn) | Mrgddb | Anaddb | AIM (Bader) | Cut3D | Optic | Mrgscr