This lesson aims at showing how to get the following physical properties, for an insulator :
This lesson should take about (to be provided) hours to be done.
Before beginning, you might consider to work in a different subdirectory as
for the other lessons. Why not "Work5" ?
The file ~ABINIT/Tutorial/t5x.files lists the file names and root names. You can copy it in the Work5 directory (and change it, as usual). Note that two pseudopotentials are mentioned in this "files" file: one for the Aluminum atom, and one for the Arsenic atom. The first to be mentioned, for Al, will define the first type of atom. The second to be mentioned, for As, will define the second type of atom. It is the first time that you encounter this situation in the tutorials.
You can also copy the file ~ABINIT/Tutorial/t51.in in Work5. This is your input file. You should edit it, read it carefully, and, because of the use of two types of atoms, have a look at the following input variables :
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 response-function calculation.You should make the run (a dozen of second on a PIII at 450 MHz), and obtain the following value for the energy, in the final echo section :
etotal -9.7626837450E+00However, we will rely later on a more accurate (more digits) value of this total energy, that can be found about a dozen of lines before this final echo :
Ewald =-8.47989583509473E+00 resulting in Etotal=-9.76268374500102E+00 hartree
The output file also mention 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 response functions will be computed. The main output of this ground-state run is the wavefunction file t51_oWFK, that you can already rename as t51_iWFK. Indeed, it will be used in the next runs, as an input file. So, in the corresponding "files" file, third line, pay attention to specify "t51_i", even if you change the root name for output files (fourth line) to "t52_o" or "t53_o" ...
We will now aim at computing 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 response function 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 now, in order to be able to perform a direct comparison with the result of a response-function 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 ~ABINIT/Tutorial/t52.in in Work5. This is your input file. You should edit it and briefly look at the two changes with respect to the file ~ABINIT/Tutorial/t51.in : the change of xred, and the reading of the wavefunction file, using the irdwfk input variable. Then, you can make the run. 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 (about one minute on a PIII 450 MHz).
From this run, it is possible to get the values of the total energy, and the value of the gradient of the total energy with respect to change of reduced coordinate:
rms dE/dt= 3.5517E-03; max dE/dt= 5.0079E-03; dE/dt below (all hartree)
1 0.005007930232 0.002526304574 0.002526304574
2 -0.005007868293 -0.002526274885 -0.002526274885
...
Ewald =-8.47988991313938E+00 resulting in Etotal=-9.76268124105590E+00 hartree
The change of reduced coordinate of the Al atom along the first axis
was rather small (1/1000), 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.76268374500102E+00 hartree) of the previous run, and the perturbed value (-9.76268124105590E+00 hartree) of the present one, is thus one half of the square of the coordinate change (1/1000) times the 2DTE. From these number, the 2DTE is 5.00791024 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.005007930232 (hartree). At first order, this quantity is the product of the 2DTE by the reduced coordinate difference. The estimate of the 2DTE is thus 5.007930232 hartree. The agreement with the other estimate is rather good (2.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 :
rms dE/dt= 7.1249E-03; max dE/dt= 1.0016E-02; dE/dt below (all hartree)
1 0.010016299779 0.005097509981 0.005097509981
2 -0.010016176675 -0.005097455174 -0.005097455174
...
Ewald =-8.47987214716789E+00 resulting in Etotal=-9.76267372899697E+00 hartree
From these data, taking into account that the perturbation
was twice stronger, the same procedure than above leads
to the values 5.00800270 hartree (from finite difference of energy)
and 5.009149889 hartree (from finite difference of forces,
the value 0.010016299779 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.0001104 Ha, which divided by three,
and subtracted from 5.00789230 hartree,
gives 5.0078555 hartree. The same higher-order procedure for
force estimates gives 5.0078552 hartree. So, the agreement
between total-energy estimate and force estimate of the 2DTE can be
observed up to the 7th digit, inclusive.
Before comparing this result with the 2DTE directly computed from the response-function 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, for example for the doubled displacement :
rms dE/dt= 7.1249E-03; max dE/dt= 1.0016E-02; dE/dt below (all hartree)
1 0.010016299779 0.005097509981 0.005097509981
2 -0.010016176675 -0.005097455174 -0.005097455174
show a small, but non-negligible difference between the
two atoms. Actually, the forces should cancel each other exactly
if the translation symmetry were 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 :
cartesian forces (hartree/bohr) at end:
1 -0.00000421123276 -0.00047199792687 -0.00047199792687
2 0.00000421123276 0.00047199792687 0.00047199792687
This effect will be seen also at the level of 2DTE. The so-called
"acoustic sum rule", imposing that the frequency
of three modes (called acoustic modes) tend
to zero with vanishing wavevector, will also be slightly broken.
In this case also, it will be rather easy to reimpose the
acoustic sum rule. In any case, taking a finer XC grid will
allow to reduce this effect.
We now compute the second derivative of the total energy with respect to the same atomic displacement, using the response-function capabilities of ABINIT.
You can copy the file ~ABINIT/Tutorial/t53.in in Work5. This is your input file. You should edit it. The changes with respect to the file ~ABINIT/Tutorial/t51.in are all gathered in the first part of this file, before
####################################################################### #Common input variables
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.
When you have understood the purpose of the input variable values specified before the "Common input variables" section, you can make the code run. It takes less than one minute on a PIII 450MHz.
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 four different files :
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. You should be able to find the place where the iterations for the minimisation (with respect to the unique perturbation) take place :
iter 2DEtotal(Ha) deltaE(Ha) residm vres2
ETOT 1 6.5156051312863 -1.464E+01 1.146E-02 1.947E+02
ETOT 2 5.0216331638978 -1.494E+00 9.267E-04 2.027E+00
ETOT 3 5.0082678671217 -1.337E-02 4.772E-06 7.929E-02
ETOT 4 5.0078677958105 -4.001E-04 1.980E-07 2.712E-03
ETOT 5 5.0078558860285 -1.191E-05 6.103E-09 5.074E-05
ETOT 6 5.0078557520180 -1.340E-07 1.258E-10 9.613E-06
ETOT 7 5.0078557017091 -5.031E-08 2.768E-11 3.841E-07
ETOT 8 5.0078557001188 -1.590E-09 2.983E-12 6.089E-09
From these data, you can see that the 2DTE determined by the
response-function technique is in excellent agreement with the
higher-order finite-difference formula for the 2DTE, determined
in the previous section :
5.0078555 hartree from the energy differences,
and 5.0078552 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 edit the t53o_DDB file, and read the corresponding section 6.4 of the respfn_help file.
Finally, the excellent agreement between the finite-difference formula and the response-function approach calls for some accuracy considerations. These can be found in section 7 of the respfn_help file.
We are now in the position to compute the full dynamical matrix at Gamma (q=0). You can copy the file ~ABINIT/Tutorial/t54.in in Work5. This is your input file. You should edit it. As for test 53, the changes with respect to the file ~ABINIT/Tutorial/t51.in are all gathered in the first part of this file. Moreover, the changes with respect to t53.in concern only the input variables rfatpol, and rfdir. Namely, all the atoms will be displaced, in all the directions.
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 section 3 of the respfn_help file.
Now, you can make the run. It is a bit longer than one minute on a PIII at 450MHz. You edit the file t54.out, and notice that the response to two perturbations were computed explicitely, while the response to the other four could be deduced by using the symmetries. Nothing mysterious : one of the two irreducible perturbations is for the Al atom, placed in a rather symmetric local site, and the other perturbation is for the As atom.
The phonon frequencies, obtained by diagonalizing the dynamical matrix (where the atomic masses have been taken into account, see amu), are given as follows :
Phonon wavevector (reduced coordinates) : 0.00000 0.00000 0.00000 Phonon energies in Hartree : 2.586632E-06 2.590723E-06 2.614440E-06 1.568560E-03 1.568560E-03 1.568560E-03 Phonon frequencies in cm-1 : - 5.677000E-01 5.685980E-01 5.738033E-01 3.442590E+02 3.442590E+02 - 3.442590E+02
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 automatic treat the comparison with respect to reference files. Except if you become a 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 is a good news about this result, and a bad news. The good news is that there are indeed three acoustic mode, 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, by P. Giannozzi, S. de Gironcoli, P. Pavone and S. Baroni, Phys. Rev. B 43, 7231 (1991) , 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 P. Giannozzi, S. de Gironcoli, P. Pavone and S. Baroni, Phys. Rev. B 43, 7231 (1991) . The reader will find in X. Gonze, Phys. Rev. B 55, 10337 (1997) and X. Gonze and C. Lee, Phys. Rev. B 55, 10355 (1997) more detailed information of this perturbation, closely related to the ABINIT implementation. There is also an extensive discussion of the Born effective charges by Ph. Ghosez, J.-P. Michenaud and X. Gonze, Phys. Rev. B 58, 6224 (1998).
In order to compute the response of solids to an homogeneous electric field, as implemented in ABINIT, the remaining sections of the respfn_help file should be read. These 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 the file ~ABINIT/Tutorial/t55.in in Work5. This is your input file. You should edit it. As for the other RF tests, the changes with respect to the file ~ABINIT/Tutorial/t51.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 run lasts about 3 minutes on a PIII 450MHz.
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 :
==> initialize data related to q vector <==
The list of irreducible perturbations for this q vector is:
1) idir= 1 ipert= 3
================================================================================
--------------------------------------------------------------------------------
Perturbation wavevector (in red.coord.) 0.000000 0.000000 0.000000
Perturbation : derivative vs k along direction 1
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.0028582985E+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 :
==> initialize data related to q vector <==
The list of irreducible perturbations for this q vector is:
1) idir= 1 ipert= 1
2) idir= 1 ipert= 2
3) idir= 1 ipert= 4
Much later, the dielectric tensor is given :
Dielectric tensor, in cartesian coordinates,
j1 j2 matrix element
dir pert dir pert real part imaginary part
1 4 1 4 9.7606048428 0.0000000000
1 4 2 4 0.0000000000 0.0000000000
1 4 3 4 0.0000000000 0.0000000000
2 4 1 4 0.0000000000 0.0000000000
2 4 2 4 9.7606048428 0.0000000000
2 4 3 4 0.0000000000 0.0000000000
3 4 1 4 0.0000000000 0.0000000000
3 4 2 4 0.0000000000 0.0000000000
3 4 3 4 9.7606048428 0.0000000000
It is diagonal and isotropic, and corresponds to a dielectric
constant of 9.7606048428 .
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 X. Gonze, Phys. Rev. B55, 10355 (1997) :
Effective charges, in cartesian coordinates, (from electric field response) ...and
Effective charges, in cartesian coordinates, (from phonon response) ...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.
Finally, the phonon frequencies are computed:
Phonon wavevector (reduced coordinates) : 0.00000 0.00000 0.00000 Phonon energies in Hartree : 2.586632E-06 2.590723E-06 2.614440E-06 1.568560E-03 1.568560E-03 1.568560E-03 Phonon frequencies in cm-1 : - 5.677000E-01 5.685980E-01 5.738033E-01 3.442590E+02 3.442590E+02 - 3.442590E+02 Phonon at Gamma, with non-analyticity in the direction (cartesian coordinates) 1.00000 0.00000 0.00000 Phonon energies in Hartree : 2.590670E-06 2.590723E-06 4.101029E-06 1.568560E-03 1.568560E-03 1.729575E-03 Phonon frequencies in cm-1 : - 5.685864E-01 5.685980E-01 9.000719E-01 3.442590E+02 3.442590E+02 - 3.795979E+02 Phonon at Gamma, with non-analyticity in the direction (cartesian coordinates) 0.00000 1.00000 0.00000 Phonon energies in Hartree : 2.586632E-06 2.614440E-06 4.088526E-06 1.568560E-03 1.568560E-03 1.729575E-03 Phonon frequencies in cm-1 : - 5.677000E-01 5.738033E-01 8.973277E-01 3.442590E+02 3.442590E+02 - 3.795979E+02 Phonon at Gamma, with non-analyticity in the direction (cartesian coordinates) 0.00000 0.00000 1.00000 Phonon energies in Hartree : 2.590723E-06 2.610339E-06 4.088540E-06 1.568560E-03 1.568560E-03 1.729575E-03 Phonon frequencies in cm-1 : - 5.685980E-01 5.729031E-01 8.973308E-01 3.442590E+02 3.442590E+02 - 3.795979E+02The 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 1cm^-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 the Gianozzi et al paper mentioned above). 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 ( X. Gonze, Ph. Ghosez and R. Godby, Phys. Rev. Lett. (1995)). 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 response-function formalism is more efficient than using the frozen-phonon technique: the use of supercell is completely avoided. For an explanation of this fact, see for example section IV of X. Gonze, Phys. Rev. B55, 10337 (1997).
You can copy the file ~ABINIT/Tutorial/t56.in in Work5. This is your input file. You should edit it. As for the other RF tests, the changes with respect to the file ~ABINIT/Tutorial/t51.in are all gathered in the first part of this file. The multi-dataset mode is used, computing from scratch the ground-state properties, then computing different dynamical matrices. The run is rather long : about 30 minutes on a PIII 450MHz. 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 95cm^-1, 216cm^-1, 337cm^-1 and 393cm^-1, while we get 92.5cm^-1, 204.6cm^-1, 313.9cm^-1 and 375.9cm^-1. With ecut=6 hartree, we get 89.7cm^-1, 212.3cm^-1, 328.5cm^-1 and 385.8cm^-1.
At L, they get 71cm^-1, 212cm^-1, 352cm^-1 and 372cm^-1, while we get 69.0cm^-1, 202.5cm^-1, 332.6cm^-1 and 352.3cm^-1. With ecut=6 hartree, we get 68.1cm^-1, 208.5cm^-1, 346.7cm^-1 and 362.6cm^-1.
At q=(0.1 0 0), we get 31.6cm^-1, 63.6cm^-1, 342.0cm^-1 and 379.7cm^-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.
This section is still to be written. You might have a look at the tests 26 to 32 of the directory ~ABINIT/Infos/Test_v2 .
The ABINIT tutorial is now finished. It might be worth to read the full list of abinis input variables. Then, proceeds to the ~ABINIT/Infos/Dirs_and_files file, to have a global view of ABINIT.