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 ($P_s$) due to lack of inversion symmetry but whether $P_s$ 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 $P_s$ is the coercive field ($E_c$). While quantitative prediction of coercive field quires larger scale simulations of domain wall motions, polarization switching barrier ($\omega_s$) at atomic scale can provide qualitative insights since $\omega_s$describe the required energy density to change $P_s$. $E_D$

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.

A general workflow to calculate polarization switching barrier

Calculating polarization switching barrier ($\omega_s$) can be broken into two general components:

1. Finding the polar (+$P_s$) and antipolar (-$P_s$) structures of a given material (with polar point group)
2. Determined the minimum energy path (MEP) between them

Once the MEP is determined, the largest barrier along the path is the $\omega_s$.

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

1. Fully relax (ionic position, cell size, and cell shape) the end images

2. Use linear interpolation (cell as well as ionic positions) between end images to get the initial intermediate images

3. 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.

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]

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.

Minimum energy path between polar and antipolar structure in $\mathrm{Al}_{1-x}\mathrm{Sc}_x\mathrm{N}$alloys. 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]