From 55ad026e7b15b077d65248ef999856eecfcc465b Mon Sep 17 00:00:00 2001
From: Sebastian Kloth <sebastian.kloth@pkm.tu-darmstadt.de>
Date: Sat, 10 Feb 2024 14:53:47 +0100
Subject: [PATCH] Worked out example for static properties

---
 examples/static_properties.py | 42 +++++++++++++++++++++++++++++++----
 1 file changed, 38 insertions(+), 4 deletions(-)

diff --git a/examples/static_properties.py b/examples/static_properties.py
index 6b0ea8e..5627e26 100644
--- a/examples/static_properties.py
+++ b/examples/static_properties.py
@@ -1,19 +1,53 @@
-bins = np.linspace(0, 1, 101)
+# %%
+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()
 
-bins = np.linspace(0, 1, 101)
+# %%
+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, segments=100, skip=0.1
+    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()
\ No newline at end of file
+plt.close()