Second tutorial on the Projector AugmentedWave (PAW) technique¶
Generation of PAW atomic datasets¶
This tutorial aims at showing how to compute atomic datasets for the Projector AugmentedWave (PAW) method.
You will learn how to generate these atomic datasets and how to control their softness and transferability. You already should know how to use ABINIT in the PAW case (see the tutorial PAW1 ).
Important
All the necessary input files to run the examples can be found in the ~abinit/tests/ directory where ~abinit is the absolute path of the abinit toplevel directory.
To execute the tutorials, you are supposed to create a working directory (Work*
) and
copy there the input files and the files file of the lesson.
The files file ending with _x (e.g. tbase1_x.files) must be edited every time you start to use a new input file. You will discover more about the files file in section 1.1 of the help file.
To make things easier, we suggest to define some handy environment variables by executing the following lines in the terminal:
export ABI_HOME=Replace_with_the_absolute_path_to_the_abinit_top_level_dir export ABI_TESTS=$ABI_HOME/tests/ export ABI_TUTORIAL=$ABI_TESTS/tutorial/ # Files for base1234, GW ... export ABI_TUTORESPFN=$ABI_TESTS/tutorespfn/ # Files specific to DFPT tutorials. export ABI_TUTOPARAL=$ABI_TESTS/tutoparal/ # Tutorials about parallel version export ABI_TUTOPLUGS=$ABI_TESTS/tutoplugs/ # Examples using external libraries. export ABI_PSPDIR=$ABI_TESTS/Psps_for_tests/ # Pseudos used in examples. export PATH=$ABI_HOME/src/98_main/:$PATH
The examples in this tutorial will use these shell variables so that one can easily copy and paste the code snippets into the terminal (remember to set ABI_HOME first!)
The last line adds the directory containing the executables to your PATH so that one can invoke the code by simply typing abinit in the terminal instead of providing the absolute path.
Finally, to run the examples in parallel with e.g. 2 MPI processes, use mpirun (mpiexec) and the syntax:
mpirun n 2 abinit < files_file > log 2> err
The standard output of the application is redirected to log
while err
collects the standard error
(runtime error messages, if any, are written here).
This tutorial should take about 2h00.
1. The PAW atomic dataset  introduction¶
The PAW method is based on the definition of spherical augmentation regions of radius r_c around the atoms of the system in which a basis of atomic partialwaves \phi_i, of pseudized partialwaves \tphi_i, and of projectors \tprj_i (dual to \tphi_i) have to be defined. This set of partialwaves and projectors functions (plus some additional atomic data) are stored in a socalled PAW dataset. A PAW dataset has to be generated for each atomic species in order to reproduce atomic behavior as accurate as possible while requiring minimal CPU and memory resources in executing ABINIT for the crystal simulations. These two constraints are obviously conflicting.
The PAW dataset generation is the purpose of this tutorial. It is done according the following procedure (all parameters that define a PAW dataset are in bold):

Choose and define the concerned chemical species (name and atomic number).

Solve the atomic allelectrons problem in a given atomic configuration. The atomic problem is solved within the DFT formalism, using an exchangecorrelation functional and either a Schrodinger (default) or scalarrelativistic approximation. This spherical problem is solved on a radial grid. The atomic problem is solved for a given electronic configuration that can be an ionized/excited one.

Choose a set of electrons that will be considered as frozen around the nucleus (core electrons). The others electrons are valence ones and will be used in the PAW basis. The core density is then deduced from the core electrons wave functions. A smooth core density equal to the core density outside a given r_{core} matching radius is computed.

Choose the size of the PAW basis (number of partialwaves and projectors). Then choose the partialwaves included in the basis. The later can be atomic eigenfunctions related to valence electrons (bound states) and/or additional atomic functions, solution of the wave equation for a given l quantum number at arbitrary reference energies (unbound states).

Generate pseudo partialwaves (smooth partialwaves build with a pseudization scheme and equal to partialwaves outside a given r_c matching radius) and associated projector functions. Pseudo partialwaves are solutions of the PAW Hamiltonian deduced from the atomic Hamiltonian by pseudizing the effective potential (a local pseudopotential is built and equal to effective potential outside a $r_{Vloc} matching radius). Projectors and partialwaves are then orthogonalized with a chosen orthogonalization scheme.

Build a compensation charge density used later in order to retrieve the total charge of the atom. This compensation charge density is located inside the PAW spheres and based on an analytical shape function (which analytic form and localization radius r_{shape} can be chosen).
The user can choose between two PAW dataset generators to produce atomic files directly readable by ABINIT.
The first one is the PAW generator ATOMPAW (originally by N. Holzwarth) and
the second one is the UltraSoft USPP generator (originally written by D. Vanderbilt). In this tutorial, we concentrate only on ATOMPAW
.
It is highly recommended to refer to the following papers to understand correctly the generation of PAW atomic datasets:

“Projector augmentedwave method”  [Bloechl1994]

“A projector Augmented Wave (PAW) code for electronic structure”  [Holzwarth2001]

“From ultrasoft pseudopotentials to the projector augmentedwave method”  [Kresse1999]

