Add all_nuclei method

Spin-Chemistry-Labs · Feb 3, 2024 · 2e4242e · 2e4242e
1 parent 22b1e6d
commit 2e4242e
Show file tree

Hide file tree

Showing 2 changed files with 68 additions and 10 deletions.
diff --git a/radicalpy/data.py b/radicalpy/data.py
@@ -422,6 +422,58 @@ def available(cls):
         paths = get_data("molecules").glob("*.json")
         return sorted([path.with_suffix("").name for path in paths])
 
+    @classmethod
+    def _check_molecule_available(cls, name):
+        if name not in cls.available():
+            lines = [f"Molecule `{name}` not found in database."]
+            lines += ["Available molecules:"]
+            lines += cls.available()
+            raise ValueError("\n".join(lines))
+
+    @classmethod
+    def all_nuclei(cls, name: str):
+        """Construct a molecule from the database with all nuclei.
+
+        Args:
+            name (str): A name of the molecule available in the
+                database (see `Molecule.available()`).
+
+        Examples:
+            >>> Molecule.all_nuclei("flavin_anion")
+            Molecule: flavin_anion
+            Nuclei:
+              14N(19337792.0, 3, 0.5141 <anisotropic available>)
+              14N(19337792.0, 3, -0.001275 <anisotropic available>)
+              14N(19337792.0, 3, -0.03654 <anisotropic available>)
+              1H(267522187.44, 2, 0.05075 <anisotropic available>)
+              1H(267522187.44, 2, -0.1371 <anisotropic available>)
+              1H(267522187.44, 2, -0.1371 <anisotropic available>)
+              1H(267522187.44, 2, -0.1371 <anisotropic available>)
+              1H(267522187.44, 2, -0.4403 <anisotropic available>)
+              1H(267522187.44, 2, 0.4546 <anisotropic available>)
+              1H(267522187.44, 2, 0.4546 <anisotropic available>)
+              1H(267522187.44, 2, 0.4546 <anisotropic available>)
+              1H(267522187.44, 2, 0.009597 <anisotropic available>)
+              1H(267522187.44, 2, 0.4263 <anisotropic available>)
+              1H(267522187.44, 2, 0.4233 <anisotropic available>)
+              1H(267522187.44, 2, -0.02004 <anisotropic available>)
+              14N(19337792.0, 3, 0.1784 <anisotropic available>)
+            Radical: E(-176085963023.0, 2, 0.0 <anisotropic not available>)
+            Info: {'units': 'mT', 'name': 'Flavin radical anion'}
+        """
+        cls._check_molecule_available(name)
+        molecule_json = cls.load_molecule_json(name)
+        info = molecule_json["info"]
+        data = molecule_json["data"]
+        nuclei_list = []
+        for nucleus in data.keys():
+            isotope = data[nucleus]["element"]
+            hfc = data[nucleus]["hfc"]
+            nuclei_list.append(Nucleus.fromisotope(isotope, hfc))
+        molecule = cls(name=name, nuclei=nuclei_list, info=info)
+        molecule.custom = False
+        return molecule
+
     @classmethod
     def fromdb(cls, name: str, nuclei: list[str] = []):
         """Construct a molecule from the database.
@@ -440,11 +492,7 @@ def fromdb(cls, name: str, nuclei: list[str] = []):
             Radical: E(-176085963023.0, 2, 0.0 <anisotropic not available>)
             Info: {'units': 'mT', 'name': 'Flavin radical anion'}
         """
-        if name not in cls.available():
-            lines = [f"Molecule `{name}` not found in database."]
-            lines += ["Available molecules:"]
-            lines += cls.available()
-            raise ValueError("\n".join(lines))
+        cls._check_molecule_available(name)
         molecule_json = cls.load_molecule_json(name)
         info = molecule_json["info"]
         data = molecule_json["data"]

diff --git a/radicalpy/semiclassical.py b/radicalpy/semiclassical.py
@@ -2,6 +2,7 @@
 import numpy as np
 
 import radicalpy as rp
+from radicalpy.data import Molecule
 
 H, N = 0.5, 1
 FAD_spin = [N, N, N, N, H, H, H, H, H, H, H, H, H, H, H, H, H, H, H, H, H]
@@ -66,11 +67,20 @@ def angle(num_samples):
         yield theta, phi
 
 
+def cmp_hfcs(mol, hfcs):
+    mol_sorted = sorted(FAD.nuclei, key=lambda k: k.hfc.isotropic)
+    hfcs_sorted = sorted(FAD_HFCs)
+    for n1, n2 in zip(mol_sorted, hfcs_sorted):
+        print(f"{n1.hfc.isotropic:10.4}, {n2:10.4}")
+
+
 if __name__ == "__main__":
     sample_number = 10
     np.random.seed(42)
-    theta, phi = semiclassical_theta_phi(sample_number)
-    np.random.seed(42)
-    for i, (t, p, (tt, pp)) in enumerate(zip(theta, phi, angle(sample_number))):
-        print(f"{t-tt=} {p-pp=} {tt=} {pp=}")
-    FAD = rp.data.Molecule.fromdb("FAD")
+    FAD = Molecule.all_nuclei("flavin_anion")
+    Trp = Molecule.all_nuclei("tryptophan_cation")
+    print(f"{len(FAD.nuclei)=}")
+    print(f"{len(FAD_HFCs)=}")
+    print(f"{len(Trp.nuclei)=}")
+    print(f"{len(Trp_HFCs)=}")
+    cmp_hfcs(FAD, FAD_HFCs)