Input parsing cleaning and PhysConsts from QCelemental

Pipfile

@@ -6,6 +6,7 @@ name = "pypi"
 [packages]
 
 [dev-packages]
+qcelemental = "<= 0.24.0"
 Pygments = "*"
 recommonmark = "*"
 Sphinx = ">=2.0"
@@ -15,3 +16,4 @@ breathe = "*"
 pyyaml = "*"
 yapf = "*"
 parselglossy = ">=0.7"
+ruamel.yaml = ">=0.17.21"

doc/users/schema_input.json

@@ -10,13 +10,15 @@
         "xyz": array[float]                  # Nuclear Cartesian coordinate
       }                                      
     ],                                       
-    "cavity_coords": array[                  # Array of solvation cavities
-      {                                      # (one entry per sphere)
-        "center": array[float],              # Cartesian coordinate of sphere center
-        "radius": float                      # Radius of cavity sphere
-      }                                      
-    ],                                       
-    "cavity_width": float                    # Width of cavity boundary
+    "cavity": {
+      "spheres": array[                      # Array of cavity spheres
+        {                                    # (one entry per sphere)
+          "center": array[float],            # Cartesian coordinate of sphere center
+          "radius": float                    # Radius of cavity sphere
+        }                                      
+      ],                                       
+      "width": float                         # Width of cavity boundary
+    }
   },                                         
   "mpi": {                                   # Section for MPI specification
     "bank_size": int,                        # Number of MPI ranks in memory bank
@@ -294,5 +296,17 @@
         }
       }
     }
+  },
+  "constants": {                             # Physical constants used throughout MRChem
+    "angstrom2bohrs": float,                 # Conversion factor from Angstrom to Bohr
+    "dipmom_au2debye": float,                # Conversion factor from atomic units to Debye
+    "electron_g_factor": float,              # Electron g factor in atomic units
+    "fine_structure_constant": float,        # Fine-structure constant in atomic units
+    "hartree2ev": float,                     # Conversion factor from Hartree to eV
+    "hartree2kcalmol": float,                # Conversion factor from Hartree to kcal/mol
+    "hartree2kjmol": float,                  # Conversion factor from Hartree to kJ/mol
+    "hartree2simagnetizability": float,      # Conversion factor from Hartree to J T^-2
+    "hartree2wavenumbers": float,            # Conversion factor from Hartree to cm^-1
+    "light_speed": float                     # Speed of light in vacuo in atomic units
   }
 }

doc/users/user_ref.rst

@@ -117,6 +117,12 @@ User input reference
     **Predicates**
       - ``50 < value < 100``
 
+   :print_constants: Print table of physical constants used by MRChem. 
+
+    **Type** ``bool``
+
+    **Default** ``False``
+
  :Plotter: Give details regarding the density and orbital plots. Three types of plots are available, line, surface and cube, and the plotting ranges are defined by three vectors (A, B and C) and an origin (O): ``line``: plots on line spanned by A, starting from O. ``surf``: plots on surface spanned by A and B, starting from O. ``cube``: plots on volume spanned by A, B and C, starting from O. 
 
   :red:`Keywords`
@@ -285,6 +291,9 @@ User input reference
 
     **Default** ``1``
 
+    **Predicates**
+      - ``value > 0``
+
    :translate: Translate coordinates such that center of mass coincides with the global gauge origin. 
 
     **Type** ``bool``
@@ -323,12 +332,6 @@ User input reference
  :ZORA: Define required parameters for the ZORA Hamiltonian. 
 
   :red:`Keywords`
-   :light_speed: Adjust speed of light. 
-
-    **Type** ``float``
-
-    **Default** ``-1.0``
-
    :include_nuclear: Include the nuclear potential ``V_nuc`` in the ZORA potential. 
 
     **Type** ``bool``
@@ -863,4 +866,67 @@ User input reference
 
         **Predicates**
           - ``value.lower() in ['exponential']``
-
+
+ :Constants: Physical and mathematical constants used by MRChem
+
+  :red:`Keywords`
+   :hartree2simagnetizability: | Conversion factor for magnetizability from atomic units to SI units  (unit: J T^-2). Affected code: Printed value of the magnetizability property.
+
+    **Type** ``float``
+
+    **Default** ``78.9451185``
+
+   :light_speed: | Speed of light in atomic units  (unit: au). Affected code: Relativistic Hamiltonians (ZORA, etc.)
+
+    **Type** ``float``
+
+    **Default** ``137.035999084``
+
+   :angstrom2bohrs: | Conversion factor for Cartesian coordinates from Angstrom to Bohr  (unit: Å^-1). Affected code: Parsing of input coordinates, printed coordinates
+
+    **Type** ``float``
+
+    **Default** ``1.8897261246257702``
+
+   :hartree2kjmol: | Conversion factor from Hartree to kJ/mol  (unit: kJ mol^-1). Affected code: Printed value of energies.
+
+    **Type** ``float``
+
+    **Default** ``2625.4996394798254``
+
+   :hartree2kcalmol: | Conversion factor from Hartree to kcal/mol  (unit: kcal mol^-1). Affected code: Printed value of energies.
+
+    **Type** ``float``
+
+    **Default** ``627.5094740630558``
+
+   :hartree2ev: | Conversion factor from Hartree to eV  (unit: ev). Affected code: Printed value of energies.
+
+    **Type** ``float``
+
+    **Default** ``27.211386245988``
+
+   :hartree2wavenumbers: | Conversion factor from Hartree to wavenumbers (unit: cm^-1). Affected code: Printed value of frequencies.
+
+    **Type** ``float``
+
+    **Default** ``219474.6313632``
+
+   :fine_structure_constant: | Fine-structure constant in atomic units (unit: au). Affected code: Certain magnetic interaction operators.
+
+    **Type** ``float``
+
+    **Default** ``0.0072973525693``
+
+   :electron_g_factor: | Electron g factor in atomic units (unit: au). Affected code: Certain magnetic interaction operators.
+
+    **Type** ``float``
+
+    **Default** ``-2.00231930436256``
+
+   :dipmom_au2debye: | Conversion factor for dipoles from atomic units to Debye (unit: ?). Affected code: Printed value of dipole moments.
+
+    **Type** ``float``
+
+    **Default** ``2.5417464739297717``
+

