This lesson aims at showing how to compute atomic data files for the projector-augmented-wave method.
You will learn how to generate the atomic data and
what the main
variables are to govern their softness and transferability.
It is supposed you already know how to use ABINIT
in the
PAW case
This lesson should take about 1h30.
Choose and define the concerned chemical species (name and atomic number).
Solve the atomic all-electrons problem in a given atomic configuration. The atomic problem is solved within the DFT formalism, using an exchange-correlation functional and either a Schrödinger (default) or scalar-relativistic approximation. It is a spherical problem and it 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 partial-waves and projectors). Then choose the partial-waves included in the basis. The later can be atomic eigen-functions 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 partial-waves (smooth partial-waves build with a pseudization scheme and equal to partial-waves outside a given r_{c} matching radius) and associated projector functions. Pseudo partial-waves 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 partial-waves 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).
It is highly recommended to refer to the following papers to understand correctly the generation of PAW atomic datasets:
Before continuing, you might consider to work in a different subdirectory as for the other lessons. Why not "Work_paw2" ?
Provided that ABINIT has been compiled with the "--with-dft-flavor="...+atompaw"" option,
the ATOMPAW code is
directly available from command line.
First, just try to type: atompaw
if "atompaw vx.y.z" message appears,
everything is fine.
Otherwise, you can try "~abinit_compilation_directory/fallbacks/exports/bin/atompaw-abinit"
In any case, in the following, we name atompaw the ATOMPAW executable.
Our test case will be NICKEL (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 ~abinit/doc/tutorial/lesson_paw2/Ni.atompaw.input1). Edit this file.This file has been built in the following way:
1-All-electrons calculation:
- A line with the maximum n quantum number for each electronic shell; here "4 4 3" means 4s, 4p, 3d.
- Definition of occupation numbers:Partial-waves basis generation:
The generated PAW
dataset (contained in Ni.atomicdata, Ni.GGA-PBE-paw.abinit
or Ni.GGA-PBE.xml
file) is a
first draft.
Several parameters have
to be adjusted, in order to get accurate results and efficient DFT
calculations.
Note that only Ni.GGA-PBE-paw.abinit file is directly usable by ABINIT.
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
You could decrease the size of the grid; by setting 400 points you should obtain:
Small grids give PAW dataset with small size (in kB) and run faster in ABINIT, but accuracy can be affected.
- Note that the final r_{PAW} value ("rc = ..." in Ni file) change 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.
- Also note that, if the results are difficult to get converged (some error produced by ATOMPAW), you should try a linear grid…
The relativistic approximation of the wave equation:
Scalar-relativistic option should give better results than non-relativistic one, but it sometimes produces difficulties for the convergence of the atomic problem (either at the all-electrons resolution step or at the PAW Hamiltonian solution step). If convergence cannot be reached, try a non-relativistic calculation (not recommended for high Z materials)
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!
Examine
the partial-waves, PS partial-waves and projectors.
These are
saved in files named
wfni, where
i ranges over the number of partial
waves used, so 6 in the present example. Each file contains 4 columns: the radius in
column 1, the partial wave φ_{i} in column 2,
the PS partial wave~φ_{i} in column 3,
and the projector~
p_{i} in column 4. Plot the three curves as a function of radius using a plotting tool of
your choice.
Here is the first s- partial wave /projector of the Ni example:
The φ_{i} should meet the~φ_{i }near or after the last maximum (or minimum). If not, it is preferable to change the value of the matching (pseudization) radius.
The maxima of the ~φ_{i }and ~p_{i} functions should have the same order of magnitude (but need not agree exactly). If not, you can try to get this in three ways:
Example: plot the wfn6 file, concerning the second d- partial wave:
This
partial wave has been
generated at E_{ref}=0 Ry and
orthogonalized with the first d-
partial wave which has an eigenenergy equal to -0.65Ry (see Ni file).
These two energies are too close and orthogonalization process produces
"high" partial waves.
Try
to
replace the reference energy for
the additional d-
partial wave. For example, put E_{ref}=1. instead of E_{ref}=0.
(line 24 of Ni.atompaw.input1
file). Run ATOMPAW again and
plot wfn6 file:
Now the PS partial wave and projector have the same order of magnitude !
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 partial waves
certainly are not judicious.
Examine the logarithmic
derivatives, i.e., derivatives of an l-state d(log(Ψ_{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 angular momentum quantum number l,
and contains three columns of data: the energy, the logarithmic
derivative of the l-state 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 two energies corresponding to the two l partial-waves. If the superimposition is not good enough, the reference energy for the second l partial-wave should be changed.
Generally a discontinuity in the logarithmic derivative curve appears at 0<=E_{0}<=4 Rydberg. 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 (as previously seen in section 5). But moving reference energies away from each other can damage accuracy of logarithmic derivatives
Summary of convergency results:
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 cut-off energy (in Hartree units) by: ecut=q_{cut}^{2}/4. These quantities are only calculated for the bound states, since the Fourier transform of an extended function is not well-defined.
Generating projectors with Blöchl’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 cut-off). Vanderbilt projectors (and even more “custom” projectors) sometimes produce instabilities during the PAW dataset generation process and/or the DFT calculations…
In most cases, after having changed the projector generation scheme, one has to restart the procedure from step 5.
Finally,
the last step is to examine carefully
the physical quantities obtained with the PAW
dataset.
Copy ~abinit/tests/tutorial/Input/tpaw2_2.in
in your working directory. Edit it, to activate the eight datasets (instead of one).
Use the
~abinit/doc/tutorial/lesson_paw2/Ni.GGA-PBE-paw.abinit.rrkj
psp file (it
has been obtained from Ni.atompaw.input2
file).
Modify
tpaw2_x.files
file according to these new files.
Run ABINIT (this may take a
while...).
ABINIT computes the converged ground state of ferromagnetic FCC Nickel for several volumes around equilibrium.
Plot the etotal vs acell curve:You should always compare results with all-electrons ones (or other PAW computations), not with experimental ones...
It can be useful to test the sensitivity of results to some ATOMPAW input parameters (see user's guide for details on keywords):
All these parameters have to be meticulously checked, especially if the PAW dataset is used for non-standard solid structures or thermodynamical domains.
Optional exercise: let's add 3s and 3p semi-core states in PAW dataset !
The idea is
quite simple: when expressing the different atomic radial functions
(φ_{i},~φ_{i}, ~p_{i})
on the plane
waves basis, the number of plane waves depends on the
"locality" of these radial functions in reciprocal space.
In the
following reference (we suggest to read it): R.D. King-Smith,
M.C. Payne, J.S. Lin, Phys. Rev. B 44,
13063 (1991)
You can try several values for g_{max} (keeping γ/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 ?