The reference directory that contains the example files for the
tutorial is no more $ABI_TUTORIAL(as for the basic tutorials and the
specialized, non-DFPT ones), but $ABI_TUTORESPFN.
This will be the case for all the DFPT based part of the tutorial.

You can copy it in the Work_rf1 directory and, as usual, change its name to trf1_1.files
and replace the occurrences of trf1_x by trf1_1.

Note that two pseudopotentials are mentioned in this files file: one
for the Aluminum atom, and one for the Arsenic atom.
The first listed in trf1_x.files (for Al) will define the first type of atom of the input file
(see input variables typat and ntypat) and the second (for As) will define the second type of atom.
It might be the first time that you encounter this situation (more than one type of atoms) in the
tutorials, at variance with the first four basic tutorials.

You can also copy the file $ABI_TUTORESPFN/Input/trf1_1.in in Work_rf1.
This is your input file. You should read it carefully.

# Crystalline AlAs : computation of the total energy##Specific to ground state calculation
tolvrs 1.0d-18 # SCF stopping criterion########################################################################Common input variables#Definition of the unit cell
acell 3*10.61 # This is equivalent to 10.61 10.61 10.61
rprim 0.0 0.5 0.5 # In tutorials 1 and 2, these primitive vectors 0.5 0.0 0.5 # (to be scaled by acell) were 1 0 0 0 1 0 0 0 1 0.5 0.5 0.0 # that is, the default.#Definition of the atom types
ntypat 2# There are two types of atom
znucl 1333# The keyword "znucl" refers to the atomic number of the # possible type(s) of atom. The pseudopotential(s) # mentioned in the "files" file must correspond# to the type(s) of atom. Here, type 1 is the Aluminum,# type 2 is the Arsenic.#Definition of the atoms
natom 2# There are two atoms
typat 12# The first is of type 1 (Al), the second is of type 2 (As).
xred # This keyword indicate that the location of the atoms# will follow, one triplet of number for each atom0.0 0.0 0.0 # Triplet giving the REDUCED coordinate of atom 1.0.25 0.25 0.25 # Triplet giving the REDUCED coordinate of atom 2.#Gives the number of band, explicitely (do not take the default)
nband 4# For an insulator (if described correctly as an insulator # by DFT), there is no need to include conduction bands # in response-function calculations#Exchange-correlation functional
ixc 1# LDA Teter Pade parametrization#Definition of the planewave basis set
ecut 3.0 # Maximal kinetic energy cut-off, in Hartree#Definition of the k-point grid
kptrlatt -4 44# In cartesian coordinates, this grid is simple cubic, and4 -4 4# actually corresponds to the so-called 8x8x8 Monkhorst-Pack grid.44 -4 # It might as well be obtained through the use of# ngkpt, nshiftk and shiftk .#Definition of the SCF procedure
nstep 15# Maximal number of SCF cycles
diemac 9.0 # Although this is not mandatory, it is worth to# precondition the SCF cycle. The model dielectric# function used as the standard preconditioner# is described in the "dielng" input variable section.# The dielectric constant of AlAs is smaller that the one of Si (=12).## After modifying the following section, one might need to regenerate the pickle database with runtests.py -r#%%<BEGIN TEST_INFO>#%% [setup]#%% executable = abinit#%% test_chain = trf1_1.in, trf1_2.in, trf1_3.in, trf1_4.in#%% [files]#%% files_to_test = #%% trf1_1.out, tolnlines= 0, tolabs= 0.000e+00, tolrel= 0.000e+00#%% psp_files = 13al.981214.fhi, 33as.pspnc#%% [paral_info]#%% max_nprocs = 2#%% [extra_info]#%% authors = X. Gonze#%% keywords = #%% description = Crystalline AlAs : computation of the total energy#%%<END TEST_INFO>

It drives a single self-consistent calculation of the total energy of AlAs to generate the corresponding
self-consistent charge density and wavefunctions, that will be used for the DFPT calculations.

Note that the value of tolvrs is rather stringent.
This is because the wavefunctions determined by the present run will be used later as starting
point of the DFPT calculation.
However, the number of steps, nstep, in this example file has been set to 15, and
you will see that this is not enough to reach the target tolvrs.
In production runs, you should choose a large enough value of nstep to reach your tolvrs target.
In the present tutorial, due to portability concerns related to automatic testing, we could
not allow a larger nstep.
This minor problem with some tutorial examples was mentioned briefly in
a side note to the answer to question 1 of tutorial 1 - just before
this section.

Danger

Do not follow blindly all examples of the tutorials: always check the convergence of your calculations while in production!

You will work at fixed ecut (3Ha) and k-point grid, defined by kptrlatt (the 8x8x8 Monkhorst-Pack grid).
In real life you should do a convergence test with respect to both parameters.
We postpone the discussion of the accuracy of these choices and
the choice of pseudopotential to the end of the
fifth section of this tutorial.
They give acceptable but not very accurate results such that the running time is reasonable for a tutorial.

You should make the run (a few seconds):

abinit < trf1_1.files > log 2> err

The resulting main output file, trf1_1.out, should be similar to the one below:

