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-2007 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.txt .
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 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 ../Input/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 ../Input/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 ../Input/tpaw1_1.in and ../Input/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
vw1_1.outAariable. 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
Using nearly the same input file as in section 1.a, you might
perform a PAW calculation for diamond.
Use the ../Input/tpaw1_x.files file as "files" file,
and ../Input/tpaw1_2.in file as "input" file. Copy them 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. Concerning the input file, note that the input variable
pawecutdg has been specified.
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.1027092138E+01
etotal2 -1.1489998109E+01
etotal3 -1.1510395783E+01
etotal4 -1.1510819813E+01
etotal5 -1.1511126450E+01
etotal6 -1.1511171822E+01
etotal7 -1.1511185750E+01
etotal8 -1.1511191275E+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 a typical
value (40 Ha, larger or equal to the different values of ecut; it
must always be greater or equal to 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.1510395779E+01
etotal2 -1.1510356584E+01
etotal3 -1.1510414462E+01
etotal4 -1.1510414645E+01
etotal5 -1.1510389843E+01
etotal6 -1.1510395783E+01
etotal7 -1.1510395793E+01
etotal8 -1.1510391819E+01
etotal9 -1.1510391826E+01
etotal10 -1.1510391782E+01
In any case, the energy is converged to 1mHa .
For pawecutdg=35 Ha (5th dataset)
or pawecutdg=40 Ha (6th
dataset), the energy is converged to 1 microHa
. 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.
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 ../Input/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.1510389843E+01
etotal2 -1.1510389843E+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 norm-conserving pseudopotentials, and 0.7 for PAW. For big unit cells or complex systems, it might be needed
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 useful 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.1510389483E+01
which is in very good agreement with the result etotal1 of
tpaw1_5.out (etotal1 -1.1510389843E+01; the
tolerance on energy was 1.10-7 Ha). For which
concerns
the atomic positions, we obtain:
xred -2.5069738915E-02 -2.5069738915E-02 -2.5069738915E-02
2.2506973892E-01 2.2506973892E-01 2.2506973892E-01
that is to say that the position of the second C atom with respect to the first 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