46 lines
1.2 KiB
Python
46 lines
1.2 KiB
Python
#%%
|
|
import numpy as np
|
|
|
|
import mdevaluate as md
|
|
import mdevaluate.extra.free_energy_landscape as fel
|
|
|
|
|
|
data_dir = "/data/skloth/python_packages/mdevaluate_examples/plots"
|
|
path_to_sim = (
|
|
"/data/skloth/sim/silica_pore/tip4p2005/D3_L6_S4.9_R0/T300_isochor/T250_nvt_long"
|
|
)
|
|
trajectory = md.open(path_to_sim, topology="run.tpr", trajectory="out/traj_full.xtc")
|
|
|
|
|
|
OW = trajectory.subset(atom_name="OW")
|
|
|
|
box = trajectory[0].box
|
|
box_voxels = (np.diag(box) // [0.05, 0.05, 0.05] + [1, 1, 1]) * [0.05, 0.05, 0.05]
|
|
occupation_matrix = fel.occupation_matrix(OW, skip=0, segments=1000)
|
|
radius_maxima = 0.05 * 3 ** (1 / 2) + 0.05 / 100
|
|
maxima_matrix = fel.find_maxima(
|
|
occupation_matrix,
|
|
box=box_voxels,
|
|
radius=radius_maxima,
|
|
pore_geometry="cylindrical"
|
|
)
|
|
maxima_matrix = fel.add_distances(maxima_matrix, "cylindrical", box / 2)
|
|
r_bins = np.arange(0, 1, 0.02)
|
|
distance_bins = np.arange(0.05, 1.55, 0.1)
|
|
energy_df = fel.distance_resolved_energies(
|
|
maxima_matrix, distance_bins, r_bins, box, "cylindrical", 225
|
|
)
|
|
result = fel.find_energy_maxima(energy_df, r_min=0.05, r_max=0.15)
|
|
|
|
|
|
#%%
|
|
print(sum(maxima_matrix["maxima"] == True))
|
|
|
|
#%%
|
|
energy_df = fel.distance_resolved_energies(
|
|
maxima_matrix, distance_bins, r_bins, box, temperature=250
|
|
)
|
|
|
|
# %%
|
|
print(energy_df)
|