This output file is not very long, so you can quickly read it entirely.
Note that one obtains the following value for the energy, in the final echo section:

etotal -9.7626837450E+00

However, we will rely later, for the purpose of doing finite differences, on a more accurate (more digits) value of this
total energy, that can be found about a dozen of lines before this final echo:

>>>>>>>>> Etotal= -9.76268374500280E+00

The output file also mentions that the forces on both atoms vanish.

The run that you just made will be considered as defining a ground-state
configuration, on top of which responses to perturbations will be computed.
The main output of this ground-state run is the wavefunction file trf1_1o_WFK,
that you can already rename as trf1_1i_WFK to use it as input wave function for the next runs.

Warning

So, in the corresponding files file for all the following runs, at third line,
pay attention to keep “trf1_1i”. By contrast, for the second run, you
should change the first line from “trf1_1” to “trf1_2”, and do similarly for the second,
fourth and fifth lines of this file.

2 Frozen-phonon calculation of a second derivative of the total energy¶

We will now compute the second derivative of the total energy with
respect to an atomic displacement by different means.
For that purpose, you must first read sections 0 and the first paragraph of section 1
of the respfn help file (the auxiliary help file, that deals specifically with the DFPT features).
We will explain later, in more detail, the signification of the different
input parameters introduced in section 1 of the respfn help file.

For the time being, in order to be able to perform a direct comparison with the result of
a DFPT calculation, we choose as a perturbation the displacement of the Al
atom along the first axis of the reduced coordinates.

You can copy the file $ABI_TUTORESPFN/Input/trf1_2.in in Work_rf1.
This is your input file. You should open it and briefly look at the two
changes with respect to trf1_1.in:
the change of xred, and the reading of the wavefunction file, using the irdwfk input variable.

Then, you can make the run, following the same command as before, with a different files file, referring to trf1_2.in.
The symmetry is lowered with respect to the ground-state geometry, so that the number of k-points
increases a lot, and of course, the CPU time.

# Crystalline AlAs : computation of the total energy# and forces in a distorted geometry
irdwfk 1# Read input wavefunctions#Specific to ground state calculation
tolvrs 1.0d-12 # SCF stopping criterion########################################################################Common input variables#Definition of the unit cell
acell 3*10.61 # This is equivalent to 10.61 10.61 10.61
rprim 0.0 0.5 0.5 # In tutorials 1 and 2, these primitive vectors 0.5 0.0 0.5 # (to be scaled by acell) were 1 0 0 0 1 0 0 0 1 0.5 0.5 0.0 # that is, the default.#Definition of the atom types
ntypat 2# There are two types of atom
znucl 1333# The keyword "znucl" refers to the atomic number of the # possible type(s) of atom. The pseudopotential(s) # mentioned in the "files" file must correspond# to the type(s) of atom. Here, type 1 is the Aluminum,# type 2 is the Arsenic.#Definition of the atoms
natom 2# There are two atoms
typat 12# The first is of type 1 (Al), the second is of type 2 (As).
xred # This keyword indicate that the location of the atoms# will follow, one triplet of number for each atom0.001 0.0 0.0 # Triplet giving the REDUCED coordinate of atom 1.0.25 0.25 0.25 # Triplet giving the REDUCED coordinate of atom 2.#Gives the number of band, explicitely (do not take the default)
nband 4# For an insulator (if described correctly as an insulator # by DFT), there is no need to include conduction bands # in response-function calculations#Exchange-correlation functional
ixc 1# LDA Teter Pade parametrization#Definition of the planewave basis set
ecut 3.0 # Maximal kinetic energy cut-off, in Hartree#Definition of the k-point grid
kptrlatt -4 44# In cartesian coordinates, this grid is simple cubic, and4 -4 4# actually corresponds to the so-called 8x8x8 Monkhorst-Pack grid.44 -4 # It might as well be obtained through the use of# ngkpt, nshiftk and shiftk .#Definition of the SCF procedure
nstep 15# Maximal number of SCF cycles
diemac 9.0 # Although this is not mandatory, it is worth to# precondition the SCF cycle. The model dielectric# function used as the standard preconditioner# is described in the "dielng" input variable section.# The dielectric constant of AlAs is smaller that the one of Si (=12).## After modifying the following section, one might need to regenerate the pickle database with runtests.py -r#%%<BEGIN TEST_INFO>#%% [setup]#%% executable = abinit#%% test_chain = trf1_1.in, trf1_2.in, trf1_3.in, trf1_4.in#%% input_prefix = trf1_1o#%% [files]#%% files_to_test = #%% trf1_2.out, tolnlines= 0, tolabs= 0.000e+00, tolrel= 0.000e+00#%% psp_files = 13al.981214.fhi, 33as.pspnc#%% [paral_info]#%% max_nprocs = 2#%% [extra_info]#%% authors = X. Gonze#%% keywords = #%% description = #%% Crystalline AlAs : computation of the total energy#%% and forces in a distorted geometry#%%<END TEST_INFO>

The change of reduced coordinate (xred) of the Al atom along the first axis was
rather small (1/1000 = 0.001), and we can make an estimate of the second derivative of
the total energy with respect to the reduced coordinate thanks to finite-difference formulas.

