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:
- 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
- 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.
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:361: 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 3 3] -> [ 5 16 16]
Set the number of paths we are looking for.
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$
dmin = 0.15
Calculate the paths¶
from gemdat.path import multiple_paths
paths = multiple_paths(
F_graph=F_graph,
start=tuple(start_point),
stop=tuple(stop_point),
n_paths=n_paths,
min_diff=dmin,
)
Plot the paths.
from gemdat import plots
fig1 = volume.plot_3d(paths=paths, structure=sites)
fig1.show()
fig2 = plots.energy_along_path(paths=paths, structure=sites, volume=volume);