Skip to content

Commit

Permalink
Turn rates and activation energies into dataframes (#250)
Browse files Browse the repository at this point in the history
  • Loading branch information
stefsmeets authored Jan 31, 2024
1 parent c0b76ce commit 05b2365
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 28 deletions.
39 changes: 19 additions & 20 deletions src/gemdat/jumps.py
Original file line number Diff line number Diff line change
Expand Up @@ -236,22 +236,18 @@ def collective(self, max_dist: float = 1) -> Collective:
)

@weak_lru_cache()
def activation_energies(
self,
n_parts: int = 10) -> dict[tuple[str, str], tuple[float, float]]:
def activation_energies(self, n_parts: int = 10) -> pd.DataFrame:
"""Calculate activation energies for jumps (UNITS?).
Returns
-------
e_act : dict[tuple[str, str], tuple[float, float]]
Dictionary with jump activation energies and standard deviations between site pairs.
n_parts : int
Number of parts to split transitions/jumps into for statistics
df : pd.DataFrame
Dataframe with jump activation energies and standard deviations between site pairs.
"""
trajectory = self.transitions.trajectory
attempt_freq, _ = SimulationMetrics(trajectory).attempt_frequency()

e_act = {}
dct = {}

temperature = trajectory.metadata['temperature']

Expand Down Expand Up @@ -285,10 +281,13 @@ def activation_energies(
e_act_arr = -np.log(eff_rate / attempt_freq) * (
Boltzmann * temperature) / elementary_charge

e_act[site_start,
site_stop] = np.mean(e_act_arr), np.std(e_act_arr, ddof=1)
dct[site_start, site_stop] = np.mean(e_act_arr), np.std(e_act_arr,
ddof=1)

return e_act
df = pd.DataFrame(dct).T
df.columns = ('energy', 'std')

return df

def jumps_counter(self) -> Counter:
"""Calculate number of jumps between sites.
Expand Down Expand Up @@ -318,18 +317,15 @@ def split(self, n_parts) -> list[Jumps]:
]

@weak_lru_cache()
def rates(self,
n_parts: int = 10) -> dict[tuple[str, str], tuple[float, float]]:
def rates(self, n_parts: int = 10) -> pd.DataFrame:
"""Calculate jump rates (total jumps / second).
Returns
-------
rates : dict[tuple[str, str], tuple[float, float]]
Dictionary with jump rates and standard deviations between site pairs
n_parts : int
Number of parts to split jumps into for statistics
df : pd.DataFrame
Dataframe with jump rates and standard deviations between site pairs
"""
rates: dict[tuple[str, str], tuple[float, float]] = {}
dct = {}

parts = [part.jumps_counter() for part in self.split(n_parts)]

Expand All @@ -342,6 +338,9 @@ def rates(self,
jump_freq_mean = np.mean(n_jumps) / denom
jump_freq_std = np.std(n_jumps, ddof=1) / denom

rates[site_pair] = float(jump_freq_mean), float(jump_freq_std)
dct[site_pair] = float(jump_freq_mean), float(jump_freq_std)

df = pd.DataFrame(dct).T
df.columns = ('rates', 'std')

return rates
return df
18 changes: 10 additions & 8 deletions tests/integration/sites_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
from math import isclose

import numpy as np
import pandas as pd
import pytest
from numpy.testing import assert_allclose

Expand Down Expand Up @@ -125,23 +126,24 @@ def test_n_jumps(self, vasp_jumps):

def test_rates(self, vasp_jumps):
rates = vasp_jumps.rates(n_parts=10)
assert isinstance(rates, dict)
assert isinstance(rates, pd.DataFrame)
assert len(rates) == 1

rates, rates_std = rates[('48h', '48h')]
assert isclose(rates, 542307692307.69226)
assert isclose(rates_std, 41893421993.683655)
row = rates.loc[('48h', '48h')]

assert isclose(row['rates'], 542307692307.69226)
assert isclose(row['std'], 41893421993.683655)

def test_activation_energies(self, vasp_jumps, vasp_sites):
activation_energies = vasp_jumps.activation_energies(n_parts=10)

assert isinstance(activation_energies, dict)
assert isinstance(activation_energies, pd.DataFrame)
assert len(activation_energies) == 1

e_act, e_act_std = activation_energies[('48h', '48h')]
row = activation_energies.loc[('48h', '48h')]

assert isclose(e_act, 0.17445, abs_tol=1e-4)
assert isclose(e_act_std, 0.004059, abs_tol=1e-6)
assert isclose(row['energy'], 0.17445, abs_tol=1e-4)
assert isclose(row['std'], 0.004059, abs_tol=1e-6)

def test_jump_diffusivity(self, vasp_jumps):
assert isclose(vasp_jumps.jump_diffusivity(3),
Expand Down

0 comments on commit 05b2365

Please sign in to comment.