Merge pull request #151 from Spin-Chemistry-Labs/149-fix-rpsimulation…

…hilbertsimluationdipolar_hamiltonian_3d 149 fix rpsimulationhilbertsimluationdipolar hamiltonian 3d
Spin-Chemistry-Labs · Apr 17, 2024 · e56c262 · e56c262
2 parents 67653f3 + 884c944
commit e56c262
Show file tree

Hide file tree

Showing 11 changed files with 358 additions and 205 deletions.
diff --git a/.github/workflows/lint.yaml b/.github/workflows/lint.yaml
@@ -6,7 +6,7 @@ jobs:
   lint:
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@master
       - uses: psf/black@stable
       - uses: isort/isort-action@master
         with:

diff --git a/examples/anisotropy_3d_polar_paper.py b/examples/anisotropy_3d_polar_paper.py
@@ -0,0 +1,67 @@
+#! /usr/bin/env python
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import radicalpy as rp
+from radicalpy import kinetics, relaxation
+from radicalpy.experiments import anisotropy
+from radicalpy.simulation import State
+
+
+def main():
+
+    fad_n5_hfc = np.array(
+        [
+            [0.280, -0.138, 0.678],
+            [-0.138, 0.043, -0.331],
+            [0.678, -0.331, 1.412],
+        ]
+    )
+
+    dipolar = rp.estimations.dipolar_interaction_anisotropic(r=22e-10)
+
+    theta = np.linspace(0, np.pi, 35)
+    phi = np.linspace(0, 2 * np.pi, 58)
+
+    flavin = rp.simulation.Molecule.fromisotopes(isotopes=["14N"], hfcs=[fad_n5_hfc])
+    Z = rp.simulation.Molecule("zorro", [])
+    sim = rp.simulation.HilbertSimulation([flavin, Z])
+
+    time = np.arange(0, 15e-6, 5e-9)
+    B0 = 0.05
+    k = 1e6
+
+    results = anisotropy(
+        sim,
+        init_state=State.SINGLET,
+        obs_state=State.SINGLET,
+        time=time,
+        theta=theta,
+        phi=phi,
+        B0=B0,
+        D=dipolar,
+        J=0,
+        kinetics=[kinetics.Exponential(k)],
+    )
+
+    Y = results["product_yield_sums"]
+    delta_phi_s, gamma_s = rp.utils.yield_anisotropy(Y, theta, phi)
+    Y_av = rp.utils.spherical_average(Y, theta, phi)
+    Y = Y - Y_av
+
+    np.save("Y", Y)
+    np.save("theta", theta)
+    np.save("phi", phi)
+    rp.plot.anisotropy_surface(theta, phi, Y)
+
+    print(f"{Y_av=}")
+    print(f"{delta_phi_s=}")
+    print(f"{gamma_s=}")
+    plt.show()
+
+    return 0
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/anisotropy_3d_polar_paper_2.py b/examples/anisotropy_3d_polar_paper_2.py
@@ -0,0 +1,84 @@
+#! /usr/bin/env python
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import radicalpy as rp
+from radicalpy import kinetics, relaxation
+from radicalpy.experiments import anisotropy
+from radicalpy.simulation import State
+
+
+def main():
+
+    fad_n5_hfc = np.array(
+        [
+            [0.280, -0.138, 0.678],
+            [-0.138, 0.043, -0.331],
+            [0.678, -0.331, 1.412],
+        ]
+    )
+
+    trp_hbeta_hfc = np.array(
+        [
+            [0.944, -0.019, 0.030],
+            [-0.019, 1.091, -0.065],
+            [0.030, -0.065, 1.070],
+        ]
+    )
+
+    dipolar = (
+        np.array(
+            [
+                [-0.225, 0.156, 0.198],
+                [0.156, 0.117, -0.082],
+                [0.198, -0.082, 0.107],
+            ]
+        )
+        * rp.data.Isotope("E").gamma_mT
+    )
+
+    theta = np.linspace(0, np.pi, 35)
+    phi = np.linspace(0, 2 * np.pi, 58)
+
+    flavin = rp.simulation.Molecule.fromisotopes(isotopes=["14N"], hfcs=[fad_n5_hfc])
+    trp = rp.simulation.Molecule.fromisotopes(isotopes=["1H"], hfcs=[trp_hbeta_hfc])
+    sim = rp.simulation.HilbertSimulation([flavin, trp])
+
+    time = np.arange(0, 15e-6, 5e-9)
+    B0 = 0.05
+    k = 1e6
+
+    results = anisotropy(
+        sim,
+        init_state=State.SINGLET,
+        obs_state=State.SINGLET,
+        time=time,
+        theta=theta,
+        phi=phi,
+        B0=B0,
+        D=dipolar,
+        J=0,
+        kinetics=[kinetics.Exponential(k)],
+    )
+
+    Y = results["product_yield_sums"]
+    delta_phi_s, gamma_s = rp.utils.yield_anisotropy(Y, theta, phi)
+    Y_av = rp.utils.spherical_average(Y, theta, phi)
+    Y = Y - Y_av
+
+    np.save("Y", Y)
+    np.save("theta", theta)
+    np.save("phi", phi)
+    rp.plot.anisotropy_surface(theta, phi, Y)
+
+    print(f"{Y_av=}")
+    print(f"{delta_phi_s=}")
+    print(f"{gamma_s=}")
+    plt.show()
+
+    return 0
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/comparison.ipynb b/examples/comparison.ipynb
@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 123,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -190,27 +190,28 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 357,
+   "execution_count": 23,
    "metadata": {},
    "outputs": [],
    "source": [
     "spins_hilbert = [\"3\", \"4\", \"5\", \"6\", \"7\"]#, \"8\"]\n",
     "hilbert = [0.047, 0.048, 0.065, 0.19, 0.847]#, 6.08]\n",
-    "spins_liouville = [\"3\", \"4\", \"5\", \"6\", \"7\"]\n",
-    "liouville = [0.019, 0.068, 0.54, 20.89, 908.59]\n",
+    "spins_liouville = [\"3\", \"4\", \"5\", \"6\"]\n",
+    "# liouville = [0.019, 0.068, 0.54, 20.89, 908.59]\n",
+    "liouville = [0.026, 0.114, 1.98, 102.32]\n",
     "semiclassical = 0.53"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 361,
+   "execution_count": 4,
    "metadata": {},
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "C:\\Users\\lewis\\AppData\\Local\\Temp\\ipykernel_15612\\2357155276.py:8: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.\n",
+      "C:\\Users\\lewis\\AppData\\Local\\Temp\\ipykernel_21480\\2357155276.py:8: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.\n",
       "  cmap = get_cmap('viridis')\n"
      ]
     }
