added mode analysed data

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

@@ -2,18 +2,9 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "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"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "import numpy as np\n",
     "import MDAnalysis as mda\n",
@@ -22,7 +13,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -34,7 +25,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -104,23 +95,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": null,
    "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"
-     ]
-    }
-   ],
+   "outputs": [],
    "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",
@@ -132,7 +112,7 @@
     "        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)"
+    "        print(n_hydrogen, T, n_intra)"
    ]
   }
  ],

examples/illustrations/bulk-water/analysis.ipynb

@@ -103,7 +103,7 @@
     "    # varying number\n",
     "    mystep = 1\n",
     "    for N in [\"N25\", \"N39\", \"N62\", \"N99\", \"N158\", \"N251\", \"N398\", \"N631\", \"N1002\", \"N1589\", \"N2521\"]:\n",
-    "        datapath = git_path + \"/nmrformd-data/bulk-water/raw-data/\"+N+\"/\"\n",
+    "        datapath = git_path + \"/data/bulk-water-tip4p/bulk-water/raw-data/\"+N+\"/\"\n",
     "        for n in np.arange(1, 10):\n",
     "            if (os.path.exists(datapath+\"prod\"+str(n)+\".xtc\")) & (os.path.exists(datapath+\"topology.data\")):\n",
     "                u = mda.Universe(datapath+\"topology.data\", datapath+\"prod\"+str(n)+\".xtc\")\n",
@@ -124,7 +124,7 @@
     "    # varying step\n",
     "    N = \"N4000\"\n",
     "    for mystep in [1, 2, 4, 8, 16, 32]:\n",
-    "        datapath = git_path + \"/nmrformd-data/bulk-water/raw-data/\"+N+\"/\"\n",
+    "        datapath = git_path + \"/data/bulk-water-tip4p/bulk-water/raw-data/\"+N+\"/\"\n",
     "        xtcs = [[datapath+\"prod1.xtc\", datapath+\"prod2.xtc\"],\n",
     "                [datapath+\"prod3.xtc\", datapath+\"prod4.xtc\"]]\n",
     "        for xtc in xtcs:\n",
@@ -139,29 +139,9 @@
     "            n_inter = save_result(inter, name = N + \"_inter_dt\" + str(dt))\n",
     "            if mystep == 1:\n",
     "                print(n_hydrogen, mystep, n_intra, dt)\n",
-    "    # high res\n",
-    "    N = \"N4000_HR\"\n",
-    "    for mystep in [1, 2, 4, 8, 16, 32]:\n",
-    "        datapath = git_path + \"/nmrformd-data/bulk-water/raw-data/\"+N+\"/\"\n",
-    "        xtc = [datapath+\"prod1.xtc\", datapath+\"prod2.xtc\",\n",
-    "               datapath+\"prod3.xtc\", datapath+\"prod4.xtc\",\n",
-    "               datapath+\"prod5.xtc\", datapath+\"prod6.xtc\",\n",
-    "               datapath+\"prod7.xtc\", datapath+\"prod8.xtc\",\n",
-    "               datapath+\"prod9.xtc\", datapath+\"prod10.xtc\"]\n",
-    "        u = mda.Universe(datapath+\"prod1.tpr\", xtc)\n",
-    "        u.transfer_to_memory(step = mystep)\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 = N + \"_intra_dt\" + str(dt))\n",
-    "        n_inter = save_result(inter, name = N + \"_inter_dt\" + str(dt))\n",
-    "        if mystep == 1:\n",
-    "            print(n_hydrogen, mystep, n_intra, dt)\n",
-    "    # iso\n",
+    "    # aniso\n",
     "    N = \"N4000\"\n",
-    "    datapath = git_path + \"/nmrformd-data/bulk-water/raw-data/\"+N+\"/\"\n",
+    "    datapath = git_path + \"/data/bulk-water-tip4p/bulk-water/raw-data/\"+N+\"/\"\n",
     "    xtcs = [[datapath+\"prod1.xtc\", datapath+\"prod2.xtc\"],\n",
     "            [datapath+\"prod3.xtc\", datapath+\"prod4.xtc\"]]\n",
     "    for xtc in xtcs:\n",
@@ -177,6 +157,39 @@
     "        n_inter = save_result(inter, name = N + \"_inter_aniso_dt\" + str(dt), aniso = True)\n",
     "        print(\"N:\", n_intra)"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "while 1<2:   \n",
