Multiple paths¶

This notebook shows how to use GEMDAT to compare multiple optimal paths between two sites.

We follow the strategy discussed in the path example, based on the free energy discretization coupled with a pathfinding algorithm of choice. In addition now, you can specify how many paths you want to compare.

As input you will need:

1. discretized grid resolution
   • float value
2. start point of the path
   • a tuple representing the (x,y,z) coordinates of the start point
3. end point of the path
   • a tuple representing the (x,y,z) coordinates of the end point
4. number of paths to compare (n_paths)

The resulting n_paths paths can be analyzed and plotted in 3d, together with the energy profile.

Loading data¶

Load simulation data and reference jump sites.

In [1]:

from gemdat import Trajectory
from gemdat.io import load_known_material
from gemdat.utils import VASPRUN

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

trajectory = Trajectory.from_vasprun(VASPRUN)
diff_trajectory = trajectory.filter('Li')

sites = load_known_material('argyrodite', supercell=(2, 1, 1))

from gemdat import Trajectory from gemdat.io import load_known_material from gemdat.utils import VASPRUN # Use your own data: # VASPRUN = 'path/to/your/vasprun.xml' trajectory = Trajectory.from_vasprun(VASPRUN) diff_trajectory = trajectory.filter('Li') sites = load_known_material('argyrodite', supercell=(2, 1, 1))

Generate free energy graph¶

The next step is to generate the free energy graph. This uses the density volume, where every grid point (voxel) becomes a node in the graph.

Resolution is the minimum size of the voxels in Ångstrom.

Notice that there are points where it is numerically impossible to find mobile atoms. In such points the free energy is infinite and it can cause numerical instabilities.

In [2]:

volume = trajectory.to_volume(resolution=0.45)

free_energy = volume.get_free_energy(temperature=trajectory.metadata['temperature'])
free_energy

volume = trajectory.to_volume(resolution=0.45) free_energy = volume.get_free_energy(temperature=trajectory.metadata['temperature']) free_energy

/home/stef/python/gemdat/src/gemdat/volume.py:372: RuntimeWarning: divide by zero encountered in log
  * np.log(prob)

Out[2]:

Data: (44, 22, 22)
Lattice, abc:  19.873726   9.919447   9.916454
      angles:  90.214114  90.859135  89.950474
         pbc:       True       True       True
Label: volume

Define path¶

Select a start and an end point for the path.

For this example, we select two peaks of the density volume. The volume is converted into a pymatgen structure to make it easier to work with. Any integer voxel coordinate (start_point, stop_point) can be used to define the path.

In [3]:

peaks = volume.to_structure()

start_site = peaks[0]
stop_site = peaks[100]

print(start_site, '->', stop_site)

start_point = volume.site_to_voxel(start_site)
stop_point = volume.site_to_voxel(stop_site)

print(start_point, '->', stop_point)

peaks = volume.to_structure() start_site = peaks[0] stop_site = peaks[100] print(start_site, '->', stop_site) start_point = volume.site_to_voxel(start_site) stop_point = volume.site_to_voxel(stop_site) print(start_point, '->', stop_point)

[10.1513076   1.57875625  1.49884936] X0+ -> [2.496801   7.49305245 7.33775091] X0+
[22  3  3] -> [ 5 16 16]

Set the number of paths we are looking for.

In [4]:

n_paths = 3

n_paths = 3

Before proceding, we have to consider that it is possible to have multiple small variations of the same path labelled as different paths. For this reason we have to introduce a difference score: $$\mathcal{D}(\mathcal{p}_1,\mathcal{p}_2) = \frac{\sum_{\mathrm{site}_i \in \mathcal{p}_1} [0\; \mathrm{ if }\;\mathrm{site}_i \in \mathcal{p}_2\; \mathrm{else }\;1 ] }{ \mathrm{Lenght}(p_1) },$$ that corresponds to the percentage of sites of $\mathcal{p}_1$ that are not in $\mathcal{p}_2$.

We say that two paths are different if $\mathcal{D}> dmin$

In [5]:

dmin = 0.15

dmin = 0.15

Calculate the paths¶

The first path is the optimal path, the other paths are small variations of the optimal path.

In [6]:

best_path, *other_paths = free_energy.optimal_n_paths(
    start=start_point,
    stop=stop_point,
    n_paths=n_paths,
    min_diff=dmin,
)
best_path

best_path, *other_paths = free_energy.optimal_n_paths( start=start_point, stop=stop_point, n_paths=n_paths, min_diff=dmin, ) best_path

Out[6]:

Path: (22, 3, 3) -> (5, 16, 16)
Steps: 28
Total energy: 14.507 eV
Dimensions: (44, 22, 22)

Plot the paths.

In [7]:

fig1 = volume.plot_3d(paths=(best_path, *other_paths), structure=sites)
fig1.show()

fig2 = best_path.plot_energy_along_path(structure=sites, other_paths=other_paths)
fig2.show()

fig1 = volume.plot_3d(paths=(best_path, *other_paths), structure=sites) fig1.show() fig2 = best_path.plot_energy_along_path(structure=sites, other_paths=other_paths) fig2.show()