“Electronic structure packages: two implementations of the Projector AugmentedWave (PAW) formalism”  [Torrent2010]

“Notes for revised form of atompaw code” (by N. Holzwarth)  PDF
2. Use of the generation code¶
Before continuing, you might consider to work in a different subdirectory as
for the other tutorials. Why not Work_paw2
?
cd $ABI_TUTORIAL/Input mkdir Work_paw2 cd Work_paw2
Provided that ABINIT has been compiled with the withdftflavor="...+atompaw"
option,
the ATOMPAW
code is directly available from command line. Try to type:
atompaw
If atompaw vx.y.z
message appears, everything is fine.
Otherwise, you can try:
$ABI_HOME/fallbacks/exports/bin/atompaw
Note
In the following, we name atompaw the ATOMPAW
executable.
How to use ATOMPAW
?
 Edit an input file in a text editor (content of input explained here)
 Run: atompaw < inputfile
Partial waves \phi_i, PS partialwaves \tphi_i and projectors \tprj_i are given in wfn.i
files.
Logarithmic derivatives from atomic Hamiltonian and PAW Hamiltonian
resolutions are given in logderiv.l
files.
A summary of the atomic allelectron computation and PAW dataset properties
can be found in the Atom_name
file (Atom_name
is the first parameter of the input file).
Resulting PAW dataset is contained in:
Atom_name.XCfunc.xml
file Normalized XML file according to the PAW XML specifications.
Atom_name.XCfuncpaw.abinit
file Proprietary legacy format for ABINIT
3. First (and basic) PAW dataset for Nickel¶
Our test case will be nickel; electronic configuration: [1s^2 2s^2 2p^6 3s^2 3p^6 3d^8 4s^2 4p^0].
In a first stage, copy a simple input file for ATOMPAW
in your working directory
(find it in $\ABI_HOME/doc/tutorial/paw2/paw2_assets/Ni.atompaw.input1).
Edit this file.
Ni 28 ! Definition of material GGAPBE scalarrelativistic loggrid 2000 ! Allelectrons calc.: GGA+scalarrelativstic  log.grid with 2000 pts 4 4 3 0 0 0 ! Max. n per angular momentum: 4s 3p 3d 3 2 8 ! Partially occupied states: 3d: occ=8 4 1 0 ! 4p: occ=0 0 0 0 ! End of occupation section c ! 1s: core state c ! 2s: core state c ! 3s: core state v ! 4s: valence state c ! 2p: core state c ! 3p: core state v ! 4p: valence state v ! 3d: valence state 2 ! Max. l for partial waves basis 2.3 ! r_PAW radius y ! Do we add an additional s partial wave ? yes 0.5 ! Reference energy for this new s partial wave (Ryd) n ! Do we add an additional s partial wave ? no y ! Do we add an additional p partial wave ? yes 0. ! Reference energy for this p new partial wave (Ryd) n ! Do we add an additional p partial wave ? no y ! Do we add an additional d partial wave ? yes 0. ! Reference energy for this new d partial wave (Ryd) n ! Do we add an additional s partial wave ? no bloechl ! Scheme for PS partial waves and projectors 3 0. troulliermartins ! Scheme for pseudopotential (l_loc=3, E_loc=0Ry) 2 ! Option for ABINIT file creation default ! All parameters set to default for ABINIT file 0 ! End of file
This file has been built in the following way:
 Allelectron calculation parameters:
 1^{st} line: define the material.
Ni 28
 2^{nd} line: choose the exchangecorrelation functional (LDAPW or GGAPBE)
and select a scalarrelativistic wave equation
(nonrelativistic or scalarrelativistic)
and a (2000 points) logarithmic grid.
GGAPBE scalarrelativistic loggrid 2000
 1^{st} line: define the material.

Electronic configuration:
How many electronic states do we need to include in the computation? Besides the fully and partially occupied states, it is recommended to add all states that could be reached by electrons in the solid. Here, for Nickel, the 4p state is concerned. So we decide to add it in the computation. 3^{rd} line: define the electronic configuration.
A line with the maximum n quantum number for each electronic shell; here
4 4 3
means4s, 4p, 3d
.4 4 3 0 0 0
 Following lines : definition of occupation numbers.
For each partially occupied shell enter the occupation number.
An excited configuration may be useful if the PAW dataset is intended for use in a context where the material is charged (such as oxides). Although, in our experience, the results are not highly dependent on the chosen electronic configuration.
We choose here the [3d^8 4s^2 4p^0] configuration. Only 3d and 4p shells are partially occupied (3 2 8
and4 1 0
lines). A0 0 0
ends the occupation section.3 2 8 4 1 0 0 0 0
 3^{rd} line: define the electronic configuration.
A line with the maximum n quantum number for each electronic shell; here

Selection of core and valence electrons.
In a first approach, select only electrons from outer shells as valence. But, if particular thermodynamical conditions are to be simulated, it is generally needed to include “semicore states” in the set of valence electrons. Semicore states are generally needed with transition metal and rareearth materials.
Note that all wave functions designated as valence electrons will be used in the partialwave basis.
Core shells are designated by a c and valence shells by a v. All s states first, then p states and finally d states.
Here:means:c c c v c c v v
1s core 2s core 3s core 4s valence 2p core 3p core 4p valence 3d valence
 Partialwaves basis generation:
 A line with l_{max} the maximum l for
