This lesson aims at showing how to get the following linear and non-linear physical properties, for an insulator :
The case of the linear response is treated independently in other tutorials. However, you will learn here how to get these quantities using either the response-function or the finite difference techniques within ABINIT. To that end, we will describe how to perform Berry phase and finite electric field calculations.
This lesson should take about 2 hours to be done.
Before beginning, you might consider to work in a different subdirectory as for the other lessons. Why not "Work_NLO" ?
In this tutorial we will assume that the ground-state properties of AlAs have been previously obtained (for more informations in this matter, you can look at Section 1 in the Tutorial concerning the dynamical and dielectric properties of AlAs. In any case, it is recommended to perform Tutorials Response-function 1 and Response-function 2 before starting this one.
All along the present tutorial we will adopt the following set of generic parameters :
acell 10.53 ixc 3 ecut 2.8 (results with ecut = 5 are also reported in the discussion) ecutsm 0.5 dilatmx 1.05 nband 4 (=number of occupied bands) nbdbuf 0 ngkpt 6 6 6 nshiftk 4 shiftk 0.5 0.5 0.5 0.5 0.0 0.0 0.0 0.5 0.0 0.0 0.0 0.5 pseudopotentials 13al.pspnc 33as.pspnc
In principle, the acell to be used should be the one corresponding to the optimized structure at the ecut, and ngkpt combined with nshiftk and shiftk, chosen for the calculations. Unfortunately, for the purpose of this tutorial, in order to limit the duration of the runs, we have to work at an unusually low cutoff of 2.8 Ha for which the optimized lattice constant is unrealistic and equal to 7.45 bohr (instead of the converged value of 10.64). In what follows, the lattice constant has been arbitrarily fixed to 10.53 bohr. For comparison, results with ecut=5 are also reported and, in that case, were obtained at the optimized lattice constant of 10.64 bohr. For those who would like to try later, convergence tests and structural optimizations can be done using the file tnlo_1.in. Before going further, you might refresh your mind concerning the other variables : ixc, ecutsm, dilatmx, nbdbuf.
Let us emphasize that, from our experience, accurate and trustable results for third order energy derivatives usually require a extremely high level of convergence. So, even if it is not always very apparent here, careful convergence studies must be explicitely performed in any practical study you could start after this tutorial.
As a theoretical support to this section of the tutorial, you might consider reading the following article:
M. Veithen, X. Gonze, and Ph. Ghosez,
Nonlinear optical susceptibilities, Raman efficiencies,
and electro-optic tensors from first-principles density functional perturbation theory
Phys. Rev. B 71, 125107 (2005).
In the first part of this tutorial, we will describe how to compute various linear and non-linear responses directly connected to second-order and third-order derivatives of the energy, using DFPT. From the (2n+1) theorem, computation of energy derivatives up to third order only requires the knowledge of the ground-state and first-order wavefunctions, see X. Gonze and J.-P. Vigneron, PRB 39, 13120 (1989) and X. Gonze, Phys. Rev. A 52, 1096 (1995). Our study will therefore include the following steps : (i) resolution of the ground-state problem, (ii) determination of the first-order wavefunctions and construction of the related databases for second and third-order energy derivatives, (iii) combination of the different databases and analysis to get the physical properties of interest.
This is closely related to what was done for the dielectric and dynamical properties except that an additional database for third-order energy derivatives will now be set up during the run. Only selected third-order derivatives are implemented at this stage and concern responses to electric field and atomic displacements:
Many different steps can be combined in one input file. For the sake of clarity we have decomposed the calculation into individual inputs that are now described.
Responses to electric fields and atomic displacements.
You can copy the file ~abinit/tests/tutorial/Input/tnlo_2.in in Work-NLO. This is your first input file. Its purpose is to build databases for second and third energy derivatives with respect to electric fields and atomic displacements. You can edit it. It is made of 5 datasets. The first four data sets are nearly the same as for a usual linear response calculation : (1) self-consistent calculation in the IBZ; (2) non self-consistent calculations to get the wave-functions over the full BZ; (3) ddk calculation, (4) derivatives with respect to electric field and atomic displacements. Some specific features must however be explicitely specified in order to prepare the non-linear response step (dataset 5). First, from dataset 2 it is mandatory to specify :
nbdbuf 0 nband 4 (= number of valence bands)
Also, in dataset 4, it is required to impose prtden, and prepanl
prtden4 1 prepanl4 1
The purpose for this is (i) to constrain kptopt = 2 even in the computation of phonons where ABINIT usually take advantages of symmetry irrespective of kptopt and (ii) compute the electric field derivatives in the 3 directions of space, irrespective of the symmetry.
You can now make the run. You can have a quick look to the output file to verify that everything is OK. It contains the values of second and third energy derivatives. It has however no immediate interest since the information is not presented in a very convenient format. The relevant information is in fact also stored in the database files (DDB) generated independently for second and third energy derivatives at the end of run steps 4 and 5. Keep these databases that will be used later for a global and convenient analysis of the results using ANADDB.
Responses to strain.
You can now copy the file ~abinit/tests/tutorial/Input/tnlo_3.in in Work-NLO. This is your second input file. Its purpose is to build databases for second energy derivatives with respect to strains. You can edit it. It is made of 4 datasets : (1) self-consistent calculation in the IBZ; (2) non self-consistant calculations to get the wave-functions over the full BZ; (3) ddk calculation, (4) strain perturbation. The ddk calculations has been included in order to access to the piezoelectric tensor.
When you are ready, you can make the run. Again, you can have a quick look to the output tnlo_3.out. It contains rigid ions elastic and piezoelectric constants as well as the internal strain coupling parameters. This information is also stored in a database file (DDB) for further convenient analysis with ANADDB.
Merge of the DDB.
At this stage, all the relevant energy derivatives have been obtained and are stored in individual databases. These must be combined with the MERGE utility in order to get a unique database tnlo_4.ddb.out. Explicitely, you should merge the files tnlo_2o_DS4_DDB, tnlo_3o_DS4_DDB, and tnlo_2o_DS5_DDB . You might have a look at the input file for MRGDDB named tnlo_4.in , and use it to perform the merge.
Analysis of the DDB.
We are now ready for the analysis of the results using ANADDB. You can copy the file ~abinit/tests/tutorial/Input/tnlo_5.in in Work-NLO. You already used ANADDB. The present input is in principle very similar to the one you have used for the analysis of dynamical and dielectric responses except that some new flags need to be activated.
For the strain perturbation you need to specify elaflag, piezoflag, and instrflag :
elaflag 3 piezoflag 3 instrflag 1
For the non-linear responses you need
nlflag 1 ramansr 1 alphon 1 prtmbm 1
nlflag=1 activates the non-linear response.
ramansr=1 will impose the sum rule on the first order change of the electronic dielectric susceptibility under atomic displacement, hereafter referred to as on dchi/dtau . It is a condition of invariance of chi under translation of the whole crystal, similar to the acoustic sum rules for phonons at Gamma or the charge neutrality sum rule on Z*.
prtmbm=1 will allow to get the mode by mode phonon contributions of the ions to the electro-optic coefficients.
alphon=1 will allow to get the previous mode by mode contribution when aligning the eigenvectors with the cartesian axis of coordinates (in the input, the principal axis should always be aligned with z for a convenient analysis of the results).
Finally, the second list of phonon, specified with nph2l and qph2l, must also be explicitely considered to obtain the Raman efficiencies of longitudinal modes (in a way similar to the computation of frequencies of longitudinal mode frequencies at Gamma):
# Wave vector list no. 2 #*********************** nph2l 1 qph2l 1.0 0.0 0.0 0.0
You can now run the code ANADDB. The results are in the file tnlo_5.out. Various interesting physical properties are now directly accessible in this output in meaningful units. You can go through the file and look in order to identify the results mention hereafter.
For comparison, we report in parenthesis (...) the values obtained with ecut = 5, and for nonlinear reponses in brackets [...] the fully converged result as reported in PRB 71, 125107 (2005).
Some physical quantities concern the dielectric and dynamical properties as discussed in Tutorial 5.
Z*_Al = 2.043399E+00 (2.105999E+00)
w_TO (cm^-1) = 3.694366E+02 (3.602635E+02) w_LO (cm^-1) = 4.031189E+02 (3.931598E+02)
Electronic dielectric tensor = 9.20199931 (9.94846084)
relaxed ion dielectric tensor = 10.95641891 (11.84823634)Some other quantities, as the piezoelectric coefficients, are related to the strain response as it is more extensively discussed in the tutorial on the strain perturbation.
clamped ion (Unit:c/m^2) = -0.65029623 (-0.69401363) relaxed ion (Unit:c/m^2) = 0.03754587 (-0.04228777)Finally, different quantities are related to non-linear responses.
d_36 (pm/V) = 21.175513 (32.772254) [fully converged :35]
Electronic EO constant (pm/V): -1.000298328 (-1.324507791) [-1.69] Full Ionic EO constant (pm/V): 0.543837627 (0.533097548) [0.64] Total EO constant (pm/V): -0.456460701 (-0.791410242) [-1.05]
alpha(TO) = -0.008489207 (-0.009114814) alpha(LO) = -0.011466205 (-0.013439375)
dchi_23/dtau_1 (Bohr^-1) of Al = -0.094488223 (-0.099889084)
a = Omega_0 * dchi/dtau = Sqrt[mu * Omega_0] alpha
a(TO) (Unit: Ang^2)= -7.7233 (-8.4222112) [-8.48] a(LO) (Unit: Ang^2)= -10.4317 (-12.418168) [-12.48]
Finite difference calculation of dchi/dtau.
For comparison with the DPFT calculation, we can compute dchi/dtau for Al from finite differences. In practice, this is achieved by computing the linear optical susceptibility for 3 different positions of Al. This is done with the file ~abinit/tests/tutorial/Input/tnlo_6.in.
The calculation is very long (more than 15 minutes). For those who want to do it you anyway, you can copy ~abinit/tests/tutorial/Input/tnlo_6.in in your working directory. You should first modify the cutoff to ecut=5 since the cutoff of 2.8 is too low and yields unrealistic results. You can have a look at this input file. It contains 12 datasets. We need to compute the linear optical susceptibility (4 datasets for SC calculation in the IBZ, NSC calculation in the full BZ, ddk, ddE) and then to repeat this for 3 sets of atomic positions : these will be referred to as tau =0 for the reference structure and tau = +0.01 and tau = -0.01 for the Al atom displaced right and left by 0.01 Bohr (look at xcart to identify the 3 cases).
You can start the run. You have now time for a Belgian Beer, why not a Gouyasse ?! ... Or you can look at the results as summarized below.
For tau = 0 :
Dielectric tensor, in cartesian coordinates, j1 j2 matrix element dir pert dir pert real part imaginary part 1 4 1 4 9.9484608457 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.9484608460 0.0000000000 2 4 3 4 0.0000000000 0.0000000000
Dielectric tensor, in cartesian coordinates, j1 j2 matrix element dir pert dir pert real part imaginary part 1 4 1 4 9.9488281596 0.0000000000 1 4 2 4 0.0000000009 0.0000000000 1 4 3 4 0.0000000009 0.0000000000 2 4 1 4 0.0000000004 0.0000000000 2 4 2 4 9.9486321634 0.0000000000 2 4 3 4 -0.0133978975 0.0000000000
Dielectric tensor, in cartesian coordinates, j1 j2 matrix element dir pert dir pert real part imaginary part 1 4 1 4 9.9488281561 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.9486321602 0.0000000000 2 4 3 4 0.0133978972 0.0000000000
dchi_23/dtau_1= (1/4 pi) (eps_23[tau=+0.01] +eps_23[tau=-0.01])/(2*tau) = (1/4 pi) (-0.0133978975 -0.0133978972)/(0.02) = -0.106617
This value is close to that obtained at ecut=5 from DFPT (-0.099889084). When convergency is reached, both approaches allow to get the right answer. You might therefore ask which approach is the most convenient and should be used in practice.
As a guide, we can mention that the finite difference approach give very similar results for a similar cutoff and k-point grid. It is however more tedious because, individual atomic displacement must be successively considered (heavy for complex crystals) and the results must then be converted into appropriate units with risk of error of manipulations.
The DFPT approach is the most convenient and avoid a lot of human work. Everything is reported together (not only dchi/dtau but also the full Raman polarizability tensors) and in appropriate units. It might therefore be considered as the best choice when available as in ABINIT.
In this section, you will learn how to perform a Berry phase calculation. As a practical problem we will try to compute the Born effective charges from finite difference of the polarization (under finite atomic displacements).
You can now copy the file ~abinit/tests/tutorial/Input/tffield_1.in in Work-NLO. You can look at it. It is made of three datasets corresponding to the reference optimized structure (tau=0) and to structure with the Al atom displaced from 0.01 bohr right and left (referred to as tau = +0.01 and tau = -0.01). This is typically the amplitude of atomic displacement to be considered in this kind of computations.
There are two implementation of the berry phase within ABINIT. One corresponds to positive values of berryopt and was implemented by NaSai. The other one corresponds to negative values of berryopt was implemented by Marek Veithen. Both are suitable to compute the polarization. Here we will focus on the implementation of Marek Veithen for two reasons. First, the results are directly provided in cartesian coordinate at the end of the run (while the implementation of NaSai report them in reduced coordinates). Second the implementation of Marek Veithen is the one to be used for the finite electric field calculation as described in the next Section. Finally, note also that Veithen's implementation works with kptopt= 1 or 2 while Na Sai implementation is restricted to kptopt=2 which is less convenient.
The file is the one of a usual self-consistent calculation. On top of usual variables, for the berry phase calculation we simply need to define berryopt and rfdir:
berryopt -1 rfdir11 1 1 1
You can make the run. You can open the output file and look for the occurrence "Berry". The output report value of the Berry phase for individual k-point strings.
Computing the polarization (Berry phase) for reciprocal vector: 0.16667 0.00000 0.00000 (in reduced coordinates) -0.01583 0.01583 0.01583 (in cartesian coordinates - atomic units) Number of strings: 144 Number of k points in string: 6 Summary of the results Electronic Berry phase -1.633637909E-02 Ionic phase -7.369993171E-01 Total phase -7.533356962E-01 Remapping in [-1,1] -7.533356962E-01 Polarization -1.569029714E-02 (a.u. of charge)/bohr^2 Polarization -8.977166249E-01 C/m^2
tau = 0
Polarization in cartesian coordinates (a.u.): (the sum of the electronic and ionic Berry phase has been fold into [-1, 1]) Electronic: -0.589330682E-03 -0.589330682E-03 -0.589330682E-03 Ionic: -0.265870611E-01 -0.265870611E-01 -0.265870611E-01 Total: -0.271763918E-01 -0.271763918E-01 -0.271763918E-01 Polarization in cartesian coordinates (C/m^2): (the sum of the electronic and ionic Berry phase has been fold into [-1, 1]) Electronic: -0.337184150E-01 -0.337184150E-01 -0.337184150E-01 Ionic: -0.152117239E+01 -0.152117239E+01 -0.152117239E+01 Total: -0.155489081E+01 -0.155489081E+01 -0.155489081E+01
tau = +0.01
Polarization in cartesian coordinates (a.u.): (the sum of the electronic and ionic Berry phase has been fold into [-1, 1]) Electronic: -0.621428977E-03 -0.589348321E-03 -0.589348321E-03 Ionic: -0.264842841E-01 -0.265870611E-01 -0.265870611E-01 Total: -0.271057131E-01 -0.271764095E-01 -0.271764095E-01 Polarization in cartesian coordinates (C/m^2): (the sum of the electronic and ionic Berry phase has been fold into [-1, 1]) Electronic: -0.355549113E-01 -0.337194243E-01 -0.337194243E-01 Ionic: -0.151529203E+01 -0.152117239E+01 -0.152117239E+01 Total: -0.155084694E+01 -0.155489181E+01 -0.155489181E+01
tau = -0.01
Polarization in cartesian coordinates (a.u.): (the sum of the electronic and ionic Berry phase has been fold into [-1, 1]) Electronic: -0.557227851E-03 -0.589301321E-03 -0.589301321E-03 Ionic: -0.266898382E-01 -0.265870611E-01 -0.265870611E-01 Total: -0.272470660E-01 -0.271763625E-01 -0.271763625E-01 Polarization in cartesian coordinates (C/m^2): (the sum of the electronic and ionic Berry phase has been fold into [-1, 1]) Electronic: -0.318816591E-01 -0.337167352E-01 -0.337167352E-01 Ionic: -0.152705275E+01 -0.152117239E+01 -0.152117239E+01 Total: -0.155893441E+01 -0.155488913E+01 -0.155488913E+01
From the previous data, we can extract the Born effective charge of Al. Values to be used are those in a.u., in order to find the charge in electron unit. It corresponds to (the volume of the primitive uni cell must be specified in Bohr too) :
Z* = \Omega_0 (P[tau=+0.01] - P[tau=-0.01]) / (2*tau) = 305.39 (-0.271057131E-01 + -0.272470660E-01) / 0.02 = 2.16
Using the Berry phase approach it is also possible to compute the piezoelectric constant from finite difference of polarization with respect to strains. This can be done considering clamped ions or relaxed ions configurations. For those who are interested, they can have a look to the files ~abinit/tests/tutorial/Input/tffield_2.in and ~abinit/tests/tutorial/Input/tffield_3.in. They should give the following final results for ecut =5 :
Proper piezoelectric constants(clamped ion)(Unit:c/m^2) = -0.6909505 Proper piezoelectric constants(relaxed ion)(Unit:c/m^2) = -0.0326305
You can now copy the file ~abinit/tests/tutorial/Input/tffield_6.in in Work-NLO. It performs a finite field calculation at clamped atomic positions. You can look at this input file to identify the features specific to the finite field calculation.
As general parameters, it is recommended to specify nband, nbdbuf and kptopt:
nband 4 nbdbuf 0 kptopt 1
As a first step (dataset 11), the code must perform a Berry phase calculation in zero field. For that purpose, it is necessary to set the values of berryopt and rfdir:
berryopt11 -1 rfdir11 1 1 1
Important warning: you cannot put berryopt to +1 for a finite field calculation. You mandatory need berryopt -1.
After that, there are different steps corresponding to various value of the electric field, thanks to efield. For those steps, it is important to take care of the following parameters here reported for dataset 21:
berryopt21 4 efield21 0.0001 0.0001 0.0001 getwfk21 11
The electric field is here applied along the [111] direction. It must be incremented step by step (it is not possible to go to high field directly). At each step, the wavefunctions of the previous step must be used. When the field gets larger, it is important to check that it does not significantly exceeds the critical field for Zener breakthrough (the drop of potential over the Born-von Karmann supercell must be smaller than the gap). In practice, the number of k-point must be enlarged to reach the convergence. However, at the same time, the critical field becomes smaller. In practice, reasonable fields can still be reached for k-point grids providing a reasonable degree of convergence. A compromize must however be found.
As these calculations are quite long, the input file has been limited to very small fields. Three cases have been selected : E = 0, E = +0.0001 and E = -0.0001. If you have time later, you can relax this constraint and perform a more exhaustive calculations for a larger set of fields.
You can now start the run. Various quantities can be extracted from the finite field calculation at clamped ions using finite difference techniques: the Z* can be extracted from the forces, the optical dielectric constant can be deduced from the polarization, the clamped ion piezoelectric tensor can be deduced from the stress tensor. As an illustration we will focus here on the computation of Z*.
You can look at your output file. For each field, the file contain various quantities that can be collected. For the present purpose, we can look at the evolution of the forces with the field and extract the following results from the output file :
E=0
cartesian forces (hartree/bohr) at end: 1 0.00000000000000 0.00000000000000 0.00000000000000 2 0.00000000000000 0.00000000000000 0.00000000000000
E = +0.0001
cartesian forces (hartree/bohr) at end: 1 0.00020614538503 0.00020614538503 0.00020614538503 2 -0.00020614538503 -0.00020614538503 -0.00020614538503
E = -0.0001
cartesian forces (hartree/bohr) at end: 1 -0.00020650900243 -0.00020650900243 -0.00020650900243 2 0.00020650900243 0.00020650900243 0.00020650900243
In a finite electric field, the force on atom A in direction i can be written as :
F_Ai = Z*_A,ii * E + \Omega_0 dchi/dtau E^2
The value for positive and negative fields above are nearly similar showing that the quadratic term is almost negligible. This clearly appears in the figure below where the field dependence of the force for a larger range of electric fields is plotted.
We can therefore extract with a good accuracy the Born effective charge as :
Z*_Al = (F_Al[E=+0.0001] - F_Al[E=-0.0001])/(2*0.0001) = (0.00020614538503 + 0.00020650900243)/0.0002 = 2.06
This value is close to the value of 2.04 reported from DFPT. If you do calculations for more electric fields, fitting them with the general expression of the force above (including the E^2 term) allow to determine the dchi/dtau. You can try to do it. At ecut=5 you should get dchi/dtau for Al = -0.06169.
Going back to the output file, you can also look at the evolution of the polarization with the field.
E=0
Polarization in cartesian coordinates (a.u.): (the sum of the electronic and ionic Berry phase has been fold into [-1, 1]) Electronic: 0.730821479E-04 0.730821479E-04 0.730821479E-04 Ionic: -0.270560574E-01 -0.270560574E-01 -0.270560574E-01 Total: -0.269829753E-01 -0.269829753E-01 -0.269829753E-01
E = +0.0001
Polarization in cartesian coordinates (a.u.): (the sum of the electronic and ionic Berry phase has been fold into [-1, 1]) Electronic: 0.129869542E-03 0.129869542E-03 0.129869542E-03 Ionic: -0.270560574E-01 -0.270560574E-01 -0.270560574E-01 Total: -0.269261879E-01 -0.269261879E-01 -0.269261879E-01
E = -0.0001
Polarization in cartesian coordinates (a.u.): (the sum of the electronic and ionic Berry phase has been fold into [-1, 1]) Electronic: 0.163575373E-04 0.163575372E-04 0.163575371E-04 Ionic: -0.270560574E-01 -0.270560574E-01 -0.270560574E-01 Total: -0.270396999E-01 -0.270396999E-01 -0.270396999E-01
In a finite electric field, the force on atom A in direction i can be written as :
P = \chi * E + \chi(2) E^2
The change of polarization for positive and negative fields above are nearly similar showing that the quadratic term is almost negligible. This clearly appears in the figure below where the field dependence of the force for a larger range of electric fields is plotted.
We can therefore extract with a good accuracy the linear optical dielectric susceptibility :
\chi(l)_11 = (P[E=+0.0001] - P[E=-0.0001])/(2*0.0001) = (-0.269261879E-01+ 0.270396999E-01)/0.0002 = 0.56756
The optical dielectric constant is :
\epsilon_11 = 1+ 4 \pi \chi(1)_11 = 8.13
This value underestimate the value of 9.20 obtained by DFPT. The agreement is not better at ecut=5. If you do calculations for more electric fields, fitting them with the general expression of the force above (including the E^2 term) allows one to determine the non-linear optical susceptibilities \chi(2). You can try to do it. At ecut=5 you should get \chi(2) = 29.769 pm/V.
Looking at the evolution of the stress (see graphic below), you should also be able to extract the piezoelectric constants. You can try to do it as an exercice. As the calculation here was at clamped ions, you will get the clamped ion proper piezoelectric constants. At ecut=5, you should obtain -0.69088 C/m^2.
Redoing the same kind of finite field calculation when allowing to the ions to relax (using the input file ~abinit/tests/tutorial/Input/tffield_7.in) you can access the relaxed ion proper piezoelectric constant. At ecut=5 the result is -0.03259.