Introduction to GEMDAT¶
This notebook as an introduction to GEMDAT.
Gemdat is a Python library for the analysis of diffusion in solid-state electrolytes from Molecular Dynamics simulations. Gemdat is built on top of Pymatgen, making it straightforward to integrate it into your Pymatgen-based workflows.
To start using GEMDAT, you will need some molecular dynamics data, e.g. from VASP. We import here some default data, but feel free to use your own:
from gemdat.utils import VASPRUN
# Use your own data:
# VASPRUN = 'path/to/your/vasprun.xml'
Trajectory¶
The entry point to GEMDAT is via a Trajectory. It is an extension of the pymatgen Trajectory class.
from gemdat import Trajectory
trajectory = Trajectory.from_vasprun(VASPRUN)
trajectory
Full Formula (Li48 P8 S40 Br8) Reduced Formula: Li6PS5Br abc : 19.873726 9.919447 9.916454 angles: 90.214114 90.859135 89.950474 pbc : True True True Constant lattice (True) Sites (104) Time steps (5000)
With a trajectory, you can plot various statistics.
trajectory.plot_displacement_per_atom()
If we are interested in a single species, such as how Lithium diffuses through the system, we can filter it from the trajectory. This creates a new Trajectory 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)
diff_trajectory.plot_displacement_per_atom()
diff_trajectory.plot_displacement_histogram()
diff_trajectory.plot_frequency_vs_occurence()
diff_trajectory.plot_vibrational_amplitudes()
Transitions and jumps¶
The trajectory can be used to find transitions between specific sites, cavities or clusters.
Sometimes the sites that the diffusing atom moves between are known in advance, but GEMDAT also includes tools to find site clusters in the trajectory. Sites can be loaded from a CIF file or from a small internal database of known materials. In the example below we load a argyrodite with a (2,1,1) supercell.
Note that sites
is a pymatgen Structure.
from gemdat.io import load_known_material
sites = load_known_material('argyrodite', supercell=(2, 1, 1))
sites[0:10]
[PeriodicSite: 48h (Li) (1.816, 1.816, 0.2382) [0.0915, 0.183, 0.024], PeriodicSite: 48h (Li) (11.74, 1.816, 0.2382) [0.5915, 0.183, 0.024], PeriodicSite: 48h (Li) (1.816, 6.778, 5.2) [0.0915, 0.683, 0.524], PeriodicSite: 48h (Li) (11.74, 6.778, 5.2) [0.5915, 0.683, 0.524], PeriodicSite: 48h (Li) (6.778, 1.816, 5.2) [0.3415, 0.183, 0.524], PeriodicSite: 48h (Li) (16.7, 1.816, 5.2) [0.8415, 0.183, 0.524], PeriodicSite: 48h (Li) (6.778, 6.778, 0.2382) [0.3415, 0.683, 0.024], PeriodicSite: 48h (Li) (16.7, 6.778, 0.2382) [0.8415, 0.683, 0.024], PeriodicSite: 48h (Li) (9.686, 1.816, 8.108) [0.488, 0.183, 0.817], PeriodicSite: 48h (Li) (19.61, 1.816, 8.108) [0.988, 0.183, 0.817]]
The first step is to find transitions between sites for the floating species, here 'Li'.
This is done using a site radius, in the example below set to 1.0
. This is an optional value, and if left out GEMDAT tries to calculate a suitable site radius based on the vibration amplitude of your floating species.
If any Li is within 1.0 of one of the given sites, it is considered to be at that site.
transitions = trajectory.transitions_between_sites(
sites=sites,
floating_specie='Li',
site_radius=1.0,
)
The transitions are collected in a series of transition events, which are stored in a DataFrame. This tracks the index the site is at, and when it leaves the site. -1 means is not near a site, e.g. transitioning between sites.
transitions.events
atom index | start site | destination site | start inner site | destination inner site | time | |
---|---|---|---|---|---|---|
0 | 0 | 94 | -1 | 94 | -1 | 431 |
1 | 0 | -1 | 94 | -1 | 94 | 442 |
2 | 0 | 94 | -1 | 94 | -1 | 473 |
3 | 0 | -1 | 94 | -1 | 94 | 495 |
4 | 0 | 94 | -1 | 94 | -1 | 502 |
... | ... | ... | ... | ... | ... | ... |
2442 | 47 | 47 | -1 | 47 | -1 | 4917 |
2443 | 47 | -1 | 55 | -1 | 55 | 4925 |
2444 | 47 | 55 | -1 | 55 | -1 | 4945 |
2445 | 47 | -1 | 55 | -1 | 55 | 4959 |
2446 | 47 | 55 | -1 | 55 | -1 | 4969 |
2447 rows × 6 columns
The transitions can be used to calculate the occupancy: the fraction of time an atom spends at a site.
occupancy = transitions.occupancy()
occupancy[0:10]
[PeriodicSite: 48h (Li:0.450) (1.816, 1.816, 0.2382) [0.0915, 0.183, 0.024], PeriodicSite: 48h (Li:0.361) (11.74, 1.816, 0.2382) [0.5915, 0.183, 0.024], PeriodicSite: 48h (Li:0.393) (1.816, 6.778, 5.2) [0.0915, 0.683, 0.524], PeriodicSite: 48h (Li:0.432) (11.74, 6.778, 5.2) [0.5915, 0.683, 0.524], PeriodicSite: 48h (Li:0.497) (6.778, 1.816, 5.2) [0.3415, 0.183, 0.524], PeriodicSite: 48h (Li:0.682) (16.7, 1.816, 5.2) [0.8415, 0.183, 0.524], PeriodicSite: 48h (Li:0.357) (6.778, 6.778, 0.2382) [0.3415, 0.683, 0.024], PeriodicSite: 48h (Li:0.393) (16.7, 6.778, 0.2382) [0.8415, 0.683, 0.024], PeriodicSite: 48h (Li:0.556) (9.686, 1.816, 8.108) [0.488, 0.183, 0.817], PeriodicSite: 48h (Li:0.309) (19.61, 1.816, 8.108) [0.988, 0.183, 0.817]]
Jumps¶
Transitions are not necessarily equal to jumps. To classify transitions as jumps, meaning when an atom has left its original site and arrived at a new site, GEMDAT uses the Jumps class.
GEMDAT uses a minimal residence time (number of time steps) that the jumping atom must stay at the destination site for the jump to be completed and counted.
jumps = transitions.jumps(minimal_residence=0)
jumps.plot_jumps_vs_time()
jumps.plot_collective_jumps()
jumps.plot_jumps_3d()