Correlation function calculations

README.md

@@ -10,7 +10,7 @@ micromamba activate soluenv
 rm -rf subprojects 
 git restore subprojects
 meson setup build --wipe
-pip install -e . --no-build-isolation
+pip install -e . #--no-build-isolation
 ```
 
 ## Usage

environment.yml

@@ -42,8 +42,12 @@ dependencies:
   - pybind11
   - meson-python
   - h5py
+  - scipy
   # Examples
   - gdown
+  - matplotlib
   - pip:
     # works for regular pip packages
     - lammps-logfile
+    # Required for plotting
+    - spirit-extras

examples/hydrogen_bond_tcf/hbond_tcf.py

@@ -1,16 +1,70 @@
 from pathlib import Path
 import gdown
+import soluanalysis as solu
+from soluanalysis.io import read_lammps_dump
+import numpy as np
 
 hbond_example_dir = Path(__file__).resolve().parent
 input_folder = hbond_example_dir / "input"
 # Create the folder if it does not exist 
 input_folder.mkdir(parents=True, exist_ok=True)
+# Create an output folder if it does not exist
+output_folder = hbond_example_dir / "output"
+output_folder.mkdir(parents=True, exist_ok=True)
+output_file =  output_folder / 'tcf.txt'
 
 # URL for the trajectory file (571 MB)
 file_id = '15sv3cc_Mx281C1RHDTzUiqMoVBjwcWNu'  
 url = f"https://drive.google.com/uc?id={file_id}"
-# Output file path 
-output_file = input_folder / 'dump-tcf.lammpstrj'
+# File path for the trajectory 
+input_file = input_folder / 'dump-tcf.lammpstrj'
 # Download the file into the 'input' folder (if it does not already exist)
-if not output_file.exists():
-    gdown.download(url, str(output_file), quiet=False)
+if not input_file.exists():
+    gdown.download(url, str(input_file), quiet=False)
+
+# Read in the trajectory
+systems, timesteps = read_lammps_dump(input_file, ':')
+print("Trajectory read.\n")
+# General system information
+fe_type = 3
+o_type = 1
+h_type = 2
+cl_type = 4
+# Cutoffs, types etc needed for hydrogen bonds
+donor_atom_types = [o_type]
+acceptor_atom_types = [cl_type, o_type]
+h_atom_types = [h_type]
+donor_acceptor_cutoff = 3.2
+max_angle_deg = 30  # in degrees 
+# List that will hold UndirectedNetwork objects corresponding to hydrogen bonds
+networks = [] 
+
+# Loop through the trajectory 
+for system in systems:
+    n_atoms = system.n_atoms()
+    network = solu.graphlib.UndirectedNetwork(n_atoms) # Create the network
+    # We will ignore hydrogens so there is no need for intramolecular hydrogen bonds
+    solu.james.add_hbonds(
+        network,
+        system,
+        donor_atom_types,
+        acceptor_atom_types,
+        h_atom_types,
+        donor_acceptor_cutoff,
+        max_angle_deg,
+        True,
+    )
+    networks.append(network)
+print("Network list generated\n")
+# Now that you have a list of UndirectedNetwork objects and the timesteps, 
+# the time autocorrelation function can be calculated
+
+# Calculate the time correlation function (using default values of start_t0 etc)
+tau_values, tcf_avg, tcf_error = solu.james.time_correlation_function(
+        networks, timesteps, 0, 1, 1, None
+    )
+print("Calculated the time correlation function\n")
+# Write these out to a CSV 
+# Write out to file  
+header_string = 'tau\ttcf\ttcf_stderr'
+np.savetxt(output_file, np.column_stack((tau_values, tcf_avg, tcf_error)), delimiter=' ', header = header_string)

examples/hydrogen_bond_tcf/lineplot.py

@@ -0,0 +1,143 @@
+from spirit_extras.plotting import Paper_Plot
+import numpy as np
+from pathlib import Path
+import matplotlib.pyplot as plt
+from PIL import Image
+import os
+from os.path import isfile, join
+import matplotlib.patheffects as path_effects
+import matplotlib.lines as mlines
+from matplotlib.legend_handler import HandlerTuple
+from soluanalysis.misc import fit_biexponential
+
+def insert_inset(image_path, axis, rel_height, 
+    rel_width, margin_x, margin_y, x_align, y_align):
+    image = pplot.open_image(image_path)  # Read the image as a numpy array
+    # image = pplot.crop_to_content(image)  # crop the image to content
+    # image = pplot.crop(image, 950,950)
+
+    image = pplot.replace_background_color(image, [0,0,0,0], background_color=None)
+
+    # Create an axis to hold the inset
+    # The inset shall be located in the top right corner
+    ax_inset = pplot.create_inset_axis(
+        containing_ax=axis,
+        rel_height=rel_height,
+        rel_width=rel_width,
+        margin_x=margin_x,
+        margin_y=margin_y,
+        x_align=x_align,
+        y_align=y_align,
+    )
+    return ax_inset, image
+
+# Figure attributes 
+hbond_example_dir = Path(__file__).resolve().parent
+input_dir = hbond_example_dir / "output"
+fname = "tcf.txt"
+
+# Lengths in inches in general 
+CM = Paper_Plot.cm  # Use this to give lengths in cm...
+
+# Params 
+params = {
+"font.size": 8,
+"font.family": ("Arial","sans-serif"),
+"mathtext.fontset": "dejavuserif",
+"xtick.labelsize": 7,
+"ytick.labelsize": 7,
+"axes.labelsize": 8
+}
+
+# For legend labels
+font = {'family': ("Arial","sans-serif"),
+        'weight': 'normal',
+        'size': 7.5,
+        }
+
+# Either set the absolute height here or set the aspect ratio
+# inside apply_absolute_margins
+
+pplot = Paper_Plot(
+    width=3.25, height=2.75,nrows=1, ncols=1, rcParams=params
+) 
+
+# Vertical margin: bottom and then top 
+# horizontal margin: left and then right
+# Golden ratio aspect ratio=1.618 
+# hspace -> height space 
+
+pplot.apply_absolute_margins(
+    aspect_ratio=None,
+    abs_horizontal_margins=[1.3 * CM, 0.12 * CM],
+    abs_vertical_margins=[1.0 * CM, 0.5 * CM],
+    abs_wspace=0.0 * CM,
+    abs_hspace=0.0 * CM,
+)
+
+print(pplot.info_string())
+
+# Get the figure and gridpsec objects
+fig = pplot.fig()
+gs = pplot.gs()
+
+# DATA 
+# Get the data (with counterions)
+infile = os.path.join(input_dir,fname)
+data = np.loadtxt(infile,skiprows=1) # Load the text file 
+tau_val = np.array(data[:,0]) # Tau values in fs
+tau_val = tau_val/1000 # Tau values in ps 
+tcf_val =  np.array(data[:,1]) # Time correlation function values  
+tcf_val_err = np.array(data[:,2]) # Error bars for the TCF  
+
+# ---------------------------------------------
+# Curve fitting 
+# initial guess A, tau1, tau2
+params, fit_t, fit_ac, lifetime = fit_biexponential(tau_val, tcf_val, [0.5, 1, 2])
+print("Lifetime is ", lifetime, "ps \n")
+A, tau1, tau2 = params
+# If A+B is around 1.01, the fit is bad even if it looks okay
+print(f"The time constants are {tau1} ps and {tau2} ps\n")
+# ---------------------------------------------
+ax1 = fig.add_subplot(gs[0,0])
+ax1.set_xlabel(r'$\tau$ (ps)')
+ax1.set_ylabel(r'$\mathrm{C_{HB}} (\tau)$',labelpad=7)  # we already handled the x-label with ax1
+
+# Data points and line for non-octahedral state 
+# (tcf_1,) = ax1.plot(tau_val, tcf_val, marker=".", markersize=6,
+#     color="orangered",markeredgecolor='black',markeredgewidth=0.6,zorder=4,
+#     label=r'Data')
+(tcf_1,) = ax1.plot(tau_val, tcf_val, marker=None,
+    color="orangered",zorder=4,
+    label=r'Data')
+(fit,) = ax1.plot(tau_val, tcf_val, marker=None,
+    color="grey",zorder=6, linestyle = '--',
+    label=r'Fit')
+# Shaded error region for non-octahedral states 
+ax1.fill_between(tau_val, tcf_val-tcf_val_err, tcf_val+tcf_val_err,linewidth=0.1,color="peachpuff", alpha=0.8,zorder=0)
+
+# Horizontal line through 0.0 
+plt.axhline(y = 0.0, color = 'black', linestyle = '--', linewidth=1) 
+
+# ax1.set_xlim([0,2.25]) # in percent 
+# xtick_vec = np.arange(0.0,2.5,0.5)
+# xtick_vec = np.append(xtick_vec, 2.25)
+# ax1.set_xticks(xtick_vec)
+# ax1.set_ylim([0.0,25])
+# ax1.set_yscale("log")
+ax1.set_xlim([0.0,10])
+
+# ---------------------------------------------------------------------
+
+# PLOT LABEL (lifetime)
+text = r'$\tau_C=$'+ "{:.2f} ps".format(lifetime)
+ax1.text(4, 0.4, text, fontsize=7.5)
+
+# LEGEND
+
+ax1.legend(handles=[tcf_1, fit], fontsize=7.4)
+
+# ---------------------------------------------------------------------
+
+# plt.show()
+fig.savefig(input_dir/'tcf.png', dpi=300)

meson.build

@@ -53,11 +53,20 @@ py.install_sources([
   subdir: 'soluanalysis'
 )
 
-# IO and extras
+# Extras
 py.install_sources([
-    'soluanalysis/io.py',
-    'soluanalysis/ion_pairs.py'
+    'soluanalysis/ion_pairs.py',
+    'soluanalysis/hdf5_io.py',
+    'soluanalysis/misc.py'
   ],
   pure: false,
   subdir: 'soluanalysis'
+)
+
+# Thirdparty ASE and IO 
+py.install_sources([
+    'soluanalysis/thirdparty/ase/io.py',
+  ],
+  pure: false,
+  subdir: 'soluanalysis/'
 )

pyproject.toml

@@ -10,7 +10,8 @@ authors = [
 dependencies = [
     "ase >= 3.22.1",
     "numpy",
-    "h5py"
+    "h5py",
+    "scipy"
     ]
 requires-python = ">=3.10"
 readme = "README.md"
@@ -19,6 +20,8 @@ license = {text = "MIT"}
 [project.optional-dependencies]
 example = [
     "gdown",
+    "spirit-extras",
+    "matplotlib"
 ]
 
 [project.urls]

python_bindings/bindings.cpp

@@ -1,6 +1,4 @@
 #pragma once
-#include "network_base.hpp"
-#include "undirected_network.hpp"
 #include <cstddef>
 #include <optional>
 #include <vector>
@@ -9,9 +7,11 @@
 #include <fstream>
 #include <iostream>
 #include <string>
+#include <tuple>
 #include <utility>
 
 // Basics
+#include "bondcorrel.hpp"
 #include "bondfinder.hpp"
 #include "directed_network.hpp"
 #include "network_base.hpp"
@@ -77,6 +77,8 @@ PYBIND11_MODULE(james, m) {
       .def_readwrite("boxLo", &James::Atoms::System::boxLo)
       .def("n_atoms", &James::Atoms::System::n_atoms)
       .def("collect_ids", &James::Atoms::System::collect_ids)
+      .def("collect_positions", &James::Atoms::System::collect_positions)
+      .def("reset_positions", &James::Atoms::System::reset_positions)
       .def("delete",
            static_cast<void (James::Atoms::System::*)(int)>(
                &James::Atoms::System::del),
@@ -140,6 +142,49 @@ PYBIND11_MODULE(james, m) {
   py::enum_<James::Path::WriteIdentifier>(m, "WriteIdentifier")
       .value("AtomID", James::Path::WriteIdentifier::AtomID)
       .value("Index", James::Path::WriteIdentifier::Index);
+  // Binding for the templated time_correlation_function, which is templated on
+  // the network type
+  // For a vector of UndirectedNetwork objects
+  m.def(
+      "time_correlation_function",
+      [](const std::vector<Graph::UndirectedNetwork<double>>
+             &network_time_series,
+         const std::vector<double> &time, int start_t0 = 0, int start_tau = 1,
+         int delta_tau = 1, std::optional<int> calc_upto_tau = std::nullopt,
+         bool continuous_bond = true) {
+        return James::Bond::Correlation::time_correlation_function<
+            Graph::UndirectedNetwork<double>>(network_time_series, time,
+                                              start_t0, start_tau, delta_tau,
+                                              calc_upto_tau, continuous_bond);
+      },
+      pybind11::arg("network_time_series"), pybind11::arg("time"),
+      pybind11::arg("start_t0") = 0, pybind11::arg("start_tau") = 1,
+      pybind11::arg("delta_tau") = 1,
+      pybind11::arg("calc_upto_tau") = std::nullopt,
+      pybind11::arg("continuous_bond") = true,
+      "Time correlation function returning tau values, the normalized "
+      "correlation function values, and the standard error in the "
+      "correlation function values.");
+  // For a vector of DirectedNetwork objects
+  m.def(
+      "time_correlation_function",
+      [](const std::vector<Graph::DirectedNetwork<double>> &network_time_series,
+         const std::vector<double> &time, int start_t0 = 0, int start_tau = 1,
+         int delta_tau = 1, std::optional<int> calc_upto_tau = std::nullopt,
+         bool continuous_bond = true) {
+        return James::Bond::Correlation::time_correlation_function<
+            Graph::DirectedNetwork<double>>(network_time_series, time, start_t0,
+                                            start_tau, delta_tau, calc_upto_tau,
+                                            continuous_bond);
+      },
+      pybind11::arg("network_time_series"), pybind11::arg("time"),
+      pybind11::arg("start_t0") = 0, pybind11::arg("start_tau") = 1,
+      pybind11::arg("delta_tau") = 1,
+      pybind11::arg("calc_upto_tau") = std::nullopt,
+      pybind11::arg("continuous_bond") = true,
+      "Time correlation function returning tau values, the normalized "
+      "correlation function values, and the standard error in the "
+      "correlation function values.");
 }
 
 PYBIND11_MODULE(graphlib, m) {