We start first from the total-energy difference.
The total energy is symmetric with respect to that perturbation, so that it has no linear term.
The difference between the ground-state value (-9.76268374500280E+00 Hartree) of
the previous run, and the perturbed value (-9.76268124105767E+00 Hartree) of
the present one, is thus one half of the square of the coordinate change
(0.001) times the second derivative of total energy (2DTE).
From these number, the 2DTE is 5.0078902589 Hartree.

Alternatively, we can start from the reduced gradients. The value of the
reduced gradient with respect to a displacement of the Al atom along the first
reduced axis is 0.005007937776 Ha. At first order, this quantity is the
product of the 2DTE by the reduced coordinate difference. The estimate of the
2DTE is thus 5.007937776 Ha. The agreement with the other estimate is rather good (4.10^-5 Hartree).

However, it is possible to do much better, thanks to the use of a higher-order finite-difference formula.
For this purpose, one can perform another calculation, in which the change of
reduced coordinate along the first axis is 0.002, instead of 0.001.
The doubling of the perturbation allows for a rather
simple higher-order estimation, as we will see later.
The results of this calculation are as follows:

From these data, taking into account that the perturbation was twice stronger,
the same procedure as above leads to the values 5.00800219 Hartree (from
finite difference of energy) and 5.008153784 Hartree (from finite difference
of forces, the value 0.010016307568 has to be multiplied by 1000/2).

The combination of these data with the previous estimate can be done thanks to an
higher-order finite-difference formula, in which the difference of estimations
(the largest perturbation minus the smallest one) is divided by three, and
then subtracted from the smallest estimation.

As far as the total-energy estimation is concerned, the difference is 0.00011194 Ha, which divided by
three and subtracted from 5.00789025 Hartree, gives 5.00785294 Hartree.
The same higher-order procedure for force estimates gives 5.0078657 Hartree. So,
the agreement between total-energy estimate and force estimate of the 2DTE can
be observed up to the 6th digit, inclusive.

Before comparing this result with the 2DTE directly computed from the DFPT
capabilities of ABINIT, a last comment is in order. One can observe that the
action-reaction law is fulfilled only approximately by the system. Indeed, the
force created on the second atom, should be exactly equal in magnitude to the
force on the first atom. The values of dE/dt, mentioned above show a small,
but non-negligible difference between the two atoms. As an example, for the
doubled perturbation, there is a difference in the absolute values of the
first component of the reduced force, 0.010016307568 and -0.010016187598.

Actually, the forces should cancel each other exactly if the translation
symmetry is perfect. This is not the case, but the breaking of this symmetry
can be shown to arise only from the presence of the exchange-correlation
grid of points. This grid does not move when atoms are displaced, and so there
is a very small variation of the total energy when the system is moved as a
whole. It is easy to restore the action-reaction law, by subtracting from
every force component the mean of the forces on all atoms. This is actually
done when the gradient with respect to reduced coordinates are transformed
into forces, and specified in cartesian coordinates, as can be seen in the
output file for the small displacement:

This effect will be seen also at the level of 2DTE. The so-called acoustic sum rule,
which imposes that the frequency of three modes (called acoustic modes)
tends to zero with vanishing wavevector, will also be slightly broken. In this
case, it will also be rather easy to reimpose the acoustic sum rule. In any
case, taking a finer XC grid will allow one to reduce this effect.

We now compute the second derivative of the total energy with respect to the
same atomic displacement through the DFPT capabilities of ABINIT.

You can copy the file $ABI_TUTORESPFN/Input/trf1_3.in in Work_rf1.
This is your input file. You should examine it. The changes with respect to
trf1_1.in are all gathered in the first part of this file, before

Accordingly, you should get familiarized with the new input variables:
rfphon, rfatpol, rfdir. Then, pay attention to the special use of
the kptopt input variable. It will be explained in more detail later.

