Exercise 5: Polarization, Born Effective Charges, and Nonanalytical 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 FHIaims.
 Learn how to compute dielectric tensor in FHIaims.
 Use nonanalytical term correction (NAC) to observe longitudinal opticaltransverse optical (LOTO) splitting in the MgO crystal.
Nonmetallic, 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 nonanalytical term correction (NAC)^{7},^{8},^{9} is calculated. We will use Phonopy following the work by Wang and coworkers^{5} 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 FHIaims code are presented in the paper by Carbogno et al^{10}. As an example, we will use the rock salt MgO.
Calculation of the dielectric tensor
To compute the dielectric tensor with DFPT using FHIaims, we need to add the following flag to the control.in
file:
DFPT dielectric
Go to the folder:
phononswithfhivibes/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 FHIaims (e.g. mpirun n 4 aims.x > aims.out
or similar).
If you don't know how to start the FHIaims calculation, please read the tutorial Basics of Running FHIaims.
After the calculations are finished, take a look at the output file. You should find the following:
DFPT for dielectric_constant:>
3.25462982052995 4.508771475482781E015 4.047716164544250E014
4.507996563003625E015 3.25462982609752 4.211307534839892E015
4.042527887519123E014 4.238589729672110E015 3.25462982605836
Here, the dielectric tensor (highfrequency 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.
where the index \(k\) runs over the atoms, the indexes \(i,j\) denote Cartesian directions, \(\Omega\) is the unit cell volume. The Born Effective Charge^{1} is determined by^{2}:
In practice, the partial derivative is approximated by a finite difference:
Phonopy needs only Z\(^{\ast}\) of symmetrically nonequivalent atoms as input. In the case of the primitive MgO unit cell, both atoms are symmetrically nonequivalent. However, one should carefully evaluate the nonequivalent atoms of the (primitive) unit cell under the study.
To proceed with the calculation of the BEC (usually labeled as \(Z^{\ast}\) in literature^{9}) for MgO, change the directory to:
phononswithfhivibes/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). FHIvibes 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 makesupercell n 8 geometry.in.primitive
FHIvibes 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 FHIaims, 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 kpoints 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 polarization^{12}. Note that the number of kpoints 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 FHIaims calculation. Please use aims.out
as the name for the FHIaims output file (or adapt postprocessing script afterwards for your aimsoutput file name). After both calculations have finished, please run the second script:
python3 postprocessing.py
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
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.644455E03 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.758697E03 (C/m^2)
 Directive 1 in direction of rec. latt. vec. 1 yields the electronic polarization (Dipole term): 2.139754E03 (C/m^2)
* Full Polarizations (Electron + Ions):
Summarizing all directives:
 Directive 1 in direction of rec. latt. vec. 1 yields the full polarization : 17.746004E03 (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 (noncubic symmetries, more noneequivalent atoms, etc.) the BEC.py
script contained in the utilities folder in root directory of FHIaims 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 FHIaims 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 NAC^{5}.
We already prepared the phonons data for the MgO (cubic supercell with 64 atoms) following the recipes we learned in exercise 2. The folder:
phononswithfhivibes/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 FHIaims. 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 nonequivalent 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 file^{3}. 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 LOTO splitting at the \(\Gamma\) point as presented in the figure below.
No NAC:
With NAC:
The obtained results could be compared to the experimental one^{13}, where you could clearly see LOTO splitting at the \(\Gamma\) point.

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 LOTO splitting between phonon modes, given by the longrange electric fields generated by the longwavelength longitudinal phonons and is an essential component for the analysis of polar materials, ferroelectricity, piezoelectricity, and others ^{11}. ↩

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. ↩

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. ↩

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

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

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

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

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

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

N. Spaldin, Journal of Sol. State. Chem., 195, 210, (2012) ↩

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