Polarization calculations with the Berry-phase approach¶
In this exercise you will calculate the macroscopic polarization in a polar material using the so-called Berry-phase approach of the modern theory of polarization (as a reference, see R.D. King-Smith and D. Vanderbilt, PRB 47, 1651 (1993))
We will obtain the behavior of the polarization as a function of the atomic displacements. With this information you can calculate the so-called dynamical or Born charges that characterize the interaction between a macroscopic electric field and the elastic deformations of the material.
Experimentally the dynamical charges can easily be probed by measuring the LO-TO splitting. This is a splitting between vibrational modes (difference of their frequencies) that, in principle, could be degenerate due to the symmetry of the crystal when q–>0 (q being the wavevector of the vibration). However, when long range electrostatic forces are introduced these modes behave in a different way when q–>0. LO modes produce accumulations of charge along the q direction and, for 3-D crystals, this translates into a macroscopic electric field when q–>0. In contrast, TO modes are not charged and not affected by the long range electrostatic part of the dynamical matrix. The LO-TO splitting is proportional to the Born or dynamical effective charges and can be extremely large for ferroelectric materials.
Note
Important note: These calculations only make sense for insulators!!!!!
Polarization as a function of the bond length in c-BN¶
Hint
Go to directory ‘c-BN.pol_vs_bond_elongation`
You will find there two SIESTA input files: ‘c-BN.init.fdf’ and ‘c-BN.elongated.fdf’.
At the end of the c-BN.init.fdf file you will notice the presence of a peculiar block:
%block PolarizationGrids
10 4 4
4 10 4
4 4 10
%endblock PolarizationGrids
This block contains the information about the k-space grids that will be used in order to compute the macroscopic polarization using the geometric Berry phase approach.
Only if this block is present, the calculation of the macroscopic polarization will be performed.
In this approach, to compute the polarization along a given direction it is necessary to perform:
a 1-D integral in k-space along the chosen direction and
ii) a 2-D integration of the result of the previous integral within the plane perpendicular (also in k-space) to the chosen direcction.
Each line in block PolarizationGrids corresponds to a grid used to calculate the polarization along the given lattice vectors: the first row specifies the grid that will be used to calculate the polarization along the direction of the first lattice vector, the second row will be used for the calculation along the direction of the second lattice vector, and the third row for the third lattice vector.
The numbers in the diagonal of the matrix specify the number of points to be used in the one-dimensional line integrals along the different directions. The other numbers specify the mesh used in the surface integrals. If the number of point in one of the grids is zero, the calculation will not be performed for this particular direction.
You can play with the numbers in this block to check for the convergence of the polarization as a function of the fineness of the grid.
Notice, however, that only changes in polarization are physically relevant. The file c-BN.elongated.fdf contains the input for a c-BN system with one elongated BN bond. As a consequence the system has lost its perfect cubic symmetry and will develop a polarization along the direction parallel to the elongated bond.
You should calculate the behavior of a polarization as a function of this distorsion. An interesting aspect is that the calculated curve might not be smooth. The polarization per unit cell can only be defined modulo a quantum of polarization: 2eR, where e is the electron charge and R and arbitrary vector of the Bravais lattice. Once this effect is taken into account you should obtain a smooth curve.
A question you should answer: We said that only the changes of polarization make sense. Is there a “right way” to set the zero of this curve?
Calculation of dynamical charges for c-BN¶
If you apply small enough displacements you can calculate the Born or dynamical charges, which are defined as:
where e is the charge of an electron and \(\Omega_0\) is the unit cell volume. The \(\alpha\beta\) component of the (tensor) charge measures how the polarization changes along the \(\alpha\) direction when the atom i is displaced along the \(\beta\) direction.
For a cubic crystal, like the c-BN in this case, the dynamical-charge tensor becomes a scalar \(Z^*_i\).
You can check that this is indeed the case performing calculations for Z(N) and Z(B) along different directions.
You can do this manually, but Siesta provides a way to automate the
process of generating the appropriate displacements involved in the
calculation of the effective charges. It uses the ‘FC’ mode typically
used for phonon calculations, and the Born-charge T option to
request the computation of the Born charges:
MD.TypeOfRun FC
MD.FCDispl 0.01 bohr
BornCharge .true.
You can run a full calculation of the Born effective charges for cubic BN in directory c-BN.BornCharges. The results appear in the c-BN.BC file.
The dynamical charges obey the so-called acoustic sum-rule: the sum of the dynamical charges for all the atoms in the cell should be zero. In our case this means that Z(N)=-Z(B). You should check to what extent this is true in your calculations, and how the fulfillment of this rule depends on the fineness of the utilized k-space grids.
Note
As a by-product, you also get a force-constants file c-BN.FC
which can be read by vibra to compute the phonons (in this case
only at the Gamma point). Just include this in the fdf file:
Eigenvectors .true.
LatticeDielectricConstant .true.
%block BandLines
1 0.0 0.0 0.0 \Gamma
%endblock BandLines
and run vibra:
vibra < c-BN.fdf
You will get the file c-BN.bands with the eigenvalues and c-BN.vectors with the eigenvectors.
Born charges for h-BN¶
Hexagonal BN is layered material. Therefore, ‘in plane’ and ‘out of plane’ dynamical charges might differ. Directories ‘h-BN.Z_parallel’ and ‘h-BN.Z_perpendicular’ contain input files to perform these calculations manually, but you can automate the process using the fdf file in directory h-BN.BornCharges. Note that we only compute the charges for the first pair of atoms.
Note
The distance between different layers in h-BN can influence the values of the dynamical charges, particularly for out of plane movements. Check this dependence performing a calculation for a larger value of the c/a ratio.
Format of the BC file¶
Siesta writes the Born effective charges in a file called SystemLabel.BC It should look like this (the lines and the x,y,z symbols have been included to guide your eyes:
BC matrix x y z x 2.7661178 0.0000107 -0.0000000 y -0.0085767 2.7538482 0.0000000 Born effective charges for B z 0.0004767 0.0000000 0.7117852 -------------------------------------------------- x -2.7471122 -0.0000108 0.0000000 y -0.0094377 -2.7590343 -0.0000000 Born effective charges for N z 0.0001423 0.0000000 -0.7309126 -------------------------------------------------- ....
(The example is for hexagonal BN, showing the charges for the first couple of atoms only.)
For each atom, the effectives charges are a tensor of dimension (3x3).
For more background, the reader is strongly recommended to read pages 106 and following of Philippe Ghosez’s thesis.