# Crystalline AlAs : computation of the second derivative of the total energy##Response-function calculation, with q=0
rfphon 1# Will consider phonon-type perturbation
rfatpol 11# Only the first atom is displaced
rfdir 100# Along the first reduced coordinate axis
nqpt 1# One wavevector is to be considered
qpt 000# This wavevector is q=0 (Gamma)
kptopt 2# Automatic generation of k points, taking# into account the time-reversal symmetry only
tolvrs 1.0d-8 # SCF stopping criterion
irdwfk 1# Read the ground-state wavefunctions########################################################################Common input variables#Definition of the unit cell
acell 3*10.61 # This is equivalent to 10.61 10.61 10.61
rprim 0.0 0.5 0.5 # In tutorials 1 and 2, these primitive vectors 0.5 0.0 0.5 # (to be scaled by acell) were 1 0 0 0 1 0 0 0 1 0.5 0.5 0.0 # that is, the default.#Definition of the atom types
ntypat 2# There are two types of atom
znucl 1333# The keyword "znucl" refers to the atomic number of the # possible type(s) of atom. The pseudopotential(s) # mentioned in the "files" file must correspond# to the type(s) of atom. Here, type 1 is the Aluminum,# type 2 is the Arsenic.#Definition of the atoms
natom 2# There are two atoms
typat 12# The first is of type 1 (Al), the second is of type 2 (As).
xred # This keyword indicate that the location of the atoms# will follow, one triplet of number for each atom0.0 0.0 0.0 # Triplet giving the REDUCED coordinate of atom 1.0.25 0.25 0.25 # Triplet giving the REDUCED coordinate of atom 2.#Gives the number of band, explicitely (do not take the default)
nband 4# For an insulator (if described correctly as an insulator # by DFT), there is no need to include conduction bands # in response-function calculations#Exchange-correlation functional
ixc 1# LDA Teter Pade parametrization#Definition of the planewave basis set
ecut 3.0 # Maximal kinetic energy cut-off, in Hartree#Definition of the k-point grid
kptrlatt -4 44# In cartesian coordinates, this grid is simple cubic, and4 -4 4# actually corresponds to the so-called 8x8x8 Monkhorst-Pack grid.44 -4 # It might as well be obtained through the use of# ngkpt, nshiftk and shiftk .#Definition of the SCF procedure
nstep 15# Maximal number of SCF cycles
diemac 9.0 # Although this is not mandatory, it is worth to# precondition the SCF cycle. The model dielectric# function used as the standard preconditioner# is described in the "dielng" input variable section.# The dielectric constant of AlAs is smaller that the one of Si (=12).## After modifying the following section, one might need to regenerate the pickle database with runtests.py -r#%%<BEGIN TEST_INFO>#%% [setup]#%% executable = abinit#%% test_chain = trf1_1.in, trf1_2.in, trf1_3.in, trf1_4.in#%% input_prefix = trf1_1o#%% [files]#%% files_to_test = #%% trf1_3.out, tolnlines= 2, tolabs= 1.100e-08, tolrel= 2.000e-4#%% psp_files = 13al.981214.fhi, 33as.pspnc#%% [paral_info]#%% max_nprocs = 2#%% [extra_info]#%% authors = X. Gonze#%% keywords = NC, DFPT#%% description = Crystalline AlAs : computation of the second derivative of the total energy#%%<END TEST_INFO>

When you have understood the purpose of the input variable values specified
before the “Common input variables” section, you can make the code run, as usual.

Then, we need to analyze the different output files. For that purpose, you should read
the content of the section 6 of the respfn help file.
Read it quickly, as we will come back to the most important points hereafter.

ABINIT has created several different files:

trf1_3.log (the log file)

trf1_3.out (the output file), possibly also trf1_3o_OUT.nc, an abridged netCDF version

trf1_3o_1WF1 (the 1st-order wavefunction file)

trf1_3o_DEN1 (the 1st-order density file)

trf1_3o_POT1 (the 1st-order potential file)

trf1_3o_DDB (the derivative database), possibly also trf1_3o_DDB.nc, its netCDF version

Let us have a look at the output file. You can follow the description provided
in the section 6.2 of the respfn help file.

From these data, you can see that the 2DTE determined by the DFPT technique is
in excellent agreement with the higher-order finite-difference values for the
2DTE, determined in the previous section: 5.007855 Hartree from the energy
differences, and 5.007852 Hartree from the force differences.

Now, you can read the remaining of the section 6.2 of the respfn help file.
Then, you should also open the trf1_3o_DDB file, and read the
corresponding section 6.5 of the respfn help file.

Finally, the excellent agreement between the finite-difference formula and the
DFPT approach calls for some accuracy considerations. These can be found in
section 7 of the respfn help file.

Tip

With AbiPy , one can easily visualize the convergence of the DFPT cycle with the abiopen script
and the syntax:

We are now in the position to compute the full dynamical matrix at the \Gamma point (q=0).
You can copy the file $ABI_TUTORESPFN/Input/trf1_4.in in Work_rf1.
This is your input file.

As for test rf1_3, the changes with respect to trf1_1.in are
all gathered in the first part of this file. Moreover, the changes with
respect to trf1_3.in concern only the input variables rfatpol, and rfdir.
Namely, all the atoms will be displaced, in all the directions.

