From daf566cf47c335175f6a0ad6e04b758f56f52527 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Wed, 10 Apr 2024 15:59:19 +0200 Subject: [PATCH 01/16] Get spherical coords directly from Orientations --- src/gemdat/plots/matplotlib/_rotations.py | 7 ++----- src/gemdat/rotations.py | 17 +++++++++++++++++ src/gemdat/utils.py | 21 +++++++++++---------- tests/integration/conftest.py | 3 +-- 4 files changed, 31 insertions(+), 17 deletions(-) diff --git a/src/gemdat/plots/matplotlib/_rotations.py b/src/gemdat/plots/matplotlib/_rotations.py index bb4db243..163d92db 100644 --- a/src/gemdat/plots/matplotlib/_rotations.py +++ b/src/gemdat/plots/matplotlib/_rotations.py @@ -10,7 +10,6 @@ calculate_spherical_areas, mean_squared_angular_displacement, ) -from gemdat.utils import cartesian_to_spherical def rectilinear_plot(*, @@ -35,7 +34,7 @@ def rectilinear_plot(*, Output figure """ # Convert the trajectory to spherical coordinates - trajectory = cartesian_to_spherical(orientations.vectors, degrees=True) + trajectory = orientations.get_vectors_spherical() az = trajectory[:, :, 0].flatten() el = trajectory[:, :, 1].flatten() @@ -92,9 +91,7 @@ def bond_length_distribution(*, fig : matplotlib.figure.Figure Output figure """ - - # Convert the trajectory to spherical coordinates - trajectory = cartesian_to_spherical(orientations.vectors, degrees=True) + trajectory = orientations.get_vectors_spherical() fig, ax = plt.subplots() diff --git a/src/gemdat/rotations.py b/src/gemdat/rotations.py index 9d5fe0bd..7d935660 100644 --- a/src/gemdat/rotations.py +++ b/src/gemdat/rotations.py @@ -241,6 +241,23 @@ def transform(self, matrix: np.ndarray) -> Orientations: return replace(self, in_vectors=vectors) + def get_vectors_spherical(self, degrees: bool = True) -> np.ndarray: + """Return vectors in spherical coordinates. + + Parameters + ---------- + degrees : bool + If True, return angles as degrees, else radians + + Returns + ------- + np.array + azimuth, elevation, length + """ + + from gemdat.utils import cartesian_to_spherical + return cartesian_to_spherical(self.vectors) + def calculate_spherical_areas(shape: tuple[int, int], radius: float = 1) -> np.ndarray: diff --git a/src/gemdat/utils.py b/src/gemdat/utils.py index c7f41556..dc2f4495 100644 --- a/src/gemdat/utils.py +++ b/src/gemdat/utils.py @@ -229,25 +229,26 @@ def _cart2sph(x: np.ndarray, y: np.ndarray, return az, el, r -def cartesian_to_spherical(direct_cart: np.ndarray, - degrees: bool) -> np.ndarray: +def cartesian_to_spherical(cart_coords: np.ndarray, + *, + degrees: bool = True) -> np.ndarray: """Trajectory from cartesian coordinates to spherical coordinates. Parameters ---------- - direct_cart : np.ndarray - Trajectory of the unit vectors in conventional coordinates + cart_coords : np.ndarray + Trajectory of the unit vectors in cartesian setting degrees : bool If true, return angles in degrees Returns ------- - direction_spherical : np.ndarray + spherical_coords : np.ndarray Trajectory of the unit vectors in spherical coordinates """ - x = direct_cart[:, :, 0] - y = direct_cart[:, :, 1] - z = direct_cart[:, :, 2] + x = cart_coords[:, :, 0] + y = cart_coords[:, :, 1] + z = cart_coords[:, :, 2] az, el, r = _cart2sph(x, y, z) @@ -255,6 +256,6 @@ def cartesian_to_spherical(direct_cart: np.ndarray, az = np.degrees(az) el = np.degrees(el) - direction_spherical = np.stack((az, el, r), axis=-1) + spherical_coords = np.stack((az, el, r), axis=-1) - return direction_spherical + return spherical_coords diff --git a/tests/integration/conftest.py b/tests/integration/conftest.py index bc621138..96da22a2 100644 --- a/tests/integration/conftest.py +++ b/tests/integration/conftest.py @@ -11,7 +11,6 @@ from gemdat.rotations import Orientations from gemdat.shape import ShapeAnalyzer from gemdat.trajectory import Trajectory -from gemdat.utils import cartesian_to_spherical from gemdat.volume import trajectory_to_volume DATA_DIR = Path(__file__).parents[1] / 'data' @@ -149,5 +148,5 @@ def vasp_orientations(vasp_traj_rotations): @pytest.fixture(scope='module') def vasp_orientations_spherical(vasp_orientations): cf = vasp_orientations.get_conventional_coordinates() - cf_spheric = cartesian_to_spherical(cf, degrees=True) + cf_spheric = cf.get_vectors_spherical() return cf_spheric From 046e193d97ca480680c6246eaf6043aabbf78545 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Thu, 11 Apr 2024 09:30:55 +0200 Subject: [PATCH 02/16] Add methods for plots --- src/gemdat/plots/__init__.py | 4 ++-- src/gemdat/plots/matplotlib/__init__.py | 4 ++-- src/gemdat/plots/matplotlib/_rotations.py | 8 ++++---- src/gemdat/rotations.py | 18 ++++++++++++++++-- tests/integration/plot_test.py | 6 +++--- 5 files changed, 27 insertions(+), 13 deletions(-) diff --git a/src/gemdat/plots/__init__.py b/src/gemdat/plots/__init__.py index 4a3bb72f..bf9da054 100644 --- a/src/gemdat/plots/__init__.py +++ b/src/gemdat/plots/__init__.py @@ -8,7 +8,7 @@ jumps_3d_animation, path_on_grid, radial_distribution, - rectilinear_plot, + rectilinear, shape, unit_vector_autocorrelation, ) @@ -45,7 +45,7 @@ 'plot_3d', 'msd_per_element', 'radial_distribution', - 'rectilinear_plot', + 'rectilinear', 'shape', 'vibrational_amplitudes', 'energy_along_path', diff --git a/src/gemdat/plots/matplotlib/__init__.py b/src/gemdat/plots/matplotlib/__init__.py index 26384cd2..c879f1a0 100644 --- a/src/gemdat/plots/matplotlib/__init__.py +++ b/src/gemdat/plots/matplotlib/__init__.py @@ -22,7 +22,7 @@ from ._rdf import radial_distribution from ._rotations import ( bond_length_distribution, - rectilinear_plot, + rectilinear, unit_vector_autocorrelation, ) from ._shape import shape @@ -46,7 +46,7 @@ 'msd_per_element', 'path_on_grid', 'radial_distribution', - 'rectilinear_plot', + 'rectilinear', 'shape', 'unit_vector_autocorrelation', 'vibrational_amplitudes', diff --git a/src/gemdat/plots/matplotlib/_rotations.py b/src/gemdat/plots/matplotlib/_rotations.py index 163d92db..d1e3e3d3 100644 --- a/src/gemdat/plots/matplotlib/_rotations.py +++ b/src/gemdat/plots/matplotlib/_rotations.py @@ -12,10 +12,10 @@ ) -def rectilinear_plot(*, - orientations: Orientations, - shape: tuple[int, int] = (90, 360), - normalize_histo: bool = True) -> plt.Figure: +def rectilinear(*, + orientations: Orientations, + shape: tuple[int, int] = (90, 360), + normalize_histo: bool = True) -> plt.Figure: """Plot a rectilinear projection of a spherical function. This function uses the transformed trajectory. diff --git a/src/gemdat/rotations.py b/src/gemdat/rotations.py index 7d935660..99f9fc51 100644 --- a/src/gemdat/rotations.py +++ b/src/gemdat/rotations.py @@ -6,6 +6,7 @@ from pymatgen.symmetry.groups import PointGroup from gemdat.trajectory import Trajectory +from gemdat.utils import cartesian_to_spherical @dataclass @@ -254,10 +255,23 @@ def get_vectors_spherical(self, degrees: bool = True) -> np.ndarray: np.array azimuth, elevation, length """ - - from gemdat.utils import cartesian_to_spherical return cartesian_to_spherical(self.vectors) + def plot_rectilinear(self, **kwargs): + """See [gemdat.plots.rectilinear][] for more info.""" + from gemdat import plots + return plots.rectilinear(orientations=self, **kwargs) + + def plot_bond_length_distribution(self, **kwargs): + """See [gemdat.plots.bond_length_distribution][] for more info.""" + from gemdat import plots + return plots.bond_length_distribution(orientations=self, **kwargs) + + def plot_unit_vector_autocorrelation(self, **kwargs): + """See [gemdat.plots.unit_vector_autocorrelation][] for more info.""" + from gemdat import plots + return plots.unit_vector_autocorrelation(orientations=self, **kwargs) + def calculate_spherical_areas(shape: tuple[int, int], radius: float = 1) -> np.ndarray: diff --git a/tests/integration/plot_test.py b/tests/integration/plot_test.py index 272cad7d..10b3987d 100644 --- a/tests/integration/plot_test.py +++ b/tests/integration/plot_test.py @@ -95,14 +95,14 @@ def test_rectilinear(vasp_orientations): [1 / 2**0.5, 1 / 6**0.5, -1 / 3**0.5], [0, 2 / 6**0.5, 1 / 3**0.5]], ) orientations = vasp_orientations.normalize().transform(matrix=matrix) - plots.rectilinear_plot(orientations=orientations, normalize_histo=False) + orientations.plot_rectilinear(normalize_histo=False) @image_comparison2(baseline_images=['bond_length_distribution']) def test_bond_length_distribution(vasp_orientations): - plots.bond_length_distribution(orientations=vasp_orientations, bins=1000) + vasp_orientations.plot_bond_length_distribution(bins=1000) @image_comparison2(baseline_images=['unit_vector_autocorrelation']) def test_unit_vector_autocorrelation(vasp_orientations): - plots.unit_vector_autocorrelation(orientations=vasp_orientations) + vasp_orientations.plot_unit_vector_autocorrelation() From 147d9bad7579f305b17e711f625241aa37753eda Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Thu, 11 Apr 2024 09:54:07 +0200 Subject: [PATCH 03/16] Rename msad -> autocorrelation --- src/gemdat/plots/__init__.py | 4 +- src/gemdat/plots/matplotlib/__init__.py | 4 +- src/gemdat/plots/matplotlib/_rotations.py | 9 +--- src/gemdat/rotations.py | 50 +++-------------------- src/gemdat/utils.py | 41 +++++++++++++++++++ tests/integration/plot_test.py | 2 +- tests/rotations_test.py | 14 +++---- 7 files changed, 61 insertions(+), 63 deletions(-) diff --git a/src/gemdat/plots/__init__.py b/src/gemdat/plots/__init__.py index bf9da054..cb219c01 100644 --- a/src/gemdat/plots/__init__.py +++ b/src/gemdat/plots/__init__.py @@ -3,6 +3,7 @@ # Matplotlib plots from .matplotlib import ( + autocorrelation, bond_length_distribution, energy_along_path, jumps_3d_animation, @@ -10,7 +11,6 @@ radial_distribution, rectilinear, shape, - unit_vector_autocorrelation, ) # Plotly plots (matplotlib version might be available) @@ -50,5 +50,5 @@ 'vibrational_amplitudes', 'energy_along_path', 'path_on_grid', - 'unit_vector_autocorrelation', + 'autocorrelation', ] diff --git a/src/gemdat/plots/matplotlib/__init__.py b/src/gemdat/plots/matplotlib/__init__.py index c879f1a0..56dba896 100644 --- a/src/gemdat/plots/matplotlib/__init__.py +++ b/src/gemdat/plots/matplotlib/__init__.py @@ -21,9 +21,9 @@ ) from ._rdf import radial_distribution from ._rotations import ( + autocorrelation, bond_length_distribution, rectilinear, - unit_vector_autocorrelation, ) from ._shape import shape from ._vibration import ( @@ -48,6 +48,6 @@ 'radial_distribution', 'rectilinear', 'shape', - 'unit_vector_autocorrelation', + 'autocorrelation', 'vibrational_amplitudes', ] diff --git a/src/gemdat/plots/matplotlib/_rotations.py b/src/gemdat/plots/matplotlib/_rotations.py index d1e3e3d3..aa23bb44 100644 --- a/src/gemdat/plots/matplotlib/_rotations.py +++ b/src/gemdat/plots/matplotlib/_rotations.py @@ -8,7 +8,6 @@ from gemdat.rotations import ( Orientations, calculate_spherical_areas, - mean_squared_angular_displacement, ) @@ -132,7 +131,7 @@ def _skewnorm_fit(x): return fig -def unit_vector_autocorrelation( +def autocorrelation( *, orientations: Orientations, ) -> plt.Figure: @@ -148,11 +147,7 @@ def unit_vector_autocorrelation( fig : matplotlib.figure.Figure Output figure """ - - # The trajectory is expected to have shape (n_times, n_particles, n_coordinates) - trajectory = orientations.vectors - - ac = mean_squared_angular_displacement(trajectory) + ac = orientations.autocorrelation() ac_std = ac.std(axis=0) ac_mean = ac.mean(axis=0) diff --git a/src/gemdat/rotations.py b/src/gemdat/rotations.py index 99f9fc51..74a5a2e0 100644 --- a/src/gemdat/rotations.py +++ b/src/gemdat/rotations.py @@ -6,7 +6,7 @@ from pymatgen.symmetry.groups import PointGroup from gemdat.trajectory import Trajectory -from gemdat.utils import cartesian_to_spherical +from gemdat.utils import autocorrelation, cartesian_to_spherical @dataclass @@ -257,6 +257,9 @@ def get_vectors_spherical(self, degrees: bool = True) -> np.ndarray: """ return cartesian_to_spherical(self.vectors) + def autocorrelation(self): + return autocorrelation(self.vectors) + def plot_rectilinear(self, **kwargs): """See [gemdat.plots.rectilinear][] for more info.""" from gemdat import plots @@ -267,10 +270,10 @@ def plot_bond_length_distribution(self, **kwargs): from gemdat import plots return plots.bond_length_distribution(orientations=self, **kwargs) - def plot_unit_vector_autocorrelation(self, **kwargs): + def plot_autocorrelation(self, **kwargs): """See [gemdat.plots.unit_vector_autocorrelation][] for more info.""" from gemdat import plots - return plots.unit_vector_autocorrelation(orientations=self, **kwargs) + return plots.autocorrelation(orientations=self, **kwargs) def calculate_spherical_areas(shape: tuple[int, int], @@ -304,44 +307,3 @@ def calculate_spherical_areas(shape: tuple[int, int], #hacky way to get rid of singularity on poles areas[0, :] = areas[-1, 0] return areas - - -def mean_squared_angular_displacement(trajectory: np.ndarray) -> np.ndarray: - """Compute the mean squared angular displacement using FFT. - - Parameters - ---------- - trajectory : np.ndarray - The input signal in direct cartesian coordinates. It is expected - to have shape (n_times, n_particles, n_coordinates) - - Returns - ------- - msad: - The mean squared angular displacement - """ - n_times, n_particles, n_coordinates = trajectory.shape - - msad = np.zeros((n_particles, n_times)) - normalization = np.arange(n_times, 0, -1) - - for c in range(n_coordinates): - signal = trajectory[:, :, c] - - # Compute the FFT of the signal - fft_signal = np.fft.rfft(signal, n=2 * n_times - 1, axis=0) - # Compute the power spectral density in-place - np.square(np.abs(fft_signal), out=fft_signal) - # Compute the inverse FFT of the power spectral density - autocorr_c = np.fft.irfft(fft_signal, axis=0) - - # Only keep the positive times - autocorr_c = autocorr_c[:n_times, :] - - msad += autocorr_c.T / normalization - - # Normalize the msad such that it starts from 1 - # (this makes the normalization independent on the dimensions) - msad = msad / msad[:, 0, np.newaxis] - - return msad diff --git a/src/gemdat/utils.py b/src/gemdat/utils.py index dc2f4495..a7d1a4f6 100644 --- a/src/gemdat/utils.py +++ b/src/gemdat/utils.py @@ -259,3 +259,44 @@ def cartesian_to_spherical(cart_coords: np.ndarray, spherical_coords = np.stack((az, el, r), axis=-1) return spherical_coords + + +def autocorrelation(coords: np.ndarray) -> np.ndarray: + """Compute the autocorrelation of the given coordinates using FFT. + + Parameters + ---------- + coords : np.ndarray + The input signal in direct cartesian coordinates. It is expected + to have shape (n_times, n_particles, n_coordinates) + + Returns + ------- + autocorrelation: np.array + The autocorrelation of the input signal + """ + n_times, n_particles, n_coordinates = coords.shape + + autocorrelation = np.zeros((n_particles, n_times)) + normalization = np.arange(n_times, 0, -1) + + for c in range(n_coordinates): + signal = coords[:, :, c] + + # Compute the FFT of the signal + fft_signal = np.fft.rfft(signal, n=2 * n_times - 1, axis=0) + # Compute the power spectral density in-place + np.square(np.abs(fft_signal), out=fft_signal) + # Compute the inverse FFT of the power spectral density + autocorr_c = np.fft.irfft(fft_signal, axis=0) + + # Only keep the positive times + autocorr_c = autocorr_c[:n_times, :] + + autocorrelation += autocorr_c.T / normalization + + # Normalize the autocorrelation such that it starts from 1 + # (this makes the normalization independent on the dimensions) + autocorrelation = autocorrelation / autocorrelation[:, 0, np.newaxis] + + return autocorrelation diff --git a/tests/integration/plot_test.py b/tests/integration/plot_test.py index 10b3987d..91e3c442 100644 --- a/tests/integration/plot_test.py +++ b/tests/integration/plot_test.py @@ -105,4 +105,4 @@ def test_bond_length_distribution(vasp_orientations): @image_comparison2(baseline_images=['unit_vector_autocorrelation']) def test_unit_vector_autocorrelation(vasp_orientations): - vasp_orientations.plot_unit_vector_autocorrelation() + vasp_orientations.plot_autocorrelation() diff --git a/tests/rotations_test.py b/tests/rotations_test.py index 6de6de6c..d3144e9f 100644 --- a/tests/rotations_test.py +++ b/tests/rotations_test.py @@ -8,8 +8,8 @@ from gemdat.rotations import ( Orientations, calculate_spherical_areas, - mean_squared_angular_displacement, ) +from gemdat.utils import autocorrelation def test_orientations_init(trajectory): @@ -90,9 +90,9 @@ def test_calculate_spherical_areas(): assert areas.shape == shape -def test_mean_squared_angular_displacement(trajectory): - msad = mean_squared_angular_displacement(trajectory.positions) - assert isinstance(msad, np.ndarray) - assert isclose(msad.mean(), 0.8142314269325723) - assert msad.shape == (trajectory.positions.shape[1], - trajectory.positions.shape[0]) +def test_autocorrelation(trajectory): + autocorr = autocorrelation(trajectory.positions) + assert isinstance(autocorr, np.ndarray) + assert isclose(autocorr.mean(), 0.8142314269325723) + assert autocorr.shape == (trajectory.positions.shape[1], + trajectory.positions.shape[0]) From e51eb979dae90ef20e4c451ade8a756d1e6226e5 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Thu, 11 Apr 2024 09:54:23 +0200 Subject: [PATCH 04/16] Fix deprecation error --- pyproject.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 4abe16ab..6fa2d4e4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -82,7 +82,7 @@ source = ["gemdat"] [tool.pytest.ini_options] testpaths = ["tests"] -[tool.ruff] +[tool.lint] # Enable Pyflakes `E` and `F` codes by default. select = [ "F", # Pyflakes @@ -93,7 +93,7 @@ select = [ line-length = 110 -[tool.ruff.isort] +[tool.lint.isort] known-first-party=["gemdat"] known-third-party = ["pymatgen"] required-imports = ["from __future__ import annotations"] From a0abf58c10d8414ac1ee78be0aa9adeef03d5c07 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Thu, 11 Apr 2024 10:14:05 +0200 Subject: [PATCH 05/16] Rename rotations -> orientations for consistency --- src/gemdat/__init__.py | 2 +- src/gemdat/{rotations.py => orientations.py} | 2 +- src/gemdat/plots/matplotlib/__init__.py | 2 +- .../{_rotations.py => _orientations.py} | 2 +- src/gemdat/utils.py | 2 +- tests/conftest.py | 2 +- tests/integration/conftest.py | 16 ++++++++-------- tests/integration/rotations_test.py | 18 +++++++++--------- tests/rotations_test.py | 2 +- 9 files changed, 24 insertions(+), 24 deletions(-) rename src/gemdat/{rotations.py => orientations.py} (99%) rename src/gemdat/plots/matplotlib/{_rotations.py => _orientations.py} (99%) diff --git a/src/gemdat/__init__.py b/src/gemdat/__init__.py index 2beb72b8..aed8be35 100644 --- a/src/gemdat/__init__.py +++ b/src/gemdat/__init__.py @@ -2,7 +2,7 @@ from .io import load_known_material, read_cif from .jumps import Jumps -from .rotations import Orientations +from .orientations import Orientations from .shape import ShapeAnalyzer from .simulation_metrics import SimulationMetrics from .trajectory import Trajectory diff --git a/src/gemdat/rotations.py b/src/gemdat/orientations.py similarity index 99% rename from src/gemdat/rotations.py rename to src/gemdat/orientations.py index 74a5a2e0..44e68e38 100644 --- a/src/gemdat/rotations.py +++ b/src/gemdat/orientations.py @@ -26,7 +26,7 @@ class Orientations: nr_central_atoms: int Number of central atoms, which corresponds to the number of cluster molecules vectors: np.ndarray - Vectors representing rotation direction + Vectors representing orientation direction """ trajectory: Trajectory center_type: str diff --git a/src/gemdat/plots/matplotlib/__init__.py b/src/gemdat/plots/matplotlib/__init__.py index 56dba896..cc90af17 100644 --- a/src/gemdat/plots/matplotlib/__init__.py +++ b/src/gemdat/plots/matplotlib/__init__.py @@ -20,7 +20,7 @@ path_on_grid, ) from ._rdf import radial_distribution -from ._rotations import ( +from ._orientations import ( autocorrelation, bond_length_distribution, rectilinear, diff --git a/src/gemdat/plots/matplotlib/_rotations.py b/src/gemdat/plots/matplotlib/_orientations.py similarity index 99% rename from src/gemdat/plots/matplotlib/_rotations.py rename to src/gemdat/plots/matplotlib/_orientations.py index aa23bb44..70373d26 100644 --- a/src/gemdat/plots/matplotlib/_rotations.py +++ b/src/gemdat/plots/matplotlib/_orientations.py @@ -5,7 +5,7 @@ from scipy.optimize import curve_fit from scipy.stats import skewnorm -from gemdat.rotations import ( +from gemdat.orientations import ( Orientations, calculate_spherical_areas, ) diff --git a/src/gemdat/utils.py b/src/gemdat/utils.py index a7d1a4f6..6b3d1790 100644 --- a/src/gemdat/utils.py +++ b/src/gemdat/utils.py @@ -13,7 +13,7 @@ # shortcut to test data VASPRUN = Path(__file__).parents[ 2] / 'tests' / 'data' / 'short_simulation' / 'vasprun.xml' -VASPCACHE_ROTATIONS = Path(__file__).parents[ +VASPCACHE_ORIENTATIONS = Path(__file__).parents[ 2] / 'tests' / 'data' / 'short_simulation' / 'vasprun_rotations.cache' DATA = files('gemdat') / 'data' diff --git a/tests/conftest.py b/tests/conftest.py index c2fa76d9..f5cf5900 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -4,7 +4,7 @@ import pytest from pymatgen.core import Species -from gemdat.rotations import Orientations +from gemdat.orientations import Orientations from gemdat.trajectory import Trajectory diff --git a/tests/integration/conftest.py b/tests/integration/conftest.py index 96da22a2..1c369f12 100644 --- a/tests/integration/conftest.py +++ b/tests/integration/conftest.py @@ -8,14 +8,14 @@ from gemdat.jumps import Jumps from gemdat.path import find_best_perc_path, free_energy_graph from gemdat.rdf import radial_distribution -from gemdat.rotations import Orientations +from gemdat.orientations import Orientations from gemdat.shape import ShapeAnalyzer from gemdat.trajectory import Trajectory from gemdat.volume import trajectory_to_volume DATA_DIR = Path(__file__).parents[1] / 'data' VASP_XML = DATA_DIR / 'short_simulation' / 'vasprun.xml' -VASP_ROTO_CACHE = DATA_DIR / 'short_simulation' / 'vasprun_rotations.cache' +VASP_ORI_CACHE = DATA_DIR / 'short_simulation' / 'vasprun_rotations.cache' def pytest_configure(): @@ -25,8 +25,8 @@ def pytest_configure(): ('Simulation data from vasprun.xml example is required for this test. ' 'Run `git submodule init`/`update`, and extract using `tar -C tests/data/short_simulation ' '-xjf tests/data/short_simulation/vasprun.xml.bz2`')) - pytest.vasprotocache_available = pytest.mark.skipif( - not VASP_ROTO_CACHE.exists(), + pytest.vasporicache_available = pytest.mark.skipif( + not VASP_ORI_CACHE.exists(), reason= ('Simulation data from vasprun_rotations.cache example is required for this test. ' 'Run `git submodule init`/`update`, and extract using `tar -C tests/data/short_simulation ' @@ -41,8 +41,8 @@ def vasp_traj(): @pytest.fixture(scope='module') -def vasp_traj_rotations(): - trajectory = Trajectory.from_cache(VASP_ROTO_CACHE) +def vasp_traj_orientations(): + trajectory = Trajectory.from_cache(VASP_ORI_CACHE) return trajectory @@ -135,11 +135,11 @@ def vasp_F_graph(vasp_path_vol): @pytest.fixture(scope='module') -def vasp_orientations(vasp_traj_rotations): +def vasp_orientations(vasp_traj_orientations): central_atoms = 'S' satellite_atoms = 'O' n_expected_neigh = 8 - orientations = Orientations(vasp_traj_rotations, central_atoms, + orientations = Orientations(vasp_traj_orientations, central_atoms, satellite_atoms, n_expected_neigh) return orientations diff --git a/tests/integration/rotations_test.py b/tests/integration/rotations_test.py index 6a64a537..8d1e3be3 100644 --- a/tests/integration/rotations_test.py +++ b/tests/integration/rotations_test.py @@ -5,14 +5,14 @@ import numpy as np import pytest -from gemdat.rotations import calculate_spherical_areas +from gemdat.orientations import calculate_spherical_areas matrix = np.array( [[1 / 2**0.5, -1 / 6**0.5, 1 / 3**0.5], [1 / 2**0.5, 1 / 6**0.5, -1 / 3**0.5], [0, 2 / 6**0.5, 1 / 3**0.5]], ) -@pytest.vasprotocache_available # type: ignore +@pytest.vasporicache_available # type: ignore def test_normalize(vasp_orientations): ret = vasp_orientations.normalize() assert isclose(ret.vectors.std(), 0.5773501929401034) @@ -21,7 +21,7 @@ def test_normalize(vasp_orientations): atol=1e-06) -@pytest.vasprotocache_available # type: ignore +@pytest.vasporicache_available # type: ignore def test_conventional(vasp_orientations): ret = vasp_orientations.transform(matrix=matrix) assert isclose(ret.vectors.std(), 0.8668510720176622) @@ -30,7 +30,7 @@ def test_conventional(vasp_orientations): atol=1e-06) -@pytest.vasprotocache_available # type: ignore +@pytest.vasporicache_available # type: ignore def test_symmetrize(vasp_orientations): ret = vasp_orientations.symmetrize(sym_group='m-3m') assert isclose(ret.vectors.std(), 0.8668511628167984) @@ -39,7 +39,7 @@ def test_symmetrize(vasp_orientations): atol=1e-06) -@pytest.vasprotocache_available # type: ignore +@pytest.vasporicache_available # type: ignore def test_normalize_symmetrize(vasp_orientations): ret = vasp_orientations.normalize().symmetrize(sym_group='m-3m') assert isclose(ret.vectors.std(), 0.577350269189626) @@ -48,7 +48,7 @@ def test_normalize_symmetrize(vasp_orientations): atol=1e-06) -@pytest.vasprotocache_available # type: ignore +@pytest.vasporicache_available # type: ignore def test_normalize_conventional_symmetrize(vasp_orientations): ret = vasp_orientations.normalize().transform(matrix=matrix).symmetrize( sym_group='m-3m') @@ -58,7 +58,7 @@ def test_normalize_conventional_symmetrize(vasp_orientations): atol=1e-06) -@pytest.vasprotocache_available # type: ignore +@pytest.vasporicache_available # type: ignore def test_conventional_normalize_symmetrize(vasp_orientations): ret = vasp_orientations.transform(matrix=matrix).normalize().symmetrize( sym_group='m-3m') @@ -68,7 +68,7 @@ def test_conventional_normalize_symmetrize(vasp_orientations): atol=1e-06) -@pytest.vasprotocache_available # type: ignore +@pytest.vasporicache_available # type: ignore def test_fractional_coordinates(vasp_orientations): direction = vasp_orientations._fractional_directions( vasp_orientations._distances) @@ -77,7 +77,7 @@ def test_fractional_coordinates(vasp_orientations): assert isclose(direction[3, 12, 2], -0.10902666) -@pytest.vasprotocache_available # type: ignore +@pytest.vasporicache_available # type: ignore def test_matching_matrix(vasp_orientations): distance = np.array([[1.0, 2.0, 1.5, 0.5, 3.0], [2.0, 1.0, 0.5, 1.5, 3.0]]) frac_coord_cent = np.array([[[0.3, 0.11, 0.78], [1.6, 1.0, 2.3]]]) diff --git a/tests/rotations_test.py b/tests/rotations_test.py index d3144e9f..039e993d 100644 --- a/tests/rotations_test.py +++ b/tests/rotations_test.py @@ -5,7 +5,7 @@ import numpy as np from pymatgen.core import Species -from gemdat.rotations import ( +from gemdat.orientations import ( Orientations, calculate_spherical_areas, ) From 49aa2be549e6012d888b5f9bd573daf098c6481b Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Thu, 11 Apr 2024 10:21:28 +0200 Subject: [PATCH 06/16] Prefix autocorrelation function to avoid name clash --- src/gemdat/orientations.py | 5 +++-- src/gemdat/utils.py | 2 +- tests/{rotations_test.py => orientations_test.py} | 6 +++--- 3 files changed, 7 insertions(+), 6 deletions(-) rename tests/{rotations_test.py => orientations_test.py} (96%) diff --git a/src/gemdat/orientations.py b/src/gemdat/orientations.py index 44e68e38..2cfc70bd 100644 --- a/src/gemdat/orientations.py +++ b/src/gemdat/orientations.py @@ -6,7 +6,7 @@ from pymatgen.symmetry.groups import PointGroup from gemdat.trajectory import Trajectory -from gemdat.utils import autocorrelation, cartesian_to_spherical +from gemdat.utils import fft_autocorrelation, cartesian_to_spherical @dataclass @@ -258,7 +258,8 @@ def get_vectors_spherical(self, degrees: bool = True) -> np.ndarray: return cartesian_to_spherical(self.vectors) def autocorrelation(self): - return autocorrelation(self.vectors) + """Compute the autocorrelation of the orientation vectors using FFT.""" + return fft_autocorrelation(self.vectors) def plot_rectilinear(self, **kwargs): """See [gemdat.plots.rectilinear][] for more info.""" diff --git a/src/gemdat/utils.py b/src/gemdat/utils.py index 6b3d1790..b1d1ad46 100644 --- a/src/gemdat/utils.py +++ b/src/gemdat/utils.py @@ -261,7 +261,7 @@ def cartesian_to_spherical(cart_coords: np.ndarray, return spherical_coords -def autocorrelation(coords: np.ndarray) -> np.ndarray: +def fft_autocorrelation(coords: np.ndarray) -> np.ndarray: """Compute the autocorrelation of the given coordinates using FFT. Parameters diff --git a/tests/rotations_test.py b/tests/orientations_test.py similarity index 96% rename from tests/rotations_test.py rename to tests/orientations_test.py index 039e993d..79303b04 100644 --- a/tests/rotations_test.py +++ b/tests/orientations_test.py @@ -9,7 +9,7 @@ Orientations, calculate_spherical_areas, ) -from gemdat.utils import autocorrelation +from gemdat.utils import fft_autocorrelation def test_orientations_init(trajectory): @@ -90,8 +90,8 @@ def test_calculate_spherical_areas(): assert areas.shape == shape -def test_autocorrelation(trajectory): - autocorr = autocorrelation(trajectory.positions) +def test_fft_autocorrelation(trajectory): + autocorr = fft_autocorrelation(trajectory.positions) assert isinstance(autocorr, np.ndarray) assert isclose(autocorr.mean(), 0.8142314269325723) assert autocorr.shape == (trajectory.positions.shape[1], From 04cc0f7252b3a127dda419408e523fd97c7434e6 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Thu, 11 Apr 2024 16:53:46 +0200 Subject: [PATCH 07/16] Minor refactor --- src/gemdat/orientations.py | 9 +++------ tests/orientations_test.py | 6 ------ 2 files changed, 3 insertions(+), 12 deletions(-) diff --git a/src/gemdat/orientations.py b/src/gemdat/orientations.py index 2cfc70bd..f2452be1 100644 --- a/src/gemdat/orientations.py +++ b/src/gemdat/orientations.py @@ -67,11 +67,6 @@ def _trajectory_sat(self) -> Trajectory: """Return trajectory of satellite atoms.""" return self.trajectory.filter(self.satellite_type) - def _fractional_coordinates(self) -> tuple[np.ndarray, np.ndarray]: - """Return fractional coordinates of central atoms and satellite - atoms.""" - return self._trajectory_cent.positions, self._trajectory_sat.positions - @property def _distances(self) -> np.ndarray: """Calculate distances between every central atom and all satellite @@ -148,7 +143,9 @@ def _fractional_directions(self, distance: np.ndarray) -> np.ndarray: direction: np.ndarray Contains the direction between central atoms and their ligands. """ - frac_coord_cent, frac_coord_sat = self._fractional_coordinates() + frac_coord_cent = self._trajectory_cent.positions + frac_coord_sat = self._trajectory_sat.positions + combinations = self._central_satellite_matrix(distance, frac_coord_cent) diff --git a/tests/orientations_test.py b/tests/orientations_test.py index 79303b04..220f4655 100644 --- a/tests/orientations_test.py +++ b/tests/orientations_test.py @@ -71,12 +71,6 @@ def test_orientations(orientations): assert orientations._trajectory_sat.species == [Species('Si')] -def test_fractional_coordinates(orientations): - frac_coord_cent, frac_coord_sat = orientations._fractional_coordinates() - assert isinstance(frac_coord_cent, np.ndarray) - assert isinstance(frac_coord_sat, np.ndarray) - - def test_distances(orientations): distances = orientations._distances assert isinstance(distances, np.ndarray) From 4d8958762c223ce1d7d47f5e46a745ed2c4306fa Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Thu, 11 Apr 2024 17:05:22 +0200 Subject: [PATCH 08/16] Remove rid of nr_central_atoms --- src/gemdat/orientations.py | 9 +++++---- tests/conftest.py | 4 +--- tests/integration/conftest.py | 3 +-- tests/orientations_test.py | 7 +------ 4 files changed, 8 insertions(+), 15 deletions(-) diff --git a/src/gemdat/orientations.py b/src/gemdat/orientations.py index f2452be1..c9d06baf 100644 --- a/src/gemdat/orientations.py +++ b/src/gemdat/orientations.py @@ -23,15 +23,12 @@ class Orientations: Type of the central atoms satellite_type: str Type of the satellite atoms - nr_central_atoms: int - Number of central atoms, which corresponds to the number of cluster molecules vectors: np.ndarray Vectors representing orientation direction """ trajectory: Trajectory center_type: str satellite_type: str - nr_central_atoms: int vectors: np.ndarray = field(init=False) in_vectors: InitVar[np.ndarray | None] = None @@ -123,7 +120,11 @@ def _central_satellite_matrix(self, distance: np.ndarray, combinations: np.ndarray Matrix of combinations between central and satellite atoms """ - index_central_atoms = np.arange(self.nr_central_atoms) + nr_central_atoms = frac_coord_cent.shape[1] + + index_central_atoms = np.arange(nr_central_atoms) + + # index_central_atoms = np.arange(self.nr_central_atoms) matching_matrix = self._matching_matrix(distance, frac_coord_cent) combinations = np.array([(i, j) for i in index_central_atoms for j in matching_matrix[i, :]]) diff --git a/tests/conftest.py b/tests/conftest.py index f5cf5900..70c1853a 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -33,6 +33,4 @@ def trajectory(): def orientations(trajectory): center_type = 'B' satellite_type = 'Si' - nr_central_atoms = 1 - return Orientations(trajectory, center_type, satellite_type, - nr_central_atoms) + return Orientations(trajectory, center_type, satellite_type) diff --git a/tests/integration/conftest.py b/tests/integration/conftest.py index 1c369f12..c6e9efba 100644 --- a/tests/integration/conftest.py +++ b/tests/integration/conftest.py @@ -138,9 +138,8 @@ def vasp_F_graph(vasp_path_vol): def vasp_orientations(vasp_traj_orientations): central_atoms = 'S' satellite_atoms = 'O' - n_expected_neigh = 8 orientations = Orientations(vasp_traj_orientations, central_atoms, - satellite_atoms, n_expected_neigh) + satellite_atoms) return orientations diff --git a/tests/orientations_test.py b/tests/orientations_test.py index 220f4655..880ab29a 100644 --- a/tests/orientations_test.py +++ b/tests/orientations_test.py @@ -15,20 +15,17 @@ def test_orientations_init(trajectory): orientations = Orientations(trajectory=trajectory, center_type='B', - satellite_type='Si', - nr_central_atoms=1) + satellite_type='Si') assert isinstance(orientations, Orientations) assert orientations.center_type == 'B' assert orientations.satellite_type == 'Si' - assert orientations.nr_central_atoms == 1 def test_normalize(trajectory): orientations = Orientations(trajectory=trajectory, center_type='B', satellite_type='Si', - nr_central_atoms=1, in_vectors=np.array([[1, 2, 2], [2, 2, 1]], dtype=float)) ret = orientations.normalize() @@ -40,7 +37,6 @@ def test_conventional(trajectory): orientations = Orientations(trajectory=trajectory, center_type='B', satellite_type='Si', - nr_central_atoms=1, in_vectors=np.array( [[1, 0, 0], [0, 1, 0], [0, 0, 1]], dtype=float)) @@ -54,7 +50,6 @@ def test_symmetrize(trajectory): orientations = Orientations(trajectory=trajectory, center_type='B', satellite_type='Si', - nr_central_atoms=1, in_vectors=np.array([[[1, 0, 0]], [[0, 1, 0]]], dtype=float)) sym_ops = np.array([[0, -1, 0], [1, 0, 0], [0, 0, -1]]) From cd8bfebc129599102ddf7dfe4bd7fe282af488fb Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Wed, 17 Apr 2024 13:22:53 +0200 Subject: [PATCH 09/16] Update src/gemdat/utils.py Co-authored-by: SCiarella <58949181+SCiarella@users.noreply.github.com> --- src/gemdat/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gemdat/utils.py b/src/gemdat/utils.py index b1d1ad46..4c887536 100644 --- a/src/gemdat/utils.py +++ b/src/gemdat/utils.py @@ -273,7 +273,7 @@ def fft_autocorrelation(coords: np.ndarray) -> np.ndarray: Returns ------- autocorrelation: np.array - The autocorrelation of the input signal + The autocorrelation of the input signal, with shape (n_particles, n_times) """ n_times, n_particles, n_coordinates = coords.shape From bb92b93495d5e1a5543bc872a6b7a5814565d947 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Wed, 17 Apr 2024 13:27:01 +0200 Subject: [PATCH 10/16] Update src/gemdat/utils.py Co-authored-by: SCiarella <58949181+SCiarella@users.noreply.github.com> --- src/gemdat/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gemdat/utils.py b/src/gemdat/utils.py index 4c887536..cbd32b94 100644 --- a/src/gemdat/utils.py +++ b/src/gemdat/utils.py @@ -296,7 +296,7 @@ def fft_autocorrelation(coords: np.ndarray) -> np.ndarray: autocorrelation += autocorr_c.T / normalization # Normalize the autocorrelation such that it starts from 1 - # (this makes the normalization independent on the dimensions) + # and make it independent of n_coordinates autocorrelation = autocorrelation / autocorrelation[:, 0, np.newaxis] return autocorrelation From 39f938f804fcd801e8b05b32597044592a2dc2a8 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Wed, 17 Apr 2024 13:27:21 +0200 Subject: [PATCH 11/16] Update tests/orientations_test.py Co-authored-by: SCiarella <58949181+SCiarella@users.noreply.github.com> --- tests/orientations_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/orientations_test.py b/tests/orientations_test.py index 880ab29a..39db8f87 100644 --- a/tests/orientations_test.py +++ b/tests/orientations_test.py @@ -20,7 +20,7 @@ def test_orientations_init(trajectory): assert isinstance(orientations, Orientations) assert orientations.center_type == 'B' assert orientations.satellite_type == 'Si' - + assert orientations.trajectory == trajectory def test_normalize(trajectory): orientations = Orientations(trajectory=trajectory, From 92981d12b86366b4157a9c0429c13d8e9b96d65c Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Wed, 17 Apr 2024 13:28:15 +0200 Subject: [PATCH 12/16] Rename notebook to orientations --- mkdocs.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mkdocs.yml b/mkdocs.yml index f9ad6e09..bd13863c 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -18,7 +18,7 @@ nav: - Pathways: notebooks/paths.ipynb - Percolation: notebooks/percolation.ipynb - Multiple paths: notebooks/multiple_paths.ipynb - - Orientation tracking: notebooks/rotations.ipynb + - Orientation tracking: notebooks/orientations.ipynb - Python API: - gemdat: api/gemdat.md - gemdat.collective: api/gemdat_collective.md From 1b6458e19b8bbd961d5947d9f3bc6f0fc2e561b0 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Wed, 17 Apr 2024 13:29:02 +0200 Subject: [PATCH 13/16] Yapf it --- tests/orientations_test.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/orientations_test.py b/tests/orientations_test.py index 39db8f87..659426b9 100644 --- a/tests/orientations_test.py +++ b/tests/orientations_test.py @@ -22,6 +22,7 @@ def test_orientations_init(trajectory): assert orientations.satellite_type == 'Si' assert orientations.trajectory == trajectory + def test_normalize(trajectory): orientations = Orientations(trajectory=trajectory, center_type='B', From ab4ee0461b6776a93a37bc1b9f369fead87b658f Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Wed, 17 Apr 2024 13:35:45 +0200 Subject: [PATCH 14/16] Update notebooks ref --- docs/notebooks | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/notebooks b/docs/notebooks index aa3caa64..7e4c0cfd 160000 --- a/docs/notebooks +++ b/docs/notebooks @@ -1 +1 @@ -Subproject commit aa3caa64339f850a3f925732c617e56c67e5aba9 +Subproject commit 7e4c0cfd0252d119e4bf5d7ad92370942e8af5ca From 3d5bac16cf8052213f05c3768f6d5ce3d615ef97 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Wed, 17 Apr 2024 13:44:06 +0200 Subject: [PATCH 15/16] Change function to property --- src/gemdat/orientations.py | 10 +++------- src/gemdat/plots/matplotlib/_orientations.py | 14 ++++++-------- 2 files changed, 9 insertions(+), 15 deletions(-) diff --git a/src/gemdat/orientations.py b/src/gemdat/orientations.py index c9d06baf..5d1dbfbd 100644 --- a/src/gemdat/orientations.py +++ b/src/gemdat/orientations.py @@ -240,13 +240,9 @@ def transform(self, matrix: np.ndarray) -> Orientations: return replace(self, in_vectors=vectors) - def get_vectors_spherical(self, degrees: bool = True) -> np.ndarray: - """Return vectors in spherical coordinates. - - Parameters - ---------- - degrees : bool - If True, return angles as degrees, else radians + @property + def vectors_spherical(self) -> np.ndarray: + """Return vectors in spherical coordinates in degrees. Returns ------- diff --git a/src/gemdat/plots/matplotlib/_orientations.py b/src/gemdat/plots/matplotlib/_orientations.py index 70373d26..e1647bea 100644 --- a/src/gemdat/plots/matplotlib/_orientations.py +++ b/src/gemdat/plots/matplotlib/_orientations.py @@ -32,11 +32,10 @@ def rectilinear(*, fig : matplotlib.figure.Figure Output figure """ - # Convert the trajectory to spherical coordinates - trajectory = orientations.get_vectors_spherical() - - az = trajectory[:, :, 0].flatten() - el = trajectory[:, :, 1].flatten() + # Convert the vectors to spherical coordinates + az, el, _ = orientations.vectors_spherical.T + az = az.flatten() + el = el.flatten() hist, xedges, yedges = np.histogram2d(el, az, shape) @@ -90,12 +89,11 @@ def bond_length_distribution(*, fig : matplotlib.figure.Figure Output figure """ - trajectory = orientations.get_vectors_spherical() + *_, bond_lengths = orientations.vectors_spherical.T + bond_lengths = bond_lengths.flatten() fig, ax = plt.subplots() - bond_lengths = trajectory[:, :, 2].flatten() - # Plot the normalized histogram hist, edges = np.histogram(bond_lengths, bins=bins, density=True) bin_centers = (edges[:-1] + edges[1:]) / 2 From a7497024810e91a3a497e64bd489d78029df55b3 Mon Sep 17 00:00:00 2001 From: Stef Smeets Date: Wed, 17 Apr 2024 13:45:50 +0200 Subject: [PATCH 16/16] Update notebook --- docs/notebooks | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/notebooks b/docs/notebooks index 7e4c0cfd..6821d690 160000 --- a/docs/notebooks +++ b/docs/notebooks @@ -1 +1 @@ -Subproject commit 7e4c0cfd0252d119e4bf5d7ad92370942e8af5ca +Subproject commit 6821d690b95d768649a3e41c50b12272fd95343b