Introduction to Polarization Calculations

This tutorial describes how to calculate spontaneous polarization using density functional theory and modern theory of polarization.

Spontaneous Polarization

Ferroelectrics play a significant role in various technological applications, including tunable capacitors, non-volatile random access memory devices, and electro-optical data storage. The phenomenon of ferroelectricity often arises due to a structural phase transition from a high-symmetry nonpolar phase to a low-symmetry polar phase as temperature decreases, leading to the spontaneous emergence of polarization.

Ferroelectrics exhibit a distinctive hysteresis loop in the polarization versus electric field curve.

The spontaneous polarization can be experimentally determined by taking half of the change in polarization at zero external fields.

The spontaneous polarization of a material is not directly observable. Instead, it is determined by measuring the change in polarization between two stable configurations of the material.

Following the modern theory of polarization, the polarization, P, of a crystal is defined as,

$$P= P_0 + \sum_{i \in \{{\rm a,b,c\}}} n_i \frac{eR_i}{\Omega}$$

where e is the charge of the electron, ni is an integer, Ri is a primitive lattice vector, and Ω is the unit cell volume. The second term on the right-hand side of the equation is the quantum of polarization, $eR_{i}/\Omega$, a consequence of translation symmetry.

Only differences in the computed polarization on the “same branch” are physically meaningful, so one needs a polar structure and a non-polar reference structure in order to perform the polarization calculations.

The following section describes how to compute spontaneous polarization values using density functional theory and the modern theory of polarization

Getting Polar and non-Polar structure

To compute the spontaneous polarization using the berry phase method implemented in VASP.

1: Relax the polar structure using VASP.

2: Create an anti-polar structure from the relaxed polar structure.

3: One can use vtst scripts to create several images in between these two endpoints, for example for AlN wurtzite structure, the interpolation path is such that the middle image looks like hexagonal boron nitride.

schematic representation of polar, non-polar and anti-polar structure

use nebmake.pl script to create 10 images between the polar and anti-polar endpoints with the following command line

nebmake.pl POSCAR_polar POSCAR_anti-polar 10 

4: Use these images to calculate the value of spontaneous polarization. Use the flags LCALPOL = .TRUE. and DIPOL = 0.5 0.5 0.5 .

ISYM = 0
NSW = 0
ISTART = 0
SYSTEM = AlN
PREC = Accurate
LCALCPOL = .TRUE.
DIPOL = 0.5 0.5 0.5
LORBIT = 10
ADDGRID = .TRUE.
IBRION = -1
ISMEAR = 0
LWAVE = .TRUE.
SIGMA = 0.01
ALGO = Normal

5: After completing this step collect all the POSCAR and OUTCAR files in a directory with proper indexing. e.g POSCAR-00 .... POSCAR-10, OUTCAR-00....OUTCAR-10

to get the polarization branch use pymatgen utility

If you are working with many Python packages, it is generally recommended you create a separate environment for each of your packages. For example:

conda create --name my_pymatgen python
source activate my_pymatgen  # OSX or Linux
activate my_pymatgen  # Windows

Now using pymatgen utility you can extract the spontaneous polarization values from the POSCAR and OUTCAR files collected in step 5.

from pymatgen.core.structure import Structure
from pymatgen.analysis.ferroelectricity.polarization import Polarization, calc_ionic, zval_dict_from_potcar, get_total_ionic_dipole
from pymatgen.io.vasp.inputs import Poscar, Potcar
from pymatgen.io.vasp.outputs import Outcar

num_structs = 10
poscars = [Structure.from_file("POSCAR-" + str(i)) for i in range(num_structs)]
outcars = [Outcar("OUTCAR-" + str(i)) for i in range(num_structs)]
potcar=Potcar.from_file("POTCAR")
for i in range(num_structs):
    outcars[i].zval_dict = zval_dict_from_potcar(potcar)


pol_from_out_struct_method = Polarization.from_outcars_and_structures(outcars,poscars,
    calc_ionic_from_zval = True)
print(pol_from_out_struct_method.get_same_branch_polarization_data(convert_to_muC_per_cm2=True))
~                                                                                                    

The spontaneous polarization can be calculated by using the formula

$$P_{S} = 0.5*(P_{polar} - P_{anti-polar})$$

Polarization curves for AlN along the $𝑐̂$ direction versus distortion from the polar to non-polar to anti-polar structure.

the static energy barrier for the polar to anti-polar switching. The middle point (non-polar) is halfway between polar and anti-polar structures.

The relative energy per formula shows the switching barrier from polar to anti-polar structures assuming a smooth homogeneous switching from polar to anti-polar structure via non-polar h-BN type intermediate. In order to get a proper switching barrier one needs to do ss-NEB calculations. One can find this useful tutorial on how to perform ss-NEB calculations