the partialwaves basis. Here l_{max}=2.
2
 A line with the r_{PAW} radius.
Select it to be slightly less than half the interatomic distance in the solid (as a first choice). Here r_{PAW}=2.3 a.u. If only one radius is input, all others pseudization radii will be equal to r_{PAW} (r_c, r_{core}, r_{Vloc} and r_{shape}).2.3
 Next lines: add additional partialwaves \phi_i if needed.
Choose to have 2 partialwaves per angular momentum in the basis (this choice is not necessarily optimal but this is the most common one; if r_{PAW} is small enough, 1 partialwave per l may suffice).
As a first guess, put all reference energies for additional partialwaves to 0 Rydberg.
For each angular momentum, first add “y” to add an additional partialwave. Then, next line, put the value in Rydberg units. Repeat this for each new partialwave and finally put “n”. Note : For each angular momentum, valence states already are included in the partialwaves basis. Here 4s, 4p and 3d states already are in the basis. In the present file:means that an additional s partialwave at E_{ref}=0.5 Ry as been added,y 0.5 n
means that an additional p partialwave at E_{ref}=0. Ry has been added,y 0. n
means that an additional d partialwave at E_{ref}=0. Ry as been added. Finally, partialwaves basis contains two s, two p and two d partialwaves.y 0. n
 Next line: definition of the generation scheme for pseudo partial
waves \tphi_i, and of projectors \tprj_i.
We begin here with a simple scheme (i.e. “Bloechl” scheme, proposed by P. Blochl [Bloechl1994]). This will probably be changed later to make the PAW dataset more efficient.bloechl
 Next line: generation scheme for local pseudopotential
V_{loc}. In order to get PS partialwaves, the atomic potential has to be “pseudized”
using an arbitrary pseudization scheme.
We choose here a “TroullierMartins” using a wave equation at l_{loc}=3 and E_{loc}=0. Ry. As a first draft, it is always recommended to put l_{loc}=1+lmax.3 0. troulliermartins
 Next two lines: a
2
(two) makesATOMPAW
generate a PAW dataset for ABINIT; The next line contains options for this ABINIT file. “default” set all parameters to their default value.2 default
 A
0
(zero) to end the file.0
 A line with l_{max} the maximum l for
the partialwaves basis. Here l_{max}=2.
At this stage, run ATOMPAW
. For this purpose, simply enter:
atompaw < Ni.atompaw.input1
Lot of files are produced. We will examine some of them.
A summary of the PAW dataset generation process has been written in a file
named Ni
(name extracted from first line of input file).
Open it. It should look like:
Atom = Ni Z = 28 Perdew \ Burke  Ernzerhof GGA Log grid  n,r0,rmax = 2000 2.2810899E04 8.0000000E+01 Scalar relativistic calculation  point nucleus allelectron results core states (zcore) = 18.0000000000000 1 1 0 2.0000000E+00 6.0358607E+02 2 2 0 2.0000000E+00 7.2163318E+01 3 3 0 2.0000000E+00 8.1627107E+00 5 2 1 6.0000000E+00 6.2083048E+01 6 3 1 6.0000000E+00 5.2469208E+00 valence states (zvale) = 10.0000000000000 4 4 0 2.0000000E+00 4.1475541E01 7 4 1 0.0000000E+00 9.0035738E02 8 3 2 8.0000000E+00 6.5223644E01 evale = 185.182300204924 selfenergy contribution = 8.13253645212050 paw parameters: lmax = 2 rc = 2.30969849741149 irc = 1445 Vloc: Normconserving TroullierMartins form; l= 3;e= 0.0000E+00 Projector method: Bloechl Sinc^2 compensation charge shape zeroed at rc Number of basis functions 6 No. n l Energy Cp coeff Occ 1 4 0 4.1475541E01 9.5091493E+00 2.0000000E+00 2 999 0 5.0000000E01 3.2926940E+00 0.0000000E+00 3 4 1 9.0035738E02 8.9594194E+00 0.0000000E+00 4 999 1 0.0000000E+00 1.0836820E+01 0.0000000E+00 5 3 2 6.5223644E01 9.1576176E+00 8.0000000E+00 6 999 2 0.0000000E+00 1.3369075E+01 0.0000000E+00 evale from matrix elements 1.85182309373359203E+02
This generated PAW dataset (contained in Ni.atomicdata, Ni.GGAPBEpaw.abinit or Ni.GGAPBE.xml file) is a first draft. Several parameters have to be adjusted, in order to get accurate results and efficient DFT calculations.
Note
Only Ni.GGAPBEpaw.abinit or Ni.GGAPBEpaw.xml files are directly usable by ABINIT.
4. Checking the sensitivity of results to some parameters¶
4.a. The radial grid¶
Try to select 700 points in the logarithmic grid and check if any noticeable
difference in the results appears.
You just have to replace 2000
by 700
in the second line of Ni.atompaw.input1 file.
Then run:
atompaw < Ni.atompaw.input1
again and look at the Ni file:
evale = 185.182300567432 evale from matrix elements 1.85182301887091256E+02
As you see, results obtained with this new grid are very close to previous ones. We can keep the 700 points grid.
You could decrease the size of the grid; by setting 400 points you should obtain:
evale = 185.182294626845 evale from matrix elements 1.85182337214119599E+02
Small grids give PAW dataset with small size (in kB) and run faster in ABINIT, but accuracy can be affected.
Note
The final r_{PAW} value (rc = ...
in Ni
file) changes with the
grid; just because r_{PAW} is adjusted in order to belong exactly to the radial grid.
By looking in ATOMPAW user’s guide, you can choose to keep it constant.
4.b. The relativistic approximation of the wave equation¶
The scalarrelativistic option should give better results than nonrelativistic one, but it sometimes produces difficulties for the convergence of the atomic problem (either at the allelectron resolution step or at the PAW Hamiltonian solution step). If convergence cannot be reached, try a nonrelativistic calculation (not recommended for high Z materials).
Note
For the following, note that you always should check the Ni
file, especially
the values of valence energy evale
. You can find the valence energy
computed for the exact atomic problem and the valence energy computed with the
PAW parameters evale from matrix elements
. These two results should be in close agreement!
5. Adjusting partialwaves and projectors¶
Examine the AE partialwaves \phi_i, PS partialwaves \tphi_i and projectors \tprj_i.
These are saved in files respectively named wfni
, where i
ranges over the number of partialwaves
used, so 6 in the present example.
Each file contains 4 columns: the radius r in column 1,
the AE partialwave \phi_i(r) in column 2, the PS partialwave \tphi_i(r) in
column 3, and the projector \tprj_i(r) in column 4.
Plot the 3 curves as a function of radius using a plotting tool of your choice.
Below the first s partialwave /projector of the Ni example:

\phi_i should meet \tphi_i near or after the last maximum (or minimum). If not, it is preferable to change the value of the matching (pseudization) radius r_c.

The maxima of \tphi_i and \tprj_i functions should roughly have the same order of magnitude. If not, you can try to get this in three ways:
 Change the matching radius for this partialwave; but this is not always possible (PAW spheres should not overlap in the solid).
 Change the pseudopotential scheme (see later).
 If there are two (or more) partialwaves for the angular momentum l under consideration, decreasing the magnitude of the projector is possible by displacing the references energies. Moving the energies away from each other generally reduces the magnitude of the projectors, but too big a difference between energies can lead to wrong logarithmic derivatives (see following section).
Example: plot the wfn6
file, related to the second d partialwave:
This partialwave has been generated at E_{ref}=0~Ry and orthogonalized with the
first d partialwave which has an eigenenergy equal to 0.65~Ry (see Ni
file).
These two energies are too close and orthogonalization process produces “high” partialwaves.
Try to replace the reference energy for the additional d partialwave.
For example, put E_{ref}=1.~Ry instead of E_{ref}=0.~Ry (line 24 of Ni.atompaw.input1
file).
Run ATOMPAW
again and plot wfn6
file:
Now the PS partialwave and projector have the same order of magnitude!
Important
Note again that you should always check the evale
values in Ni
file and make
sure they are as close as possible.
If not, choices for projectors and/or partialwaves are certainly not judicious.
6. Examine the logarithmic derivatives¶
Examine the logarithmic derivatives, i.e., derivatives of an lstate
\frac{d(log(\Psi_l(E))}{dE} computed for the exact atomic problem and with the PAW dataset.
They are printed in the logderiv.l
files. Each logderiv.l
file corresponds to an
angular momentum l and contains three columns of data: the
energy, the logarithmic derivative of the lstate of the exact atomic problem
and of the pseudized problem. In our Ni example, l=0, 1 or 2.
The logarithmic derivatives should have the following properties:

The 2 curves should be superimposed as much as possible. By construction, they are superimposed at the 2 energies corresponding to the 2 l partialwaves. If the superimposition is not good enough, the reference energy for the second l partialwave should be changed.

Generally a discontinuity in the logarithmic derivative curve appears at 0~Ry<=E_0<=4~Ry. A reasonable choice is to choose the 2 reference energies so that E_0 is in between.

Too close reference energies produce “hard” projector functions (see section 5). But moving reference energies away from each other can damage accuracy of logarithmic derivatives
Here are the three logarithmic derivative curves for the current dataset:
As you can see, except for l=2, exact and PAW logarithmic derivatives do not match!
According to the previous remarks, try other values for the references
energies of the s and p additional partialwaves.
First, edit again the Ni.atompaw.input1
file and put E_{ref}=3~Ry for the
additional s state (line 18); run ATOMPAW
again. Plot the logderiv.0
file.
You should get:
Then put E_{ref}=4~Ry for the second p state (line 21); run ATOMPAW
again.
Plot again the logderiv.1
file. You should get:
Now, all PAW logarithmic derivatives match with the exact ones in a reasonable interval.
Note
It is possible to change the interval of energies used to plot logarithmic
derivatives (default is [5;5]) and also to compute them at more points
(default is 200). Just add the following keywords at the end of the SECOND
LINE of the input file if you want ATOMPAW
to output logarithmic derivatives
for energies in [10;10] at 500 points:
logderivrange 10 10 500
Additional information concerning logarithmic derivatives: ghost states
Another possible issue could be the presence of a discontinuity in the PAW logarithmic derivative curve at an energy where the exact logarithmic derivative is continuous. This generally shows the presence of a ghost state.
 First, try to change to value of reference energies; this sometimes can make the ghost state disappear.
 If not, it can be useful to change the pseudopotential scheme. Normconserving pseudopotentials are
sometimes too attractive near r=0.
 A 1^{st} solution is to change the quantum number used to generate the normconserving pseudopotential. But this is generally not sufficient.
 A 2^{nd} solution is to select a
ultrasoft
pseudopotential, freeing the norm conservation constraint (simply replacetroulliermartins
byultrasoft
in the input file).  A 3^{rd} solution is to select a simple
bessel
pseudopotential (replacetroulliermartins
bybessel
in the input file). But, in that case, one has to noticeably decrease the matching radius r_{Vloc} if one wants to keep reasonable physical results. Selecting a value of r_{Vloc} between 0.6~r_{PAW} and 0.8~r_{PAW} is a good choice.
To change the value of r_{Vloc}, one has to explicitely put all matching radii: r_{PAW}, r_{shape}, r_{Vloc} and r_{core}; see user’s guide.
 Last solution : try to change the matching radius r_c for one (or both) l partialwave(s). In some cases, changing r_c can remove ghost states.
In most cases (changing pseudopotential or matching radius), one has to restart the procedure from step 5.
To see an example of ghost state, use the
$ABI_HOME/doc/tutorial/paw2_assets/Ni.ghost.atompaw.input file and run it with ATOMPAW
.
Look at the l=1 logarithmic derivatives (logderiv.1
file). They look like:
Now, edit the Ni.ghost.atompaw.input file and replace troulliermartins
by
ultrasoft
. Run ATOMPAW
again… and look at logderiv.1
file.
The ghost state has moved!
Edit again the file and replace ultrasoft
by bessel
; then change the 17^{th}
line 2.0 2.0 2.0 2.0
by 2.0 2.0 1.8 2.0
(decreasing the r_{Vloc} radius from 2.0 to 1.8).
Run ATOMPAW
: the ghost state disappears!
Start from the original state of Ni.ghost.atompaw.input file and put 1.8
for
the matching radius of p states (put 1.8
on lines 31 and 32).
Run ATOMPAW
: the ghost state disappears!
7. Testing the “efficiency” of a PAW dataset¶
Let’s use again our Ni.atompaw.input1 file for Nickel (with all our modifications). You get a file Ni.GGAPBEpaw.abinit containing the PAW dataset designated for ABINIT.
To test the efficiency of the generated PAW dataset, we finally will use ABINIT! You are about to run a DFT computation and determine the size of the plane wave basis needed to reach a given accuracy. If the cutoff energy defining the plane waves basis is too high (higher than 20 Hartree), some changes have to be made in the input file.
Copy $ABI_TUTORIAL/Input/tpaw2_x.files and $ABI_TUTORIAL/Input/tpaw2_1.in in your working directory. Edit tpaw2_1.in, and activate the 8 datasets. Run ABINIT with them.
tpaw2_1.in tpaw2_1.out tpaw2_i tpaw2_o tpaw2 Ni.GGAPBEpaw.abinit
# Nickel ferromagnetic fcc structure for testing ecut convergence # ndtset 8 # Uncomment this line to recover the results for the tutorial ndtset 1 # Comment this line to recover the results for the tutorial ecut: 8. ecut+ 2. pawecutdg 40. nstep 50 tolvrs 1.0d9 occopt 7 tsmear 0.0075 ntypat 1 acell 3*3.52 angstrom rprim 0.0 1/2 1/2 1/2 0.0 1/2 1/2 1/2 0.0 natom 1 typat 1 xred 0.0 0.0 0.0 znucl 28 nband 14 nsppol 2 spinat 0. 0. 4. ngkpt 6 6 6 nshiftk 4 shiftk 1/2 1/2 1/2 1/2 0.0 0.0 0.0 1/2 0.0 0.0 0.0 1/2 getwfk 1 prteig 0 prtden 0 ## After modifying the following section, one might need to regenerate the pickle database with runtests.py r #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% [files] #%% files_to_test = #%% tpaw2_1.out, tolnlines= 4, tolabs= 1.1e4, tolrel= 2.0e3, fld_options=medium #%% psp_files = Ni.GGAPBEpaw.abinit.bloechl #%% [paral_info] #%% max_nprocs = 6 #%% [extra_info] #%% authors = M. Torrent #%% keywords = PAW #%% description = Nickel ferromagnetic fcc structure for testing ecut convergence #%%<END TEST_INFO>
ABINIT computes the total energy of ferromagnetic FCC Nickel for several values of ecut.
At the end of output file, you get this:
ecut1 8.00000000E+00 Hartree ecut2 1.00000000E+01 Hartree ecut3 1.20000000E+01 Hartree ecut4 1.40000000E+01 Hartree ecut5 1.60000000E+01 Hartree ecut6 1.80000000E+01 Hartree ecut7 2.00000000E+01 Hartree ecut8 2.20000000E+01 Hartree etotal1 3.9300291581E+01 etotal2 3.9503638785E+01 etotal3 3.9583278145E+01 etotal4 3.9613946329E+01 etotal5 3.9623543087E+01 etotal6 3.9626889070E+01 etotal7 3.9628094989E+01 etotal8 3.9628458879E+01
etotal
convergence (at 1 mHartree) is achieve for 18<=e_{cut}<=20 Hartree
etotal
convergence (at 0,1 mHartree) is achieve for e_{cut}>22 Hartree
This is not a good result for a PAW dataset; let’s try to optimize it.

1^{st} possibility: use
vanderbilt
projectors instead ofbloechl
ones.
Vanderbilt’s projectors generally are more localized in reciprocal space than Bloechl’s ones .
Keywordbloechl
has to be replaced byvanderbilt
in theATOMPAW
input file and r_c values have to be added at the end of the file (one for each PS partialwave). See this input file: $ABI_HOME/doc/tutorial/paw2_assets/Ni.atompaw.input.vanderbilt. 
2^{nd} possibility: use
RRKJ
pseudization scheme for projectors.
Use this input file forATOMPAW
: $ABI_HOME/doc/tutorial/paw2_assets/Ni.atompaw.input2. As you can seebloechl
has been changed bycustom rrkj
and 6 r_c values have been added at the end of the file, each one corresponding to the matching radius of one PS partialwave. Repeat the entire procedure (ATOMPAW
+ABINIT
)… and get a new ABINIT output file.
Note: You have check again at log derivatives.
Ni 28 ! Definition of material GGAPBE scalarrelativistic loggrid 700 ! Allelectrons calc.: GGA+scalarrelativstic  log.grid with 700 pts 4 4 3 0 0 0 ! Max. n per angular momentum: 4s 3p 3d 3 2 8 ! Partially occupied states: 3d: occ=8 4 1 0 ! 4p: occ=0 0 0 0 ! End of occupation section c ! 1s: core state c ! 2s: core state c ! 3s: core state v ! 4s: valence state c ! 2p: core state c ! 3p: core state v ! 4p: valence state v ! 3d: valence state 2 ! Max. l for partial waves basis 2.3 ! r_PAW radius y ! Do we add an additional s partial wave ? yes 3. ! Reference energy for this new s partial wave (Ryd) n ! Do we add an additional s partial wave ? no y ! Do we add an additional p partial wave ? yes 4. ! Reference energy for this p new partial wave (Ryd) n ! Do we add an additional p partial wave ? no y ! Do we add an additional d partial wave ? yes 1. ! Reference energy for this new d partial wave (Ryd) n ! Do we add an additional s partial wave ? no custom rrkj ! Scheme for PS partial waves and projectors 3 0. troulliermartins ! Scheme for pseudopotential (l_loc=3, E_loc=0Ry) 2.3 ! r_c matching radius for first s partial wave 2.3 ! r_c matching radius for second s partial wave 2.3 ! r_c matching radius for first p partial wave 2.3 ! r_c matching radius for second p partial wave 2.3 ! r_c matching radius for first d partial wave 2.3 ! r_c matching radius for second d partial wave 2 ! Option for ABINIT file creation default ! All parameters set to default for ABINIT file 0 ! End of file
ecut1 8.00000000E+00 Hartree ecut2 1.00000000E+01 Hartree ecut3 1.20000000E+01 Hartree ecut4 1.40000000E+01 Hartree ecut5 1.60000000E+01 Hartree ecut6 1.80000000E+01 Hartree ecut7 2.00000000E+01 Hartree ecut8 2.20000000E+01 Hartree etotal1 3.9600401638E+01 etotal2 3.9627563690E+01 etotal3 3.9627901781E+01 etotal4 3.9628482371E+01 etotal5 3.9628946655E+01 etotal6 3.9629072497E+01 etotal7 3.9629079826E+01 etotal8 3.9629097793E+01
etotal
convergence (at 1 mHartree) is achieve for 12 <= e_{cut} <= 14 Hartree
etotal
convergence (at 0,1 mHartree) is achieve for 16 <= e_{cut} <= 18 Hartree
This is a reasonable result for a PAW dataset!
 3^{rd} possibility: use enhanced polynomial pseudization scheme for projectors.
Edit $ABI_HOME/doc/tutorial/paw2_assets/Ni.atompaw.input2 and replacecustom rrkj
bycustom polynom2 7 10
. Repeat the entire procedure (ATOMPAW
+ABINIT
)… and look at theecut
convergence.
Optional_exercise¶
Let’s go back to Vanderbilt projectors.
Repeat the procedure (ATOMPAW
+ ABINIT
) with the
$ABI_HOME/doc/tutorial/paw2_assets/Ni.atompaw.input.vanderbilt file.
Ni 28 ! Definition of material GGAPBE scalarrelativistic loggrid 700 ! Allelectrons calc.: GGA+scalarrelativstic  log.grid with 700 pts 4 4 3 0 0 0 ! Max. n per angular momentum: 4s 3p 3d 3 2 8 ! Partially occupied states: 3d: occ=8 4 1 0 ! 4p: occ=0 0 0 0 ! End of occupation section c ! 1s: core state c ! 2s: core state c ! 3s: core state v ! 4s: valence state c ! 2p: core state c ! 3p: core state v ! 4p: valence state v ! 3d: valence state 2 ! Max. l for partial waves basis 2.3 ! r_PAW radius y ! Do we add an additional s partial wave ? yes 3. ! Reference energy for this new s partial wave (Ryd) n ! Do we add an additional s partial wave ? no y ! Do we add an additional p partial wave ? yes 4. ! Reference energy for this p new partial wave (Ryd) n ! Do we add an additional p partial wave ? no y ! Do we add an additional d partial wave ? yes 1. ! Reference energy for this new d partial wave (Ryd) n ! Do we add an additional s partial wave ? no vanderbilt ! Scheme for PS partial waves and projectors 3 0. troulliermartins ! Scheme for pseudopotential (l_loc=3, E_loc=0Ry) 2.3 ! r_c matching radius for first s partial wave 2.3 ! r_c matching radius for second s partial wave 2.3 ! r_c matching radius for first p partial wave 2.3 ! r_c matching radius for second p partial wave 2.3 ! r_c matching radius for first d partial wave 2.3 ! r_c matching radius for second d partial wave 2 ! Option for ABINIT file creation default ! All parameters set to default for ABINIT file 0 ! End of file
As you can see ABINIT convergence cannot be achieved!
You can try whatever you want with radii and/or references energies in the
ATOMPAW
input file: ABINIT always diverges!
The solution here is to change the pseudization scheme for the local pseudopotential.
Try to replace the troulliermartins
keyword by ultrasoft
.
Repeat the procedure (ATOMPAW
+ ABINIT
).
ABINIT can now reach convergence!
Results are below:
ecut1 8.00000000E+00 Hartree ecut2 1.00000000E+01 Hartree ecut3 1.20000000E+01 Hartree ecut4 1.40000000E+01 Hartree ecut5 1.60000000E+01 Hartree ecut6 1.80000000E+01 Hartree ecut7 2.00000000E+01 Hartree ecut8 2.20000000E+01 Hartree etotal1 3.9609714395E+01 etotal2 3.9615187859E+01 etotal3 3.9618367959E+01 etotal4 3.9622476129E+01 etotal5 3.9624707476E+01 etotal6 3.9625234480E+01 etotal7 3.9625282524E+01 etotal8 3.9625330757E+01
etotal
convergence (at 1 mHartree) is achieved for 14 <= e_{cut} <= 16 Hartree
etotal
convergence (at 0,1 mHartree) is achieved for 20 <= e_{cut} <= 22 Hartree
Note
You could have tried the bessel
keyword instead of ultrasoft
one.
Summary of convergence results
Final_remarks

The localization of projectors in reciprocal space can (generally) be predicted by a look at
tprod.i
files. Such a file contains the curve of as a function of q (reciprocal space variable). q is given in Bohr^{1} units; it can be connected to ABINIT plane waves cutoff energy (in Hartree units) by: e_{cut}=\frac{q_{cut}^2}{4}. These quantities are only calculated for the bound states, since the Fourier transform of an extended function is not welldefined. 
Generating projectors with Blochl’s scheme often gives the guaranty to have stable calculations.
ATOMPAW
ends without any convergence problem and DFT calculations run without any divergence (but they need high plane wave cutoff). Vanderbilt projectors (and even morecustom
projectors) sometimes produce instabilities during the PAW dataset generation process and/or the DFT calculations but are more efficient. 
In most cases, after having changed the projector generation scheme, one has to restart the procedure from step 5.
8 Testing against physical quantities¶
The last step is to examine carefully the physical quantities obtained with our PAW dataset.
Copy $ABI_TUTORIAL/Input/tpaw2_2.in in your working directory. Edit it, activate the 8 datasets, change tpaw2_x.files to use $ABI_HOME/doc/tutorial/paw2_assets/Ni.GGAPBEpaw.abinit.rrkj psp file (obtained from Ni.atompaw.input2 file). Run ABINIT (this may take a while…).
# Nickel ferromagnetic fcc structure for testing ecut convergence # ndtset 7 # Uncomment this line to reproduce the results of the tutorial ndtset 1 # Comment this line to reproduce the results of the tutorial acell: 3*3.5150 angstrom acell+ 3*0.0025 angstrom ecut 15. pawecutdg 40. ecutsm 0.5 nstep 30 tolvrs 1.0d12 occopt 7 tsmear 0.0075 ntypat 1 rprim 0.0 1/2 1/2 1/2 0.0 1/2 1/2 1/2 0.0 natom 1 typat 1 xred 0.0 0.0 0.0 znucl 28 nband 12 nsppol 2 spinat 0. 0. 4. ngkpt 8 8 8 nshiftk 4 shiftk 1/2 1/2 1/2 1/2 0.0 0.0 0.0 1/2 0.0 0.0 0.0 1/2 getwfk 1 prteig 0 prtden 0 ## After modifying the following section, one might need to regenerate the pickle database with runtests.py r #%%<BEGIN TEST_INFO> #%% [setup] #%% executable = abinit #%% [files] #%% files_to_test = #%% tpaw2_2.out, tolnlines= 15, tolabs= 3.2e4, tolrel= 2.0e1, fld_options=easy #%% psp_files = Ni.GGAPBEpaw.abinit.rrkj #%% [paral_info] #%% max_nprocs = 6 #%% [extra_info] #%% authors = M. Torrent #%% keywords = PAW #%% description = Nickel ferromagnetic fcc structure for testing ecut convergence #%%<END TEST_INFO>
ABINIT computes the converged ground state of ferromagnetic FCC Nickel for several volumes around equilibrium.
Plot the etotal
vs acell
curve:
From this graph and output file, you can extract some physical quantities:
Equilibrium cell parameter: a0 = 3.523 angstrom Bulk modulus: B = 190 GPa Magnetic moment at equilibrium: mu = 0.60
Compare these results with published results:
 allelectron GGAFLAPW from [Kresse1999]:
a0 = 3.52 angstrom B = 200 GPa mu = 0.60
 GGAPAW with VASP code from [Kresse1999]:
a0 = 3.52 angstrom B = 194 GPa mu = 0.61
 Experimental results from from [Dewaele2008]:
a0 = 3.52 angstrom B = 183 GPa
You should always compare results with allelectron ones (or other PAW computations), not with experimental ones
Additional remark:
It can be useful to test the sensitivity of results to some ATOMPAW
input parameters
(see user’s guide for details on keywords):
 The analytical form and the cutoff radius r_{shape} of the shape function used in
compensation charge density definition,
By default a
sinc
function is used but agaussian
shape can have an influence on results.Bessel
shapes are efficient and generally need a smaller cutoff radius (r_{shape}=0.8~r_{PAW}).  The matching radius $r_{core} used to generate the pseudo core density from atomic core density,
 The inclusion of additional (“semicore”) states in the set of valence electrons,
 The pseudization scheme used to get pseudopotential Vloc(r).
All these parameters have to be meticulously checked, especially if the PAW dataset is used for nonstandard solid structures or thermodynamical domains.
Optional_exercise
Let’s add 3s and 3p semicore states in PAW dataset!
Repeat the procedure (ATOMPAW
+ ABINIT
) with $ABI_HOME/doc/tutorial/paw2_assets/Ni.atompaw.input.semicore
file. The execution time is a bit longer as more electrons have to be treated by ABINIT.
Look at a_0, B or \mu variation.
Note: this new PAW dataset has a smaller r_{PAW} radius (because semicore states are localized).
a0 = 3.519 angstrom B = 194 GPa mu = 0.60
9 The Real Space Optimization (RSO)  experienced users¶
In this section, an additional optimization of the atomic data is proposed which can contribute, in some cases, to an acceleration of the convergence on ecut. This optimization is not essential to produce efficient PAW datasets but can be useful. We advise experienced users to try it.
The idea is quite simple: when expressing the different atomic radial functions (\phi_i, \tphi_i, \tprj_i) on the plane wave basis, the number of plane waves depends on the “locality” of these radial functions in reciprocal space.
In this paper a method to enforce the locality (in reciprocal space) of projectors \tprj_i is presented; the projectors expressed in reciprocal space \tprj_i(g) are modified according to the following scheme: The reciprocal space is divided in 3 regions:

If g < g_{max}, \tprj_i(g) is unchanged

If g > \gamma, \tprj_i(g) is set to zero

If $ g_{max} < g < \gamma$, \tprj_i(g) is modified so that the contribution of \tprj_i(r) is conserved with an error W (as small as possible).
The above transformation of \tprj_i(g) is only possible if \tprj_i(r) is defined outside
the spherical augmentation region up to a radius R_0, with R_0 > r_c.
In practice we have to:
 Impose an error W (W is the maximum error admitted on total energy)
 Adjust g_{max} according to E_{cut} (g_{max} <= E_{cut})
 Choose \gamma so that 2 g_{max} < \gamma < 3 g_{max}
and let the ATOMPAW
code apply the transformation to \tprj_i and deduce R_0 radius.
You can test it now. In your working directory, reuse the dataset $ABI_HOME/doc/tutorial/paw2_assets/Ni.atompaw.input3 (Bloechl’s projectors). Replace the ABINIT options (penultimate line):
2 defaults
rsoptim 8. 2 0.0001
Run ATOMPAW. You get a new psp file for ABINIT. Run ABINIT with it using the $ABI_TUTORIAL/Input/tpaw2_1.in file. Compare the results with those obtained in section 7.
You can try several values for g_{max} (keeping \frac{\gamma}{g_{max}} and W constant) and compare the efficiency of the atomic data; do not forget to test physical properties again.
How to choose the RSO parameters?
\frac{\gamma}{g_{max}} = 2 and 0.0001 < W < 0.001 is a good choice. g_{max} has to be adjusted. The lower g_{max} the faster the convergence is but too low g_{max} can produce unphysical results.