78 lines
2.2 KiB
Python
78 lines
2.2 KiB
Python
from functools import partial
|
|
|
|
import numpy as np
|
|
import matplotlib.pyplot as plt
|
|
from matplotlib import cm
|
|
|
|
import mdevaluate as md
|
|
|
|
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_short"
|
|
)
|
|
trajectory = md.open(path_to_sim, topology="run.tpr", trajectory="out/traj_full.xtc")
|
|
|
|
oxygen_water = trajectory.subset(atom_name="OW", residue_name="SOL")
|
|
|
|
time, result_all = md.correlation.shifted_correlation(
|
|
partial(md.correlation.isf, q=22.7), oxygen_water, segments=100, skip=0.1
|
|
)
|
|
time, result_wall = md.correlation.shifted_correlation(
|
|
partial(md.correlation.isf, q=22.7),
|
|
oxygen_water,
|
|
selector=partial(md.coordinates.selector_radial_cylindrical, r_min=1.0, r_max=1.5),
|
|
segments=100,
|
|
skip=0.1,
|
|
)
|
|
time, result_center = md.correlation.shifted_correlation(
|
|
partial(md.correlation.isf, q=22.7),
|
|
oxygen_water,
|
|
selector=partial(md.coordinates.selector_radial_cylindrical, r_min=0.0, r_max=0.5),
|
|
segments=100,
|
|
skip=0.1,
|
|
)
|
|
|
|
plt.figure()
|
|
plt.plot(time, result_all, "k-", label="all")
|
|
plt.plot(time, result_wall, "r.", label="wall")
|
|
plt.plot(time, result_center, "b.", label="center")
|
|
plt.legend()
|
|
plt.xscale("log")
|
|
plt.xlabel(r"$t$ / ps")
|
|
plt.ylabel(r"S_q(t)")
|
|
plt.savefig(f"{data_dir}/selector.png", dpi=300, bbox_inches="tight")
|
|
plt.show()
|
|
|
|
|
|
def multi_radial_selector(atoms, bins):
|
|
indices = []
|
|
for i in range(len(bins) - 1):
|
|
index = md.coordinates.selector_radial_cylindrical(
|
|
atoms, r_min=bins[i], r_max=bins[i + 1]
|
|
)
|
|
indices.append(index)
|
|
return indices
|
|
|
|
|
|
bins = np.arange(0.0, 1.6, 0.1)
|
|
r = (bins[:-1] + bins[1:]) / 2
|
|
|
|
time, results = md.correlation.shifted_correlation(
|
|
partial(md.correlation.isf, q=22.7),
|
|
oxygen_water,
|
|
selector=partial(multi_radial_selector, bins=bins),
|
|
segments=100,
|
|
skip=0.1,
|
|
)
|
|
|
|
c = [cm.plasma(i) for i in np.linspace(0, 1, len(r))]
|
|
plt.figure()
|
|
for i, result in enumerate(results):
|
|
plt.plot(time, result, "-", c=c[i], label=round(r[i], 2))
|
|
plt.legend(title=r"$r$ / nm", ncols=2)
|
|
plt.xscale("log")
|
|
plt.xlabel(r"$t$ / ps")
|
|
plt.ylabel(r"S_q(t)")
|
|
plt.savefig(f"{data_dir}/multi_selector.png", dpi=300, bbox_inches="tight")
|
|
plt.show()
|