NEB / CC Diagram Path Generation
Here we will generate the a set of linearly-interpolated structures along the transformation path from one defect to another (specifically the -1 and -2 charged selenium vacancies in trigonal selenium), which can be used to calculate the configuration coordinate (CC) diagram (i.e. potential energy surfaces of the two defect states) or as the initial guess for Nudged Elastic Band (NEB) calculations.
CC diagrams can then be used to compute non-radiative recombination rates with nonrad and/or CarrierCapture.jl (Alkauskas et al. PRB 2014), to predict optical transitions and their luminescence lineshapes (Alkauskas et al. PRL 2012), to predict radiative recombination rates, and more.
Tip
You can run this notebook interactively through Google Colab or Binder – just click the links!
If running on Colab, then you’ll need to run !pip install doped in a cell to install the package, and !git clone https://github.com/SMTG-Bham/doped to download the example data (and update paths in the code cells accordingly).
CCD for Point Defects
Here we will generate the linearly-interpolated path between the -1 and -2 charged selenium vacancies in trigonal selenium, as was done to generate the CC diagram and compute their electron/hole charge capture rates (and thus contribution to non-radiative electron-hole recombination) in this study.
from doped.thermodynamics import DefectThermodynamics
Se_intrinsic_thermo = DefectThermodynamics.from_json("Se/Se_Intrinsic_Thermo.json.gz")
V_Se_m1_supercell = Se_intrinsic_thermo["vac_1_Se_-1"].defect_supercell
V_Se_m2_supercell = Se_intrinsic_thermo["vac_1_Se_-2"].defect_supercell
Something that we need to be careful with when obtaining these interpolated structures, for both CC diagrams and NEB paths, is that they should be appropriately oriented relative to eachother, so that they are actually along the closest linear path between them. E.g. for a vacancy, the two vacancy supercells should have the vacancy in the same position in the supercell (so that we don’t unintentionally model a long migration pathway!), and rotated so that the local environments most closely match eachother (essentially minimising the root-mean-square displacement between the geometries). A rough sketch of this point, showing the correct symmetry-equivalent orientation of \(V^0\) to give the closest linear path to the \(V^{+1}\) structure is shown below:
To ensure we have the correct matching orientation, we can use the orient_s2_like_s1 function, which works similarly to the StructureMatcher.get_s2_like_s1 method in pymatgen, but accelerated and extended to ensure matching atomic indices and lattice vector definitions.
from doped.utils.configurations import orient_s2_like_s1
# set verbose to True to show the difference in RMS displacement (assuming matched site indices)
# between the original and re-oriented structures
V_Se_m2_like_m1 = orient_s2_like_s1(V_Se_m1_supercell, V_Se_m2_supercell, verbose=True)
# if we don't set verbose to True, there is no printed output here, just the re-oriented structure is returned
ΔQ(s1/s2) = 9.91 amu^(1/2)Å
ΔQ(s2_like_s1/s2) = 4.88 amu^(1/2)Å
ΔQ(s1/s2_like_s1) = 8.63 amu^(1/2)Å
Here we can see that the RMS displacement (assuming matched site indices) between the \(V_{Se}^{-1}\) and \(V_{Se}^{-2}\) structures has been reduced by reorienting (but not changing) the \(V_{Se}^{-2}\) geometry, now corresponding to the closest linear path between symmetry-equivalent geometries of the two defects. We can confirm this by checking the structure symmetry and minimum interatomic distance before and after re-orientation:
from doped.utils.symmetry import point_symmetry_from_structure
from doped.utils.supercells import min_dist
# get the point group symmetry and minimum interatomic distances of the structures:
V_Se_m1_symm = point_symmetry_from_structure(V_Se_m1_supercell)
V_Se_m1_min_dist = min_dist(V_Se_m1_supercell)
V_Se_m2_symm = point_symmetry_from_structure(V_Se_m2_supercell)
V_Se_m2_min_dist = min_dist(V_Se_m2_supercell)
V_Se_m2_like_m1_symm = point_symmetry_from_structure(V_Se_m2_like_m1)
V_Se_m2_like_m1_min_dist = min_dist(V_Se_m2_like_m1)
for name, symm, min_dist in zip(
["V_Se^-1", "V_Se^-2 Original", "V_Se^-2 Reoriented"],
[V_Se_m1_symm, V_Se_m2_symm, V_Se_m2_like_m1_symm],
[V_Se_m1_min_dist, V_Se_m2_min_dist, V_Se_m2_like_m1_min_dist]
):
print(f"{name} point group symmetry: {symm}, minimum interatomic distance: {min_dist:.2f} Å")
V_Se^-1 point group symmetry: C2, minimum interatomic distance: 2.30 Å
V_Se^-2 Original point group symmetry: C2, minimum interatomic distance: 2.33 Å
V_Se^-2 Reoriented point group symmetry: C2, minimum interatomic distance: 2.33 Å
Now we choose our set of structures along the linearly-interpolated path between \(V_{Se}^{-1}\) and \(V_{Se}^{-2}\) to compute, by specifying the fractional displacements along this path to generate the structures for. Here we will generate structures at the fractional displacements [-1.5, -1.2, -1.0, -0.8, -0.6, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0, 1.2, 1.5], in units of ΔQ (the displacement from structure to structure 2):
import numpy as np
positive_displacements = np.array([0.1, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0, 1.2, 1.5])
displacements = np.concatenate([-positive_displacements[::-1], np.array([0.0]), positive_displacements])
print(displacements)
[-1.5 -1.2 -1. -0.8 -0.6 -0.4 -0.3 -0.2 -0.1 0. 0.1 0.2 0.3 0.4
0.6 0.8 1. 1.2 1.5]
Optional: We can generate and inspect the interpolated structures with get_path_structures:
from doped.utils.configurations import get_path_structures
disp_dict_m1, disp_dict_m2 = get_path_structures(
V_Se_m1_supercell, V_Se_m2_like_m1, displacements=displacements)
disp_dict_m1.keys()
dict_keys(['delQ_-1.5', 'delQ_-1.2', 'delQ_-1.0', 'delQ_-0.8', 'delQ_-0.6', 'delQ_-0.4', 'delQ_-0.3', 'delQ_-0.2', 'delQ_-0.1', 'delQ_0.0', 'delQ_0.1', 'delQ_0.2', 'delQ_0.3', 'delQ_0.4', 'delQ_0.6', 'delQ_0.8', 'delQ_1.0', 'delQ_1.2', 'delQ_1.5'])
disp_dict_m1["delQ_0.6"]
Structure Summary
Lattice
abc : 13.032718 13.03271811400469 14.890713999999997
angles : 90.0 90.0 119.99999971063251
volume : 2190.3633066388306
A : np.float64(13.032718) np.float64(0.0) np.float64(0.0)
B : np.float64(-6.516359) np.float64(11.286664999999998) np.float64(0.0)
C : np.float64(0.0) np.float64(0.0) np.float64(14.890713999999997)
pbc : True True True
PeriodicSite: Se (0.981, 0.04976, 6.643) [0.07748, 0.004408, 0.4461]
PeriodicSite: Se (-5.543, 11.24, 11.57) [0.07242, 0.9955, 0.7772]
PeriodicSite: Se (-1.25, 3.787, 1.663) [0.0719, 0.3356, 0.1117]
PeriodicSite: Se (-1.259, 3.769, 6.614) [0.07039, 0.334, 0.4442]
PeriodicSite: Se (-1.281, 3.837, 11.59) [0.07167, 0.3399, 0.7782]
PeriodicSite: Se (-3.457, 7.603, 1.674) [0.07153, 0.6736, 0.1124]
PeriodicSite: Se (-3.435, 7.61, 6.606) [0.07351, 0.6742, 0.4436]
PeriodicSite: Se (-3.512, 7.665, 11.63) [0.07007, 0.6791, 0.781]
PeriodicSite: Se (-1.138, 11.28, 1.649) [0.4122, 0.9991, 0.1108]
PeriodicSite: Se (5.398, -0.01405, 6.581) [0.4135, -0.001245, 0.442]
PeriodicSite: Se (-1.122, 11.29, 11.61) [0.4139, 1.0, 0.7798]
PeriodicSite: Se (3.053, 3.692, 1.63) [0.3978, 0.3271, 0.1095]
PeriodicSite: Se (2.997, 3.631, 6.566) [0.3908, 0.3217, 0.4409]
PeriodicSite: Se (3.08, 3.69, 11.59) [0.3998, 0.3269, 0.7783]
PeriodicSite: Se (0.9517, 7.493, 1.675) [0.405, 0.6639, 0.1125]
PeriodicSite: Se (0.9511, 7.48, 6.632) [0.4043, 0.6627, 0.4454]
PeriodicSite: Se (0.9462, 7.515, 11.62) [0.4055, 0.6659, 0.7805]
PeriodicSite: Se (10.14, 0.00654, 1.654) [0.7783, 0.0005794, 0.1111]
PeriodicSite: Se (3.302, 11.26, 6.612) [0.752, 0.9973, 0.4441]
PeriodicSite: Se (9.817, 0.0377, 11.58) [0.7549, 0.00334, 0.7779]
PeriodicSite: Se (7.469, 3.795, 1.649) [0.7412, 0.3362, 0.1107]
PeriodicSite: Se (7.462, 3.782, 6.592) [0.7401, 0.3351, 0.4427]
PeriodicSite: Se (7.469, 3.809, 11.58) [0.7419, 0.3375, 0.7777]
PeriodicSite: Se (5.265, 7.489, 1.637) [0.7358, 0.6636, 0.11]
PeriodicSite: Se (5.227, 7.438, 6.601) [0.7306, 0.659, 0.4433]
PeriodicSite: Se (5.254, 7.494, 11.57) [0.7351, 0.664, 0.7771]
PeriodicSite: Se (-0.2509, 0.4703, 3.063) [0.001582, 0.04167, 0.2057]
PeriodicSite: Se (-0.4911, 0.8559, 8.296) [0.0002323, 0.07583, 0.5571]
PeriodicSite: Se (-0.4503, 0.8526, 13.25) [0.003222, 0.07554, 0.8898]
PeriodicSite: Se (10.32, 4.613, 3.316) [0.9963, 0.4087, 0.2227]
PeriodicSite: Se (10.32, 4.661, 8.252) [0.9986, 0.413, 0.5542]
PeriodicSite: Se (10.3, 4.668, 13.25) [0.9975, 0.4136, 0.8898]
PeriodicSite: Se (-4.919, 8.528, 3.301) [0.0003415, 0.7556, 0.2217]
PeriodicSite: Se (-4.888, 8.479, 8.262) [0.0005476, 0.7512, 0.5549]
PeriodicSite: Se (8.046, 8.5, 13.3) [0.9939, 0.7531, 0.8929]
PeriodicSite: Se (3.873, 0.822, 3.264) [0.3336, 0.07283, 0.2192]
PeriodicSite: Se (3.97, 0.8335, 8.258) [0.3416, 0.07385, 0.5546]
PeriodicSite: Se (3.889, 0.824, 13.24) [0.3349, 0.07301, 0.8894]
PeriodicSite: Se (1.601, 4.502, 3.288) [0.3223, 0.3989, 0.2208]
PeriodicSite: Se (1.615, 4.483, 8.266) [0.3225, 0.3972, 0.5551]
PeriodicSite: Se (1.62, 4.468, 13.26) [0.3222, 0.3958, 0.8905]
PeriodicSite: Se (-0.4883, 8.354, 3.332) [0.3326, 0.7402, 0.2237]
PeriodicSite: Se (-0.4735, 8.372, 8.29) [0.3346, 0.7418, 0.5567]
PeriodicSite: Se (-0.4934, 8.376, 13.28) [0.3332, 0.7421, 0.8919]
PeriodicSite: Se (8.469, 0.7873, 3.256) [0.6847, 0.06976, 0.2187]
PeriodicSite: Se (8.318, 0.8128, 8.252) [0.6743, 0.07202, 0.5541]
PeriodicSite: Se (8.355, 0.8235, 13.27) [0.6776, 0.07297, 0.8911]
PeriodicSite: Se (5.987, 4.608, 3.296) [0.6635, 0.4083, 0.2214]
PeriodicSite: Se (5.984, 4.585, 8.253) [0.6623, 0.4062, 0.5542]
PeriodicSite: Se (5.986, 4.617, 13.23) [0.6639, 0.4091, 0.8886]
PeriodicSite: Se (3.843, 8.346, 3.305) [0.6646, 0.7395, 0.2219]
PeriodicSite: Se (3.807, 8.305, 8.269) [0.66, 0.7359, 0.5553]
PeriodicSite: Se (3.852, 8.388, 13.25) [0.6672, 0.7432, 0.8896]
PeriodicSite: Se (1.598, 2.759, 0.00044) [0.2448, 0.2445, 2.955e-05]
PeriodicSite: Se (1.535, 2.768, 4.899) [0.2404, 0.2452, 0.329]
PeriodicSite: Se (1.635, 2.802, 9.936) [0.2495, 0.2483, 0.6673]
PeriodicSite: Se (-0.5257, 6.676, 0.0256) [0.2554, 0.5915, 0.001719]
PeriodicSite: Se (-0.5288, 6.666, 4.981) [0.2547, 0.5906, 0.3345]
PeriodicSite: Se (-0.5287, 6.704, 9.962) [0.2564, 0.594, 0.669]
PeriodicSite: Se (-2.661, 10.47, 0.0375) [0.2596, 0.9277, 0.002518]
PeriodicSite: Se (-2.636, 10.47, 4.95) [0.2616, 0.9276, 0.3324]
PeriodicSite: Se (-2.561, 10.46, 9.934) [0.2669, 0.9269, 0.6671]
PeriodicSite: Se (6.033, 2.929, 14.88) [0.5927, 0.2595, 0.9994]
PeriodicSite: Se (6.026, 2.912, 4.935) [0.5914, 0.258, 0.3314]
PeriodicSite: Se (6.047, 2.915, 9.922) [0.5931, 0.2583, 0.6663]
PeriodicSite: Se (3.794, 6.681, 14.88) [0.5871, 0.5919, 0.999]
PeriodicSite: Se (3.77, 6.629, 4.938) [0.583, 0.5874, 0.3316]
PeriodicSite: Se (3.785, 6.631, 9.936) [0.5842, 0.5875, 0.6672]
PeriodicSite: Se (1.975, 10.49, 0.05121) [0.6163, 0.9295, 0.003439]
PeriodicSite: Se (1.846, 10.45, 4.93) [0.6047, 0.9261, 0.3311]
PeriodicSite: Se (1.807, 10.47, 9.943) [0.6025, 0.9276, 0.6678]
PeriodicSite: Se (10.35, 2.952, 14.89) [0.9249, 0.2616, 0.9997]
PeriodicSite: Se (10.35, 2.902, 4.942) [0.9229, 0.2571, 0.3319]
PeriodicSite: Se (10.32, 2.99, 9.921) [0.9242, 0.265, 0.6662]
PeriodicSite: Se (8.132, 6.769, 0.01877) [0.9238, 0.5997, 0.001261]
PeriodicSite: Se (8.145, 6.822, 4.934) [0.9272, 0.6044, 0.3313]
PeriodicSite: Se (8.149, 6.795, 9.93) [0.9263, 0.6021, 0.6669]
PeriodicSite: Se (6.256, 10.83, 0.2637) [0.9599, 0.9599, 0.01771]
PeriodicSite: Se (6.078, 10.44, 4.958) [0.9288, 0.9248, 0.333]
PeriodicSite: Se (6.025, 10.43, 9.915) [0.9242, 0.9237, 0.6659]
We can directly write the interpolated structures to file, to use for VASP calculations to generate our CC diagram / NEB initial guess (in which case we’d set n_images rather than displacements for get_path_structures/write_path_structures) using write_path_structures:
from doped.utils.configurations import write_path_structures
write_path_structures(
V_Se_m1_supercell, V_Se_m2_like_m1, displacements=displacements, output_dir="V_Se_-1_to_-2"
)
!ls V_Se_-1_to_-2
PES_1 PES_2
!ls V_Se_-1_to_-2/PES_1
delQ_-0.1 delQ_-0.4 delQ_-1.0 delQ_0.0 delQ_0.3 delQ_0.8 delQ_1.5
delQ_-0.2 delQ_-0.6 delQ_-1.2 delQ_0.1 delQ_0.4 delQ_1.0
delQ_-0.3 delQ_-0.8 delQ_-1.5 delQ_0.2 delQ_0.6 delQ_1.2
!ls V_Se_-1_to_-2/PES_1/delQ_0.6
POSCAR
Note
The same workflow applies for generating initial guess structures for NEB paths, but the n_images argument should then be used with get_path_structures/write_path_structures.
See the nonrad and/or CarrierCapture.jl tutorials for examples of plotting and analysis of the results of these calculations.
NEB Between Crystal Polymorphs
These functions can also be very useful when one wants to perform NEB calculations (or simply singlepoint calculations along an interpolated path) between different crystal structures – as performed in this study on AgBiS₂ – TODO
In this example, we want to calculate the minimum energy pathway along the structural PES from the \(R\bar{3}m\) structure (previously thought to be the ordered ground-state of AgBiS₂) and the recently-discovered \(P3_221\) phase:
from pymatgen.core.structure import Structure
RPA_R3m = Structure.from_file("AgBiS2/R-3m_RPA_PBE_Min_Energy_Relaxed_CONTCAR")
RPA_P3_221 = Structure.from_file("AgBiS2/P3_221_RPA_PBE_Min_Energy_Relaxed_CONTCAR")
RPA_R3m.lattice
Lattice
abc : 6.741702332938467 6.741702275281131 6.741702747042641
angles : 34.537890604070625 34.53789369373768 34.53788801879601
volume : 87.87241365397952
A : np.float64(3.8313124954623725) np.float64(0.0067957894543121) np.float64(5.5472108960546604)
B : np.float64(1.7344427349110751) np.float64(3.416241975347925) np.float64(5.547210896812548)
C : np.float64(0.01104923956544) np.float64(0.0067957900781565) np.float64(6.7416902673605925)
pbc : True True True
RPA_P3_221.lattice
Lattice
abc : 4.0672456500258996 4.0672456500258996 19.281437325056643
angles : 90.0 90.0 119.99999999999999
volume : 276.22999999999973
A : np.float64(4.0672456500258996) np.float64(1.2e-15) np.float64(1e-16)
B : np.float64(-2.0336228250129498) np.float64(3.522338056354181) np.float64(3e-16)
C : np.float64(-4e-16) np.float64(2e-16) np.float64(19.281437325056643)
pbc : True True True
As we can see, the unit cell sizes and lattice vectors differ for these two polymorphs, so we need to determine commensurate supercell sizes to perform our NEB calculation. There is no visually-obvious supercell conversion here, so let’s see what the conventional cell definitions look like:
from doped.utils.symmetry import get_BCS_conventional_structure
R3m_conv_struct = get_BCS_conventional_structure(RPA_R3m)[0]
R3m_conv_struct.lattice
Lattice
abc : 4.002647591940581 4.002647591940581 18.999773676
angles : 90.0 90.0 119.99999999901789
volume : 263.6172159502596
A : np.float64(2.001323796) np.float64(-3.466394497) np.float64(0.0)
B : np.float64(2.001323796) np.float64(3.466394497) np.float64(0.0)
C : np.float64(0.0) np.float64(0.0) np.float64(18.999773676)
pbc : True True True
Great, in this case the conventional cell definition for \(R\bar{3}m\) matches that of \(P3_221\), so we can use this. Alternatively, we could also have tried using the Lattice.find_all_mappings() method from pymatgen, with adjusted ltol and atol, to help find commensurate cell definitions.
As above, we need to ensure matching orientation and atomic indices between our structures, so that our NEB initial guess corresponds to the closest linear path between these structures:
from doped.utils.configurations import orient_s2_like_s1
RPA_R3m_like_P3_221 = orient_s2_like_s1(RPA_P3_221, R3m_conv_struct)
/Users/kavanase/Packages/doped/doped/utils/configurations.py:78: UserWarning: The lattices of the two input structures have been detected to be (symmetry-)inequivalent. This is usually not desirable for defect NEBs/CC diagrams, but may be the case for e.g. NEBs between polymorphs.
warnings.warn(
/Users/kavanase/Packages/doped/doped/utils/configurations.py:84: UserWarning: Note that the lattice definitions may differ between the output structure and ``struct1``. See the NEB/CC diagram tutorial for details.
warnings.warn(
Here we get a warning about these being inequivalent lattices, which is typically not desired for defect CCD/NEB calculations, but is expected in this case. We also get a warning that the lattice definitions for the two structures might still differ in this case, due to difficulties introduced with the symmetry-inequivalent lattices. Let’s check:
RPA_R3m_like_P3_221.lattice
Lattice
abc : 4.002647591940581 4.002647592 18.999773676
angles : 90.0 90.0 120.00000000049107
volume : 263.6172159502596
A : np.float64(-2.001323796) np.float64(-3.466394497) np.float64(0.0)
B : np.float64(4.002647592) np.float64(0.0) np.float64(0.0)
C : np.float64(0.0) np.float64(0.0) np.float64(-18.999773676)
pbc : True True True
RPA_P3_221.lattice
Lattice
abc : 4.0672456500258996 4.0672456500258996 19.281437325056643
angles : 90.0 90.0 119.99999999999999
volume : 276.22999999999973
A : np.float64(4.0672456500258996) np.float64(1.2e-15) np.float64(1e-16)
B : np.float64(-2.0336228250129498) np.float64(3.522338056354181) np.float64(3e-16)
C : np.float64(-4e-16) np.float64(2e-16) np.float64(19.281437325056643)
pbc : True True True
Here we can see that the lattice vectors are equivalent, but with different definitions (e.g. the B vector for RPA_R3m_like_P3_221 matches the A vector for RPA_P3_221, and the y and z directions of the lattice vectors have opposite signs). We can redefine our lattice for RPA_R3m_like_P3_221 to give the closest match to that of RPA_P3_221 (note that we avoid just directly using the lattice of RPA_P3_221 here, as this corresponds to a different space group (despite being similar)). We can do this a couple of different equivalent ways:
Using
dopedandpymatgenstructure manipulation functions:
from doped.utils.symmetry import swap_axes
# swap A and B lattice vectors:
RPA_R3m_like_P3_221_AB_swapped = swap_axes(RPA_R3m_like_P3_221, axes=[1,0,2])
# rotate lattice to flip the signs of the y and z lattice vector components:
from pymatgen.core.operations import SymmOp
import numpy as np
flip_y_z = np.array([ # rotation matrix to flip the y and z lattice vector components
[1, 0, 0], # no change to x
[0, -1, 0], # flip y
[0, 0, -1] # flip z
])
# create symmetry operation with this rotataion and no translation
op = SymmOp.from_rotation_and_translation(flip_y_z, [0, 0, 0])
RPA_R3m_like_P3_221_w_matching_lattice = RPA_R3m_like_P3_221_AB_swapped.copy()
RPA_R3m_like_P3_221_w_matching_lattice.apply_operation(op)
RPA_R3m_like_P3_221_w_matching_lattice.lattice # matches RPA_P3_221!
Lattice
abc : 4.002647592 4.002647591940581 18.999773676
angles : 90.0 90.0 120.00000000049107
volume : 263.6172159502596
A : np.float64(4.002647592) np.float64(0.0) np.float64(0.0)
B : np.float64(-2.001323796) np.float64(3.466394497) np.float64(0.0)
C : np.float64(0.0) np.float64(0.0) np.float64(18.999773676)
pbc : True True True
Alternatively, achieve the same transformation by directly redefining the lattice matrix: (Note that manual redefinition like this is only valid for symmetry-preserving transformations, which can be checked using
Structure.get_space_group_info()after the transformation).
from pymatgen.core.lattice import Lattice
from pymatgen.core.sites import PeriodicSite
import numpy as np
RPA_R3m_like_P3_221_lattice = Lattice(np.array([
RPA_R3m_like_P3_221.lattice.matrix[1], # B -> A
[
RPA_R3m_like_P3_221.lattice.matrix[0][0],
-RPA_R3m_like_P3_221.lattice.matrix[0][1],
RPA_R3m_like_P3_221.lattice.matrix[0][2]], # flip sign of A y coordinate to match P3_221 definition, and set to B vector
-RPA_R3m_like_P3_221.lattice.matrix[2] # -C -> C
])
) # Note that these are fully equivalent lattice definitions
RPA_R3m_like_P3_221_w_matching_lattice = Structure.from_sites(
[
PeriodicSite(
site.specie,
site.frac_coords,
RPA_R3m_like_P3_221_lattice,
properties=site.properties,
to_unit_cell=True,
)
for site in RPA_R3m_like_P3_221.sites
]
)
# confirm our space group is correct:
RPA_R3m_like_P3_221_w_matching_lattice.get_space_group_info()
('R-3m', 166)
Tip
If the appropriate lattice transformation was not obvious, we could also tried the Lattice.find_all_mappings() method with adjusted ltol and atol, to help find commensurate cell definitions.
Visually confirm matching atomic indices (similar fractional coordinates for site 1, 2, 3 etc in both structures):
RPA_R3m_like_P3_221_w_matching_lattice
Structure Summary
Lattice
abc : 4.002647592 4.002647591940581 18.999773676
angles : 90.0 90.0 120.00000000049107
volume : 263.6172159502596
A : np.float64(4.002647592) np.float64(0.0) np.float64(0.0)
B : np.float64(-2.001323796) np.float64(3.466394497) np.float64(0.0)
C : np.float64(-0.0) np.float64(-0.0) np.float64(18.999773676)
pbc : True True True
PeriodicSite: Ag (0.6671, 3.466, 3.167) [0.6667, 1.0, 0.1667]
PeriodicSite: Ag (0.6671, 1.155, 9.5) [0.3333, 0.3333, 0.5]
PeriodicSite: Ag (-1.334, 2.311, 15.83) [0.0, 0.6667, 0.8333]
PeriodicSite: Bi (0.6671, 3.466, 12.67) [0.6667, 1.0, 0.6667]
PeriodicSite: Bi (0.6671, 1.155, 1.647e-07) [0.3333, 0.3333, 8.667e-09]
PeriodicSite: Bi (-1.334, 2.311, 6.333) [0.0, 0.6667, 0.3333]
PeriodicSite: S (0.6671, 1.155, 14.3) [0.3333, 0.3333, 0.7524]
PeriodicSite: S (-1.334, 2.311, 1.629) [0.0, 0.6667, 0.08576]
PeriodicSite: S (0.6671, 3.466, 7.963) [0.6667, 1.0, 0.4191]
PeriodicSite: S (0.6671, 1.155, 4.704) [0.3333, 0.3333, 0.2476]
PeriodicSite: S (-1.334, 2.311, 11.04) [0.0, 0.6667, 0.5809]
PeriodicSite: S (0.6671, 3.466, 17.37) [0.6667, 1.0, 0.9142]
RPA_P3_221
Structure Summary
Lattice
abc : 4.0672456500258996 4.0672456500258996 19.281437325056643
angles : 90.0 90.0 119.99999999999999
volume : 276.22999999999973
A : np.float64(4.0672456500258996) np.float64(1.2e-15) np.float64(1e-16)
B : np.float64(-2.0336228250129498) np.float64(3.522338056354181) np.float64(3e-16)
C : np.float64(-4e-16) np.float64(2e-16) np.float64(19.281437325056643)
pbc : True True True
PeriodicSite: Ag (2.874, 8.812e-16, 3.214) [0.7066, -0.0, 0.1667]
PeriodicSite: Ag (0.5967, 1.033, 9.641) [0.2934, 0.2934, 0.5]
PeriodicSite: Ag (-1.437, 2.489, 16.07) [-0.0, 0.7066, 0.8333]
PeriodicSite: Bi (3.047, 1.032e-15, 12.85) [0.7491, 0.0, 0.6667]
PeriodicSite: Bi (0.5102, 0.8837, 1.003e-16) [0.2509, 0.2509, -0.0]
PeriodicSite: Bi (-1.523, 2.639, 6.427) [-0.0, 0.7491, 0.3333]
PeriodicSite: S (1.264, 1.192, 14.46) [0.4801, 0.3385, 0.7498]
PeriodicSite: S (2.434, 1.831, 1.602) [0.8584, 0.5199, 0.0831]
PeriodicSite: S (2.403, 0.4988, 8.029) [0.6615, 0.1416, 0.4164]
PeriodicSite: S (0.4003, 1.691, 4.825) [0.3385, 0.4801, 0.2502]
PeriodicSite: S (-0.7693, 2.33, 11.25) [0.1416, 0.6615, 0.5836]
PeriodicSite: S (0.369, 3.023, 17.68) [0.5199, 0.8584, 0.9169]
We can also numerically confirm matching atomic positions and lattices:
# checking matching atomic fractional coordinates (still some mismatch as they are different structures,
# but far less than with original orientations)
from doped.utils.configurations import get_dQ
for name, struct in zip(
["R3m_conv_struct", "RPA_R3m_like_P3_221", "RPA_R3m_like_P3_221_w_matching_lattice"],
[R3m_conv_struct, RPA_R3m_like_P3_221, RPA_R3m_like_P3_221_w_matching_lattice]
):
dQ = get_dQ(struct, RPA_P3_221)
print(f"{name}: ΔQ = {dQ:.2f} amu^(1/2)Å " + ("(matching!)" if dQ < 25 else "(not matching)"))
R3m_conv_struct: ΔQ = 204.04 amu^(1/2)Å (not matching)
RPA_R3m_like_P3_221: ΔQ = 11.88 amu^(1/2)Å (matching!)
RPA_R3m_like_P3_221_w_matching_lattice: ΔQ = 11.88 amu^(1/2)Å (matching!)
# checking matching lattices, using Frobenius norm (still some mismatch as they are different structures,
# but far less than with original orientations)
import numpy as np
for name, struct in zip(
["R3m_conv_struct", "RPA_R3m_like_P3_221", "RPA_R3m_like_P3_221_w_matching_lattice"],
[R3m_conv_struct, RPA_R3m_like_P3_221, RPA_R3m_like_P3_221_w_matching_lattice]
):
fro_norm = np.linalg.norm(struct.lattice.matrix - RPA_P3_221.lattice.matrix, ord="fro")
print(f"{name}: Frobenius norm of lattice matrix difference = {fro_norm:.2f} Å")
R3m_conv_struct: Frobenius norm of lattice matrix difference = 5.71 Å
RPA_R3m_like_P3_221: Frobenius norm of lattice matrix difference = 39.54 Å
RPA_R3m_like_P3_221_w_matching_lattice: Frobenius norm of lattice matrix difference = 0.30 Å
Now we can write our interpolated structures to file and use them for NEB calculations:
from doped.utils.configurations import write_path_structures
write_path_structures(
RPA_R3m_like_P3_221_w_matching_lattice, RPA_P3_221, n_images=6, output_dir="AgBiS2/R-3m_to_P3_221"
) # here we set n_images to generate one set of evenly-spaced structures along the linear path, to use
# for NEB calculations
!ls AgBiS2/R-3m_to_P3_221
00 01 02 03 04 05 06
!ls AgBiS2/R-3m_to_P3_221/03
POSCAR
Results:
# plot final energies from NEB calculations:
import numpy as np
import matplotlib.pyplot as plt
plt.style.use("~/Packages/publication_style.mplstyle")
energies = np.array([
-56.62423270,
-56.64091892,
-56.67399708,
-56.72550242,
-56.78033792,
-56.82650156,
-56.85690367,
-56.86827542,
])
energies -= np.min(energies)
energies_per_atom = (energies/len(RPA_R3m_like_P3_221_w_matching_lattice))*1000 # meV/atom
# mirror around 0, because we have two equivalent R-3m -> P3_221 distortion paths
energies_per_atom = np.concatenate([energies_per_atom[::-1], energies_per_atom[1:]])
# interpolate:
from scipy.interpolate import CubicSpline
spl = CubicSpline(np.arange(len(energies_per_atom)), energies_per_atom)
x_interp = np.linspace(0, len(energies_per_atom)-1, 1000)
y_interp = spl(x_interp)
plt.plot(x_interp, y_interp, alpha=0.5)
plt.scatter(np.arange(len(energies_per_atom)), energies_per_atom, marker='o')
plt.ylabel("Energy (meV/atom)")
# set x-labels as space group symbols:
plt.xticks(
[0, len(energies_per_atom)//2, len(energies_per_atom)-1], [r"$P3_221$", r"$R\bar{3}m$", r"$P3_221$"]
)
([<matplotlib.axis.XTick at 0x17dc12f90>,
<matplotlib.axis.XTick at 0x17f844690>,
<matplotlib.axis.XTick at 0x17f844e10>],
[Text(0, 0, '$P3_221$'), Text(7, 0, '$R\\bar{3}m$'), Text(14, 0, '$P3_221$')])
Advanced Manipulations: \(C2/c\) -> \(P3_221\)
In this example, we want to calculate the minimum energy pathway along the structural PES from the \(C2/c\) to \(P3_221\) candidate ground-state ordered structures of AgBiS₂.
from pymatgen.core.structure import Structure
RPA_C2c = Structure.from_file("AgBiS2/C2c_RPA_PBE_Min_Energy_Relaxed_CONTCAR")
RPA_P3_221 = Structure.from_file("AgBiS2/P3_221_RPA_PBE_Min_Energy_Relaxed_CONTCAR")
RPA_C2c.lattice
Lattice
abc : 4.039854199079283 4.0398541990792705 12.92447481944885
angles : 84.14385813351586 95.85614186648412 118.88312418556708
volume : 183.39319781406442
A : np.float64(3.478163280494634) np.float64(2.053930162897532) np.float64(0.0645989936225031)
B : np.float64(-3.4781632804946274) np.float64(2.053930162897519) np.float64(-0.0645989936225019)
C : np.float64(-1.769438330034499) np.float64(-2.6e-15) np.float64(12.80277850135556)
pbc : True True True
RPA_P3_221.lattice
Lattice
abc : 4.0672456500258996 4.0672456500258996 19.281437325056643
angles : 90.0 90.0 119.99999999999999
volume : 276.22999999999973
A : np.float64(4.0672456500258996) np.float64(1.2e-15) np.float64(1e-16)
B : np.float64(-2.0336228250129498) np.float64(3.522338056354181) np.float64(3e-16)
C : np.float64(-4e-16) np.float64(2e-16) np.float64(19.281437325056643)
pbc : True True True
from doped.utils.symmetry import get_BCS_conventional_structure
get_BCS_conventional_structure(RPA_P3_221)[0].lattice # better match
Lattice
abc : 4.067245649712695 4.067245649712695 19.281437325
angles : 90.0 90.0 119.99999999532659
volume : 276.2299999696536
A : np.float64(2.033622825) np.float64(-3.522338056) np.float64(0.0)
B : np.float64(2.033622825) np.float64(3.522338056) np.float64(0.0)
C : np.float64(0.0) np.float64(0.0) np.float64(19.281437325)
pbc : True True True
Here the conventional cell definition for \(P3_221\) has a closer match to the \(C2/c\) lattice (though with the x and y axes swapped, and with different c-axis lengths), so let’s use this.
len(RPA_C2c)
8
len(get_BCS_conventional_structure(RPA_P3_221)[0])
12
The \(C2/c\) cell has 8 atoms (2 formula units), while the \(P3_221\) cell has 12 atoms (3 formula units), and we can see that the difference comes from the \(c\) lattice vector which is 1.5 times longer in \(P3_221\) than in \(C2/c\). So we can make supercell expansions along \(c\) to get matching cell definitions:
RPA_C2c_matched_to_112_P3_221 = RPA_C2c * [1,1,3]
RPA_P3_221_matched_to_113_C2c = get_BCS_conventional_structure(RPA_P3_221)[0] * [1, 1, 2]
from doped.utils.configurations import orient_s2_like_s1
RPA_P3_221_matched_to_113_C2c_oriented = orient_s2_like_s1(
RPA_C2c_matched_to_112_P3_221, RPA_P3_221_matched_to_113_C2c
)
RPA_P3_221_matched_to_113_C2c_oriented.lattice
/Users/kavanase/Packages/doped/doped/utils/configurations.py:78: UserWarning: The lattices of the two input structures have been detected to be (symmetry-)inequivalent. This is usually not desirable for defect NEBs/CC diagrams, but may be the case for e.g. NEBs between polymorphs.
warnings.warn(
/Users/kavanase/Packages/doped/doped/utils/configurations.py:84: UserWarning: Note that the lattice definitions may differ between the output structure and ``struct1``. See the NEB/CC diagram tutorial for details.
warnings.warn(
Lattice
abc : 4.067245649712695 4.067245649712695 39.20105563367645
angles : 81.04665188307624 98.95334811692376 119.99999999532659
volume : 552.4599999393072
A : np.float64(-2.033622825) np.float64(3.522338056) np.float64(0.0)
B : np.float64(-2.033622825) np.float64(-3.522338056) np.float64(0.0)
C : np.float64(0.0) np.float64(-7.044676112) np.float64(38.56287465)
pbc : True True True
As above, we get a warning about inequivalent lattices as expected for this specific case of NEB between different polymorphs. We also get a warning that the lattice definitions for the two structures might still differ in this case, due to difficulties with the symmetry-inequivalent lattices.
Let’s compare the lattice definition of RPA_P3_221_matched_to_113_C2c_oriented (output above) with that of RPA_C2c_matched_to_112_P3_221:
RPA_C2c_matched_to_112_P3_221.lattice
Lattice
abc : 4.039854199079283 4.0398541990792705 38.77342445834655
angles : 84.14385813351586 95.85614186648412 118.88312418556708
volume : 550.1795934421933
A : np.float64(3.478163280494634) np.float64(2.053930162897532) np.float64(0.0645989936225031)
B : np.float64(-3.4781632804946274) np.float64(2.053930162897519) np.float64(-0.0645989936225019)
C : np.float64(-5.308314990103497) np.float64(-7.8e-15) np.float64(38.40833550406668)
pbc : True True True
So we can see that the x and y lattice vector components are swapped between the two structures, and there is a sign change between the x component of the \(P3_221\) lattice vector and the y component of the \(C2/c\) lattice vector. We can re-orient this to match using the structure manipulation functions in doped and pymatgen:
from doped.utils.symmetry import swap_axes
# rotate lattice to swap the x and y lattice vector components:
from pymatgen.core.operations import SymmOp
import numpy as np
# rotation matrix to swap the x and y lattice vector components, and flip the x component of the P3_221 lattice vector
swap_x_y = np.array([
[0, 1, 0], # swap x and y
[-1, 0, 0], # swap x and y, and flip x sign
[0, 0, 1] # no change to z
])
# create symmetry operation with this rotataion and no translation
op = SymmOp.from_rotation_and_translation(swap_x_y, [0, 0, 0])
RPA_P3_221_matched_to_113_C2c_oriented_w_matching_lattice = RPA_P3_221_matched_to_113_C2c_oriented.copy()
RPA_P3_221_matched_to_113_C2c_oriented_w_matching_lattice.apply_operation(op)
RPA_P3_221_matched_to_113_C2c_oriented_w_matching_lattice.lattice # matches RPA_C2c!
Lattice
abc : 4.067245649712695 4.067245649712695 39.20105563367645
angles : 81.04665188307624 98.95334811692376 119.99999999532659
volume : 552.4599999393072
A : np.float64(3.522338056) np.float64(2.033622825) np.float64(0.0)
B : np.float64(-3.522338056) np.float64(2.033622825) np.float64(0.0)
C : np.float64(-7.044676112) np.float64(0.0) np.float64(38.56287465)
pbc : True True True
RPA_P3_221_matched_to_113_C2c_oriented_w_matching_lattice.get_space_group_info() # confirms symmetry-preserving transformation
('P3_221', 154)
Numerically confirm matching atomic positions and lattices:
# checking matching atomic fractional coordinates:
from doped.utils.configurations import get_dQ
for name, struct in zip(
["RPA_P3_221_matched_to_113_C2c", "RPA_P3_221_matched_to_113_C2c_oriented_w_matching_lattice"],
[RPA_P3_221_matched_to_113_C2c, RPA_P3_221_matched_to_113_C2c_oriented_w_matching_lattice]
):
dQ = get_dQ(struct, RPA_C2c_matched_to_112_P3_221)
print(f"{name}: ΔQ = {dQ:.2f} amu^(1/2)Å " + ("(matching!)" if dQ < 25 else "(not matching)"))
RPA_P3_221_matched_to_113_C2c: ΔQ = 453.30 amu^(1/2)Å (not matching)
RPA_P3_221_matched_to_113_C2c_oriented_w_matching_lattice: ΔQ = 20.02 amu^(1/2)Å (matching!)
# checking matching lattices, using Frobenius norm:
import numpy as np
for name, struct in zip(
["RPA_P3_221_matched_to_113_C2c", "RPA_P3_221_matched_to_113_C2c_oriented_w_matching_lattice"],
[RPA_P3_221_matched_to_113_C2c, RPA_P3_221_matched_to_113_C2c_oriented_w_matching_lattice]
):
fro_norm = np.linalg.norm(struct.lattice.matrix - RPA_C2c_matched_to_112_P3_221.lattice.matrix, ord="fro")
print(f"{name}: Frobenius norm of lattice matrix difference = {fro_norm:.2f} Å")
RPA_P3_221_matched_to_113_C2c: Frobenius norm of lattice matrix difference = 9.69 Å
RPA_P3_221_matched_to_113_C2c_oriented_w_matching_lattice: Frobenius norm of lattice matrix difference = 1.75 Å
Now we can write our interpolated structures to file and use them for NEB calculations:
from doped.utils.configurations import write_path_structures
write_path_structures(
RPA_C2c_matched_to_112_P3_221, RPA_P3_221_matched_to_113_C2c_oriented_w_matching_lattice,
n_images=6, output_dir="AgBiS2/C2c_to_P3_221"
) # here we set n_images to generate one set of evenly-spaced structures along the linear path, to use
# for NEB calculations
# tutorial cleanup:
!rm -r AgBiS2/R-3m_to_P3_221
!rm -r AgBiS2/C2c_to_P3_221
!rm -r V_Se_-1_to_-2
rm: AgBiS2/R-3m_to_P3_221: No such file or directory
rm: V_Se_-1_to_-2: No such file or directory