# Crystalline AlAs : computation of the dynamical matrix at Gamma##Response-function calculation, with q=0
rfphon 1# Will consider phonon-type perturbation
rfatpol 12# All the atoms will be displaced
rfdir 111# Along all reduced coordinate axis
nqpt 1# One wavevector is to be considered
qpt 000# This wavevector is q=0 (Gamma)
kptopt 2# Automatic generation of k points, taking# into account the time-reversal symmetry only
tolvrs 1.0d-8 # SCF stopping criterion
irdwfk 1# Read the ground-state wavefunctions########################################################################Common input variables#Definition of the unit cell
acell 3*10.61 # This is equivalent to 10.61 10.61 10.61
rprim 0.0 0.5 0.5 # In tutorials 1 and 2, these primitive vectors 0.5 0.0 0.5 # (to be scaled by acell) were 1 0 0 0 1 0 0 0 1 0.5 0.5 0.0 # that is, the default.#Definition of the atom types
ntypat 2# There are two types of atom
znucl 1333# The keyword "znucl" refers to the atomic number of the # possible type(s) of atom. The pseudopotential(s) # mentioned in the "files" file must correspond# to the type(s) of atom. Here, type 1 is the Aluminum,# type 2 is the Arsenic.#Definition of the atoms
natom 2# There are two atoms
typat 12# The first is of type 1 (Al), the second is of type 2 (As).
xred # This keyword indicate that the location of the atoms# will follow, one triplet of number for each atom0.0 0.0 0.0 # Triplet giving the REDUCED coordinate of atom 1.0.25 0.25 0.25 # Triplet giving the REDUCED coordinate of atom 2.#Gives the number of band, explicitely (do not take the default)
nband 4# For an insulator (if described correctly as an insulator # by DFT), there is no need to include conduction bands # in response-function calculations#Exchange-correlation functional
ixc 1# LDA Teter Pade parametrization#Definition of the planewave basis set
ecut 3.0 # Maximal kinetic energy cut-off, in Hartree#Definition of the k-point grid
kptrlatt -4 44# In cartesian coordinates, this grid is simple cubic, and4 -4 4# actually corresponds to the so-called 8x8x8 Monkhorst-Pack grid.44 -4 # It might as well be obtained through the use of# ngkpt, nshiftk and shiftk .#Definition of the SCF procedure
nstep 15# Maximal number of SCF cycles
diemac 9.0 # Although this is not mandatory, it is worth to# precondition the SCF cycle. The model dielectric# function used as the standard preconditioner# is described in the "dielng" input variable section.# The dielectric constant of AlAs is smaller that the one of Si (=12).## After modifying the following section, one might need to regenerate the pickle database with runtests.py -r#%%<BEGIN TEST_INFO>#%% [setup]#%% executable = abinit#%% test_chain = trf1_1.in, trf1_2.in, trf1_3.in, trf1_4.in#%% input_prefix = trf1_1o#%% [files]#%% files_to_test = #%% trf1_4.out, tolnlines= 1, tolabs= 4.0e-07, tolrel= 5.000e-08, fld_options = -medium#%% psp_files = 13al.981214.fhi, 33as.pspnc#%% [paral_info]#%% max_nprocs = 2#%% [extra_info]#%% authors = X. Gonze#%% keywords = NC, DFPT#%% description = Crystalline AlAs : computation of the dynamical matrix at Gamma#%%<END TEST_INFO>

There are six perturbations to consider.
So, one might think that the CPU time will raise accordingly.
This is not true, as ABINIT is able to determine which perturbations are the symmetric of another perturbation,
see section 3 of the respfn help file.

Now, you can make the run. You open the file trf1_4.out, and notice that the
response to two perturbations were computed explicitly, while the response to
the other four could be deduced from the two first by using the symmetries.

You might wonder about the dash sign present in the first column of the two
lines giving the frequencies in cm^{-1}. The first column of the main ABINIT
output files is always dedicated to signs needed to automatically treat the
comparison with respect to reference files. Except if you become an ABINIT
developer, you should ignore these signs. In the present case, they should not
be interpreted as a minus sign for the floating numbers that follow them…

There are a good and a bad news about this result. The good news is that
there are indeed three acoustic modes, with frequency rather close to zero
(less than 1 cm^{-1}, which is rather good!). The bad news comes when the three
other frequencies are compared with experimental results, or other theoretical
results. Indeed, in the present run, one obtains three degenerate modes, while
there should be a (2+1) splitting. This can be seen in the paper Ab initio
calculation of phonon dispersions in semiconductors [Giannozzi1991], especially Fig. 2.

Actually, we have forgotten to take into account the coupling between atomic
displacements and the homogeneous electric field, that exists in the case of
polar insulators, for so-called “Longitudinal Optic (LO) modes”. A splitting
appears between these modes and the “Transverse Optic (TO) modes”. This
splitting (Lyddane-Sachs-Teller LO-TO splitting) is presented in simple terms
in standard textbooks, and should not be forgotten in doing Ab initio
calculations of phonon frequencies.

Thus we have now to treat correctly the homogeneous electric field type perturbation.

The treatment of the homogeneous electric field perturbation is formally much
more complex than the treatment of atomic displacements.
This is primarily because the change of potential associated with an homogeneous electric field
is not periodic, and thus does not satisfy the Born-von Karman periodic boundary conditions.

For the purpose of the present tutorial, one should read the section II.C of
the above-mentioned paper [Giannozzi1991].
The reader will find in [Gonze1997] and [Gonze1997a]
more detailed information about this perturbation, closely related to the ABINIT implementation.
There is also an extensive discussion of the Born
effective charges by [Ghosez1998].

In order to compute the response of solids to an homogeneous electric field in ABINIT,
the remaining sections of the respfn help file
should be read. These sections also present the information needed to compute
phonons with non-zero q wavevector, which will be the subject of the next
section of the present tutorial. The sections to be read are:

You are now in the position to compute the full dynamical matrix at Gamma (q=0),
including the coupling with an homogeneous electric field.
You can copy $ABI_TUTORESPFN/Input/trf1_5.in in Work_rf1.
This is your input file.

