Sketched simple phys constants class. not very pretty...

python/mrchem/CUBEparser.py

@@ -28,7 +28,9 @@
 from json import dump
 
 from .input_parser.plumbing import pyparsing as pp
-from.physical_constants import PhysicalConstants as PC
+from.physical_constants import PhysicalConstants
+
+pc = PhysicalConstants()
 
 
 def write_cube_dict(file_dict, world_unit):
@@ -173,7 +175,7 @@ def count(s, l, t):
 
     # get cube origin data
     N_atoms = parsed_cube["NATOMS"]
-    origin = parsed_cube["ORIGIN"] if (world_unit == "bohr") else [p*PC.ANGSTROM_2_BOHR for p in parsed_cube["ORIGIN"]]
+    origin = parsed_cube["ORIGIN"] if (world_unit == "bohr") else [p * pc.CONVERT_ANGSTROM_2_BOHR for p in parsed_cube["ORIGIN"]]
 
     # Set the amount of values depending on if the DSET_IDs were present or not
     if (len(parsed_cube["DSET_IDS"]) != 0):
@@ -198,9 +200,9 @@ def count(s, l, t):
     if (world_unit == "bohr"):
         Voxel_axes = [parsed_cube["XAXIS"]["VECTOR"], parsed_cube["YAXIS"]["VECTOR"], parsed_cube["ZAXIS"]["VECTOR"]]
     else:
-        X_voxel = [p*PC.ANGSTROM_2_BOHR for p in parsed_cube["XAXIS"]["VECTOR"]]
-        Y_voxel = [p*PC.ANGSTROM_2_BOHR for p in parsed_cube["YAXIS"]["VECTOR"]]
-        Z_voxel = [p*PC.ANGSTROM_2_BOHR for p in parsed_cube["ZAXIS"]["VECTOR"]]
+        X_voxel = [p * pc.CONVERT_ANGSTROM_2_BOHR for p in parsed_cube["XAXIS"]["VECTOR"]]
+        Y_voxel = [p * pc.CONVERT_ANGSTROM_2_BOHR for p in parsed_cube["YAXIS"]["VECTOR"]]
+        Z_voxel = [p * pc.CONVERT_ANGSTROM_2_BOHR for p in parsed_cube["ZAXIS"]["VECTOR"]]
         Voxel_axes = [X_voxel, Y_voxel, Z_voxel]
 
     # get the atom coordinates
@@ -209,7 +211,7 @@ def count(s, l, t):
 
     Z_n = [atom["ATOMIC_NUMBER"] for atom in parsed_cube["GEOM"]]
     atom_charges = [atom["CHARGE"] for atom in parsed_cube["GEOM"]]
-    atom_coords = [atom["POSITION"] if (world_unit == "bohr") else [p*PC.ANGSTROM_2_BOHR for p in atom["POSITION"]] for atom in parsed_cube["GEOM"]]
+    atom_coords = [atom["POSITION"] if (world_unit == "bohr") else [p * pc.CONVERT_ANGSTROM_2_BOHR for p in atom["POSITION"]] for atom in parsed_cube["GEOM"]]
 
     # construct the CUBE vector. Indexing is CUBE_vector[MO_ID][i*N_vals[1]*N_vals[2] + j*N_vals[2] + k] where i, j and k correspond to steps in the X, Y and Z voxel axes directions respectively.
     CUBE_vector = []

python/mrchem/api.py

@@ -34,16 +34,18 @@
     write_rsp_calc,
     parse_wf_method
 )
-from .physical_constants import PhysicalConstants as PC
+from .physical_constants import PhysicalConstants
 from .periodictable import PeriodicTable as PT, PeriodicTableByZ as PT_Z
 from .validators import MoleculeValidator
 
+pc = PhysicalConstants()
+
 
 def translate_input(user_dict):
     # get the origin in the desired units of measure
     origin = user_dict["world_origin"]
     if user_dict["world_unit"] == "angstrom":
-        origin = [PC.ANGSTROM_2_BOHR * r for r in origin]
+        origin = [pc.CONVERT_ANGSTROM_2_BOHR * r for r in origin]
 
     # prepare bits and pieces
     mol_dict = write_molecule(user_dict, origin)

python/mrchem/helpers.py

@@ -23,9 +23,11 @@
 # <https://mrchem.readthedocs.io/>
 #
 
-from .physical_constants import PhysicalConstants as PC
+from .physical_constants import PhysicalConstants
 from .CUBEparser import write_cube_dict
 