@@ -241,14 +242,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 362,
+   "execution_count": 24,
    "metadata": {},
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "C:\\Users\\lewis\\AppData\\Local\\Temp\\ipykernel_15612\\2157149891.py:8: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.\n",
+      "C:\\Users\\lewis\\AppData\\Local\\Temp\\ipykernel_21480\\715634988.py:8: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.\n",
       "  cmap = get_cmap('viridis')\n"
      ]
     }
@@ -282,23 +283,23 @@
     "ax.set_yscale('log')\n",
     "plt.tick_params(labelsize=18)\n",
     "plt.gcf().set_size_inches(10, 5)\n",
-    "path = \"comparison\" + f\"_{5}.png\"\n",
+    "path = \"comparison\" + f\"_{15}.png\"\n",
     "plt.savefig(path, dpi=300, bbox_inches=\"tight\")\n",
     "plt.close()"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 309,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "C:\\Users\\lewis\\AppData\\Local\\Temp\\ipykernel_15612\\369739269.py:5: RuntimeWarning: divide by zero encountered in true_divide\n",
+      "C:\\Users\\lewis\\AppData\\Local\\Temp\\ipykernel_21480\\369739269.py:5: RuntimeWarning: divide by zero encountered in true_divide\n",
       "  one_nuc_mary = one_nuc[\"MARY\"] / one_nuc[\"MARY\"].max()\n",
-      "C:\\Users\\lewis\\AppData\\Local\\Temp\\ipykernel_15612\\369739269.py:5: RuntimeWarning: invalid value encountered in true_divide\n",
+      "C:\\Users\\lewis\\AppData\\Local\\Temp\\ipykernel_21480\\369739269.py:5: RuntimeWarning: invalid value encountered in true_divide\n",
       "  one_nuc_mary = one_nuc[\"MARY\"] / one_nuc[\"MARY\"].max()\n"
      ]
     }
