From e8b670db2f0cbb601bf4b0501a1c4fcfd49f4d70 Mon Sep 17 00:00:00 2001 From: lucy whalley Date: Tue, 21 Feb 2023 11:48:56 +0000 Subject: [PATCH] updates and fix to script --- Kabsch_interpolation.ipynb | 70 +++++++++++++++----------------------- kabsch_interpolation.py | 14 ++++---- 2 files changed, 35 insertions(+), 49 deletions(-) diff --git a/Kabsch_interpolation.ipynb b/Kabsch_interpolation.ipynb index d76939e..cd78fdd 100644 --- a/Kabsch_interpolation.ipynb +++ b/Kabsch_interpolation.ipynb @@ -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. Too large or too small values for `mic_cutoff` can lead to unexpected output!\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." ] }, { @@ -45,7 +26,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 6, "id": "2e78162e", "metadata": {}, "outputs": [], @@ -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", @@ -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" ] }, { @@ -78,7 +61,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 7, "id": "e5b0c808", "metadata": {}, "outputs": [], @@ -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", + " " ] }, { @@ -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)" ] }, @@ -352,7 +336,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 5, "id": "baf88ab5", "metadata": {}, "outputs": [ @@ -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" ] } @@ -431,7 +415,7 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -445,7 +429,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.13" + "version": "3.9.6" } }, "nbformat": 4, diff --git a/kabsch_interpolation.py b/kabsch_interpolation.py index 78a8d6f..cdb041c 100644 --- a/kabsch_interpolation.py +++ b/kabsch_interpolation.py @@ -112,7 +112,10 @@ 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: @@ -120,14 +123,11 @@ def interpolate_structures(start_atoms, end_atoms, molecular_formulas=None, # 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( @@ -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)