Calculating Polarization Switching Barrier using SS-NEB
This tutorial is based on VASP compiled with VTST package
Written by Cheng-Wei Lee (clee2 [at] mines [dot] edu)
Ferroelectric or not?
Pyroelectric materials by definition have spontaneous polarization (Ps) due to lack of inversion symmetry but whether Ps can be experimentally reversed by an external electric field before dielectric breakdown happens determines if a material is ferroelectric. The critical electric field to reverse Ps is the coercive field (Ec). While quantitative prediction of coercive field quires larger scale simulations of domain wall motions, polarization switching barrier (ωs) at atomic scale can provide qualitative insights since ωsdescribe the required energy density to change Ps. ED
A general workflow to calculate polarization switching barrier
Determined the minimum energy path (MEP) between them
Finding antipolar structure
We can determine the antipolar structure once the nonpolar structure is determined. The nonpolar parent structure can be determined using Group-super group relations[1] and one examplar online tool to find the pseudosymmetry is PSEUDO. Here, a python script to find antipolar structure in wurtzite-type is provided. The script assumes the layered hexagonal structure as the nonpolar structure for wurtzite-type polar structures.
from pylada.crystal import read,write,neighbors
import numpy as np
def myneighbors(strc=None, howmany=12, pos=None, tol=0.2):
"""
A function that replaces the original pylada
neighbors function that seems to be broken.
It still uses the original one, just makes sure
the output is correct. Also, it returns only the
first coordination shell.
The tolerance is in fraction of the distance to
the closest atom and is 20% by default (test first)
"""
nghs=[]
# Loop ove pylada neighbors (arbitrary small tolerance)
for ng in neighbors(strc,howmany,pos,0.1):
# If empty just add ng and the first one should have the shortest distance
if len(nghs)==0:
nghs.append(ng)
continue
# Make sure the distance is within the tolerance
elif nghs[0][-1]<=ng[-1]<=(1.+tol)*nghs[0][-1]:
nghs.append(ng)
return nghs
'''
This script assumes that it is the cations that are moving during polarization switching
'''
# Anion list; make sure all the anions in the dataset are added
anions = ['N','O','S','Se']
# read polar structure
strc = read.poscar('POSCAR_p')
# find polar direction
### first find all the cation-to-anion vectors
vectot = []
for atom in strc:
if atom.type not in anions:
nghs = myneighbors(strc,4,atom.pos,0.2) # based on our knowledge for wurtzite-like structure
for ngh in nghs:
vectot.append(ngh[1])
### The average should reveal the polar direction for wurtzite-like structure
tmp_dir = np.average(vectot,axis=0)
### Get rid of noise
polar_dir = np.round(tmp_dir/max(tmp_dir,key=abs),0)
#print(polar_dir)
shifts=[]
# find distance from the hexagonal plane along polar direction
for atom in strc:
if atom.type not in anions: # only look at cations
nghs = myneighbors(strc,4,atom.pos,0.2)
vecs = []
for ngh in nghs:
if np.linalg.norm(ngh[1]*polar_dir/np.linalg.norm(ngh[1])) < 0.9: # i.e. not in polarization direction, only the neighbors in the polar direction can get large number
vecs.append(ngh[1])
assert len(vecs) == 3, "The code only works for wurtzite like structure" # assert function to catch unexpected situtations
shifts.append(np.average(vecs,axis=0)) # the average of distance vectors from anion to cation in non-polar direction gives the distance from cation to hexagonal plane
shift = np.average(shifts,axis=0) # get the average value out of all the cations
for atom in strc:
if atom.type not in anions: # only move the cations
atom.pos = atom.pos + shift*2.0 # factor of 2.0 is to get to antipolar structure; factor of 1.0 is roughly hexagonal structure
with open("POSCAR_a", "w") as file: write.poscar(strc, file, vasp5=True)// Some code
Calculating MEP using SS-NEB
Once the end images (polar and antipolar structures) are determined and created, here are the steps to perform the SS-NEB calculations
Fully relax (ionic position, cell size, and cell shape) the end images
Use linear interpolation (cell as well as ionic positions) between end images to get the initial intermediate images
Run SS-NEB calculations till convergence is reached.
nebmake.pl (provided within the VTST package) can be used to do the linear interpolation
Here's an exemplar VASP input file for the SS-NEB calculations:
ADDGRID = .TRUE.
ISPIN = 2 # depends on the system of interest
ISTART = 1
ICHARG = 0
ISYM = 0
LMAXMIX = 4
LVHAR = .FALSE.
LORBIT = 10
LWAVE = .TRUE.
LCHARG = .TRUE.
LVTOT = .FALSE.
LPEAD = .FALSE.
LPLANE = .TRUE.
NELMIN = 5
NPAR = 3 # depends on the total number of requested cores and number of images
NELECT = 288.0 # for a supercell with charged defects
ALGO = Normal
EDIFF = 1e-06
EDIFFG = -0.01
ENCUT = 340.0
PREC = Accurate
NSW = 10 # small number since we want to restart the calculations to have updated FFT grid
IBRION = 3
ISIF = 3
ISMEAR = 0
SIGMA = 0.01
POTIM = 0
ICHAIN = 0
IOPT = 3
SPRING = -5 #default value generally works fine
LCLIMB = .FALSE.
LNEBCELL = .TRUE. # True is for SS-NEB; False is for NEB
IMAGES = 36 # number of intermediate images
A large number of intermediate images are usually needed to find the SS-NEB pathway between polar and antipolar structures of multinary wurtzite-type materials. Since the current VTST implementation is limited to 99 intermediate images and to speed up the process, segmentation of the pathway is generally needed (see below). We recommend to have 2-3x the number cations within the simulation cell as the initial number of intermediate images and run roughly 100-200 SS-NEB steps to capture the the unconverged MPE. We then separate the path into multiple segments based on the local minima and fully relaxed them to be the new end images. Further segmentation may be needed and we converge each segment independently and connect them at the end.
Calculating polarization switching barrier
Once the minimum energy path is determined by the SS-NEB method, we extract the energy density with respect to the polar structure. It is calculated by the energy difference from the polar structure and normalized by the volume of the polar structure. The polarization barrier is then determined by the largest energy density peak along the switching pathway (see the image below)
While the unit of eV/f.u. is commonly used in the literature, it is not a good choice when comparing materials with different number of cations[2] since the choice leads to spuriously lower barrier for materials with higher number of cations.
Calculating polarization switching barrier (ωs) can be broken into two general components:
Finding the polar (+Ps) and antipolar (-Ps) structures of a given material (with polar point group)
Once the MEP is determined, the largest barrier along the path is the ωs.
Schematics for polarization-electric field (P-E) curve and polarization switching barrier. Ps and Pr are spontaneous and remnant polarizations, respectively. Ec is the coercive field.
Segmentations of the MEP between polar and nonpolar structures. Structures at Local minima are fully relaxed and are the end images of each separate SS-NEB calculations. The figure was taken from [3]
Minimum energy path between polar and antipolar structure in Al1−xScxNalloys. The barrier is represented by the first peak for x=0.36. There is only one peak for x=0.22. Figure is taken from [3]
