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.
import plotly.io as pio
# Render plotly figures as static HTML so they display on the docs website
pio.renderers.default = 'plotly_mimetype+notebook_connected'
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.
# Create a new trajectory object that only consists of Lithium
diff_trajectory = trajectory.filter('Li')
diff_trajectory
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.
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.
sites = vol.to_structure(specie='Li', background_level=0.1)
sites[0:10]
[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 (16.94, 6.331, 0.2842) [0.8526, 0.6377, 0.0427], PeriodicSite: Li (17.62, 7.61, -0.1168) [0.8864, 0.7667, 0.00303], 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], PeriodicSite: Li (12.79, 5.601, 8.527) [0.647, 0.5657, 0.8707], PeriodicSite: Li (12.77, 6.666, 9.699) [0.6461, 0.6732, 0.9892], PeriodicSite: Li (16.78, 1.992, 4.979) [0.8464, 0.2011, 0.5152], PeriodicSite: Li (17.15, 0.9282, 5.55) [0.8652, 0.09394, 0.5727]]
Verify the generated sites by plotting them on the volume.
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.
# Write it to a vasp file for inspection in external tools
vol.to_vasp_volume(sites, filename='/tmp/example.vasp')
<pymatgen.io.vasp.outputs.VolumetricData at 0x7b0f72c18440>
Transitions¶
Using these sites, we can calculate the transitions and jumps of Lithium.
transitions = trajectory.transitions_between_sites(sites=sites, floating_specie='Li')
Note that this gives a warning that two sites (26, 27) are too close, so we delete one of them.
del sites[26]
transitions = trajectory.transitions_between_sites(sites=sites, floating_specie='Li')
Analyze and visualize jumps.
jumps = transitions.jumps(minimal_residence=0)
jumps.plot_jumps_vs_time(n_parts=5)
jumps.plot_jumps_3d()