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:

In [1]:

from gemdat.utils import VASPRUN

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

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.

In [2]:

from gemdat import Trajectory

trajectory = Trajectory.from_vasprun(VASPRUN)
trajectory

from gemdat import Trajectory trajectory = Trajectory.from_vasprun(VASPRUN) trajectory

Out[2]:

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.

In [3]:

trajectory.plot_displacement_per_atom()

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.

In [4]:

diff_trajectory = trajectory.filter('Li')
diff_trajectory

diff_trajectory = trajectory.filter('Li') diff_trajectory

Out[4]:

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)

In [5]:

diff_trajectory.plot_displacement_per_atom()

diff_trajectory.plot_displacement_per_atom()

In [6]:

diff_trajectory.plot_displacement_histogram()

diff_trajectory.plot_displacement_histogram()

In [7]:

diff_trajectory.plot_frequency_vs_occurence()

diff_trajectory.plot_frequency_vs_occurence()

In [8]:

diff_trajectory.plot_vibrational_amplitudes()

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.

In [9]:

from gemdat.io import load_known_material

sites = load_known_material('argyrodite', supercell=(2, 1, 1))
sites[0:10]

from gemdat.io import load_known_material sites = load_known_material('argyrodite', supercell=(2, 1, 1)) sites[0:10]

Out[9]:

[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.

In [10]:

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

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.

In [11]:

transitions.events

transitions.events

Out[11]:

	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.

In [12]:

occupancy = transitions.occupancy()
occupancy[0:10]

occupancy = transitions.occupancy() occupancy[0:10]

Out[12]:

[PeriodicSite: 48h (Li:0.45) (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.

In [13]:

jumps = transitions.jumps(minimal_residence=0)

jumps = transitions.jumps(minimal_residence=0)

In [14]:

jumps.plot_jumps_vs_time()

jumps.plot_jumps_vs_time()

In [15]:

jumps.plot_collective_jumps()

jumps.plot_collective_jumps()

In [16]:

jumps.plot_jumps_3d()

jumps.plot_jumps_3d()

Radial distribution function¶

The cell below shows how to calculate the radial distribution function using GEMDAT.

In [17]:

rdf_data = transitions.radial_distribution(
    floating_specie='Li',
)

for rdfs in rdf_data.values():
    fig = rdfs.plot()
    fig.show()

rdf_data = transitions.radial_distribution( floating_specie='Li', ) for rdfs in rdf_data.values(): fig = rdfs.plot() fig.show()

Output()