+    "    # vary T\n",
+    "    for T in [280, 290, 310, 320]:\n",
+    "        datapath = git_path + \"/data/bulk-water-tip4p/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",
+    "        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)\n",
+    "    # vary co\n",
+    "    for co in [6, 8, 10, 12]:\n",
+    "        datapath = git_path + \"/data/bulk-water-tip4p/bulk-water/raw-data/N4000_co\" + str(co) + \"/\"\n",
+    "        xtc = [datapath+\"prod1.xtc\", datapath+\"prod2.xtc\"]\n",
+    "        u = mda.Universe(datapath+\"prod1.tpr\", xtc)\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_co\" + str(co))\n",
+    "        n_inter = save_result(inter, name = \"N4000_inter_co\" + str(co))\n",
+    "        print(n_hydrogen, co, n_intra)"
+   ]
   }
  ],
  "metadata": {

examples/illustrations/bulk-water/analysis.py

@@ -0,0 +1,119 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# In[ ]:
+
+
+import numpy as np
+import MDAnalysis as mda
+import nmrformd
+
+
+# In[ ]:
+
+
+import sys, os, git
+current_path = os.getcwd()
+git_repo = git.Repo(current_path, search_parent_directories=True)
+git_path = git_repo.git.rev_parse("--show-toplevel")
+
+
+# In[ ]:
+
+
+def save_result(data, name = "intra_H2O", aniso = False):
+    """Save the correlation functions in dictionary"""
+    if not os.path.exists("raw_data/"):
+        os.makedirs("raw_data/")
+    saving_file = "raw_data/" + name + ".npy"
+    t = data.t
+    f = data.f
+    if aniso:
+        C1 = data.gij[0]
+        C2 = data.gij[1]
+        C3 = data.gij[2]
+    else:
+        C = data.gij[0]
+    R1 = data.R1
+    R2 = data.R2
+    N = data.group_j.atoms.n_atoms
+    try:
+        previous_dictionary = np.load(saving_file, allow_pickle=True)
+        t_prev = np.real(previous_dictionary.item()["t"])
+        assert len(t_prev) == len(t)
+        if aniso:
+            C1_prev = np.real(previous_dictionary.item()["C1"])
+            C2_prev = np.real(previous_dictionary.item()["C2"])
+            C3_prev = np.real(previous_dictionary.item()["C3"])
+        else:
+            C_prev = np.real(previous_dictionary.item()["C"])
+        R1_prev = np.real(previous_dictionary.item()["R1"])
+        R2_prev = np.real(previous_dictionary.item()["R2"])
+        N_prev = np.real(previous_dictionary.item()["N"])
+        if aniso:
+            C1 = (C1*N + C1_prev*N_prev) / (N_prev + N)
+            C2 = (C2*N + C2_prev*N_prev) / (N_prev + N)
+            C3 = (C3*N + C3_prev*N_prev) / (N_prev + N)
+        else:
+            C = (C*N + C_prev*N_prev) / (N_prev + N)
+        R1 = (R1*N + R1_prev*N_prev) / (N_prev + N)
+        R2 = (R2*N + R2_prev*N_prev) / (N_prev + N)
+        N += N_prev
+    except:
+        pass
+    if aniso:
+        dictionary = {
+        "t": t,
+        "f": f,
+        "C1": C1,
+        "C2": C2,
+        "C3": C3,
+        "N": N,
+        "R1": R1,
+        "R2": R2,
+        }
+    else:
+        dictionary = {
+        "t": t,
+        "f": f,
+        "C": C,
+        "N": N,
+        "R1": R1,
+        "R2": R2,
+        }
+    np.save(saving_file, dictionary)
+    return N
+
+
+# In[ ]:
+
+
+# In[8]:
+
+
+while 1<2:   
+    # vary T
+    for T in [280, 290, 310, 320]:
+        datapath = git_path + "/data/bulk-water-tip4p/bulk-water/raw-data/N4000_" + str(T) + "K/"
+        xtc = [datapath+"prod1.xtc", datapath+"prod2.xtc"]
+        u = mda.Universe(datapath+"prod1.tpr", xtc)
+        hydrogen = u.select_atoms("name HW1 HW2")
+        n_hydrogen = hydrogen.n_atoms
+        intra = nmrformd.NMR(u, atom_group = hydrogen, neighbor_group = hydrogen, type_analysis = "intra_molecular", number_i=5)
+        inter = nmrformd.NMR(u, atom_group = hydrogen, neighbor_group = hydrogen, type_analysis = "inter_molecular", number_i=1)
+        n_intra = save_result(intra, name = "N4000_intra_T" + str(T) + "K")
+        n_inter = save_result(inter, name = "N4000_inter_T" + str(T) + "K")
+        print(n_hydrogen, T, n_intra)
+    # vary co
+    for co in [6, 8, 10, 12]:
+        datapath = git_path + "/data/bulk-water-tip4p/bulk-water/raw-data/N4000_co" + str(co) + "/"
+        xtc = [datapath+"prod1.xtc", datapath+"prod2.xtc"]
+        u = mda.Universe(datapath+"prod1.tpr", xtc)
+        hydrogen = u.select_atoms("name HW1 HW2")
+        n_hydrogen = hydrogen.n_atoms
+        intra = nmrformd.NMR(u, atom_group = hydrogen, neighbor_group = hydrogen, type_analysis = "intra_molecular", number_i=5)
+        inter = nmrformd.NMR(u, atom_group = hydrogen, neighbor_group = hydrogen, type_analysis = "inter_molecular", number_i=1)
+        n_intra = save_result(intra, name = "N4000_intra_co" + str(co))
+        n_inter = save_result(inter, name = "N4000_inter_co" + str(co))
+        print(n_hydrogen, co, n_intra)
+