@@ -314,7 +315,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 310,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -327,7 +328,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 311,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -340,7 +341,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 312,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -353,70 +354,59 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 313,
+   "execution_count": 14,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "C:\\Users\\lewis\\AppData\\Local\\Temp\\ipykernel_15612\\3517896710.py:5: RuntimeWarning: divide by zero encountered in true_divide\n",
-      "  one_nuc_liouv_mary = one_nuc_liouv[\"MARY\"] / one_nuc_liouv[\"MARY\"].max()\n",
-      "C:\\Users\\lewis\\AppData\\Local\\Temp\\ipykernel_15612\\3517896710.py:5: RuntimeWarning: invalid value encountered in true_divide\n",
-      "  one_nuc_liouv_mary = one_nuc_liouv[\"MARY\"] / one_nuc_liouv[\"MARY\"].max()\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "one_nuc_liouv = np.load(\n",
-    "        \"./data/fad_mary/results_1nuc_liouville.npy\", allow_pickle=True\n",
+    "        \"./data/fad_mary/results_1nuc_liouville_relaxation.npy\", allow_pickle=True\n",
     "    ).item()\n",
     "\n",
     "one_nuc_liouv_mary = one_nuc_liouv[\"MARY\"] / one_nuc_liouv[\"MARY\"].max()"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 314,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [],
    "source": [
     "two_nuc_liouv = np.load(\n",
-    "        \"./data/fad_mary/results_2nuc_liouville.npy\", allow_pickle=True\n",
+    "        \"./data/fad_mary/results_2nuc_liouville_relaxation.npy\", allow_pickle=True\n",
     "    ).item()\n",
     "\n",
     "two_nuc_liouv_mary = two_nuc_liouv[\"MARY\"] / two_nuc_liouv[\"MARY\"].max()"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 315,
+   "execution_count": 16,
    "metadata": {},
    "outputs": [],
    "source": [
     "three_nuc_liouv = np.load(\n",
-    "        \"./data/fad_mary/results_3nuc_liouville.npy\", allow_pickle=True\n",
+    "        \"./data/fad_mary/results_3nuc_liouville_relaxation.npy\", allow_pickle=True\n",
     "    ).item()\n",
     "\n",
     "three_nuc_liouv_mary = three_nuc_liouv[\"MARY\"] / three_nuc_liouv[\"MARY\"].max()"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 316,
+   "execution_count": 21,
    "metadata": {},
    "outputs": [],
    "source": [
     "four_nuc_liouv = np.load(\n",
-    "        \"./data/fad_mary/results_4nuc_liouville.npy\", allow_pickle=True\n",
+    "        \"./data/fad_mary/results_4nuc_liouville_relaxation.npy\", allow_pickle=True\n",
     "    ).item()\n",
     "\n",
     "four_nuc_liouv_mary = four_nuc_liouv[\"MARY\"] / four_nuc_liouv[\"MARY\"].max()"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 318,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -427,7 +417,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 330,
+   "execution_count": 19,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -466,7 +456,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 363,
+   "execution_count": 22,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -489,7 +479,7 @@
     "            fontsize=16, loc=\"center right\", bbox_to_anchor = (1.4, 0.5))#, ncol=3)\n",
     "plt.tick_params(labelsize=18)\n",
     "plt.gcf().set_size_inches(10, 5)\n",
-    "path = \"comparison\" + f\"_{9}.png\"\n",
+    "path = \"comparison\" + f\"_{12}.png\"\n",
     "plt.savefig(path, dpi=300, bbox_inches=\"tight\")\n",
     "plt.close()"
    ]

diff --git a/examples/kinetics_fad_semiclassical_wavelength.py b/examples/kinetics_fad_semiclassical_wavelength.py
@@ -245,7 +245,7 @@ def main():
         relaxations=[RandomFields(kr), SingletTripletDephasing(kstd)],
     )
 
-    np.save("./examples/data/fad_mary/results_new.npy", results)
+    # np.save("./examples/data/fad_mary/results_new.npy", results)
 
     total_yield = np.zeros((len(time), len(Bs), len(wavelength)), dtype=complex)
     zero_field = np.zeros((len(time), len(Bs), len(wavelength)), dtype=complex)