Miscellaneous Tips & Tricks#
If certain relaxations are not converging after multiple continuation calculations (i.e. if
resubmitting certain relaxations), this is likely due to an error in the underlying calculation, extreme forces and/or
small residual forces which the structure optimisation algorithm is struggling to relax. In most of these cases,
ShakeNBreak will automatically detect and handle these calculations with
snb-run, but some tricky cases
may require manual tuning from the user:
VASP, a common culprit is
EDWAVin the output file, which can typically be avoided by reducing
KPAR. Other errors related to forces / unreasonable interatomic distances (like
ZPOTRF), as well as cases with positive energies, are automatically handled by
snb-run(folder is renamed to
Bond_Distortion_X_High_Energyand subsequently ignored).
If some relaxations are still not converging after multiple continuations, you should check the calculation output files to see if this requires fixing. Often this may require changing a specific input file setting (e.g. in the
VASP), and using the updated setting(s) for any other relaxations which are struggling to converge.
For codes other than
VASP, forces errors or high / positive energies are not automatically handled, and so in the rare cases where this occurs, you should rename the folder(s) to
ShakeNBreakwill subsequently ignore them.
If the calculation outputs show that the relaxation is proceeding fine, without any errors, just not converging to completion (i.e. residual forces), then
ShakeNBreakwill consider the calculation converged if the energy is changing by <2 meV with >50 ionic steps. Alternatively, convergence of the forces can be aided by:
Switching the ionic relaxation algorithm (e.g. change
Reducing the ionic step width (e.g. change
Switching the electronic minimisation algorithm (e.g. change
VASP), if electronic concergence seems to be causing issues.
Tightening/reducing the electronic convergence criterion (e.g. change
Rattling the structure slightly, using the
rattlefunction. An example of using this fucntion is shown in the ‘Bulk Phase Transformations’ section below).
In the other rare case where all distortions yield high energies, relative to the
Unperturbed structure, this is
typically indicative of an unreasonable defect charge state (with the extreme excess charge inducing many false local
minima on the PES).
ShakeNBreak will print a warning in these cases, with advice on how to proceed if this is
not the case and the charge state is reasonable (see below).
The default bond distortion range of -60% to +60%, and rattling standard deviation (
stdev = 10% of the bulk bond
length) are reasonable choices for most materials, typically giving best performance, but in some rare cases these may
need to be adjusted. This can be the case with extremely hard/ionic/oxide materials which typically yield larger forces
in response to bond distortion. If this issue occurs, it will manifest as:
If the rattle standard deviation is too large, it may result in high energies for each distorted & rattled structure (consistently higher energy than the unperturbed structure). As mentioned in Tricky Relaxations above,
ShakeNBreakwill print a warning in these cases, and often it is the result of unreasonable defect charge states. If the calculations have finished ok and the defect charge states are reasonable, then you likely need to reduce the rattle standard deviation (
stdev) to 7.5% of the bulk bond length (or 5% if this still causes higher energies) to avoid this – if you’re unsure of the bulk bond length for your material, just look at the previous info messages or output
ShakeNBreak, for which the default
stdevwill be equal to 10% of the bulk bond length. If this occurs, it often indicates a complex PES with multiple energy minima, thus energy-lowering distortions particularly likely, so important to test these cases with reduced
stdev! Typically the largest rattle standard deviation for which the relaxations run without issue is best for performance in terms of finding groundstate structures.
Note that strongly-correlated / magnetic materials in particular can be extremely sensitive to large structural noise, and so these typically require rattle standard deviations (
stdev) ≤ 0.05 Å.
High energies / non-converging calculations for the ±60% endpoints. As mentioned in Tricky Relaxations above, these are automatically handled by
VASPand so no changes are required, but for other codes you should rename the folder(s) to
ShakeNBreakwill subsequently ignore them.
If you are unsure but suspect this could be an issue for your material, the best strategy is typically to begin the
calculations with the default settings for one defect, then parse the results and
ShakeNBreak will warn you if
this issue is occurring – if not, all good!
For polar/ionic systems, defects are often associated with polarons (region of localised charge formed when charge
carriers interact with the lattice ions). When using
ShakeNBreak to identify these localised solutions, we recommend
Not specify initial magnetic moments for each atom (
VASP). Generally, the lattice distortion is enough to localise the charge in the defect vicinity, with the advantage of less bias as the user does not have to specify the atoms where the charge is localised. If the defect is surrounded by two types of ions, and the charge is expected to preferentially localise in one of them, you can use the
neighbour_elementskeyword to only distort the specified element (see next section).
Specify the total magnetic moment (
VASP). We recommend to use a wider distortion mesh (
delta = 0.2) and test the main
NUPDOWNpossibilities, e.g. if there are two extra/missing electrons, trial
NUPDOWN = 0(anti-parallel arrangement) and
NUPDOWN = 2(parallel arrangement).
If your host system is magnetic (e.g. (anti-)ferromagnetic), you should initialise the magnetic moments
(i.e. set the
MAGMOM values in the
INCAR files using
Distortions.write_vasp_files()) to match this.
In these cases, the defect ground-state often involves a polaronic localisation which alters the local
Often the defect formation breaks the (magnetic) symmetry of the host, and so the choice of spin up/down
for an anti-ferromagnetic (AFM) system is no longer arbitrary within a
VASP calculation, as the
initialisation favours spin-up as the majority spin of the system (by only allowing positive
values). To account for this magnetic symmetry-breaking and thus differing localisation solutions for
different choices of AFM orderings, it is recommended to test both equivalent‡ AFM orderings for
each potential polaronic defect state (i.e. when
NUPDOWN is not zero).
This can be done by ‘flipping’ the
MAGMOM values (to flip the AFM ordering) by changing the sign of
MAGMOM value for the supercell. If doing this, it is recommended to use a wider distortion mesh
delta = 0.2) to reduce the number of calculations required, as the testing of both AFM orderings
should allow sufficient coverage of the defect PES with this reduced sampling.
‡ = equivalent for the pristine host system, but not for the defect supercell with
NUPDOWN ≠ 0.
neighbour_elements Use Cases#
When generating atomic distortions with
neighbour_elements optional parameter can be
particularly useful in certain cases if:
You have a complex multi-cation / multi-anion system, and believe that the most likely distorting species about certain defect sites are not the nearest neighbour atoms. For example, in a rare case you might have two cations (A and B), where the nearest neighbours of cation A are cations B, but it is the (second-nearest-neighbour) anions which are likely most prone to rearrange upon formation of a A cation vacancy.
You have ‘spectator’ ions (e.g. the A-site cation in ABX3 perovskites) that are nearest neighbours to the defect, but unlikely to distort or rebond. This has been seen in studies in our research groups (reference to be added when preprinted).
Bulk Phase Transformations#
If you perform
ShakeNBreak calculations with a supercell structure for which a lower energy polymorph exists
(i.e. by using a bulk structure which has imaginary phonon modes), often the symmetry-breaking introduced by
ShakeNBreak can cause your defect supercell to relax to this lower energy bulk structure. If this is the case,
it will typically cause all distortion calculations to yield a significantly lower energy structure than the
Unperturbed relaxation, for most (if not all) defects calculated; as any significant symmetry-breaking in the
supercell is allowing the structure to rearrange to the lower energy polymorph.
ShakeNBreak will print a cautionary
warning in these cases, and one way to manually check this is to visually compare the relaxed structure from one of the
low energy distortion calculations, with the relaxed
Unperturbed structure, and the bulk supercell, and see how the
regions away from the defect site compare.
Often this is useful information, as it may reveal a previously-unknown low-energy polymorph for your host system.
However, it also means that your original higher energy bulk structure is no longer an appropriate reference structure
for calculating your final defect formation energies, and so you should instead obtain the bulk supercell corresponding
to this lower energy polymorph, and use this as the reference structure for calculating defect formation energies.
The recommended workflow for doing this is to firstly try obtaining this lower energy polymorph in the defect-free host
system using similar atom rattling to that used in the
ShakeNBreak distorted structure generation:
from shakenbreak.distortions import rattle, Structure bulk = Structure.from_file("path/to/Bulk_Supercell_POSCAR") rattled_supercell = rattle(bulk) rattled_supercell.to(fmt="POSCAR", filename="Rattled_Bulk_Supercell_POSCAR")
Then calculate the energy of this rattled supercell and compare to the original high-symmetry bulk supercell.
If it’s significantly lower energy (similar to the energy difference between your
Unperturbed and distorted defect
relaxations), then this is likely the lower energy polymorph for your host system.
If this does not yield a significantly lower energy polymorph, then it’s recommended to calculate the phonon dispersion of your host material, and check if there are any imaginary phonon modes (indicating the presence of a nearby lower-symmetry lower-energy polymorph). If this is the case, then you can try to obtain this lower energy polymorph using a code like ModeMap or similar, to generate the distorted structure corresponding to this imaginary mode.
This workflow also serves to explicitly test if indeed a phase transformation is occurring in your defect supercell(s).
If this does indeed reveal a significantly lower energy polymorph for your host material, depending on how different this
structure is, you might want to regenerate the defects in this new supercell and possibly re-run
these – particularly for any defects that did not show an energy-lowering relative to
Unperturbed in the original
supercell (suggesting that they remained ‘stuck’ in the original higher-energy arrangement).
In certain cases, it’s possible that you witness behaviour similar to the above, despite no lower energy polymorph
existing for your host material. This can be the case in certain high-symmetry materials (e.g. we have witnessed this
in incipient ferroelectrics), where defects introduce strong local symmetry-breaking (i.e. a ‘local phase
transformation’) which is missed by the standard approach with unperturbed defect relaxations. These are cases where
ShakeNBreak is particularly important for your material, you can continue your calculations as normal, ignoring
ShakeNBreak is optimised to work with isolated point defects. However, it can also be used with
complex defects (and has been found to be important in these cases as well, e.g. this Chem Sci paper), but requires some
workarounds as the
pymatgen defect functions are not natively built for this.
This involves generating the defect complex as a
Defect object using one of the point
defects as the ‘bulk’ structure and the other as the ‘defect’, then feeding this to
ShakeNBreak in order to
generate the distortions. If you are trying to do this and are running into issues, you can contact the developers and
we can share some guidance for this (until improved
pymatgen-based functionality comes about for complex defects).
ShakeNBreak workflow is primarily geared toward ground-state structure identification, it can also be
applicable to finding metastable states, as described in the method paper.
For this, you can use the optional
metastable argument for
see docs here.
A current known issue with
pymatgenis that it might give an error similar to this:
ValueError: numpy.ndarray size changed, may indicate binary incompatibility. Expected 88 from C header, got 80 from PyObject
This is due to a recent change in the
numpyC API, see here. It should be fixed by reinstalling
pymatgen, so that it is rebuilt with the new
pip uninstall pymatgen pip install pymatgen
Have any tips for users from using ShakeNBreak? Please share it with the developers and we’ll add them here!