# Crystalline AlAs : computation of the response to homogeneous# electric field and atomic displacements, at q=0#
ndtset 3#Ground state calculation
kptopt1 1# Automatic generation of k points, taking# into account the symmetry
tolvrs1 1.0d-18 # SCF stopping criterion#Response Function calculation : d/dk
rfelfd2 2# Activate the calculation of the d/dk perturbation
rfdir2 100# Need to consider the perturbation in the x-direction only# This is rather specific, due to the high symmetry of the AlAs crystal# In general, just use rfdir 1 1 1# In the present version of ABINIT (v4.6), symmetry cannot be used# to reduce the number of ddk perturbations
nqpt2 1
qpt2 0.0 0.0 0.0 # This is a calculation at the Gamma point
getwfk2 -1 # Uses as input the output wf of the previous dataset
kptopt2 2# Automatic generation of k points,# using only the time-reversal symmetry to decrease# the size of the k point set.
iscf2 -3 # The d/dk perturbation must be treated # in a non-self-consistent way
tolwfr2 1.0d-22 # Must use tolwfr for non-self-consistent calculations# Here, the value of tolwfr is very low.#Response Function calculation : electric field perturbation and phonons
rfphon3 1# Activate the calculation of the atomic dispacement perturbations
rfatpol3 12# All the atoms will be displaced
rfelfd3 3# Activate the calculation of the electric field perturbation
rfdir3 111# All directions are selected. However, symmetries will be used to decrease# the number of perturbations, so only the x electric field is needed# (and this explains why in the second dataset, rfdir was set to 1 0 0).
nqpt3 1
qpt3 0.0 0.0 0.0 # This is a calculation at the Gamma point
getwfk3 -2 # Uses as input wfs the output wfs of the dataset 1
getddk3 -1 # Uses as input ddk wfs the output of the dataset 2
kptopt3 2# Automatic generation of k points,# using only the time-reversal symmetry to decrease# the size of the k point set.
tolvrs3 1.0d-8
########################################################################Common input variables#Definition of the unit cell
acell 3*10.61 # This is equivalent to 10.61 10.61 10.61
rprim 0.0 0.5 0.5 # In tutorials 1 and 2, these primitive vectors 0.5 0.0 0.5 # (to be scaled by acell) were 1 0 0 0 1 0 0 0 1 0.5 0.5 0.0 # that is, the default.#Definition of the atom types
ntypat 2# There are two types of atom
znucl 1333# The keyword "znucl" refers to the atomic number of the # possible type(s) of atom. The pseudopotential(s) # mentioned in the "files" file must correspond# to the type(s) of atom. Here, type 1 is the Aluminum,# type 2 is the Arsenic.#Definition of the atoms
natom 2# There are two atoms
typat 12# The first is of type 1 (Al), the second is of type 2 (As).
xred # This keyword indicate that the location of the atoms# will follow, one triplet of number for each atom0.0 0.0 0.0 # Triplet giving the REDUCED coordinate of atom 1.0.25 0.25 0.25 # Triplet giving the REDUCED coordinate of atom 2.#Gives the number of band, explicitely (do not take the default)
nband 4# For an insulator (if described correctly as an insulator # by DFT), there is no need to include conduction bands # in response-function calculations#Exchange-correlation functional
ixc 1# LDA Teter Pade parametrization#Definition of the planewave basis set
ecut 3.0 # Maximal kinetic energy cut-off, in Hartree#Definition of the k-point grid
kptrlatt -4 44# In cartesian coordinates, this grid is simple cubic, and4 -4 4# actually corresponds to the so-called 8x8x8 Monkhorst-Pack grid.44 -4 # It might as well be obtained through the use of# ngkpt, nshiftk and shiftk .#Definition of the SCF procedure
nstep 15# Maximal number of SCF cycles
diemac 9.0 # Although this is not mandatory, it is worth to# precondition the SCF cycle. The model dielectric# function used as the standard preconditioner# is described in the "dielng" input variable section.# The dielectric constant of AlAs is smaller that the one of Si (=12).## 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 = #%% trf1_5.out, tolnlines= 1, tolabs= 6.000e-07, tolrel= 5.000e-08, fld_options = -medium#%% psp_files = 13al.981214.fhi, 33as.pspnc#%% [paral_info]#%% max_nprocs = 2#%% [extra_info]#%% authors = X. Gonze#%% keywords = NC, DFPT#%% description = #%% Crystalline AlAs : computation of the response to homogeneous#%% electric field and atomic displacements, at q=0#%%<END TEST_INFO>

As for the other DFPT tests, the changes with respect to the trf1_1.in are all gathered
in the first part of this file.
Unlike the other tests, however, the multi-dataset mode was used, computing from scratch
the ground-state properties, then computing the effect of the ddk perturbation, then the effect of all
other perturbations (electric field as well as atomic displacements).

The analysis of the output file is even more cumbersome than the previous ones.
Let us skip the first dataset. In the dataset 2 section, one perturbation is correctly selected:

The analysis of the output for this particular perturbation is not
particularly interesting, except for the f-sum rule ratio

loper3 : ek2= 1.6833336546E+01
f-sum rule ratio= 1.0028582975E+00