+pc = PhysicalConstants()
+
 # yapf: disable
 SHORTHAND_FUNCTIONALS = [
     'svwn3',
@@ -258,7 +260,7 @@ def write_scf_plot(user_dict):
         plot_dict["plotter"] = user_dict["Plotter"]
         if user_dict["world_unit"] == "angstrom":
             plot_dict["plotter"] = {
-                k: [PC.ANGSTROM_2_BOHR * r for r in plot_dict["plotter"][k]]
+                k: [pc.CONVERT_ANGSTROM_2_BOHR * r for r in plot_dict["plotter"][k]]
                 for k in plot_dict["plotter"].keys()
             }
     return plot_dict

python/mrchem/physical_constants.py

@@ -23,6 +23,66 @@
 # <https://mrchem.readthedocs.io/>
 #
 
+from pathlib import Path
+
 class PhysicalConstants:
-    BOHR_2_METER = 5.29177210903e-11
-    ANGSTROM_2_BOHR = 1.0e-10 / BOHR_2_METER
+    """Based on CODATA2018 from https://physics.nist.gov/cuu/Constants/Table/allascii.txt"""
+
+    DATA_FILE = Path(__file__).parent.parent.parent.joinpath('share/codata_2018.txt')
+
+    def __init__(self):
+        self.pc = self.parse_data_file()
+
+        # Non-derived constants
+        self.CONSTANT_JOULE_2_KCAL = 4.184
+        self.CONSTANT_PI = 3.1415926535897932384
+        self.CONSTANT_SQRT_PI = 1.7724538509055160273
+
+        # Conversions
+        self.CONVERT_AU_2_KJ_PER_MOLE = self.pc['hartree-joule relationship']['value'] * self.pc['Avogadro constant']['value'] / 1000
+        self.CONVERT_AU_2_KCAL_PER_MOLE = self.AU_2_KJ / self.JOULE_2_KCAL
+        self.CONVERT_AU_2_EV = self.pc['hartree-ev relationship']['value']
+        self.CONVERT_AU_2_JT_m2 = self.pc['atomic unit of magnetizability']['value']
+        self.CONVERT_AU_2_DEBYE = 2.54174623105  # How derive this one?
+        self.CONVERT_AU_2_CM_m1 = self.pc['hartree-inverse meter relationship']['value'] / 100
+        self.CONVERT_ANGSTROM_2_BOHR = self.pc['atomic unit of length']['value'] * 1.0e-10
+        self.CONVERT_BOHR_2_ANGSTROM = 1.0 / self.ANGSTROM_2_BOHR
+
+        # Constants
+        self.PHYS_FINE_STRUCTURE_CONSTANT = self.pc['fine-structure constant']['value']
+        self.PHYS_SPEED_OF_LIGHT = self.pc['inverse fine-structure constant']['value']
+        self.PHYS_FREE_ELECTRON_G_VALUE = abs(self.pc['electron g factor']['value'])
+        self.PHYS_BOHR_MAGNETON = self.pc['Bohr magneton']['value'] / self.pc['atomic unit of mag. dipole mom.']['value']
+        self.PHYS_NUCLEAR_MAGNETON = self.pc['nuclear magneton']['value'] / self.pc['atomic unit of mag. dipole mom.']['value']
+
+    def parse_data_file(self):
+        data_starts_at = 12 - 1
+        fields = {
+            'quan': {'start': 0, 'stop': 60},
+            'valu': {'start': 60, 'stop': 85},
+            'unce': {'start': 85, 'stop': 110},
+            'unit': {'start': 110, 'stop': -1}
+        }
+
+        with self.DATA_FILE.open() as f:
+            lines = f.readlines()
+
+        constants = {}
+        for line in lines[data_starts_at:]:
+            quan = line[fields['quan']['start']:fields['quan']['stop']].strip()
+            valu = line[fields['valu']['start']:fields['valu']['stop']].strip()
+            # unce = line[fields['unce']['start']:fields['unce']['stop']].strip()
+            unit = line[fields['unit']['start']:fields['unit']['stop']].strip()
+
+            # Remove whitespace and "..." from values
+            valu = valu.replace(' ', '').replace('...', '')
+
+            try:
+                valu = float(valu)
+            except ValueError:
+                raise RuntimeError(f'{quan} could not be converted to float: {valu}')
+
+            # Store values
+            constants[quan] = {'value': valu, 'unit': unit}
+
+        return constants

python/mrchem/validators.py

@@ -27,9 +27,11 @@
 import itertools
 import re
 
-from .physical_constants import PhysicalConstants as PC
+from .physical_constants import PhysicalConstants
 from .periodictable import PeriodicTable, PeriodicTableByZ
 
+pc = PhysicalConstants()
+
 
 class MoleculeValidator:
     """Sanity check routines for the user's Molecule section."""
@@ -127,7 +129,7 @@ def __init__(self, user_dict, origin):
             self.atomic_coords = self.ang2bohr_array(self.atomic_coords)
             self.cavity_coords = self.ang2bohr_array(self.cavity_coords)
             self.cavity_radii = self.ang2bohr_vector(self.cavity_radii)
-            self.cavity_width *= PC.ANGSTROM_2_BOHR
+            self.cavity_width *= pc.CONVERT_ANGSTROM_2_BOHR
 
     def get_coords_in_program_syntax(self):
         """Convert nuclear coordinates from JSON syntax to program syntax."""
@@ -328,10 +330,10 @@ def euclidian_distance(a, b):
     def ang2bohr_array(coords):
         """Helper function. Convert List[List[float]] from angstrom to bohr."""
         return [
-            [c * PC.ANGSTROM_2_BOHR for c in element] for element in coords
+            [c * pc.CONVERT_ANGSTROM_2_BOHR for c in element] for element in coords
         ]
 
     @staticmethod
     def ang2bohr_vector(vec):
         """Helper function. Convert List[float] from angstrom to bohr"""
-        return [el * PC.ANGSTROM_2_BOHR for el in vec]
+        return [el * pc.CONVERT_ANGSTROM_2_BOHR for el in vec]