1. Multi-valuedness of polarization and polarization branch¶
Polarization in a periodic crystal exhibits multi-valuedness, depending on the choice of unit cells used to define the electric dipoles. The key point is that one should focus on changes in polarization rather than absolute values, as these changes are uniquely defined and directly correspond to experimentally measurable quantities such as the switching current.
As illustrated in Fig. 1, polarization values \(\mathbf{P_1}\) and \(\mathbf{P_2}\) obtained from different choices of unit cells can even point in opposite directions. Their difference, \(\mathbf{P_q}=\mathbf{P_1}-\mathbf{P_2}\), is referred to as a polarization quanta. For a three-dimensional crystal, the polarization quanta is given by
\[
\mathbf{P}_{\mathbf{q},i} = \frac{1}{\Omega} e \mathbf{R}_i,
\]
where \(\Omega\) is the cell volume, \(e\) the elementary charge, \(\mathbf{R}_i\) the lattice basis vector in the \(i\) direction.
Figure 1. Multi-valued polarization in a one-dimensional chain of alternating anions and cations. Dashed rectangles indicate unit cells with lattice parameter $a$. Extracted from Ref. [1].
To determine the spontaneous polarization of a polar structure, one should specify a nonpolar reference structure and establish a continuous polarization branch connecting the nonpolar and polar states. As shown in Fig. 2, the polarization evolves monotonically along the continuous branch; therefore, the correct polarization corresponds to \(\mathbf{P_1}\) rather than \(\mathbf{P_2}\).
Figure 2. Polarization as a function of percentage ferroelectric distortion. Extracted from Ref. [1].
In the absence of a priori knowledge of the polarization value for a polar structure, it is necessary to perform a series of Berry-phase calculations for a set of intermediate structures connecting a specified nonpolar reference to the target polar structure. Here, I would like to share a bash script that aligns the polarization branch using the OUTCAR and POSCAR files obtained from these calculations.
In this script, the raw polarization values are obtained by dividing the total dipole moments (ionic + electronic) by the corresponding cell volumes. The target polarization components (Pa_tmp, Pb_tmp, Pc_tmp) for the first, nonpolar reference structure are set to zero. For subsequent structures, the polarization branch is aligned such that the change in polarization between two adjacent structures remains smaller than half of the polarization quantum PQ, i.e., within the range (-PQ/2, PQ/2). Consequently, a sufficient number of intermediate structures must be included, particularly for large supercells where the polarization quanta is typically smaller than that of a unit cell.
Note that the dipole moments reported in the OUTCAR files are expressed in Cartesian coordinates (x, y, z), whereas the polarization quanta are defined along the lattice basis vectors. Therefore, the script includes a coordinate transformation between the Cartesian and lattice basis (a, b, c) representations. This step is particularly important for non-orthogonal lattices.
#!/bin/bash# This script calculates the polarization variation from nonpolar phase to polar phase, and prints the polarization quanta# To avoid ambiguity, it is recommended to set the lattice cell of the nonpolar phase the same as that of the polar phase# The polarization change between two adjacent structures must be smaller than half of the polarization quanta, i.e., within (-PQ/2, PQ/2)# This script involves transformation between Cartesian coordinates (x,y,z) and lattice basis coordinates (a,b,c)n=5# number of structuresconst=1602.176565# polarization unit: e/A^2 --> uC/cm^2Pa_tmp=0;Pb_tmp=0;Pc_tmp=0echo"P(uC/cm2)""Pa""Pb""Pc""Ptot""PQa""PQb""PQc"|awk'{printf "%10s %6s %6s %6s %6s %6s %6s %6s\n",$1,$2,$3,$4,$5,$6,$7,$8}'>P.dat
for((i=1;i<=n;i++));doDix=$(grep"Ionic dipole moment"$i/OUTCAR|awk'{print $5}')# dipole: (x,y,z)Diy=$(grep"Ionic dipole moment"$i/OUTCAR|awk'{print $6}')Diz=$(grep"Ionic dipole moment"$i/OUTCAR|awk'{print $7}')Dex=$(grep"Total electronic dipole moment"$i/OUTCAR|awk'{print $6}')Dey=$(grep"Total electronic dipole moment"$i/OUTCAR|awk'{print $7}')Dez=$(grep"Total electronic dipole moment"$i/OUTCAR|awk'{print $8}')Dx=$(echo"(-1)*($Dix+$Dex)"|bc-l)# dipole: electrons Angst --> eADy=$(echo"(-1)*($Diy+$Dey)"|bc-l)Dz=$(echo"(-1)*($Diz+$Dez)"|bc-l)a=($(sed-n'3 p'$i/POSCAR))b=($(sed-n'4 p'$i/POSCAR))c=($(sed-n'5 p'$i/POSCAR))La=$(echo"sqrt(${a[0]}^2+${a[1]}^2+${a[2]}^2)"|bc-l)Lb=$(echo"sqrt(${b[0]}^2+${b[1]}^2+${b[2]}^2)"|bc-l)Lc=$(echo"sqrt(${c[0]}^2+${c[1]}^2+${c[2]}^2)"|bc-l)V=$(echo"${a[0]}*(${b[1]}*${c[2]}+(-1)*${b[2]}*${c[1]})+${a[1]}*(${b[2]}*${c[0]}+(-1)*${b[0]}*${c[2]})+${a[2]}*(${b[0]}*${c[1]}+(-1)*${b[1]}*${c[0]})"|bc-l)for((j=0;j<=2;j++));doa_norm[$j]=$(echo"${a[j]}/$La"|bc-l)b_norm[$j]=$(echo"${b[j]}/$Lb"|bc-l)c_norm[$j]=$(echo"${c[j]}/$Lc"|bc-l)donedeterminant=$(echo"${a_norm[0]}*(${b_norm[1]}*${c_norm[2]}+(-1)*${b_norm[2]}*${c_norm[1]}) \ +${a_norm[1]}*(${b_norm[2]}*${c_norm[0]}+(-1)*${b_norm[0]}*${c_norm[2]}) \ +${a_norm[2]}*(${b_norm[0]}*${c_norm[1]}+(-1)*${b_norm[1]}*${c_norm[0]})"|bc-l)a_norm_inv[0]=$(echo"(${b_norm[1]}*${c_norm[2]}+(-1)*${b_norm[2]}*${c_norm[1]})/$determinant"|bc-l)a_norm_inv[1]=$(echo"(${c_norm[1]}*${a_norm[2]}+(-1)*${c_norm[2]}*${a_norm[1]})/$determinant"|bc-l)a_norm_inv[2]=$(echo"(${a_norm[1]}*${b_norm[2]}+(-1)*${a_norm[2]}*${b_norm[1]})/$determinant"|bc-l)b_norm_inv[0]=$(echo"(${b_norm[2]}*${c_norm[0]}+(-1)*${b_norm[0]}*${c_norm[2]})/$determinant"|bc-l)b_norm_inv[1]=$(echo"(${c_norm[2]}*${a_norm[0]}+(-1)*${c_norm[0]}*${a_norm[2]})/$determinant"|bc-l)b_norm_inv[2]=$(echo"(${a_norm[2]}*${b_norm[0]}+(-1)*${a_norm[0]}*${b_norm[2]})/$determinant"|bc-l)c_norm_inv[0]=$(echo"(${b_norm[0]}*${c_norm[1]}+(-1)*${b_norm[1]}*${c_norm[0]})/$determinant"|bc-l)c_norm_inv[1]=$(echo"(${c_norm[0]}*${a_norm[1]}+(-1)*${c_norm[1]}*${a_norm[0]})/$determinant"|bc-l)c_norm_inv[2]=$(echo"(${a_norm[0]}*${b_norm[1]}+(-1)*${a_norm[1]}*${b_norm[0]})/$determinant"|bc-l)Da=$(echo"${a_norm_inv[0]}*$Dx+${b_norm_inv[0]}*$Dy+${c_norm_inv[0]}*$Dz"|bc-l)# dipole: (x,y,z) --> (a,b,c)Db=$(echo"${a_norm_inv[1]}*$Dx+${b_norm_inv[1]}*$Dy+${c_norm_inv[1]}*$Dz"|bc-l)Dc=$(echo"${a_norm_inv[2]}*$Dx+${b_norm_inv[2]}*$Dy+${c_norm_inv[2]}*$Dz"|bc-l)Pa=$(echo"$const*$Da/$V"|bc-l)# polarization: (a,b,c)Pb=$(echo"$const*$Db/$V"|bc-l)Pc=$(echo"$const*$Dc/$V"|bc-l)PQa=$(echo"$const*$La/$V"|bc-l)# polarization quanta PQPQb=$(echo"$const*$Lb/$V"|bc-l)PQc=$(echo"$const*$Lc/$V"|bc-l)dPa=$(echo"$Pa+(-1)*$Pa_tmp"|bc-l)# polarization change between two adjacent structuresdPb=$(echo"$Pb+(-1)*$Pb_tmp"|bc-l)dPc=$(echo"$Pc+(-1)*$Pc_tmp"|bc-l)dPa=$(echo"$dPa % $PQa"|bc)# constrain dP within (-PQ/2, PQ/2)if[`echo"$dPa/$PQa > 0.5"|bc-l`-eq1];thendPa=$(echo"$dPa+(-1)*$PQa"|bc-l)elif[`echo"$dPa/$PQa < -0.5"|bc-l`-eq1];thendPa=$(echo"$dPa+$PQa"|bc-l)fidPb=$(echo"$dPb % $PQb"|bc)if[`echo"$dPb/$PQb > 0.5"|bc-l`-eq1];thendPb=$(echo"$dPb+(-1)*$PQb"|bc-l)elif[`echo"$dPb/$PQb < -0.5"|bc-l`-eq1];thendPb=$(echo"$dPb+$PQb"|bc-l)fidPc=$(echo"$dPc % $PQc"|bc)if[`echo"$dPc/$PQc > 0.5"|bc-l`-eq1];thendPc=$(echo"$dPc+(-1)*$PQc"|bc-l)elif[`echo"$dPc/$PQc < -0.5"|bc-l`-eq1];thendPc=$(echo"$dPc+$PQc"|bc-l)fiPa=$(echo"$dPa+$Pa_tmp"|bc-l)Pb=$(echo"$dPb+$Pb_tmp"|bc-l)Pc=$(echo"$dPc+$Pc_tmp"|bc-l)Pa_tmp=$Pa;Pb_tmp=$Pb;Pc_tmp=$PcPx=$(echo"${a_norm[0]}*$Pa+${b_norm[0]}*$Pb+${c_norm[0]}*$Pc"|bc-l)# polarization: (a,b,c) --> (x,y,z)Py=$(echo"${a_norm[1]}*$Pa+${b_norm[1]}*$Pb+${c_norm[1]}*$Pc"|bc-l)Pz=$(echo"${a_norm[2]}*$Pa+${b_norm[2]}*$Pb+${c_norm[2]}*$Pc"|bc-l)Ptot=$(echo"sqrt($Px^2+$Py^2+$Pz^2)"|bc-l)echo$i$Pa$Pb$Pc$Ptot$PQa$PQb$PQc|awk'{printf "%10s %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n",$1,$2,$3,$4,$5,$6,$7,$8}'>>P.dat
donecatP.dat
As an example, the aligned polarization results are shown for the three polar phases of \(\mathrm{BaTiO}_3\): P4mm, Amm2, and R3m, which exhibit spontaneous polarizations along the [001], [110], and [111] directions, respectively. All polarization values are expressed in units of \(\mathrm{\mu C/cm^2}\). The components (Pa, Pb, Pc) represent the polarization projected along the lattice basis vectors (a, b, c), whereas (PQa, PQb, PQc) denote the corresponding polarization quanta along these directions. The outputs correctly reproduce the expected polarization orientations for the different polar phases. Moreover, the polarization varies monotonically along structures 1–5, where structure 1 corresponds to the nonpolar reference and structure 5 to the polar phase.