that should be close to 1, and becomes closer to it when ecut is
increased, and the sampling of k points is improved. (In the present status of
ABINIT, the f-rule ratio is not computed correctly when ecutsm/=0)

In the third dataset section, three irreducible perturbations are considered:

It is diagonal and isotropic, and corresponds to a dielectric constant of 9.7606052146.

Then, the Born effective charges are given, either computed from the
derivative of the wavefunctions with respect to the electric field, or
computed from the derivative of the wavefunctions with respect to an atomic
displacement, as explained in section II of [Gonze1997a]:

Namely, the Born effective charge of the Al atom is 2.104, and the one of the
As atom is -2.127 . The charge neutrality sum rule is not fulfilled exactly.
When ecut is increased, and the sampling of k points is improved, the sum
of the two charges goes closer to zero.

The first few lines discard any effect of the homogeneous electric field,
while the next sections consider it along the three cartesian coordinates.

In the present material, the directionality of the electric field has no
influence. We note that there are still three acoustic mode, below 1 cm^{-1},
while the optic modes have the correct degeneracies: two TO modes at 344.3
cm^{-1}, and one LO mode at 379.6 cm^{-1}.

These values can be compared to experimental (361 cm^{-1}, 402 cm^{-1}) as well
as theoretical (363 cm^{-1}, 400 cm^{-1}) values (again [Giannozzi1991]).
Most of the discrepancy comes from the too low value of
ecut. Using ABINIT with ecut=6 Hartree gives (358.8 cm^{-1}, 389.8 cm^{-1}).
The remaining of the discrepancy may come partly from the
pseudopotentials, that are particularly soft.

The comparison of Born effective charges is also interesting. After imposition
of the neutrality sum rule, the Al Born effective charge is 2.116. The value
from Gianozzi et al is 2.17, the experimental value is 2.18.
Increasing ecut to 6 Hartree in ABINIT gives 2.168.

For the dielectric tensor, it is more delicate. The value from Gianozzi et al
is 9.2, while the experimental value is 8.2 . The agreement is not very good,
a fact that can be attributed to the LDA lack of polarization-dependence [Gonze1995a].
Still, the agreement of our calculation with the theoretical result is not very good. With
ecut = 3 Hartree, we have 9.76. Changing it to 6 Hartree gives 10.40 . A
better k point sampling (8x8x8), with ecut = 6 Hartree, reduces the value to 9.89.
Changing pseudopotentials finally improves the agreement: with the
much harder 13al.pspgth and 33as.psphgh pseudopotentials with adequate
ecut = 16 Hartree and 8x8x8 Monkhorst-Pack sampling, we reach a value of 9.37.
This illustrates that the dielectric tensor is a much more sensitive quantity than the others.

The computation of phonon frequencies at non-zero q is actually simpler than the one at Gamma.
One must distinguish two cases. Either the q wavevector connects k points that belong
to the same grid, or the wavevector q is general.
In any case, the computation within the reciprocal space DFPT formalism is more
efficient than the real space frozen-phonon technique since the use of supercells is
completely avoided with DFPT. For an explanation of this fact, see for example section IV of [Gonze1997].

You can copy the file $ABI_TUTORESPFN/Input/trf1_6.in in Work_rf1.
This is your input file.

