ShakeNBreak for Polarons#
ShakeNBreak applied to polarons in CuSbSe₂#
In this tutorial, we use ShakeNBreak (SnB) to apply bond distortions and rattling to aid location of small polaron geometries in CuSbSe₂ (a candidate solar cell material). In this example case, we will:
Apply rattling and bond distortions to atoms to break symmetry and favour polaronic localisation.
Add an electron/hole (using
NELECT) and an unpaired spin (NUPDOWN=1) to the system.Generate a set of
VASPinput files to relax and post-process.
Reference:
Fu‡ & Lohan‡ et al. Structural and electronic features enabling delocalized charge-carriers in CuSbSe₂ Nature Communications 2025
Note
Searching for stable distorted polaron geometries in this way is very similar to the original bond distortion method introduced by Deskins et al. (upon which much of the defect structure-searching methodology in ShakeNBreak was inspired), with some additional benefits from partially-constrained Monte Carlo atomic rattling to break symmetry and further favour lower-energy localised solutions.
Please cite the original bond distortion method if using this approach to search for polaron structures!
Rationale for bond distortion workflow with polarons#
Polarons, much like point defects, often result in atom-centred structural distortions along with a change in the local charge density. As there may be a finite energy barrier between the undistorted free-electron case and distorted polaronic groundstate, it is usually necessary to perturb the system in some way before it will fall into a localised groundstate (if one exists).
from pymatgen.core.structure import Structure
bulk_supercell = Structure.from_file("CuSbSe2_POSCAR") # load our bulk supercell structure
Simple Approach: Rattling#
For the most simple first step, we can apply atomic rattling (random perturbation of all atomic positions using a bond-length-dependent Gaussian distribution with a Monte Carlo algorithm; see docstring) to break the symmetry of the system and disrupt the long-range lattice potential, which often can be sufficient to identify a localised polaronic groundstate – particularly in cases of strong trapping.
from shakenbreak.distortions import rattle
rattled_supercell = rattle(bulk_supercell)
# note there are many options we can use with ``rattle()`` to alter the rattling process
# (see https://shakenbreak.readthedocs.io/en/latest/shakenbreak.distortions.html#shakenbreak.distortions.rattle; e.g. ``stdev``)
# but the default settings are usually good for the vast majority of cases
Bond Distortions#
In many cases however, rattling will not be sufficient on its own to locate polaronic geometries, and so we can apply targeted bond distortions to specific lattice sites to encourage localisation, using the distortions code in ShakeNBreak.
First, let’s determine the different symmetry-inequivalent sites in our material. We can do this using the SpacegroupAnalyzer class in pymatgen:
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
sga = SpacegroupAnalyzer(bulk_supercell)
sga.get_symmetrized_structure() # 1 symmetry-inequivalent Cu site, 1 Sb site, 2 Se sites
SymmetrizedStructure
Full Formula (Cu16 Sb16 Se32)
Reduced Formula: CuSbSe2
Spacegroup: Pnma (62)
abc : 10.334277 10.334277 16.265031
angles: 75.639685 75.639685 102.667344
Sites (64)
# SP a b c Wyckoff
--- ---- -------- -------- -------- ---------
0 Cu 0.603183 0.978183 0.670123 16c
1 Sb 0.727246 0.102246 0.442623 16c
2 Se 0.287447 0.412448 0.175795 16c
3 Se 0.047292 0.422292 0.403673 16c
Or we can quickly use doped to do this as well, which generates vacancies/antisites at each inequivalent site (and thus shows 1 symmetry-inequivalent Cu site, 1 Sb site and 2 Se sites also):
from doped.generation import DefectsGenerator
defect_gen = DefectsGenerator( # no interstitials, don't generate supercell
bulk_supercell, interstitial_gen_kwargs=False, generate_supercell=False
)
Generating DefectEntry objects: 100.0%|██████████| [00:04, 21.82it/s]
Vacancies Guessed Charges Conv. Cell Coords Wyckoff
-------------- ----------------- ------------------- ---------
v_Cu [+1,0,-1] [0.248,0.250,0.170] 4c
v_Sb [+1,0,-1,-2,-3] [0.272,0.750,0.443] 4c
v_Se_Cs_Cu2.41 [+2,+1,0,-1] [0.376,0.250,0.324] 4c
v_Se_Cs_Cu2.44 [+2,+1,0,-1] [0.127,0.250,0.596] 4c
Substitutions Guessed Charges Conv. Cell Coords Wyckoff
--------------- --------------------------- ------------------- ---------
Cu_Sb [0,-1,-2] [0.272,0.750,0.443] 4c
Cu_Se_Cs_Cu2.41 [+4,+3,+2,+1,0] [0.376,0.250,0.324] 4c
Cu_Se_Cs_Cu2.44 [+4,+3,+2,+1,0] [0.127,0.250,0.596] 4c
Sb_Cu [+4,+3,+2,+1,0,-1,-2,-3,-4] [0.248,0.250,0.170] 4c
Sb_Se_Cs_Cu2.41 [+7,+6,+5,+4,+3,+2,+1,0,-1] [0.376,0.250,0.324] 4c
Sb_Se_Cs_Cu2.44 [+7,+6,+5,+4,+3,+2,+1,0,-1] [0.127,0.250,0.596] 4c
Se_Cu [+3,+2,+1,0,-1,-2,-3] [0.248,0.250,0.170] 4c
Se_Sb [+1,0,-1,-2,-3,-4,-5] [0.272,0.750,0.443] 4c
The number in the Wyckoff label is the site multiplicity/degeneracy of that defect in the conventional ('conv.') unit cell, which comprises 4 formula unit(s) of CuSbSe2.
Hole/electron polarons will typically localise in the orbitals which make up the VBM/CBM of the material, as these are (initially) the lowest energy states available for these charge carriers.
In CuSbSe₂, the VBM is primarily composed of Se p and Cu d orbitals, so we will target Se-centred and Cu-centred distortions for a hole polaron example here.
First we get the atomic site we’re going to distort around:
Cu_site_idx = bulk_supercell.sites.index(defect_gen["v_Cu_0"].defect_supercell_site)
Then we can generate the distorted structures. Let’s try +/-40% bond length distortions (i.e. expansion and contraction) around this lattice site:
from shakenbreak.distortions import distort_and_rattle
trial_hole_polaron_supercells = {f"Cu_{distortion_factor:+.1%}": distort_and_rattle(
structure=bulk_supercell,
site_index=Cu_site_idx,
num_nearest_neighbours=4, # let's distort the 4 nearest neighbours
distortion_factor=distortion_factor,
)["distorted_structure"] for distortion_factor in [-0.4, 0.4]}
And for Se-centred distortions:
# Se sites: v_Se_Cs_Cu2.41 & v_Se_Cs_Cu2.44
Se_Cs_Cu2pt41_site_idx = bulk_supercell.sites.index(defect_gen["v_Se_Cs_Cu2.41_0"].defect_supercell_site)
Se_Cs_Cu2pt44_site_idx = bulk_supercell.sites.index(defect_gen["v_Se_Cs_Cu2.44_0"].defect_supercell_site)
trial_hole_polaron_supercells.update({f"Se_Cs_Cu2.41_{distortion_factor:+.1%}": distort_and_rattle(
structure=bulk_supercell,
site_index=Se_Cs_Cu2pt41_site_idx,
num_nearest_neighbours=4, # let's distort the 4 nearest neighbours
distortion_factor=distortion_factor,
)["distorted_structure"] for distortion_factor in [-0.4, 0.4]})
trial_hole_polaron_supercells.update({f"Se_Cs_Cu2.44_{distortion_factor:+.1%}": distort_and_rattle(
structure=bulk_supercell,
site_index=Se_Cs_Cu2pt44_site_idx,
num_nearest_neighbours=4, # let's distort the 4 nearest neighbours
distortion_factor=distortion_factor,
)["distorted_structure"] for distortion_factor in [-0.4, 0.4]})
We can then write these structures to file to run VASP polaron-searching calculations on them:
for name, structure in trial_hole_polaron_supercells.items():
structure.to(fmt="POSCAR", filename=f"{name}_POSCAR")
!ls *_POSCAR # check the POSCARs have been written
CuSbSe2_POSCAR Se_Cs_Cu2.41_-40.0%_POSCAR
Cu_+40.0%_POSCAR Se_Cs_Cu2.44_+40.0%_POSCAR
Cu_-40.0%_POSCAR Se_Cs_Cu2.44_-40.0%_POSCAR
Se_Cs_Cu2.41_+40.0%_POSCAR
Input Files#
With our trial distorted supercells, we can then add an electron/hole to the system by increasing/decreasing NELECT in the INCAR by 1 (relative to a neutral cell calculation), and setting NUPDOWN=1 to constrain the unpaired spin.
The DopedDictSet/DefectDictSet classes from doped.vasp can be used to do this through python, or of course you can manually generate & edit INCAR, KPOINTS and POTCAR files for this.
The NELECT value for a neutral cell calculation can be obtained from the OUTCAR of a previous neutral cell / bulk calculation (or by momentarily running VASP without specifying NELECT and looking for NELECT in the OUTCAR), or by summing the number of valence electrons of each element in the POTCAR times their atom count in the supercell.
We would then run VASP relaxations for these structures, to see if localised polarons are obtained. For this, we would typically compare the final energies to that of a single-shot calculation (no ionic relaxation) of the same supercell with the same NELECT (which should give the energy of a delocalised carrier in this supercell), with this energy difference giving an estimate of the polaron binding energy. Of course, the binding energy calculated this way can be sensitive to the supercell size, and so it is recommended to either perform these calculations over multiple supercell sizes and check their convergence / extrapolate to the value at infinite cell size (e.g. see Figs S14,S15 in https://pubs.acs.org/doi/suppl/10.1021/acs.jpclett.2c02436/suppl_file/jz2c02436_si_001.pdf), or employ a charge correction scheme for the localised polaron case (such as the FNV scheme implemented in doped’s codebase or sxdefectalign or that of Falletta et al.).
Further Analysis#
After comparing the energies of these relaxed supercells, you may want to plot the energies vs distortion and perform further analysis of the electronic densities (e.g. by plotting orbital densities with PARCHG outputs, for which the Wavecar.get_parchg function from pymatgen is useful, by plotting magnetisation densities using CHGCAR files (e.g. in VESTA or CrystalMaker), by analysing the electronic structure (e.g. using get_eigenvalue_analysis in doped) etc) – e.g. see Fig. S10 in Structural and electronic features enabling delocalized charge-carriers in CuSbSe₂.
The hole-finder.py script provided in the ShakeNBreak docs directory can be useful for identifying hole states in VASP calculations, for PARCHG analysis.