python/README.md

@@ -1,13 +1,48 @@
-# How to update the input parser
+# How to update the input parser...
 
-Run:
+## ...when you have updated the `Constants` input section
+The `Constants` section require a bit of special treatment, because the constants are derived from NIST data via an interface
+to `qcelemental`. 
+All constants should be defined in the class `MRChemPhysConstants` located in `python/mrchem/physical_constants.py`.
+The utility script `python/mrchem/update_input_parser.py` will automatically replace the current `Constants`
+section in `python/template.yml` with a new one generated by instantiating `MRChemPhysConstants` and
+reading the subset of constants defined therein.
 
-``` bash
+To run the script:
+
+```bash
+$ cd python/mrchem
+$ python update_input_parser.py
+```
+
+This does three things:
+
+1. Update `template.yml` with the new `Constants` section
+2. Run `parselglossy` to update the input parser
+3. Copy the user reference file to `doc/users/user_ref.rst`
+
+Each of these steps can be deactivated by appropriate flags. For help on running the script, run
+
+```bash
 $ cd python
-$ parselglossy generate --template template.yml --docfile user_ref.rst --doc-header="User input reference" --target="mrchem/input_parser"
+$ python update_input_parser.py -h
 ```
-Remember to also update the documentation:
 
-``` bash
-cp python/mrchem/input_parser/docs/user_ref.rst doc/users/user_ref.rst
+## ...when you have updated any other input section
+If the `Constants` input section is unchanged, you can update the input parser by performing steps 2 and 3 above.
+The utility script can also do this, by specifying the flag `--skip-template`:
+
+```bash
+$ cd python
+$ python update_input_parser --skip-template
+```
+
+This will perform steps 2-3, updating the input parser and the user reference.
+
+You can also perform these steps manually without using the utility script:
+
+```bash
+$ cd python
+$ parselglossy generate --template template.yml --docfile user_ref.rst --doc-header="User input reference" --target="mrchem/input_parser"
+$ cp mrchem/input_parser/docs/user_ref.rst ../doc/users/user_ref.rst
 ```

python/mrchem/CUBEparser.py

@@ -23,17 +23,20 @@
 # For information on the complete list of contributors to MRChem, see:
 # <https://mrchem.readthedocs.io/>
 #
-from .input_parser.plumbing import pyparsing as pp
+
 import os
 from json import dump
 
-BOHR_2_METER = 5.29177210903e-11
-"""Conversion from atomic units of length (Bohr) to meter (CODATA 2018)"""
-ANGSTROM_2_BOHR = 1e-10 / BOHR_2_METER
-#ANGSTROM_2_BOHR = 1.889725989
-"""Conversion factor from Angstrom to Bohr"""
+from .input_parser.plumbing import pyparsing as pp
+
+global pc
+
+
+def write_cube_dict(user_dict):
+    file_dict = user_dict['Files']
+    world_unit = user_dict['world_unit']
+    pc = user_dict['Constants']
 
-def write_cube_dict(file_dict, world_unit):
     all_path_list = []
     all_path_list.append(sort_paths(file_dict["guess_cube_p"]))
     all_path_list.append(sort_paths(file_dict["guess_cube_a"]))
@@ -175,7 +178,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*ANGSTROM_2_BOHR for p in parsed_cube["ORIGIN"]]
+    origin = parsed_cube["ORIGIN"] if (world_unit == "bohr") else [p * pc['angstrom2bohrs'] 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):
@@ -200,9 +203,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*ANGSTROM_2_BOHR for p in parsed_cube["XAXIS"]["VECTOR"]]
-        Y_voxel = [p*ANGSTROM_2_BOHR for p in parsed_cube["YAXIS"]["VECTOR"]]
-        Z_voxel = [p*ANGSTROM_2_BOHR for p in parsed_cube["ZAXIS"]["VECTOR"]]
+        X_voxel = [p * pc['angstrom2bohrs'] for p in parsed_cube["XAXIS"]["VECTOR"]]
+        Y_voxel = [p * pc['angstrom2bohrs'] for p in parsed_cube["YAXIS"]["VECTOR"]]
+        Z_voxel = [p * pc['angstrom2bohrs'] for p in parsed_cube["ZAXIS"]["VECTOR"]]
         Voxel_axes = [X_voxel, Y_voxel, Z_voxel]
 
     # get the atom coordinates
@@ -211,7 +214,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*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['angstrom2bohrs'] 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 = []