Find jump sites from Trajectory¶

This notebook shows how to find the jump sites from the trajectory directly.

The idea is the diffusing species naturally cluster on (meta-)stable sites in the simulation. These sites can be used as the jump sites to calculate transitions.

The procedure is as follows:

Load trajectory
Convert trajectory into density volume
Find peaks in volume
Convert peaks into pymatgen structure

First, we load the data from a VASP simulation.

In [1]:

from gemdat import Trajectory
from gemdat.utils import VASPRUN

# Use your own data:
# VASPRUN = 'path/to/your/vasprun.xml'

trajectory = Trajectory.from_vasprun(VASPRUN)

Create a new trajectory object that only consists of the diffusing species, Lithium.

In [2]:

# Create a new trajectory object that only consists of Lithium
diff_trajectory = trajectory.filter('Li')
diff_trajectory

Out[2]:

Full Formula (Li48)
Reduced Formula: Li
abc   :  19.873726   9.919447   9.916454
angles:  90.214114  90.859135  89.950474
pbc   :       True       True       True
Constant lattice (True)
Sites (48)
Time steps (5000)

Converting a Trajectory to a Volume¶

To convert a Trajectory to a Volume object we use Trajectory.to_volume() method provided by gemdat. A volume in this context is a discretized grid version of the trajectory. Each grid voxel acts as a bin, and each coordinate in the trajectory is assigned to a bin. The number of of points in each bin is counted to create the density volume.

The resolution is the minimum size of the voxels in Ångstrom.

Note that a voxel does not necessarily represent a cubic or even orthorhombic volume, and depends on the lattice parameters of the trajectory.

In [3]:

vol = diff_trajectory.to_volume(resolution=0.3)

Converting a Volume to Sites¶

We can directly convert a volume to a structure with the Volume.to_structure() method. This uses a watershed algorithm to find the sites.

The background_level sets the sensitivity of the peak finder. Set it to a higher value if you get too many spurious peaks. Sometimes the generated sites need to be cleaned up a bit manually.

In [4]:

sites = vol.to_structure(specie='Li', background_level=0.1)
sites[0:10]

Out[4]:

[PeriodicSite: Li (7.214, 5.229, 2.009) [0.3636, 0.5273, 0.2091],
 PeriodicSite: Li (8.248, 5.08, 3.038) [0.4161, 0.5123, 0.3137],
 PeriodicSite: Li (3.125, 8.8, 5.397) [0.1591, 0.8879, 0.5485],
 PeriodicSite: Li (4.325, 8.024, 5.877) [0.2197, 0.8098, 0.5976],
 PeriodicSite: Li (3.117, 9.76, 6.566) [0.1591, 0.9848, 0.6667],
 PeriodicSite: Li (1.977, 9.701, 7.965) [0.1023, 0.9792, 0.8068],
 PeriodicSite: Li (16.76, 2.665, 9.515) [0.847, 0.2697, 0.9727],
 PeriodicSite: Li (17.97, 2.203, 9.503) [0.9078, 0.2231, 0.9723],
 PeriodicSite: Li (11.9, 9.459, 8.221) [0.6019, 0.9545, 0.8401],
 PeriodicSite: Li (12.21, 8.683, 9.274) [0.6178, 0.8765, 0.9463]]

Verify the generated sites by plotting them on the volume.

In [5]:

vol.plot_3d(structure=sites)

Helper functions are available on the Volume object, for example, finding peaks or converting a volume to a vasp file. See gemdat.Volume for more info.

This is useful to read the volume in a dedicated viewer like VESTA.

In [6]:

# Write it to a vasp file for inspection in external tools
vol.to_vasp_volume(sites, filename='/tmp/example.vasp')

Transitions¶

Using these sites, we can calculate the transitions and jumps of Lithium.

In [7]:

transitions = trajectory.transitions_between_sites(sites=sites, floating_specie='Li')

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[7], line 1
----> 1 transitions = trajectory.transitions_between_sites(sites=sites, floating_specie='Li')

File ~/local_projects/GEMDAT/src/gemdat/trajectory.py:535, in Trajectory.transitions_between_sites(self, sites, floating_specie, site_radius, site_inner_fraction)
    513 """Compute transitions between given sites for floating specie.
    514 
    515 Parameters
   (...)
    532 transitions: Transitions
    533 """
    534 from gemdat.transitions import Transitions
--> 535 return Transitions.from_trajectory(
    536     trajectory=self,
    537     sites=sites,
    538     floating_specie=floating_specie,
    539     site_radius=site_radius,
    540     site_inner_fraction=site_inner_fraction,
    541 )

File ~/local_projects/GEMDAT/src/gemdat/transitions.py:112, in Transitions.from_trajectory(cls, trajectory, sites, floating_specie, site_radius, site_inner_fraction)
    108 if site_radius is None:
    109     vibration_amplitude = SimulationMetrics(
    110         diff_trajectory).vibration_amplitude()
--> 112     site_radius = _compute_site_radius(
    113         trajectory=trajectory,
    114         sites=sites,
    115         vibration_amplitude=vibration_amplitude)
    117 if isinstance(site_radius, float):
    118     site_radius = {'': site_radius}

File ~/local_projects/GEMDAT/src/gemdat/transitions.py:422, in _compute_site_radius(trajectory, sites, vibration_amplitude)
    418             lines.append(f'{site_j.specie.name}({j}) {site_j.frac_coords}')
    420         msg = ''.join(lines)
--> 422         raise ValueError(
    423             f'Crystallographic sites are too close together (expected: >{site_radius*2:.4f}, '
    424             f'got: {min_dist:.4f} for {msg}')
    426 return site_radius

ValueError: Crystallographic sites are too close together (expected: >0.3255, got: 0.3355 for 
Too close:Li(26) [0.88636364 0.77272727 0.01515152] - Li(27) [0.88636364 0.75757576 0.98484848]
Too close:Li(27) [0.88636364 0.75757576 0.98484848] - Li(26) [0.88636364 0.77272727 0.01515152]

Note that this gives a warning that two sites (26, 27) are too close, so we delete one of them.

In [8]:

del sites[26]

transitions = trajectory.transitions_between_sites(sites=sites, floating_specie='Li')

Analyze and visualize jumps.

In [9]:

jumps = transitions.jumps(minimal_residence=0)
jumps.plot_jumps_vs_time(n_parts=5)

In [10]:

jumps.plot_jumps_3d()

In [ ]: