Finished data management example

examples/data_management.py

@@ -1,15 +1,74 @@
-DATA_DIR = "/data/skloth/results/workshop_data"
-np.save(f"{DATA_DIR}/msd.npy", result_msd_nojump)
-print(np.load(f"{DATA_DIR}/msd.npy"))
+# %%
+import numpy as np
+import matplotlib.pyplot as plt
 
+import mdevaluate as md
+
+DATA_DIR = "/data/skloth/python_packages/mdevaluate_examples/data"
+path_to_sim = "/data/skloth/sim/IL_water/C4MIM_BF4_x50/T300_nvt"
+trajectory = md.open(path_to_sim, topology="run.tpr", trajectory="out/traj_full.xtc")
+
+oxygen_water = trajectory.subset(atom_name="OW", residue_name="SOL")
+cation_indices = trajectory.subset(residue_name="amim").atom_subset.indices
+cation_com = md.coordinates.center_of_masses(trajectory, atom_indices=cation_indices)
+anion_com = trajectory.subset(atom_name="B", residue_name="BF4")
+
+
+# %%
+time, msd_oxygen_water = md.correlation.shifted_correlation(
+    md.correlation.msd, oxygen_water.nojump, segments=100, skip=0.1, average=True
+)
+time, msd_cation_com = md.correlation.shifted_correlation(
+    md.correlation.msd, cation_com.nojump, segments=100, skip=0.1, average=True
+)
+time, msd_anion_com = md.correlation.shifted_correlation(
+    md.correlation.msd, anion_com.nojump, segments=100, skip=0.1, average=True
+)
+
+plt.figure()
+plt.plot(time, msd_oxygen_water, ".", label="Water")
+plt.plot(time, msd_cation_com, ".", label="Cation")
+plt.plot(time, msd_anion_com, ".", label="Anion")
+plt.legend()
+plt.xscale("log")
+plt.yscale("log")
+plt.xlabel(r"$t$ in ps", fontsize=16)
+plt.ylabel(r"$\langle r^2\rangle(t)$ in nm$^2$", fontsize=16)
+plt.show()
+plt.close()
+
+# %%
 import pandas as pd
-msd_df = pd.DataFrame({"t": result_msd[0], "msd_water": result_msd[1]})
+msd_df = pd.DataFrame(
+    {"t": time,
+     "msd_water": msd_oxygen_water,
+     "msd_cation": msd_cation_com,
+     "msd_anion": msd_anion_com
+     }
+)
 print(msd_df)
-msd_df["T"] = 300
 msd_df["cation"] = "C4mim"
 msd_df["anion"] = "DCA"
+msd_df["x_w"] = 50
+msd_df["T"] = 300
 print(msd_df)
+
+
+# %%
+def add_description(dataframe, cation, anion, x_w, T):
+    dataframe["cation"] = cation
+    dataframe["anion"] = anion
+    dataframe["x_w"] = x_w
+    dataframe["T"] = T
+    return dataframe
+
+
+msd_df = add_description(msd_df, "C4mim", "DCA", 50, 300)
+print(msd_df)
+
+#%%
 msd_df.to_hdf(f"{DATA_DIR}/result.h5", key="msd")
+
 msd_df = pd.read_hdf(f"{DATA_DIR}/result.h5", key="msd")
-print(msd_df.info())
-md.utils.cleanup_h5(f"{DATA_DIR}/result.h5")
+
+md.utils.cleanup_h5(f"{DATA_DIR}/result.h5")

examples/dynamic_properties.py

@@ -0,0 +1,57 @@
+# %%
+from functools import partial
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import mdevaluate as md
+
+data_dir = "/data/skloth/python_packages/mdevaluate_examples/plots"
+path_to_sim = "/data/skloth/sim/IL_water/C4MIM_BF4_x50/T300_nvt"
+trajectory = md.open(path_to_sim, topology="run.tpr", trajectory="out/traj_full.xtc")
+
+oxygen_water = trajectory.subset(atom_name="OW", residue_name="SOL")
+
+# %%
+print(md.correlation.msd(trajectory[0], trajectory[-1]))
+
+result_msd_pbc = md.correlation.shifted_correlation(
+    md.correlation.msd, oxygen_water, segments=100, skip=0.1, average=True
+)
+result_msd_nojump = md.correlation.shifted_correlation(
+    md.correlation.msd, oxygen_water.nojump, segments=100, skip=0.1, average=True
+)
+
+plt.figure()
+plt.plot(result_msd_pbc[0], result_msd_pbc[1], ".", label="pbc")
+plt.plot(result_msd_nojump[0], result_msd_nojump[1], ".", label="nojump")
+plt.legend()
+plt.xscale("log")
+plt.yscale("log")
+plt.xlabel(r"$t$ in ps", fontsize=16)
+plt.ylabel(r"$\langle r^2\rangle(t)$ in nm$^2$", fontsize=16)
+plt.show()
+plt.close()
+
+# %%
+cation = trajectory.subset(residue_name="amim")
+cation_indices = cation.atom_subset.indices
+cation_com = md.coordinates.center_of_masses(trajectory, atom_indices=cation_indices)
+
+result_msd_cation = md.correlation.shifted_correlation(
+    md.correlation.msd, cation.nojump, segments=100, skip=0.1, average=True
+)
+result_msd_cation_com = md.correlation.shifted_correlation(
+    md.correlation.msd, cation_com.nojump, segments=100, skip=0.1, average=True
+)
+
+plt.figure()
+plt.plot(result_msd_cation[0], result_msd_cation[1], ".", label="cation all atoms")
+plt.plot(result_msd_cation_com[0], result_msd_cation_com[1], ".", label="cation com")
+plt.legend()
+plt.xscale("log")
+plt.yscale("log")
+plt.xlabel(r"$t$ in ps", fontsize=16)
+plt.ylabel(r"$\langle r^2\rangle(t)$ in nm$^2$", fontsize=16)
+plt.show()
+plt.close()

examples/free_energy_landscape.py

@@ -1,3 +1,4 @@
+#%%
 import numpy as np
 
 import mdevaluate as md
@@ -13,14 +14,32 @@ trajectory = md.open(path_to_sim, topology="run.tpr", trajectory="out/traj_full.
 
 OW = trajectory.subset(atom_name="OW")
 
-box = np.diag(trajectory[0].box)
-box_voxels = (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=len(OW))
-maxima_matrix = fel.find_maxima(occupation_matrix, box=box_voxels, edge_length=0.05)
+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, 2, 0.02)
+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, 225
+    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)

examples/static_properties.py

@@ -0,0 +1,53 @@
+# %%
+from functools import partial
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+import mdevaluate as md
+
+data_dir = "/data/skloth/python_packages/mdevaluate_examples/plots"
+path_to_sim = "/data/skloth/sim/IL_water/C4MIM_BF4_x50/T300_nvt"
+trajectory = md.open(path_to_sim, topology="run.tpr", trajectory="out/traj_full.xtc")
+
+oxygen_water = trajectory.subset(atom_name="OW", residue_name="SOL")
+
+# %%
+bins = np.arange(0, 1, 0.01)
+result_rdf = md.distribution.rdf(oxygen_water[-1], bins=bins)
+plt.figure()
+plt.plot(bins[:-1], result_rdf)
+plt.xlabel(r"$r$ in nm", fontsize=16)
+plt.ylabel(r"$g(r)$", fontsize=16)
+plt.savefig(f"{data_dir}/rdf.png", dpi=300, bbox_inches="tight")
+plt.show()
+plt.close()
+
+# %%
+bins = np.arange(0, 1, 0.01)
+result_rdf = md.distribution.time_average(
+    partial(md.distribution.rdf, bins=bins), oxygen_water, segments=1000, skip=0.1
+)
+plt.figure()
+plt.plot(bins[:-1], result_rdf)
+plt.xlabel(r"$r$ in nm", fontsize=16)
+plt.ylabel(r"$g(r)$", fontsize=16)
+plt.show()
+plt.close()
+
+# %%
+boron_anion = trajectory.subset(atom_name="B", residue_name="BF4")
+bins = np.arange(0, 1, 0.01)
+result_rdf = md.distribution.time_average(
+    partial(md.distribution.rdf, bins=bins),
+    oxygen_water,
+    coordinates_b=boron_anion,
+    segments=1000,
+    skip=0.1,
+)
+plt.figure()
+plt.plot(bins[:-1], result_rdf)
+plt.xlabel(r"$r$ in nm", fontsize=16)
+plt.ylabel(r"$g(r)$", fontsize=16)
+plt.show()
+plt.close()

examples/subsets.py

@@ -0,0 +1,25 @@
+import mdevaluate as md
+import numpy as np
+
+path_to_sim = "/data/skloth/sim/IL_water/C4MIM_BF4_x50/T300_nvt"
+trajectory = md.open(path_to_sim, topology="run.tpr", trajectory="out/traj_full.xtc")
+
+print(trajectory[0].atom_names)
+print(trajectory[0].residue_names)
+
+print(np.unique(trajectory[0].atom_names))
+print(np.unique(trajectory[0].residue_names))
+
+print(np.unique(trajectory[0].atom_names[trajectory[0].residue_names == "amim"]))
+print(np.unique(trajectory[0].atom_names[trajectory[0].residue_names == "BF4"]))
+print(np.unique(trajectory[0].atom_names[trajectory[0].residue_names == "SOL"]))
+
+carbon_atoms = trajectory.subset(atom_name="C")
+hydrogen_atoms = trajectory.subset(atom_name="H")
+cation_chain_end_atoms = trajectory.subset(atom_name="CT")
+
+cation_atoms = trajectory.subset(residue_name="amim")
+anion_atoms = trajectory.subset(residue_name="BF4")
+water_atoms = trajectory.subset(residue_name="SOL")
+
+hydrogen_water = trajectory.subset(atom_name="HW", residue_name="SOL")

examples/vectors.py

@@ -0,0 +1,59 @@
+import mdevaluate as md
+import numpy as np
+
+path_to_sim = "/data/skloth/sim/IL_water/C4MIM_BF4_x50/T300_nvt"
+trajectory = md.open(path_to_sim, topology="run.tpr", trajectory="out/traj_full.xtc")
+
+atom_indices_a = trajectory.subset(
+    atom_name="HW.", residue_name="SOL"
+).atom_subset.indices
+atom_indices_b = trajectory.subset(
+    atom_name="OW", residue_name="SOL"
+).atom_subset.indices
+atom_indices_b = np.vstack([atom_indices_b] * 2).T.reshape((-1,))
+OH_vectors = md.coordinates.vectors(
+    trajectory,
+    atom_indices_a=atom_indices_a,
+    atom_indices_b=atom_indices_b,
+    normed=True,
+)
+print(OH_vectors[-1])
+
+atom_indices_a = trajectory.subset(
+    atom_name="B", residue_name="BF4"
+).atom_subset.indices
+atom_indices_b = trajectory.subset(
+    atom_name="F.", residue_name="BF4"
+).atom_subset.indices
+atom_indices_a = np.vstack([atom_indices_a] * 4).T.reshape((-1,))
+BF_vectors = md.coordinates.vectors(
+    trajectory,
+    atom_indices_a=atom_indices_b,
+    atom_indices_b=atom_indices_a,
+    normed=True,
+)
+print(BF_vectors[-1])
+
+atom_indices = trajectory.subset(residue_name="SOL").atom_subset.indices
+water_dipole_vectors = md.coordinates.dipole_vector(
+    trajectory, atom_indices=atom_indices, normed=True
+)
+print(water_dipole_vectors[-1])
+
+atom_indices_a = trajectory.subset(
+    atom_name="NA1", residue_name="amim"
+).atom_subset.indices
+atom_indices_b = trajectory.subset(
+    atom_name="CR", residue_name="amim"
+).atom_subset.indices
+atom_indices_c = trajectory.subset(
+    atom_name="CW1", residue_name="amim"
+).atom_subset.indices
+cation_ring_normal = md.coordinates.normal_vectors(
+    trajectory,
+    atom_indices_a=atom_indices_a,
+    atom_indices_b=atom_indices_b,
+    atom_indices_c=atom_indices_c,
+    normed=True,
+)
+print(cation_ring_normal[-1])