added tip3p data

.gitmodules

@@ -19,3 +19,6 @@
 [submodule "data/HEWL-in-water"]
 	path = data/HEWL-in-water
 	url = [email protected]:simongravelle/HEWL-in-water.git
+[submodule "data/bulk-water-tip3p"]
+	path = data/bulk-water-tip3p
+	url = [email protected]:simongravelle/bulk-water-tip3p.git

examples/illustrations/bulk-water-tip3p/analysis.ipynb

@@ -0,0 +1,160 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/simon/.local/lib/python3.11/site-packages/MDAnalysis/topology/TPRParser.py:161: DeprecationWarning: 'xdrlib' is deprecated and slated for removal in Python 3.13\n",
+      "  import xdrlib\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "import MDAnalysis as mda\n",
+    "import nmrformd"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys, os, git\n",
+    "current_path = os.getcwd()\n",
+    "git_repo = git.Repo(current_path, search_parent_directories=True)\n",
+    "git_path = git_repo.git.rev_parse(\"--show-toplevel\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def save_result(data, name, aniso = False):\n",
+    "    \"\"\"Save the correlation functions in dictionary\"\"\"\n",
+    "    if not os.path.exists(\"raw_data/\"):\n",
+    "        os.makedirs(\"raw_data/\")\n",
+    "    saving_file = \"raw_data/\" + name + \".npy\"\n",
+    "    t = data.t\n",
+    "    f = data.f\n",
+    "    if aniso:\n",
+    "        C1 = data.gij[0]\n",
+    "        C2 = data.gij[1]\n",
+    "        C3 = data.gij[2]\n",
+    "    else:\n",
+    "        C = data.gij[0]\n",
+    "    R1 = data.R1\n",
+    "    R2 = data.R2\n",
+    "    N = data.group_j.atoms.n_atoms\n",
+    "    try:\n",
+    "        previous_dictionary = np.load(saving_file, allow_pickle=True)\n",
+    "        t_prev = np.real(previous_dictionary.item()[\"t\"])\n",
+    "        assert len(t_prev) == len(t)\n",
+    "        if aniso:\n",
+    "            C1_prev = np.real(previous_dictionary.item()[\"C1\"])\n",
+    "            C2_prev = np.real(previous_dictionary.item()[\"C2\"])\n",
+    "            C3_prev = np.real(previous_dictionary.item()[\"C3\"])\n",
+    "        else:\n",
+    "            C_prev = np.real(previous_dictionary.item()[\"C\"])\n",
+    "        R1_prev = np.real(previous_dictionary.item()[\"R1\"])\n",
+    "        R2_prev = np.real(previous_dictionary.item()[\"R2\"])\n",
+    "        N_prev = np.real(previous_dictionary.item()[\"N\"])\n",
+    "        if aniso:\n",
+    "            C1 = (C1*N + C1_prev*N_prev) / (N_prev + N)\n",
+    "            C2 = (C2*N + C2_prev*N_prev) / (N_prev + N)\n",
+    "            C3 = (C3*N + C3_prev*N_prev) / (N_prev + N)\n",
+    "        else:\n",
+    "            C = (C*N + C_prev*N_prev) / (N_prev + N)\n",
+    "        R1 = (R1*N + R1_prev*N_prev) / (N_prev + N)\n",
+    "        R2 = (R2*N + R2_prev*N_prev) / (N_prev + N)\n",
+    "        N += N_prev\n",
+    "    except:\n",
+    "        pass\n",
+    "    if aniso:\n",
+    "        dictionary = {\n",
+    "        \"t\": t,\n",
+    "        \"f\": f,\n",
+    "        \"C1\": C1,\n",
+    "        \"C2\": C2,\n",
+    "        \"C3\": C3,\n",
+    "        \"N\": N,\n",
+    "        \"R1\": R1,\n",
+    "        \"R2\": R2,\n",
+    "        }\n",
+    "    else:\n",
+    "        dictionary = {\n",
+    "        \"t\": t,\n",
+    "        \"f\": f,\n",
+    "        \"C\": C,\n",
+    "        \"N\": N,\n",
+    "        \"R1\": R1,\n",
+    "        \"R2\": R2,\n",
+    "        }\n",
+    "    np.save(saving_file, dictionary)\n",
+    "    return N"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "8000 280 1 1.0\n",
+      "8000 290 1 1.0\n",
+      "8000 300 1 1.0\n"
+     ]
+    }
+   ],
+   "source": [
+    "while 1<2:   \n",
+    "    # high res\n",
+    "    N = \"N4000_HR\"\n",
+    "    for T in [280, 290, 300, 310, 320]:\n",
+    "        datapath = git_path + \"/data/bulk-water-tip3p/bulk-water/raw-data/N4000_\" + str(T) + \"K/\"\n",
+    "        xtc = [datapath+\"prod1.xtc\", datapath+\"prod2.xtc\"]\n",
+    "        u = mda.Universe(datapath+\"prod1.tpr\", xtc)\n",
+    "        dt = np.round(u.trajectory.dt,2)\n",
+    "        hydrogen = u.select_atoms(\"name HW1 HW2\")\n",
+    "        n_hydrogen = hydrogen.n_atoms\n",
+    "        intra = nmrformd.NMR(u, atom_group = hydrogen, neighbor_group = hydrogen, type_analysis = \"intra_molecular\", number_i=5)\n",
+    "        inter = nmrformd.NMR(u, atom_group = hydrogen, neighbor_group = hydrogen, type_analysis = \"inter_molecular\", number_i=1)\n",
+    "        n_intra = save_result(intra, name = \"N4000_intra_T\" + str(T) + \"K\")\n",
+    "        n_inter = save_result(inter, name = \"N4000_inter_T\" + str(T) + \"K\")\n",
+    "        print(n_hydrogen, T, n_intra, dt)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}