improved experimental comparison

examples/illustrations/bulk-water-tip4p-vs-temperature/analysis.py

@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+import numpy as np
+import MDAnalysis as mda
+import nmrformd
+
+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")
+
+sys.path.append(git_path + '/examples/shared/')
+from utilities import save_result, measure_N, measure_R1
+
+path = "/data/nmrformd-data/bulk-water-tip4p-vs-temperature/raw_data/"
+run = True
+R1_intra = 1e-5
+R1_inter = 1e-5
+while run:
+    for T in np.arange(280, 330, 10):
+        # Read current status
+        name_intra = "N4000_intra_T" + str(T) + "K"
+        name_inter = "N4000_inter_T" + str(T) + "K"
+        n_intra = measure_N(name_intra)
+        n_inter = measure_N(name_inter)
+        # Import MDA universe
+        datapath = path + "N4000_" + str(T) + "K/"
+        u = mda.Universe(datapath+"prod.tpr", datapath+"prod.xtc")
+        hydrogen = u.select_atoms("name HW1 HW2")
+        n_hydrogen = hydrogen.n_atoms
+        if n_intra < n_hydrogen:
+            # Calculated NMR properties
+            intra = nmrformd.NMR(u, atom_group = hydrogen,
+                                neighbor_group = hydrogen,
+                                type_analysis = "intra_molecular",
+                                number_i = 1)
+            save_result(intra, name = name_intra)
+            if np.random.random() < 0.2:
+                inter = nmrformd.NMR(u, atom_group = hydrogen,
+                                    neighbor_group = hydrogen,
+                                    type_analysis = "inter_molecular",
+                                    number_i = 1)
+                save_result(inter, name = "N4000_inter_T" + str(T) + "K")
+        else:
+            run = False
+        # Print information
+        R1_intra = measure_R1(name_intra)
+        R1_inter = measure_R1(name_inter)
+        print("T =", T, "K --", n_intra, n_hydrogen, 
+            "R1:", np.round(R1_intra,2), np.round(R1_inter,2), 
+            "T1:", np.round(1/(R1_intra+R1_inter),2))
+    print(" ")

examples/shared/utilities.py

@@ -0,0 +1,84 @@
+import numpy as np
+import os
+
+def save_result(data, name, 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)
+
+def measure_N(name):
+    """measure_N"""
+    try:
+        saving_file = "raw_data/" + name + ".npy"
+        previous_dictionary = np.load(saving_file, allow_pickle=True)
+        N = np.real(previous_dictionary.item()["N"])
+    except:
+        N = 0
+    return N
+
+def measure_R1(name):
+    """measure_N"""
+    try:
+        saving_file = "raw_data/" + name + ".npy"
+        previous_dictionary = np.load(saving_file, allow_pickle=True)
+        R1 = np.real(previous_dictionary.item()["R1"])
+        return R1[0]
+    except:
+        return 1e-5