From 8352542066f8711b0d402d7f5a537f4754aabca2 Mon Sep 17 00:00:00 2001 From: Harry-Rich Date: Wed, 22 Jan 2025 13:44:43 +0000 Subject: [PATCH 1/4] Added improved method for Centre of mass calculation --- kinisi/parser.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/kinisi/parser.py b/kinisi/parser.py index f480434..f545d1b 100644 --- a/kinisi/parser.py +++ b/kinisi/parser.py @@ -207,7 +207,8 @@ def get_disps(self, elif self.sampling == 'multi-origin': disp = np.concatenate([ drift_corrected[self.indices, np.newaxis, time_interval - 1], - np.subtract(drift_corrected[self.indices, time_interval:], drift_corrected[self.indices, :-time_interval]) + np.subtract(drift_corrected[self.indices, time_interval:], + drift_corrected[self.indices, :-time_interval]) ], axis=1) if np.multiply(*disp[:, ::time_interval].shape[:2]) <= 1: @@ -567,6 +568,7 @@ def __init__(self, structure, coords, latt, volume = self.get_structure_coords_latt(universe, sub_sample_atoms, sub_sample_traj, progress) + if specie is not None: indices = self.get_indices(structure, specie, framework_indices) elif isinstance(specie_indices, (list, tuple)): @@ -673,7 +675,8 @@ def _get_molecules(structure: "ase.atoms.Atoms" or "pymatgen.core.structure.Stru framework_indices) -> Tuple[np.ndarray, np.ndarray, Tuple[np.ndarray, np.ndarray]]: """ Determine framework and non-framework indices for an :py:mod:`ase` or :py:mod:`pymatgen` or :py:mod:`MDAnalysis` compatible file when specie_indices are provided and contain multiple molecules. Warning: This function changes the structure without changing the object. - + Calculates Centre of mass of provided particles groups as per the method described within: Link paper + :param structure: Initial structure. :param coords: fractional coordinates for all atoms. :param indices: indices for the atoms in the molecules in the trajectory used in the calculation @@ -716,10 +719,15 @@ def _get_molecules(structure: "ase.atoms.Atoms" or "pymatgen.core.structure.Stru zeta = np.sin(theta) xi_bar = np.average(xi, axis=-2, weights=weights) zeta_bar = np.average(zeta, axis=-2, weights=weights) - theta_bar = np.arctan2(zeta_bar, xi_bar) + theta_bar = np.arctan2(-zeta_bar, -xi_bar) + np.pi new_s_coords = theta_bar / (2 * np.pi) - new_coords = np.concatenate((new_s_coords, sq_coords[:, drift_indices]), axis=1) + # Implementation of improved method for center of mass calculation as described with: XXXX + pseudo_com_recentering = ((s_coords - (new_s_coords + 0.5)[:, :, np.newaxis]) % 1) + com_pseudo_space = np.average(pseudo_com_recentering, weights=masses, axis=2) + corrected_com = (com_pseudo_space + (new_s_coords + 0.5)) % 1 + + new_coords = np.concatenate((corrected_com, sq_coords[:, drift_indices]), axis=1) new_indices = list(range(n_molecules)) new_drift_indices = list(range(n_molecules, n_molecules + len(drift_indices))) From fe805a5b5bccb10b554fec6a1023d34c39167f3c Mon Sep 17 00:00:00 2001 From: Andrew McCluskey Date: Fri, 24 Jan 2025 14:18:02 +0000 Subject: [PATCH 2/4] Update kinisi/parser.py --- kinisi/parser.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/kinisi/parser.py b/kinisi/parser.py index f545d1b..a4b3f92 100644 --- a/kinisi/parser.py +++ b/kinisi/parser.py @@ -675,7 +675,7 @@ def _get_molecules(structure: "ase.atoms.Atoms" or "pymatgen.core.structure.Stru framework_indices) -> Tuple[np.ndarray, np.ndarray, Tuple[np.ndarray, np.ndarray]]: """ Determine framework and non-framework indices for an :py:mod:`ase` or :py:mod:`pymatgen` or :py:mod:`MDAnalysis` compatible file when specie_indices are provided and contain multiple molecules. Warning: This function changes the structure without changing the object. - Calculates Centre of mass of provided particles groups as per the method described within: Link paper + Calculates the centre of mass of provided particle groups using the pseudo-centre of mass approach (paper to come). :param structure: Initial structure. :param coords: fractional coordinates for all atoms. From b063fb67d3672aec7abd36a0180edd4f737c59c8 Mon Sep 17 00:00:00 2001 From: Andrew McCluskey Date: Fri, 24 Jan 2025 14:18:08 +0000 Subject: [PATCH 3/4] Update kinisi/parser.py --- kinisi/parser.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/kinisi/parser.py b/kinisi/parser.py index a4b3f92..c6800e7 100644 --- a/kinisi/parser.py +++ b/kinisi/parser.py @@ -722,7 +722,7 @@ def _get_molecules(structure: "ase.atoms.Atoms" or "pymatgen.core.structure.Stru theta_bar = np.arctan2(-zeta_bar, -xi_bar) + np.pi new_s_coords = theta_bar / (2 * np.pi) - # Implementation of improved method for center of mass calculation as described with: XXXX + # Implementation of pseudo-centre of mass approach to centre of mass calculation (paper to come). pseudo_com_recentering = ((s_coords - (new_s_coords + 0.5)[:, :, np.newaxis]) % 1) com_pseudo_space = np.average(pseudo_com_recentering, weights=masses, axis=2) corrected_com = (com_pseudo_space + (new_s_coords + 0.5)) % 1 From 17785c96f6d2b8fb83ec76ecab5760e3e6e384f5 Mon Sep 17 00:00:00 2001 From: Harry-Rich Date: Fri, 24 Jan 2025 14:32:48 +0000 Subject: [PATCH 4/4] Update test_parser.py Updated COM and COG tests --- kinisi/tests/test_parser.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/kinisi/tests/test_parser.py b/kinisi/tests/test_parser.py index 08b9f98..b64a4b2 100644 --- a/kinisi/tests/test_parser.py +++ b/kinisi/tests/test_parser.py @@ -131,7 +131,7 @@ def test_pymatgen_init_with_COG(self): assert_almost_equal(data.time_step, 1) assert_almost_equal(data.step_skip, 1) assert_equal(data.indices, [0]) - assert_almost_equal(data.coords_check, [[[0.269634905, 0.262183827, 0.2]]]) + assert_almost_equal(data.coords_check, [[[0.2733333, 0.2666667, 0.2 ]]]) def test_pymatgen_init_with_COM(self): xd = Xdatcar(os.path.join(os.path.dirname(kinisi.__file__), 'tests/inputs/example_center.XDATCAR')) @@ -146,7 +146,7 @@ def test_pymatgen_init_with_COM(self): assert_almost_equal(data.time_step, 1) assert_almost_equal(data.step_skip, 1) assert_equal(data.indices, [0]) - assert_almost_equal(data.coords_check, [[[0.382421597, 0.2087361, 0.2]]]) + assert_almost_equal(data.coords_check, [[[0.3788889, 0.2111111, 0.2 ]]]) def test_pymatgen_init_with_framwork_indices(self): xd = Xdatcar(os.path.join(os.path.dirname(kinisi.__file__), 'tests/inputs/example_drift.XDATCAR')) @@ -200,7 +200,7 @@ def test_ase_init_with_COG(self): assert_almost_equal(data.time_step, 1) assert_almost_equal(data.step_skip, 1) assert_equal(data.indices, [0]) - assert_almost_equal(data.coords_check, [[[0.269634905, 0.262183827, 0.2]]]) + assert_almost_equal(data.coords_check, [[[0.2733333, 0.2666667, 0.2]]]) def test_ase_init_with_COM(self): traj = Trajectory(os.path.join(os.path.dirname(kinisi.__file__), 'tests/inputs/example_ase_center.traj')) @@ -215,7 +215,7 @@ def test_ase_init_with_COM(self): assert_almost_equal(data.time_step, 1) assert_almost_equal(data.step_skip, 1) assert_equal(data.indices, [0]) - assert_almost_equal(data.coords_check, [[[0.382421597, 0.2087361, 0.2]]]) + assert_almost_equal(data.coords_check, [[[0.3788889, 0.2111111, 0.2]]]) def test_ase_init_with_framwork_indices(self): traj = Trajectory(os.path.join(os.path.dirname(kinisi.__file__), 'tests/inputs/example_ase_drift.traj')) @@ -266,7 +266,7 @@ def test_mda_init_with_molecules(self): data = parser.MDAnalysisParser(xd, **da_params) assert_almost_equal(data.time_step, 0.005) assert_almost_equal(data.step_skip, 250) - assert_almost_equal(data.coords_check[0, 0], [0.46465184, 0.03090944, 0.4023758]) + assert_almost_equal(data.coords_check[0, 0], [0.5169497, 0.1174514, 0.3637794]) assert_equal(data.indices, list(range(len(molecules)))) def test_mda_init_with_COG(self): @@ -279,7 +279,7 @@ def test_mda_init_with_COG(self): assert_almost_equal(data.time_step, 1) assert_almost_equal(data.step_skip, 1) assert_equal(data.indices, [0]) - assert_almost_equal(data.coords_check, [[[0.269634905, 0.262183827, 0.2]]]) + assert_almost_equal(data.coords_check, [[[0.2733333, 0.2666667, 0.1999999]]]) def test_mda_init_with_COM(self): xd = mda.Universe(os.path.join(os.path.dirname(kinisi.__file__), 'tests/inputs/example_LAMMPS_center.data'), @@ -297,7 +297,7 @@ def test_mda_init_with_COM(self): assert_almost_equal(data.time_step, 1) assert_almost_equal(data.step_skip, 1) assert_equal(data.indices, [0]) - assert_almost_equal(data.coords_check, [[[0.382421597, 0.2087361, 0.2]]]) + assert_almost_equal(data.coords_check, [[[0.3788889, 0.2111111, 0.2]]]) def test_mda_init_with_framwork_indices(self): xd = mda.Universe(os.path.join(os.path.dirname(kinisi.__file__), 'tests/inputs/example_LAMMPS_drift.data'), @@ -319,4 +319,4 @@ def test_mda_init_with_framwork_indices(self): assert_equal(data_2.drift_indices, list(range(1, 9))) print(data_2.dc[0]) disp_array = [[0.0, 0.0, 0.0], [0.2, 0.2, 0.2], [0.4, 0.4, 0.4], [1, 1, 1]] - assert_almost_equal(data_2.dc[0], disp_array, decimal=6) \ No newline at end of file + assert_almost_equal(data_2.dc[0], disp_array, decimal=6)