Optimal path between sites¶
This notebook shows how to use GEMDAT to compute the optimal path between two sites.
It first computes the free energy on a discretized grid from the trajectory, then it finds the nearest sites of the structure to each point of the grid to map the space. Now the user can specify the start and end points of the path, and the metrics to compute the best path.
As input you will need:
- discretized grid resolution
- float value
- start point of the path
- a tuple representing the (x,y,z) coordinates of the start point
- end point of the path
- a tuple representing the (x,y,z) coordinates of the end point
- pathfinding algorithm, that can be one of the following:
- Dijkstra algorithm through
networkx.short_path
- Bellman-Ford algorithm through
networkx.short_path
- Minimal maximum energy
- Dijkstra but with exponentially scaled energy
- Dijkstra algorithm through
The resulting path can be analyzed and plotted in 3d, together with the energy profile along the path.
Loading data¶
Load simulation data and reference jump sites.
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))
/home/stef/python/gemdat/.venv/lib64/python3.12/site-packages/MDAnalysis/topology/TPRParser.py:161: DeprecationWarning: 'xdrlib' is deprecated and slated for removal in Python 3.13 import xdrlib
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.
from gemdat.path import free_energy_graph
volume = trajectory.to_volume(resolution=0.45)
F = volume.get_free_energy(temperature=trajectory.metadata['temperature'])
F_graph = free_energy_graph(F, max_energy_threshold=1e7, diagonal=True)
/home/stef/python/gemdat/src/gemdat/volume.py:329: RuntimeWarning: divide by zero encountered in log 'Boltzmann constant in eV/K'][0] * np.log(prob)
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.
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 4 4] -> [ 6 17 16]
0. Unweighted path¶
First we define a path as the direct (unweighted) path between the start and the stop point. This result is unrealistic and has low probability. It gives a baseline energy and shows the contrast with the energy optimized path.
from gemdat.path import optimal_path
path0 = optimal_path(
F_graph,
tuple(start_point),
tuple(stop_point),
method='simple',
)
path0
Path: (22, 4, 4) -> (6, 17, 16) Steps: 27 Total energy: 15.110 eV
from gemdat import plots
fig1 = volume.plot_3d(structure=sites, paths=path0)
fig1.show()
fig2 = plots.energy_along_path(paths=path0, structure=sites, volume=volume);
1. Dijkstra¶
path1 = optimal_path(
F_graph,
tuple(start_point),
tuple(stop_point),
method='dijkstra',
)
path1
Path: (22, 4, 4) -> (6, 17, 16) Steps: 27 Total energy: 14.295 eV
fig1 = volume.plot_3d(structure=sites, paths=path1)
fig1.show()
fig2 = plots.energy_along_path(paths=path1, structure=sites, volume=volume);
2. Bellman-Ford¶
path2 = optimal_path(
F_graph,
tuple(start_point),
tuple(stop_point),
method='bellman-ford',
)
path2
Path: (22, 4, 4) -> (6, 17, 16) Steps: 27 Total energy: 14.295 eV
fig1 = volume.plot_3d(structure=sites, paths=path2)
fig1.show()
fig2 = plots.energy_along_path(paths=path2, structure=sites, volume=volume);
3. Minimum maximum energy¶
This method tries to find the path with the lowest overall maximum energy. This can result in an path that is longer, but overall more efficient.
path3 = optimal_path(
F_graph,
tuple(start_point),
tuple(stop_point),
method='minmax-energy',
)
path3
Path: (22, 4, 4) -> (6, 17, 16) Steps: 27 Total energy: 14.295 eV
fig1 = volume.plot_3d(structure=sites, paths=path3)
fig1.show()
fig2 = plots.energy_along_path(paths=path3, structure=sites, volume=volume);
4. exponential Dijkstra¶
path4 = optimal_path(
F_graph,
tuple(start_point),
tuple(stop_point),
method='dijkstra-exp',
)
path4
Path: (22, 4, 4) -> (6, 17, 16) Steps: 27 Total energy: 14.295 eV
fig1 = volume.plot_3d(structure=sites, paths=path4)
fig1.show()
fig2 = plots.energy_along_path(paths=path4, structure=sites, volume=volume);
Comparison¶
As a final comparison, plot the energy of the paths identified with the different algorithms.
Notice from the figure above that dijkstra
, dijksta-exp
and bellman-ford
identify the same path, while minmax-energy
is able to find a path with a smaller highest-barrier, but sometimes a bit longer.
import matplotlib.pyplot as plt
fig, axs = plt.subplots(1, 1, figsize=(8, 6))
axs.plot(
path0.energy,
label=f'Direct: E={path0.total_energy:.3f} eV',
linestyle='solid',
marker='+',
)
axs.plot(
path1.energy,
label=f'Dijkstra: E={path1.total_energy:.3f} eV',
linestyle='-',
marker='o',
)
axs.plot(
path2.energy,
label=f'Bellman-Ford: E={path2.total_energy:.3f} eV',
linestyle='--',
marker='s',
)
axs.plot(
path3.energy,
label=f'Minmax: E={path3.total_energy:.3f} eV',
linestyle='-.',
marker='^',
)
axs.plot(
path4.energy,
label=f'Dijkstra-exp: E={path4.total_energy:.3f} eV',
linestyle=':',
marker='d',
)
axs.set_xlabel('Path steps')
axs.set_ylabel('Energy (eV)')
axs.legend()
axs.set_title(f'Energy along paths\nfrom{start_site}\nto {stop_site}');