# Crystalline AlAs : computation of the response to homogeneous# electric field and atomic displacements, at q=X, q=L and # an acoustic mode close to q=Gamma
ndtset 5#Ground state calculation
kptopt1 1# Automatic generation of k points, taking# into account the symmetry
tolvrs1 1.0d-18 # SCF stopping criterion
prtden1 1# Will be needed for dataset 4#Response Function calculation : phonon at X
rfphon2 1# Activate the calculation of the atomic dispacement perturbations
rfatpol2 12# All the atoms will be displaced
rfdir2 111# Need to consider the perturbation in the x-direction only
nqpt2 1
qpt2 0.0 0.5 0.5 # This is a calculation at the X point
getwfk2 1# Uses as input wfs the output wfs of the dataset 1
kptopt2 3# Automatic generation of k points,# no use of symmetries to decrease # the size of the k point set.
tolvrs2 1.0d-8
#Response Function calculation : phonon at L
rfphon3 1# Activate the calculation of the atomic dispacement perturbations
rfatpol3 12# All the atoms will be displaced
rfdir3 111# Need to consider the perturbation in the x-direction only
nqpt3 1
qpt3 0.5 0.0 0.0 # This is a calculation at the L point
getwfk3 1# Uses as input wfs the output wfs of the dataset 1
kptopt3 3# Automatic generation of k points,# no use of symmetries to decrease # the size of the k point set.
tolvrs3 1.0d-8
#Non-self consistent ground-state calculation, with q=(0.1 0 0)
nqpt4 1
qpt4 0.1 0.0 0.0 # This is a calculation along the Gamma-L line
getwfk4 1# Uses as input wfs the output wfs of the dataset 1
getden4 1# Uses as input density the output density of the dataset 1
kptopt4 3# Automatic generation of k points,# no use of symmetries to decrease # the size of the k point set.
tolwfr4 1.0d-18
iscf4 -2 # Non-self-consistent calculation#Response Function calculation : phonon close to Gamma
rfphon5 1# Activate the calculation of the atomic dispacement perturbations
rfatpol5 12# All the atoms will be displaced
rfdir5 111# Need to consider the perturbation in the x-direction only
nqpt5 1
qpt5 0.1 0.0 0.0 # This is a calculation along the Gamma-L line
getwfk5 1# Uses as input wfs the output wfs of the dataset 1
getwfq5 4# Uses as input ddk wfs the output of the dataset 4
kptopt5 3# Automatic generation of k points,# no use of symmetries to decrease # the size of the k point set.
tolvrs5 1.0d-8
########################################################################Common input variables#Definition of the unit cell
acell 3*10.61 # This is equivalent to 10.61 10.61 10.61
rprim 0.0 0.5 0.5 # In tutorials 1 and 2, these primitive vectors 0.5 0.0 0.5 # (to be scaled by acell) were 1 0 0 0 1 0 0 0 1 0.5 0.5 0.0 # that is, the default.#Definition of the atom types
ntypat 2# There are two types of atom
znucl 1333# The keyword "znucl" refers to the atomic number of the # possible type(s) of atom. The pseudopotential(s) # mentioned in the "files" file must correspond# to the type(s) of atom. Here, type 1 is the Aluminum,# type 2 is the Arsenic.#Definition of the atoms
natom 2# There are two atoms
typat 12# The first is of type 1 (Al), the second is of type 2 (As).
xred # This keyword indicate that the location of the atoms# will follow, one triplet of number for each atom0.0 0.0 0.0 # Triplet giving the REDUCED coordinate of atom 1.0.25 0.25 0.25 # Triplet giving the REDUCED coordinate of atom 2.#Gives the number of band, explicitely (do not take the default)
nband 4# For an insulator (if described correctly as an insulator # by DFT), there is no need to include conduction bands # in response-function calculations#Exchange-correlation functional
ixc 1# LDA Teter Pade parametrization#Definition of the planewave basis set
ecut 3.0 # Maximal kinetic energy cut-off, in Hartree#Definition of the k-point grid
kptrlatt -4 44# In cartesian coordinates, this grid is simple cubic, and4 -4 4# actually corresponds to the so-called 8x8x8 Monkhorst-Pack grid.44 -4 # It might as well be obtained through the use of# ngkpt, nshiftk and shiftk .#Definition of the SCF procedure
nstep 15# Maximal number of SCF cycles
diemac 9.0 # Although this is not mandatory, it is worth to# precondition the SCF cycle. The model dielectric# function used as the standard preconditioner# is described in the "dielng" input variable section.# The dielectric constant of AlAs is smalle that the one of Si (=12).## 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 = #%% trf1_6.out, tolnlines= 11, tolabs= 1.1e-03, tolrel= 2.000e-04, fld_options = -medium#%% psp_files = 13al.981214.fhi, 33as.pspnc#%% [paral_info]#%% max_nprocs = 2#%% [extra_info]#%% authors = X. Gonze#%% keywords = NC, DFPT#%% description = #%% Crystalline AlAs : computation of the response to homogeneous#%% electric field and atomic displacements, at q=X, q=L and #%% an acoustic mode close to q=Gamma#%%<END TEST_INFO>

As for the other RF tests, the changes with respect to trf1_1.in are
all gathered in the first part of this file.
The multi-dataset mode is used, computing from scratch the ground-state wave functions, then computing different
dynamical matrices with DFPT.
The run is about 2 minutes on a 2.8 GHz machine.
So, you would better leave your computer running, and either read more of the ABINIT documentation
(why not the mrgddb_help and the anaddb_help), or make a walk).

The results of this simulation can be compared to those provided in the Gianozzi et al paper.
The agreement is rather good, despite the low cut-off energy, and different pseudopotentials.

At X, they get 95 cm^{-1}, 216 cm^{-1}, 337 cm^{-1} and 393 cm^{-1}, while we get
92.5 cm^{-1}, 204.6 cm^{-1}, 313.9 cm^{-1} and 375.9 cm^{-1}.
With ecut=6 Hartree, we get 89.7 cm^{-1}, 212.3 cm^{-1}, 328.5 cm^{-1} and 385.8 cm^{-1}.

At L, they get 71 cm^{-1}, 212 cm^{-1}, 352 cm^{-1} and 372 cm^{-1}, while we get
69.0 cm^{-1}, 202.5 cm^{-1}, 332.6 cm^{-1} and 352.3 cm^{-1}. With ecut=6 Hartree, we
get 68.1 cm^{-1}, 208.5 cm^{-1}, 346.7 cm^{-1} and 362.6 cm^{-1}.

At q=(0.1 0 0), we get 31.6 cm^{-1}, 63.6 cm^{-1}, 342.0 cm^{-1} and 379.7 cm^{-1}. The
acoustic modes tends (nearly-)linearly to zero, while the optic modes are
close to their values at Gamma : 344.3 cm^{-1} and 379.6 cm^{-1}.

Note

This ABINIT tutorial is now finished.
You are advised to go through the second tutorial on DFPT to
make some post-processing analysis (phonon dispersions, thermodynamical properties, etc)