Skip to content

Exercise 5: Polarization, Born Effective Charges, and Non-analytical term correction

Estimated total CPU time: 15 min

Warning

In the following exercises, computational settings including the reciprocal space grid (tag k_grid), the basis set, and supercell's size, have been chosen to allow for a rapid computation of the exercises in the limited time and within the CPU resources available during the tutorial session. Without loss of generality, these settings allow to demonstrate trends of the lattice dynamics of materials. In the production calculation, all computational parameters should be converged.

Important

The general theory needed for this exercise was described during the lecture. In the following exercise we do not describe the theoretical part in details, but provide a brief summary. To deeper familiarize yourself with the theory behind, explore the Phonopy Formulations and books on the fundamental of lattice dynamics 7,8,9 or at least read this Physics stachexchange issue.

In this exercise, you will:

  • Learn how to compute Born Effective Charges (BEC) tensor in FHI-aims.
  • Learn how to compute dielectric tensor in FHI-aims.
  • Use non-analytical term correction (NAC) to observe longitudinal optical-transverse optical (LO-TO) splitting in the MgO crystal.

Non-metallic, polar crystals are polarized by atomic displacements and the generated macroscopic field changes the force constants near the \(\Gamma\) point. To account for this contribution the non-analytical term correction (NAC)7,8,9 is calculated. We will use Phonopy following the work by Wang and co-workers5 to compute the NAC. In practice, Phonopy requires two input quantities for the calculation of the NAC:

  • The dielectric tensor and
  • Born Effective Charges (BEC) tensor.

In this exercise, Density Functional Perturbation Theory (DFPT)4,6 is used for the calculation of the dielectric tensor. The general theory and details of the polarization and BEC calculations in the FHI-aims code are presented in the paper by Carbogno et al10. As an example, we will use the rock salt MgO.

Calculation of the dielectric tensor

To compute the dielectric tensor with DFPT using FHI-aims, we need to add the following flag to the control.in file: 

DFPT dielectric

Go to the folder: 

phonons-with-fhi-vibes/Tutorial/phonons/5_BEC/dielectric

to inspect the control.in file for our MgO calculation. In the same folder you will also find the geometry.in with the primitive unit cell of MgO. You can use both input files to start the calculation with FHI-aims (e.g. mpirun -n 4 aims.x > aims.out or similar). If you don't know how to start the FHI-aims calculation, please read the tutorial Basics of Running FHI-aims.

After the calculations are finished, take a look at the output file. You should find the following: 

DFPT for dielectric_constant:--->
   3.25462982052995       4.508771475482781E-015  4.047716164544250E-014
  4.507996563003625E-015   3.25462982609752      -4.211307534839892E-015
  4.042527887519123E-014 -4.238589729672110E-015   3.25462982605836

Here, the dielectric tensor (high-frequency one) is printed. Note that your results could be slightly different, because of the numerical accuracy. In the case of MgO, the dielectric tensor is diagonal with identical numbers on the diagonal (approximately 3.2546).

Calculation of Polarization and the Born Effective Charges

As a next step, we calculate the Born effective charges (BEC). The concept of BEC formalizes the relationship between the change in polarization \(\Delta P\) once with atomic displacement, i.e.

\[ \Delta P_{i} = \frac{e}{\Omega}\sum_{k,j} Z^{\ast}_{k,ij} \Delta R_{k,j} \]

where the index \(k\) runs over the atoms, the indexes \(i,j\) denote Cartesian directions, \(\Omega\) is the unit cell volume. The Born Effective Charge1 is determined by2:

\[ Z^{\ast}_{k,ij} = \frac{\Omega}{e} \frac{\partial P_{i}}{\partial R_{k,j}} \]

In practice, the partial derivative is approximated by a finite difference:

\[ Z^{\ast}_{k,ij}=\frac{\Omega}{e} \frac{\Delta P_{i}}{\Delta R_{k,j}}. \qquad (1) \]

Phonopy needs only Z\(^{\ast}\) of symmetrically non-equivalent atoms as input. In the case of the primitive MgO unit cell, both atoms are symmetrically non-equivalent. However, one should carefully evaluate the non-equivalent atoms of the (primitive) unit cell under the study.

To proceed with the calculation of the BEC (usually labeled as \(Z^{\ast}\) in literature9) for MgO, change the directory to: 

phonons-with-fhi-vibes/Tutorial/phonons/5_BEC/born_charges

As a first step, we need to create a supercell (we generate a supercell, just to illustrate that some atoms are equivalent and in practice the calculations should be done with just a primitive cell). FHI-vibes can automatically suggest a cubic supercell (or at least as close as possible to the cubic shape) for a certain number of atoms. Type:

vibes utils make-supercell -n 8 geometry.in.primitive

FHI-vibes then produces cubic conventional cell of MgO with 8 atoms. Rename the file of the supercell to geometry.in.

To trigger the calculations of polarization in FHI-aims, you need to specify the following line in the control.in file: 

output polarization 1 20 10 10

The first number specifies that polarization is computed along the first lattice vector and the last three integers corresponds to the number of k-points taken along each of the reciprocal lattice vectors for calculations of the wannier center of charge evolution. To familiarize yourself with the theory behind, you could inspect the beginner's guide to the modern theory of polarization12. Note that the number of k-points is a parameter and needs to be tested for convergence.

There are two python scripts, which help with the automation of the BEC calculations, namely:

  • preprocessing.py
  • postprocessing.py

The code inside them is self explanatory -- please, inspect it. First, run the first script by:

python3 preprocessing.py

which creates one directory for each geometry with one displaced atom (for MgO, there will be two folders). The corresponding control.in file is copied to each of the directories, too. As a next step, start the FHI-aims calculation. Please use aims.out as the name for the FHI-aims output file (or adapt postprocessing script afterwards for your aims-output file name). After both calculations have finished, please run the second script:

python3 postprocessing.py
Then you would obtain BEC's for Mg and O atoms:

The displacement is 1% from the lattice vector lenghts: 0.04226199420000001 AA
Volume of the cell is: 75.48314005099871 AA^3
Z* for Mg is 1.9785090242786756
Z* for O is -1.9792159848836068
Note that because of the cubic symmetry the BEC and dielectric tensors are diagonal. Thus, in the calculations of BECs we have only displaced atom along one direction. This script extracted the necessary information from the aims.out files. The polarization output in the aims.out file looks similar to the following:

------------------------------------------------------------

  Summary for all Polarization Directives:

  * Ionic Polarization along rec. latt. vecs. is (C/m^2):       107.644455E-03        0.000000E+00        0.000000E+00
  * Electronic Contributions:
    Summarizing all directives:
    - Directive    1 in direction of rec. latt. vec.  1 yields the electronic polarization (Berry term):         -87.758697E-03 (C/m^2)
    - Directive    1 in direction of rec. latt. vec.  1 yields the electronic polarization (Dipole term):         -2.139754E-03 (C/m^2)
  * Full Polarizations (Electron + Ions):
    Summarizing all directives:
    - Directive    1 in direction of rec. latt. vec.  1 yields the full polarization      :        17.746004E-03 (C/m^2)

------------------------------------------------------------

The script postprocessing.py extracted the value for the full polarization and calculates the BEC according to the Eq. (1). Note that BEC's of Mg and O in principle should sum up to zero -- slight deviation from this rule comes from the accuracy of our calculations.

For more general cases (non-cubic symmetries, more none-equivalent atoms, etc.) the BEC.py script contained in the utilities folder in root directory of FHI-aims distribution can be used to simplify the work. An example on how to use the BEC.py script for the MgO system discussed here can be found in the solutions folder of this exercise. More Information about the usage of the script is given the FHI-aims manual, see Sec. How to calculate Born Effective Charges?.

Calculation of phonon band structure including NAC

Now, we can generate the BORN file, which Phonopy can read and use to apply the NAC5. We already prepared the phonons data for the MgO (cubic supercell with 64 atoms) following the recipes we learned in exercise 2. The folder:

phonons-with-fhi-vibes/Tutorial/phonons/5_BEC/MgO_phonons

contains the necessary data. Change directory to that folder and make a file called BORN, which should contain:

14.400
3.254629 0 0 0 3.254629 0 0 0 3.254629
1.9786 0 0 0 1.9786 0 0 0 1.9786
-1.9793 0 0 0 -1.9793 0 0 0 -1.9793

The first line is a conversion factor used by Phonopy to utilize units from FHI-aims. The second line is a dielectric tensor (9 components) written in one line (the resuts, which were obtained above, but slightly rounded for convenience). The next two lines are the BEC tensor (again 9 components) for Mg and O atoms, respectively. Since our conventional unit cell of MgO contains only two non-equivalent atoms BEC should be presented only for them. Once the born file is created, we can calculate the phonon band structure including the NAC. In order to make comparison between the band structures with and without NAC, we will plot both. First, we compute the band structure without NAC This is done as before, by typing:

 vibes output phonopy phonopy/trajectory.son -bs

Then, the folder phonopy contains a folder output containing the band structure, plotted and stored in pdf file3. Rename the folder output to output_without_nac.

In order to account for the NAC, type:

 vibes output phonopy phonopy/trajectory.son --born BORN -bs

This command would again produce output folder, but this time including the NAC. For this second plot of the band structure, you find the LO-TO splitting at the \(\Gamma\) point as presented in the figure below.

No NAC:

With NAC:

The obtained results could be compared to the experimental one13, where you could clearly see LO-TO splitting at the \(\Gamma\) point.

  1. The BEC can differ considerably from the formal (static) charge of a given atom. It is often referred to as a dynamic charge, because it embodies the dynamical contribution of the atomic lattice to the polarization. This dynamic quantity enters the description of the LO-TO splitting between phonon modes, given by the long-range electric fields generated by the long-wavelength longitudinal phonons and is an essential component for the analysis of polar materials, ferroelectricity, piezoelectricity, and others 11. ↩

  2. Analogously, this charge can be thought of as the response of the atomic forces upon an applied electric field, but here we will stick to the above definition. ↩

  3. Note that you would probably again see small (of the order of 0.001) imaginary (negative) frequencies at the \(\Gamma\) point -- this is a numerical artifact. ↩

  4. S. Baroni, et.al.,Rev. Mod. Phys. 73, 515 (2001) ↩

  5. Y. Wang,et.al.,J. Phys. Condens. Matter. 22, 202201 (2010) ↩↩

  6. H. Shang, C. Carbogno, P. Rinke, M. Scheffler, Comp. Phys. Com. 215. (2017) ↩

  7. R. M. Pick, M. H. Cohen, and R. M. Martin, Phys. Rev. B 1, 910, (1970). ↩↩

  8. M. Dove, Introduction to Lattice Dynamics, Cambridge University Press (1993). ↩↩

  9. [M. Born, and K. Huang, Dynamical Theory of Crystal Lattices, Oxford University Press (1962).] ↩↩↩

  10. [C. Carbogno, et.al, in preparation (2022).] ↩

  11. X. Gonze, C. Lee, Phys. Rev. B, 55, 10355-10368 (1997) ↩

  12. N. Spaldin, Journal of Sol. State. Chem., 195, 2-10, (2012) ↩

  13. M. J. L. Sangster et. al., J. Phys. C: Solid State Phys., 3, 1026, (1970) ↩