Skip to content

Commit

Permalink
updates and fix to script
Browse files Browse the repository at this point in the history
  • Loading branch information
@lucydot
lucydot committed Feb 21, 2023
1 parent 3234457 commit e8b670d
Show file tree
Hide file tree
Showing 2 changed files with 35 additions and 49 deletions.
70 changes: 27 additions & 43 deletions Kabsch_interpolation.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,26 +13,7 @@
"- [Citing the associated paper](https://arxiv.org/abs/2302.08412) (currently under review)\n",
"- [Citing the Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/faq.html#how-should-i-cite-ase)\n",
"\n",
"### What does this code do?\n",
"- This code interpolates between two crystal structures\n",
"- It is aimed at hybrid organic-inorganic materials where distortions to the inorganic framework can be accurately described using translational interpolation, whereas molecular species require a description of translation and rotation.\n",
"- For isolated atoms it linearly interpolates along the translation vector mapping between the atom start position and atom end position. \n",
"- For molecules it linearly interpolates along the translation vector mapping between the molecule start centre-of-mass and molecule end centre-of-mass, and combines this with a linear interpolation along the rotation vector mapping between the molecule start orientation and molecule end orientation. To do so it uses the [Kabsch algorithm](https://en.wikipedia.org/wiki/Kabsch_algorithm). \n",
"- It makes extensive use of the [Atomic Simulation Environment](https://wiki.fysik.dtu.dk/ase/index.html) and so can be used to read-in/generate structures for various atomistic simulation codes.\n",
"\n",
"### What does this code assume?\n",
"- any molecular distortion can be ignored; the molecules are treated as rigid bodies.\n",
"- molecules can be identified using natural cutoffs. Please see the [following link](https://wiki.fysik.dtu.dk/ase/ase/neighborlist.html#ase.neighborlist.natural_cutoffs) for more information. Alternatively, a neighbours list can be explicitly provided.\n",
"- atoms in the start Atoms object and end Atoms object are listed in the same order. \n",
"\n",
"### How do I use this code for my own research?\n",
"- you need to have the [ASE](https://wiki.fysik.dtu.dk/ase/index.html), [scipy](https://scipy.org/) and [numpy](https://numpy.org/doc/stable/index.html) Python libraries installed. If you are new to Python, [Anaconda](https://www.anaconda.com/products/distribution) is recommended.\n",
"- read the docstring for the `interpolate_structures` function. This is the function you need to call to generate the interpolated structures (there is an example call to this function in the penultimate notebook cell). \n",
"- You can use this notebook, or the script `kabsch_interpolation.py`, to generate the interpolated structures.\n",
"- If your molecules are not interpolating as expected you may need to vary the `mic_cutoff` parameter. <mark>Too large or too small values for `mic_cutoff` can lead to unexpected output!</mark>\n",
"- ASE can also be used to identify unusual changes in bond lengths. There is example code for this in the final cell of this notebook.\n",
"- If you have any questions about this code please contact l.whalley@northumbria.ac.uk or (preferably) raise an issue on the [repository issue tracker](https://github.com/NU-CEM/Kabsch_interpolation/issues). \n",
"\n"
"Please see https://github.com/NU-CEM/Kabsch_interpolation for more information."
]
},
{
Expand All @@ -45,7 +26,7 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 6,
"id": "2e78162e",
"metadata": {},
"outputs": [],
Expand All @@ -57,6 +38,7 @@
"from ase import Atoms\n",
"from ase import neighborlist\n",
"\n",
"from IPython.display import HTML\n",
"import numpy as np\n",
"from numpy.linalg import norm \n",
"from scipy.spatial.transform import Rotation as R\n",
Expand All @@ -65,7 +47,8 @@
"\n",
"import itertools\n",
"from collections import Counter\n",
"from copy import deepcopy"
"from copy import deepcopy\n",
"from tempfile import NamedTemporaryFile"
]
},
{
Expand All @@ -78,7 +61,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 7,
"id": "e5b0c808",
"metadata": {},
"outputs": [],
Expand Down Expand Up @@ -293,7 +276,8 @@
" if end_atoms[i].position[j] - position > start_atoms.cell.cellpar()[j]/2:\n",
" end_atoms[i].position[j] = (end_atoms[i].position[j] - \n",
" start_atoms.cell.cellpar()[j])\n",
" return start_atoms, end_atoms"
" return start_atoms, end_atoms\n",
" "
]
},
{
Expand Down Expand Up @@ -335,8 +319,8 @@
}
],
"source": [
"start_atoms = ase.io.read(\"./structures/POSCAR_MixedNegOP2-14.vasp\")\n",
"end_atoms = ase.io.read(\"./structures/POSCAR_MixedNeutOP2-14.vasp\")\n",
"start_atoms = ase.io.read(\"./POSCAR_start.vasp\")\n",
"end_atoms = ase.io.read(\"./POSCAR_end.vasp\")\n",
"interpolate_structures(start_atoms, end_atoms, molecular_formulas = [\"CNH6\"], number_intermediates=9, fformat=\"vasp\", reverse=True, molecular_indices=None, translation_species = [\"Cs\",\"Pb\",\"I\"], mic_cutoff=0.5)"
]
},
Expand All @@ -352,7 +336,7 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 5,
"id": "baf88ab5",
"metadata": {},
"outputs": [
Expand All @@ -362,39 +346,39 @@
"text": [
"~~~~POSCAR_000.vasp~~~~~\n",
"There are 94 Pb-I bonds in a supercell of 16 12-atom primitive cells.\n",
"The average Pb-I bond length is 3.1913613806222125.\n",
"The standard deviation is 0.06097140338001101.\n",
"The maximum bond is: 3.411827858717294\n",
"The minimum bond is: 3.072450446327065\n",
"The average Pb-I bond length is 3.191361380622213.\n",
"The standard deviation is 0.060971403380010926.\n",
"The maximum bond is: 3.4118278587172934\n",
"The minimum bond is: 3.0724504463270645\n",
"~~~~POSCAR_002.vasp~~~~~\n",
"There are 95 Pb-I bonds in a supercell of 16 12-atom primitive cells.\n",
"The average Pb-I bond length is 3.191447340531054.\n",
"The standard deviation is 0.06050571815005562.\n",
"The maximum bond is: 3.41812658395874\n",
"The average Pb-I bond length is 3.1914473405310537.\n",
"The standard deviation is 0.06050571815005559.\n",
"The maximum bond is: 3.418126583958739\n",
"The minimum bond is: 3.0759907227815773\n",
"~~~~POSCAR_004.vasp~~~~~\n",
"There are 96 Pb-I bonds in a supercell of 16 12-atom primitive cells.\n",
"The average Pb-I bond length is 3.1919412052555316.\n",
"The standard deviation is 0.06098370188758939.\n",
"The maximum bond is: 3.42566730700801\n",
"The standard deviation is 0.060983701887589385.\n",
"The maximum bond is: 3.4256673070080095\n",
"The minimum bond is: 3.0702041001002494\n",
"~~~~POSCAR_006.vasp~~~~~\n",
"There are 96 Pb-I bonds in a supercell of 16 12-atom primitive cells.\n",
"The average Pb-I bond length is 3.190709941597333.\n",
"The standard deviation is 0.05913309672396905.\n",
"The average Pb-I bond length is 3.1907099415973335.\n",
"The standard deviation is 0.05913309672396904.\n",
"The maximum bond is: 3.4344418470158464\n",
"The minimum bond is: 3.0507873566948955\n",
"~~~~POSCAR_008.vasp~~~~~\n",
"There are 96 Pb-I bonds in a supercell of 16 12-atom primitive cells.\n",
"The average Pb-I bond length is 3.190510977396766.\n",
"The standard deviation is 0.06044411627898002.\n",
"The standard deviation is 0.06044411627897989.\n",
"The maximum bond is: 3.444440774734817\n",
"The minimum bond is: 3.0170944032241285\n",
"~~~~POSCAR_010.vasp~~~~~\n",
"There are 94 Pb-I bonds in a supercell of 16 12-atom primitive cells.\n",
"The average Pb-I bond length is 3.1852200661633847.\n",
"The standard deviation is 0.04973975081201891.\n",
"The maximum bond is: 3.311637603070428\n",
"The standard deviation is 0.049739750812018854.\n",
"The maximum bond is: 3.3116376030704275\n",
"The minimum bond is: 3.0019002233926653\n"
]
}
Expand Down Expand Up @@ -431,7 +415,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -445,7 +429,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.13"
"version": "3.9.6"
}
},
"nbformat": 4,
Expand Down
14 changes: 8 additions & 6 deletions kabsch_interpolation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,22 +112,22 @@ def interpolate_structures(start_atoms, end_atoms, molecular_formulas=None,
atom_translation = (end_atoms[atom_index].position -
start_atoms[atom_index].position)
beta = atom_translation*interval
positions[atom_index] = start_atoms[atom_index].position + beta
atoms[atom_index].position = start_atoms[atom_index].position + beta

# update positions with the interpolated positions
positions = atoms.positions

# interpolation along translation and rotation vectors
for molecule_indices in molecular_indices_list:

# ASE atoms object to describe particular molecule
start_molecule = start_atoms[molecule_indices]
end_molecule = end_atoms[molecule_indices]

# deepcopy to avoid overwriting start_atoms
molecule = deepcopy(start_molecule)

# need to apply minimum image convention to any molecule that may
# move between neighbouring unit cells during relaxation
if ((np.abs(start_molecule.positions) < mic_cutoff).any() or
(np.abs(end_molecule.positions) < mic_cutoff).any()):
if (((np.abs(start_molecule.positions) < mic_cutoff).any()) or
((np.abs(end_molecule.positions) < mic_cutoff).any())):
start_molecule.set_positions(ase.geometry.geometry.find_mic(
start_molecule.positions, start_atoms.cell)[0])
end_molecule.set_positions(ase.geometry.geometry.find_mic(
Expand All @@ -146,6 +146,8 @@ def interpolate_structures(start_atoms, end_atoms, molecular_formulas=None,
alpha = angle * interval

# apply the translation and rotation
# deepcopy to avoid overwriting start_atoms
molecule = deepcopy(start_molecule)
molecule.rotate(alpha, axis, center="COM")
molecule.translate(delta)

Expand Down

0 comments on commit e8b670d

Please sign in to comment.