Implement energy landscapes and pathways from simulation

[SCiarella]

fix #144

This branch implements the analysis of optimal pathways between sites.
Only site-to-site pathways are implemented for now, but later it will be extended to identify percolating pathways and reproduce literature results.
The notebook 'paths.ipynb' shows the addition of this branch.

[stefsmeets]

Nice work! Let me know when this is finished, and I'll happily review this code.

Couple of things so far:

• Please put any documenting notebooks here to avoid ballooning this repo.
• Do we need a custom Dijksta implementation? Why not use a well-tested one from scipy or networkx?
• Can you add some tests for these new functions?
• What is this direct path for?

[SCiarella]

The branch has been updated and it is ready for review. In particular, the specific comments by @stefsmeets have been addressed in the following way:

• Please put any documenting notebooks here to avoid ballooning this repo.

Done

• Do we need a custom Dijksta implementation? Why not use a well-tested one from scipy or networkx?

Switched to networkx

• Can you add some tests for these new functions?

I am not an expert with testing, so I would like to discuss next time about the 'best practices' in this direction

• What is this direct path for?

It is used to show that it does not work, as mentioned by @tfamprikis in #144. I have now moved direct_path from the main gemdat package to path.ipynb, since its use is only pedagogical.

[v1kko]

Looks good to me if you add a test for the perculating pathway

[stefsmeets]

Hi @SCiarella , Nice work! Just made a pass through the code, some small comments from my side for now. I will have a look at the notebook tomorrow and the PR as a whole! 😃

[stefsmeets]

Maybe not for this PR, but for the long term we could consider that all methods like this return an instance of Volume with some flag that tells that it is a probability or free energy, and have some .units attribute to keep track of the units.

[stefsmeets]

I just checked the notebook, and I'm wondering if you can think a bit of how you want people to use this code. I find some parts not very obvious and I think it can be streamlined quite a bit.

Feel free to tackle this in a follow-up PR if that makes more sense.

Let me know what you think!

Free energy

This is the current api in this PR:

from gemdat.transitions import optimal_path, free_energy_graph

diff_volume = trajectory_to_volume(diff_trajectory, resolution=0.3)

F = diff_volume.get_free_energy(kBT=trajectory.metadata['temperature'])

F_graph = free_energy_graph(F, max_energy_threshold=1e7, diagonal=True)

start_point = (49, 26, 10)
end_point = (10, 0, 8)

path, path_energy = optimal_path(
    F_graph,
    start_point,
    end_point,
)

fig1 = plots.path_on_landscape(diff_volume, path, structure)
fig1.show()

fig2 = plots.energy_along_path(energy_path=path_energy)
fig2.show()

Playing around a little bit, this is what I would suggest we aim for:

• we can grab the temperature from the trajectory metadata
• F is a an instance of Volume or a subclass of Volume with energy related methods
• no need for any extra gemdat imports
• does not hide the networkx graph, so users can do whatever they want with it

import networkx as nx

diff_volume = trajectory.to_volume(resolution=0.3)

F = diff_volume.to_free_energy()
G = F.to_graph(max_energy_threshold=1e7, diagonal=True)

source = (49, 26, 10)
target = (10, 0, 8)

optimal_path = nx.shortest_path(G,
                                source=source,
                                target=target,
                                weight='weight')

fig1 = plots.path_on_landscape(diff_volume, optimal_path, structure)
fig1.show()

energy = [G.nodes[node]['energy'] for node in optimal_path]

fig2 = plots.energy_along_path(energy_path=energy)
fig2.show()

Percolation

Same goes for the percolation, my suggestion would be for an api that looks something like this:

• no need to import anything = more discoverable
• returning a dataclass of some sort makes it easier to extend the code
• hides the wrapping in a method

peaks = F.find_peaks()
best_path = F.find_best_perc_path(peaks, perc_xyz=(True, False, False))

print(f"Total Energy required: {best_path.total_energy_cost}")
print(f"Starting Point: {best_path.starting_point}")
print(f"Best Path: {best_path.best_path}")
print(f"Best Path Energy: {best_path.best_path_energy}")

fig1 = plots.path_on_landscape(diff_volume, structure, path=best_path)
fig1.show()

fig2 = plots.energy_along_path(path=best_path)
fig2.show()

[v1kko]

Looks good to me, I suggested a few small changes, but otherwise I think its good to go, any remaining comments can be converted into issues 🚀

[v1kko]

I think that a function for this single statement is not necessary, and the function can be ommited

[SCiarella]

Ok, I have removed the function

[SCiarella]

I have updated the branch as suggested. The remaining comments about API and usage can be discussed in another issue

[v1kko]

I have updated the branch as suggested. The remaining comments about API and usage can be discussed in another issue

Perfect, after you create the issue, let's merge this 🚀