Polarization calculations with the Berry-phase approach¶
In this exercise you will calculate the macroscopic polarization in a polar material using the Berry-phase approach of the modern theory of polarization.
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. Note that these calculations only make sense for insulators!
Note
You can follow these slides from Javier Junquera as an introduction.
For further references, see R.D. King-Smith and D. Vanderbilt, PRB 47, 1651 (1993) and Philippe Ghosez’s thesis.
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 then proportional to the Born or dynamical effective charges and can be extremely large for ferroelectric materials.
Calculating the Polarization¶
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
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 the 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.
Once your calculation is done, you will find the following lines in your SIESTA output file:
siesta: Macroscopic polarization per unit cell (Debye):
siesta: Along the lattice vectors 15.331110 15.331110 15.331110
siesta: Along cartesian directions 21.681464 21.681464 21.681464
You can play with the numbers in the PolarizationGrids block to check for
the convergence of the polarization as a function of the fineness of the grid.
Notice, however, that the absolute value of the polarization lacks phyisical
meaning in this context. Only changes in polarization are physically relevant.
Polarization as a function of the bond length in c-BN¶
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; however, an interesting aspect is that the calculated curve might not be smooth. In truth, the polarization per unit cell will be defined in terms of multiples of a quantum of polarization, 2eR, where e is the electron charge and R an arbitrary vector of the Bravais lattice. This gap in the polarization cannot be predicted, but you can verify after the fact that is indeed a multiple of 2eR. Afterwards, you can subtract this gap from the rest of the function and therefore obtain a smooth curve.
A question you may be asking yourself is: We said that only the changes in polarization make sense, not the absolute values. So, 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\).
We 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 .true. 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 (see at the bottom of this tutorial for the specifics).
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 (we have added lines and the x,y,z symbols as a guide):
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).