diff --git a/Pipfile b/Pipfile
index e3426f3d4..9369a1805 100644
--- a/Pipfile
+++ b/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"
\ No newline at end of file
diff --git a/doc/users/schema_input.json b/doc/users/schema_input.json
index 9abc3d538..61bbc4ffb 100644
--- a/doc/users/schema_input.json
+++ b/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
}
}
diff --git a/doc/users/user_ref.rst b/doc/users/user_ref.rst
index f7c3a80b4..e9c458687 100644
--- a/doc/users/user_ref.rst
+++ b/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']``
-
\ No newline at end of file
+
+ :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``
+
\ No newline at end of file
diff --git a/python/README.md b/python/README.md
index 59f4e5b62..94144711d 100644
--- a/python/README.md
+++ b/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
```
diff --git a/python/mrchem/CUBEparser.py b/python/mrchem/CUBEparser.py
index e8b1574b8..0597dfd5e 100644
--- a/python/mrchem/CUBEparser.py
+++ b/python/mrchem/CUBEparser.py
@@ -23,17 +23,20 @@
# For information on the complete list of contributors to MRChem, see:
#
#
-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 = []
diff --git a/python/mrchem/api.py b/python/mrchem/api.py
index f10146e0b..d0b4176b3 100644
--- a/python/mrchem/api.py
+++ b/python/mrchem/api.py
@@ -25,16 +25,25 @@
import math
-from .helpers import *
-from .CUBEparser import *
+from .helpers import (
+ write_scf_fock,
+ write_scf_guess,
+ write_scf_solver,
+ write_scf_properties,
+ write_scf_plot,
+ write_rsp_calc,
+ parse_wf_method
+)
from .periodictable import PeriodicTable as PT, PeriodicTableByZ as PT_Z
+from .validators import MoleculeValidator
def translate_input(user_dict):
# get the origin in the desired units of measure
origin = user_dict["world_origin"]
+ pc = user_dict['Constants']
if user_dict["world_unit"] == "angstrom":
- origin = [ANGSTROM_2_BOHR * r for r in origin]
+ origin = [pc['angstrom2bohrs'] * r for r in origin]
# prepare bits and pieces
mol_dict = write_molecule(user_dict, origin)
@@ -53,6 +62,7 @@ def translate_input(user_dict):
"molecule": mol_dict,
"scf_calculation": scf_dict,
"rsp_calculations": rsp_dict,
+ "constants": user_dict["Constants"]
}
return program_dict
@@ -111,119 +121,17 @@ def write_mra(user_dict, mol_dict):
def write_molecule(user_dict, origin):
# Translate into program syntax
- coords_raw = user_dict["Molecule"]["coords"]
- coords_dict = []
- for line in coords_raw.split("\n"):
- sp = line.split()
- if len(sp) > 0:
-
- try:
- int(sp[0])
- except:
- atom = sp[0].lower()
- else:
- atom = PT_Z[int(sp[0])].symbol.lower()
-
- xyz = list(map(float, sp[1:]))
- if len(xyz) != 3:
- raise RuntimeError(f"Invalid coordinate: {atom.upper()} {str(xyz)}")
- coords_dict.append({"atom": atom, "xyz": xyz})
-
- # Convert angstrom -> bohr
- if user_dict["world_unit"] == "angstrom":
- for coord in coords_dict:
- coord["xyz"] = [ANGSTROM_2_BOHR * r for r in coord["xyz"]]
-
- # Check for singularity
- for a in range(len(coords_dict)):
- for b in range(a + 1, len(coords_dict)):
- A = coords_dict[a]
- B = coords_dict[b]
- R = math.sqrt(sum([(a - b) ** 2 for a, b in zip(A["xyz"], B["xyz"])]))
- if R < 1.0e-6:
- msg = f"ABORT: Atoms are too close\n {A['atom']}: {str(A['xyz'])}\n {B['atom']}: {str(B['xyz'])}"
- raise RuntimeError(msg)
- elif R < 1.0e-3:
- msg = f"WARNING: Atoms are too close\n {A['atom']}: {str(A['xyz'])}\n {B['atom']}: {str(B['xyz'])}"
- print(msg)
-
- # Translate center of mass to origin
- if user_dict["Molecule"]["translate"]:
- # Calc center of mass
- M = 0.0
- CoM = [0.0, 0.0, 0.0]
- for coord in coords_dict:
- m = PT[coord["atom"]].mass
- M += m
- CoM = [m * r + o for r, o in zip(coord["xyz"], CoM)]
- CoM = [x / M - o for x, o in zip(CoM, origin)]
-
- # Translate coords
- for coord in coords_dict:
- coord["xyz"] = [r - o for r, o in zip(coord["xyz"], CoM)]
-
- # initialize the cavity
- cav_coords_dict = []
- if len(user_dict["Environment"]["Cavity"]["spheres"]) > 0:
- cav_coords_raw = user_dict["Environment"]["Cavity"]["spheres"]
- for line in cav_coords_raw.split('\n'):
- sp = line.split()
- if len(sp) > 0:
- center = list(map(float, sp[:-1]))
- radius = float(sp[-1])
- if len(center) != 3:
- raise RuntimeError(f"Invalid coordinate: {center}")
- cav_coords_dict.append({
- "center": center,
- "radius": radius
- })
- else:
- for atom in coords_dict:
- radius = float(PT[atom["atom"]].radius)
- cav_coords_dict.append({
- "center": atom["xyz"],
- "radius": radius
- })
-
- if user_dict["world_unit"] == "angstrom" :
- for coord in cav_coords_dict:
- coord["center"] = [ANGSTROM_2_BOHR * r for r in coord["center"]]
- coord["radius"] *= ANGSTROM_2_BOHR
- cavity_width = user_dict["Environment"]["Cavity"]["cavity_width"]*ANGSTROM_2_BOHR
- else:
- cavity_width = user_dict["Environment"]["Cavity"]["cavity_width"]
-
- # Make sure the specified charge and multiplicity make sense
- mult = user_dict["Molecule"]["multiplicity"]
- charge = user_dict["Molecule"]["charge"]
- restricted = user_dict["WaveFunction"]["restricted"]
-
- # Collect number of electrons
- n_electrons = sum([PT[atom["atom"]].Z for atom in coords_dict]) - charge
- n_unpaired = mult - 1
-
- # Helper function
- parity = lambda n: 'even' if n % 2 == 0 else 'odd'
-
- # Check for unphysical multiplicity (not likely to be much of a problem, but still...)
- if n_unpaired > n_electrons:
- raise RuntimeError(f"The specified multiplicity requires more unpaired electrons ({mult - 1}) than are available ({n_electrons}))")
-
- # Check for restricted open-shell
- elif restricted and n_unpaired > 0:
- raise RuntimeError("Restricted open-shell calculations are not available")
-
- # Check for invalid spin multiplicity
- elif parity(n_electrons) == parity(mult):
- raise RuntimeError(f"The specified multiplicity ({parity(mult)}) is not compatible with the number of electrons ({parity(n_electrons)})")
-
+ mol = MoleculeValidator(user_dict, origin)
mol_dict = {
- "multiplicity": mult,
- "charge": charge,
- "coords": coords_dict,
- "cavity_coords": cav_coords_dict,
- "cavity_width": cavity_width
+ "multiplicity": mol.mult,
+ "charge": mol.charge,
+ "coords": mol.get_coords_in_program_syntax(),
}
+ if user_dict["Environment"]["run_environment"]:
+ mol_dict["cavity"] = {
+ "spheres": mol.get_cavity_in_program_syntax(),
+ "width": mol.cavity_width,
+ }
return mol_dict
diff --git a/python/mrchem/helpers.py b/python/mrchem/helpers.py
index 56b2b0104..8eadda3cb 100644
--- a/python/mrchem/helpers.py
+++ b/python/mrchem/helpers.py
@@ -23,13 +23,8 @@
#
#
-from .CUBEparser import *
+from .CUBEparser import write_cube_dict
-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"""
# yapf: disable
SHORTHAND_FUNCTIONALS = [
'svwn3',
@@ -58,7 +53,6 @@ def write_scf_fock(user_dict, wf_dict, origin):
# ZORA
if user_dict["WaveFunction"]["relativity"].lower() == "zora":
fock_dict["zora_operator"] = {
- "light_speed": user_dict["ZORA"]["light_speed"],
"include_nuclear": user_dict["ZORA"]["include_nuclear"],
"include_coulomb": user_dict["ZORA"]["include_coulomb"],
"include_xc": user_dict["ZORA"]["include_xc"]
@@ -164,7 +158,7 @@ def write_scf_guess(user_dict, wf_dict):
file_dict = user_dict["Files"]
if (guess_type == 'cube'):
- write_cube_dict(file_dict, user_dict["world_unit"])
+ write_cube_dict(user_dict)
vector_dir = file_dict["cube_vectors"]
guess_dict = {
@@ -254,7 +248,7 @@ def write_scf_plot(user_dict):
plot_dict["plotter"] = user_dict["Plotter"]
if user_dict["world_unit"] == "angstrom":
plot_dict["plotter"] = {
- k: [ANGSTROM_2_BOHR * r for r in plot_dict["plotter"][k]]
+ k: [user_dict['Constants']['angstrom2bohrs'] * r for r in plot_dict["plotter"][k]]
for k in plot_dict["plotter"].keys()
}
return plot_dict
diff --git a/python/mrchem/input_parser/README.md b/python/mrchem/input_parser/README.md
index 53ac9c1e9..83ddb199d 100644
--- a/python/mrchem/input_parser/README.md
+++ b/python/mrchem/input_parser/README.md
@@ -1,2 +1,2 @@
-This file was automatically generated by parselglossy on 2022-03-18
+This file was automatically generated by parselglossy on 2022-04-25
Editing is *STRONGLY DISCOURAGED*
diff --git a/python/mrchem/input_parser/__init__.py b/python/mrchem/input_parser/__init__.py
index 7266e970e..3a5186327 100644
--- a/python/mrchem/input_parser/__init__.py
+++ b/python/mrchem/input_parser/__init__.py
@@ -1,4 +1,4 @@
# -*- coding: utf-8 -*-
-# This file was automatically generated by parselglossy on 2022-03-18
+# This file was automatically generated by parselglossy on 2022-04-25
# Editing is *STRONGLY DISCOURAGED*
diff --git a/python/mrchem/input_parser/api.py b/python/mrchem/input_parser/api.py
index 7b98cccac..a41ae6be6 100644
--- a/python/mrchem/input_parser/api.py
+++ b/python/mrchem/input_parser/api.py
@@ -1,6 +1,6 @@
# -*- coding: utf-8 -*-
-# This file was automatically generated by parselglossy on 2022-03-18
+# This file was automatically generated by parselglossy on 2022-04-25
# Editing is *STRONGLY DISCOURAGED*
from copy import deepcopy
@@ -159,7 +159,10 @@ def stencil() -> JSONDict:
{ 'default': 75,
'name': 'print_width',
'predicates': ['50 < value < 100'],
- 'type': 'int'}],
+ 'type': 'int'},
+ { 'default': False,
+ 'name': 'print_constants',
+ 'type': 'bool'}],
'name': 'Printer'},
{ 'keywords': [ { 'default': 'plots',
'name': 'path',
@@ -255,6 +258,7 @@ def stencil() -> JSONDict:
'type': 'int'},
{ 'default': 1,
'name': 'multiplicity',
+ 'predicates': ['value > 0'],
'type': 'int'},
{ 'default': False,
'name': 'translate',
@@ -293,10 +297,7 @@ def stencil() -> JSONDict:
"'nzora']"],
'type': 'str'}],
'name': 'WaveFunction'},
- { 'keywords': [ { 'default': -1.0,
- 'name': 'light_speed',
- 'type': 'float'},
- { 'default': True,
+ { 'keywords': [ { 'default': True,
'name': 'include_nuclear',
'type': 'bool'},
{ 'default': True,
@@ -579,4 +580,35 @@ def stencil() -> JSONDict:
'in '
"['exponential']"],
'type': 'str'}],
- 'name': 'Permittivity'}]}]}
+ 'name': 'Permittivity'}]},
+ { 'keywords': [ { 'default': 78.9451185,
+ 'name': 'hartree2simagnetizability',
+ 'type': 'float'},
+ { 'default': 137.035999084,
+ 'name': 'light_speed',
+ 'type': 'float'},
+ { 'default': 1.8897261246257702,
+ 'name': 'angstrom2bohrs',
+ 'type': 'float'},
+ { 'default': 2625.4996394798254,
+ 'name': 'hartree2kjmol',
+ 'type': 'float'},
+ { 'default': 627.5094740630558,
+ 'name': 'hartree2kcalmol',
+ 'type': 'float'},
+ { 'default': 27.211386245988,
+ 'name': 'hartree2ev',
+ 'type': 'float'},
+ { 'default': 219474.6313632,
+ 'name': 'hartree2wavenumbers',
+ 'type': 'float'},
+ { 'default': 0.0072973525693,
+ 'name': 'fine_structure_constant',
+ 'type': 'float'},
+ { 'default': -2.00231930436256,
+ 'name': 'electron_g_factor',
+ 'type': 'float'},
+ { 'default': 2.5417464739297717,
+ 'name': 'dipmom_au2debye',
+ 'type': 'float'}],
+ 'name': 'Constants'}]}
diff --git a/python/mrchem/input_parser/cli.py b/python/mrchem/input_parser/cli.py
index 1f2671e87..1ea10a97b 100644
--- a/python/mrchem/input_parser/cli.py
+++ b/python/mrchem/input_parser/cli.py
@@ -1,6 +1,6 @@
# -*- coding: utf-8 -*-
-# This file was automatically generated by parselglossy on 2022-03-18
+# This file was automatically generated by parselglossy on 2022-04-25
# Editing is *STRONGLY DISCOURAGED*
import argparse
diff --git a/python/mrchem/input_parser/docs/user_ref.rst b/python/mrchem/input_parser/docs/user_ref.rst
index f7c3a80b4..e9c458687 100644
--- a/python/mrchem/input_parser/docs/user_ref.rst
+++ b/python/mrchem/input_parser/docs/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']``
-
\ No newline at end of file
+
+ :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``
+
\ No newline at end of file
diff --git a/python/mrchem/input_parser/plumbing/lexer.py b/python/mrchem/input_parser/plumbing/lexer.py
index 95bb2a29f..c52aa63c7 100644
--- a/python/mrchem/input_parser/plumbing/lexer.py
+++ b/python/mrchem/input_parser/plumbing/lexer.py
@@ -1,6 +1,6 @@
# -*- coding: utf-8 -*-
-# This file was automatically generated by parselglossy on 2022-03-18
+# This file was automatically generated by parselglossy on 2022-04-25
# Editing is *STRONGLY DISCOURAGED*
import json
diff --git a/python/mrchem/physical_constants.py b/python/mrchem/physical_constants.py
new file mode 100644
index 000000000..c71f68d70
--- /dev/null
+++ b/python/mrchem/physical_constants.py
@@ -0,0 +1,146 @@
+#
+# MRChem, a numerical real-space code for molecular electronic structure
+# calculations within the self-consistent field (SCF) approximations of quantum
+# chemistry (Hartree-Fock and Density Functional Theory).
+# Copyright (C) 2022 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors.
+#
+# This file is part of MRChem.
+#
+# MRChem is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# MRChem is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with MRChem. If not, see .
+#
+# For information on the complete list of contributors to MRChem, see:
+#
+#
+
+import math
+from pydoc import doc
+from types import SimpleNamespace
+
+from qcelemental.physical_constants.context import PhysicalConstantsContext
+
+
+class MRChemPhysConstants:
+ """Class holding physical and math constants relevant to MRChem.
+ The constants are fetched from QCElemental, and new constants are
+ derived from those in QCElemental (as far as possible).
+
+ Source in ascii:
+ https://physics.nist.gov/cuu/Constants/Table/allascii.txt
+ """
+
+ # Explicitly defined constants here
+ # (not present in qcelemental or cannot be derived)
+ PI = 3.1415926535897932384
+ HARTREE2SIMAGNETIZABILITY = 78.9451185 # TODO: We should be able to derive this one
+
+ def __init__(self):
+ """Here we extract those constants that we need to run any MRChem calculation.
+ Those that are not directly available in QCElemental, we compute by via the existing
+ constants in QCElementlal.
+ """
+ self.context = "CODATA2018"
+ self.nist = "https://physics.nist.gov/cuu/Constants/Table/allascii.txt"
+ self.qce = PhysicalConstantsContext(context=self.context)
+ self.data = []
+
+ ########################
+ # Define our constants
+ ########################
+
+ # Convert from au to SI units for magnetizability
+ name = 'hartree2simagnetizability'
+ unit = 'J T^-2'
+ value = self.HARTREE2SIMAGNETIZABILITY
+ docstring = f"| Conversion factor for magnetizability from atomic units to SI units (unit: {unit}). Affected code: Printed value of the magnetizability property."
+ self.add_constant(name, unit, value, docstring)
+
+ # Speed of light in atomic units
+ name = 'light_speed'
+ unit = 'au'
+ value = self.qce.c_au
+ docstring = f"| Speed of light in atomic units (unit: {unit}). Affected code: Relativistic Hamiltonians (ZORA, etc.)"
+ self.add_constant(name, unit, value, docstring)
+
+ # Convert from Angstrom to Bohr
+ name = 'angstrom2bohrs'
+ unit = 'Å^-1'
+ value = 1.0 / self.qce.bohr2angstroms
+ docstring = f"| Conversion factor for Cartesian coordinates from Angstrom to Bohr (unit: {unit}). Affected code: Parsing of input coordinates, printed coordinates"
+ self.add_constant(name, unit, value, docstring)
+
+ # Convert from Hartree to kJ/mol
+ name = 'hartree2kjmol'
+ unit = 'kJ mol^-1'
+ value = self.qce.hartree2kJmol
+ docstring = f"| Conversion factor from Hartree to kJ/mol (unit: {unit}). Affected code: Printed value of energies."
+ self.add_constant(name, unit, value, docstring)
+
+ # Convert from Hartree to kcal/mol
+ name = 'hartree2kcalmol'
+ unit = 'kcal mol^-1'
+ value = self.qce.hartree2kcalmol
+ docstring = f"| Conversion factor from Hartree to kcal/mol (unit: {unit}). Affected code: Printed value of energies."
+ self.add_constant(name, unit, value, docstring)
+
+ # Convert from Hartree to eV
+ name = 'hartree2ev'
+ unit = 'ev'
+ value = self.qce.hartree2ev
+ docstring = f"| Conversion factor from Hartree to eV (unit: {unit}). Affected code: Printed value of energies."""
+ self.add_constant(name, unit, value, docstring)
+
+ # Convert from Hartree to cm-1
+ name = 'hartree2wavenumbers'
+ unit = 'cm^-1'
+ value = self.qce.hartree2wavenumbers
+ docstring = f"| Conversion factor from Hartree to wavenumbers (unit: {unit}). Affected code: Printed value of frequencies."
+ self.add_constant(name, unit, value, docstring)
+
+ # Fine-structure constant in atomic units
+ name = 'fine_structure_constant'
+ unit = 'au'
+ value = self.qce.fine_structure_constant
+ docstring = f"| Fine-structure constant in atomic units (unit: {unit}). Affected code: Certain magnetic interaction operators."
+ self.add_constant(name, unit, value, docstring)
+
+ # Electron g factor in atomic units
+ name = 'electron_g_factor'
+ unit = 'au'
+ value = self.qce.electron_g_factor
+ docstring = f"| Electron g factor in atomic units (unit: {unit}). Affected code: Certain magnetic interaction operators."
+ self.add_constant(name, unit, value, docstring)
+
+ # Convert from atomic units to Debye
+ name = 'dipmom_au2debye'
+ unit = '?'
+ value = self.qce.dipmom_au2debye
+ docstring = f"| Conversion factor for dipoles from atomic units to Debye (unit: {unit}). Affected code: Printed value of dipole moments."
+ self.add_constant(name, unit, value, docstring)
+
+ # Set our constants to instance attributes for easy access
+ for name, _, value, _ in self.data:
+ self.__setattr__(name.lower(), float(value))
+
+ def add_constant(self, name, unit, value, docstring):
+ """Add constant to `data` instance attribute."""
+ self.data.append((name, unit, value, docstring))
+
+ def print_constants_for_tests(self, varname='testConstants'):
+ """Helper function for printing constants for copy/pasting into the c++ code.
+ We need to store the constants internally for the unit tests to pass."""
+ for name, _, value, _ in sorted(self.data, key=lambda x: x[0]):
+ print(f'{{\"{name}\", {value}}},')
+
+if __name__ == '__main__':
+ MRChemPhysConstants().print_constants_for_tests()
diff --git a/python/mrchem/update_input_parser.py b/python/mrchem/update_input_parser.py
new file mode 100644
index 000000000..5207ab2a0
--- /dev/null
+++ b/python/mrchem/update_input_parser.py
@@ -0,0 +1,117 @@
+#!/usr/bin/env python
+
+#
+# MRChem, a numerical real-space code for molecular electronic structure
+# calculations within the self-consistent field (SCF) approximations of quantum
+# chemistry (Hartree-Fock and Density Functional Theory).
+# Copyright (C) 2022 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors.
+#
+# This file is part of MRChem.
+#
+# MRChem is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# MRChem is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with MRChem. If not, see .
+#
+# For information on the complete list of contributors to MRChem, see:
+#
+#
+
+
+from pathlib import Path
+from ruamel.yaml import YAML
+import subprocess
+import os
+import shutil
+import argparse
+
+from physical_constants import MRChemPhysConstants
+
+root = Path.cwd().parent
+target = root.joinpath('mrchem', 'input_parser')
+
+yaml = YAML()
+yaml.indent(sequence=4, mapping=2, offset=2)
+
+def update_constants():
+ pc = MRChemPhysConstants()
+ f_template = root.joinpath("template.yml")
+ template = yaml.load(f_template)
+
+ new = {
+ "keywords": template["keywords"],
+ "sections": [section for section in template['sections'] if section['name'] != "Constants"]}
+
+ # Build new constants section
+ constants = {
+ "name": "Constants",
+ "docstring": "Physical and mathematical constants used by MRChem",
+ "keywords": [
+ {
+ "name": name,
+ "default": value,
+ "type": "float",
+ "docstring": docstring
+ } for name, _, value, docstring in pc.data
+ ]
+ }
+
+ # Construct the full template file and dump
+ new["sections"].append(constants)
+ yaml.dump(new, f_template)
+
+
+def run_parselglossy():
+ os.chdir(root)
+ cmd = [
+ 'parselglossy',
+ 'generate',
+ '--template', 'template.yml',
+ '--docfile', 'user_ref.rst',
+ '--doc-header', 'User input reference',
+ '--target', target
+ ]
+
+ subprocess.call(cmd)
+
+def update_doc():
+ src = target.joinpath('docs/user_ref.rst')
+ dst = root.parent.joinpath('doc/users/user_ref.rst')
+ shutil.copyfile(src, dst)
+
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser(description="CLI for synchronizing the template with current constants, plus updating the input parser and user ref.")
+ parser.add_argument('-sT', '--skip-template', action='store_true', help="Do not update constants section in template.yml")
+ parser.add_argument('-sP', '--skip-parselglossy', action='store_true', help='Do not update input parser with parselglossy')
+ parser.add_argument('-sD', '--skip-doc', action='store_true', help='Do not update the user reference file')
+ parser.add_argument('-D', '--dryrun', action='store_true', help='Skip all actions')
+ args = parser.parse_args()
+
+ if args.dryrun:
+ args.skip_template = True
+ args.skip_parselglossy = True
+ args.skip_doc = True
+
+ if not args.skip_template:
+ print(f'{"Updating template":20} ... ', end='')
+ update_constants()
+ print('done')
+
+ if not args.skip_parselglossy:
+ print(f'{"Running parselglossy":20} ... ', end='')
+ run_parselglossy()
+ print('done')
+
+ if not args.skip_doc:
+ print(f'{"Updating user ref":20} ... ', end='')
+ update_doc()
+ print('done')
diff --git a/python/mrchem/validators.py b/python/mrchem/validators.py
new file mode 100644
index 000000000..0979e8b89
--- /dev/null
+++ b/python/mrchem/validators.py
@@ -0,0 +1,335 @@
+#
+# MRChem, a numerical real-space code for molecular electronic structure
+# calculations within the self-consistent field (SCF) approximations of quantum
+# chemistry (Hartree-Fock and Density Functional Theory).
+# Copyright (C) 2022 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors.
+#
+# This file is part of MRChem.
+#
+# MRChem is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# MRChem is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with MRChem. If not, see .
+#
+# For information on the complete list of contributors to MRChem, see:
+#
+#
+
+import math
+import itertools
+import re
+
+from .periodictable import PeriodicTable, PeriodicTableByZ
+
+
+class MoleculeValidator:
+ """Sanity check routines for the user's Molecule section."""
+
+ # Thresholds
+ THRESHOLD_NUCLEAR_SINGULARITY_ERROR = 1.0e-6
+ THRESHOLD_NUCLEAR_SINGULARITY_WARNING = 1.0e-3
+
+ # Unit strings
+ UNIT_ANGSTROM = 'angstrom'
+ UNIT_BOHR = 'bohr'
+
+ # Error/warning messages
+ ERROR_MESSAGE_ATOMIC_COORDINATES = lambda self, details: f'ABORT: INVALID ATOMIC COORDINATES: {details}'
+ ERROR_MESSAGE_ATOMIC_SYMBOLS = lambda self, details: f'ABORT: INVALID ATOMIC SYMBOLS: {details}'
+ ERROR_MESSAGE_CAVITY_COORDINATES = lambda self, details: f'ABORT: INVALID CAVITY COORDINATES: {details}'
+ ERROR_MESSAGE_CAVITY_RADII = lambda self, details: f'ABORT: INVALID CAVITY RADII: {details}'
+
+ ERROR_MESSAGE_NUCLEAR_SINGULARITY = lambda self, details: f'ABORT: SOME ATOMS TOO CLOSE (norm < {MoleculeValidator.THRESHOLD_NUCLEAR_SINGULARITY_ERROR}):\n{details}'
+ WARNING_MESSAGE_NUCLEAR_SINGULARITY = lambda self, details: f'WARNING: SOME ATOMS VERY CLOSE (norm < {MoleculeValidator.THRESHOLD_NUCLEAR_SINGULARITY_WARNING}):\n{details}'
+
+ ERROR_INCOMPATIBLE_MULTIPLICITY = lambda self, details: f'ABORT: INCOMPATIBLE MULTIPLICITY: {details}'
+ ERROR_UNPHYSICAL_MULTIPLICITY = lambda self, details: f'ABORT: UNPHYSICAL MULTIPLICITY: {details}'
+ ERROR_UNPHYSICAL_CHARGE = lambda self, details: f'ABORT: UNPHYSICAL CHARGE: {details}'
+
+ ERROR_RESTRICTED_OPEN_SHELL = 'ABORT: Restricted open-shell not implemented'
+
+
+ def __init__(self, user_dict, origin):
+ """
+ Raises RunTimeError with helpful messages when invalid format
+ or unphysical input value are detected:
+
+ Atomic coordinates
+ - Correct XYZ format checked with regexes (both atomic symbols
+ and numbers are valid atom identifiers, and can be used
+ interchancably in the same input)
+ - Nuclear singularities
+ - Atomic symbols checked against periodic table
+
+ Cavity spheres
+ - Correct format checked with regexes
+ - Negative radii not allowed
+
+ Electrons, charge and multiplicity
+ - charge > Z not allowed
+ - n_unpaired > n_electrons not allowed
+ - multiplicity and n_electrons both odd/even not allowed
+ - restricted open-shell not allowed
+
+ Unit conversions
+ - convert to bohr if user specified angstrom
+
+ Parameters
+ ----------
+ user_dict: dict, dictionary of user input
+ origin: List[float], origin to be used in the calcuation
+ """
+ self.user_dict = user_dict
+ self.origin = origin
+ self.unit = user_dict['world_unit']
+ self.pc = user_dict['Constants']
+
+ # Molecule related data
+ self.user_mol = user_dict['Molecule']
+ self.charge = self.user_mol['charge']
+ self.mult = self.user_mol['multiplicity']
+ self.do_translate = self.user_mol['translate']
+ self.coords_raw = self.user_mol['coords']
+
+ # Cavity related data
+ self.cavity_dict = user_dict['Environment']['Cavity']
+ self.spheres_raw = self.cavity_dict['spheres']
+ self.cavity_width = self.cavity_dict['cavity_width']
+ self.has_sphere_input = len(self.spheres_raw.strip().splitlines()) > 0
+
+ # Validate atomic coordinates
+ self.atomic_symbols, self.atomic_coords = self.validate_atomic_coordinates()
+
+ # Translate center of mass if requested
+ # We must test for translation before validating the cavity,
+ # in case the nuclear coordinates are to be used for the
+ # sphere centers
+ if self.do_translate:
+ self.atomic_coords = self.translate_com_to_origin()
+
+ # Validate cavity spheres
+ self.cavity_radii, self.cavity_coords = self.validate_cavity()
+
+ # Perform some sanity checks
+ self.check_for_nuclear_singularities()
+ self.check_for_invalid_electronic_configuration()
+
+ # Convert to bohrs if user gave angstroms
+ if self.unit == self.UNIT_ANGSTROM:
+ 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 *= self.pc['angstrom2bohrs']
+
+ def get_coords_in_program_syntax(self):
+ """Convert nuclear coordinates from JSON syntax to program syntax."""
+ return [
+ {
+ "atom": label,
+ "xyz": coord
+ } for label, coord in zip(self.atomic_symbols, self.atomic_coords)
+ ]
+
+ def get_cavity_in_program_syntax(self):
+ """Convert cavity spheres from JSON syntax to program syntax."""
+ # Use sphere coordinates and radii if given
+ if self.has_sphere_input:
+ return [
+ {
+ "center": center,
+ "radius": radius
+ } for center, radius in zip(self.cavity_coords, self.cavity_radii)
+ ]
+ # If not build spheres from nuclear coordinates and default radii
+ else:
+ return [
+ {
+ "center": coord,
+ "radius": PeriodicTable[label.lower()].radius
+ } for label, coord in zip(self.atomic_symbols, self.atomic_coords)
+ ]
+
+ def validate_atomic_coordinates(self):
+ """Parse the $coords block and ensure correct formatting."""
+ # Regex components
+ line_start = r'^'
+ line_end = r'$'
+ symbol = r'[a-zA-Z]{1,3}'
+ decimal = r'[+-]?([0-9]+\.?[0-9]*|\.[0-9]+)'
+ integer = r'[0-9]+'
+ one_or_more_whitespace = r'[\s]+'
+ zero_or_more_whitespace = r'[\s]*'
+
+ # Build regex
+ atom_with_symbol = line_start + zero_or_more_whitespace + symbol + (one_or_more_whitespace + decimal)*3 + zero_or_more_whitespace + line_end
+ atom_with_number = line_start + zero_or_more_whitespace + integer + (one_or_more_whitespace + decimal)*3 + zero_or_more_whitespace + line_end
+
+ p_with_symbol = re.compile(atom_with_symbol)
+ p_with_number = re.compile(atom_with_number)
+
+ # Parse coordinates
+ coords = []
+ labels = []
+ bad_atoms = []
+ for atom in self.coords_raw.strip().splitlines():
+ match_symbol = p_with_symbol.match(atom)
+ match_number = p_with_number.match(atom)
+ if match_symbol:
+ g = match_symbol.group()
+ labels.append(g.split()[0].strip().lower())
+ coords.append([float(c.strip()) for c in g.split()[1:]])
+ elif match_number:
+ g = match_number.group()
+ Z = int(g.split()[0].strip())
+ labels.append(PeriodicTableByZ[Z].symbol.lower())
+ coords.append([float(c.strip()) for c in g.split()[1:]])
+ else:
+ bad_atoms.append(atom)
+
+ if bad_atoms:
+ newline = '\n'
+ raise RuntimeError(self.ERROR_MESSAGE_ATOMIC_COORDINATES(
+ f'One or more atomic coordinates had an invalid input format:\n{newline.join(bad_atoms)}'
+ ))
+
+ # Check that the atomic symbols represent valid elements
+ fltr = filter(lambda x: x not in PeriodicTable, labels)
+ if any(list(fltr)):
+ newline = '\n'
+ raise RuntimeError(self.ERROR_MESSAGE_ATOMIC_SYMBOLS(
+ f'One or more invalid atomic symbols:\n{newline.join(set(fltr))}'
+ ))
+
+ return labels, coords
+
+
+ def validate_cavity(self):
+ """Parse the $spheres block and ensure correct formatting."""
+ # Regex components
+ line_start = r'^'
+ line_end = r'$'
+ decimal = r'[+-]?([0-9]+\.?[0-9]*|\.[0-9]+)'
+ one_or_more_whitespace = r'[\s]+'
+ zero_or_more_whitespace = r'[\s]*'
+
+ # Build regex
+ valid_sphere = line_start + zero_or_more_whitespace + decimal + (one_or_more_whitespace + decimal)*3 + zero_or_more_whitespace + line_end
+ p = re.compile(valid_sphere)
+
+ # Parse spheres
+ coords = []
+ radii = []
+ bad_spheres = []
+
+ for sphere in self.spheres_raw.strip().splitlines():
+ match = p.match(sphere)
+ if match:
+ g = match.group()
+ radii.append(float(g.split()[-1].strip()))
+ coords.append([float(c.strip()) for c in g.split()[:-1]])
+ else:
+ bad_spheres.append(sphere)
+
+ if bad_spheres:
+ newline = '\n'
+ raise RuntimeError(self.ERROR_MESSAGE_CAVITY_COORDINATES(
+ f'One or more cavity spheres had an invalid input format:\n{newline.join(bad_spheres)}'
+ ))
+
+ # Check for negative radii
+ if any([r < 0 for r in radii]):
+ raise RuntimeError(self.ERROR_MESSAGE_CAVITY_RADII(
+ 'Cavity radii cannot be negative'
+ ))
+
+ return radii, coords
+
+ def check_for_nuclear_singularities(self):
+ """Check for singularities in the nuclear positions."""
+ # Bad pairs will be stored here
+ error_pairs = []
+ warning_pairs = []
+
+ # Loop over all unique atom pairs and compute euclidian distance
+ for (ca, la), (cb, lb) in itertools.combinations(zip(self.atomic_coords, self.atomic_symbols), 2):
+ pair_label = f'{la}: {ca}\n{lb}: {cb}'
+ R = self.euclidian_distance(ca, cb)
+
+ # Compare distance to internal thresholds
+ if R < self.THRESHOLD_NUCLEAR_SINGULARITY_ERROR:
+ error_pairs.append(pair_label)
+ elif R < self.THRESHOLD_NUCLEAR_SINGULARITY_WARNING:
+ warning_pairs.append(pair_label)
+
+ # Print warnings and raise exception if necessary
+ if warning_pairs:
+ msg = self.WARNING_MESSAGE_NUCLEAR_SINGULARITY('\n\n'.join(warning_pairs))
+ print(msg)
+
+ if error_pairs:
+ msg = self.ERROR_MESSAGE_NUCLEAR_SINGULARITY('\n\n'.join(error_pairs))
+ raise RuntimeError(msg)
+
+ def check_for_invalid_electronic_configuration(self):
+ """Check that the number of electrons and spin multiplicity are compatible.
+ Also check for restricted open-shell calculation."""
+ restricted = self.user_dict['WaveFunction']['restricted']
+ Z = sum([PeriodicTable[atom.lower()].Z for atom in self.atomic_symbols])
+ n_electrons = Z - self.charge
+ n_unpaired = self.mult - 1
+
+ # Helper function
+ parity = lambda n: 'even' if n % 2 == 0 else 'odd'
+
+ # Check for impossible charge
+ if self.charge > Z:
+ raise RuntimeError(self.ERROR_UNPHYSICAL_CHARGE(
+ f'The specified charge ({self.charge}) cannot be larger than the nuclear charge ({Z})'
+ ))
+
+ # Check for unphysical multiplicity
+ elif n_unpaired > n_electrons:
+ raise RuntimeError(self.ERROR_UNPHYSICAL_MULTIPLICITY(
+ f"The specified multiplicity requires more unpaired electrons ({self.mult - 1}) than are available ({n_electrons}))."
+ ))
+
+ # Check for invalid spin multiplicity
+ elif parity(n_electrons) == parity(self.mult):
+ raise RuntimeError(self.ERROR_INCOMPATIBLE_MULTIPLICITY(
+ f"The specified multiplicity ({parity(self.mult)}) is not compatible with the number of electrons ({parity(n_electrons)})"
+ ))
+
+ # Check for restricted open-shell
+ elif restricted and n_unpaired > 0:
+ raise RuntimeError(self.ERROR_RESTRICTED_OPEN_SHELL)
+
+ def translate_com_to_origin(self):
+ """Translate center of mass to the origin."""
+ masses = [PeriodicTable[label.lower()].mass for label in self.atomic_symbols]
+ M = sum(masses)
+ return [masses[i] * (self.atomic_coords[i] - self.origin[i]) / M for i in range(3)]
+
+ @staticmethod
+ def euclidian_distance(a, b):
+ """Helper function for the nuclear singularies validation.
+ Computes the euclidian distance between two vectors, a and b."""
+ squared_deviations = [(a[i] - b[i])**2 for i in range(3)]
+ return math.sqrt(sum(squared_deviations))
+
+ def ang2bohr_array(self, coords):
+ """Convert List[List[float]] from angstrom to bohr."""
+ return [
+ [c * self.pc['angstrom2bohrs'] for c in element] for element in coords
+ ]
+
+ def ang2bohr_vector(self, vec):
+ """Convert List[float] from angstrom to bohr"""
+ return [el * self.pc['angstrom2bohrs'] for el in vec]
diff --git a/python/template.yml b/python/template.yml
index 98486bd5a..7e8845508 100644
--- a/python/template.yml
+++ b/python/template.yml
@@ -2,23 +2,23 @@ keywords:
- name: world_prec
type: float
predicates:
- - '1.0e-10 < value < 1.0'
+ - 1.0e-10 < value < 1.0
docstring: |
Overall relative precision in the calculation.
- name: world_size
type: int
default: -1
predicates:
- - 'value <= 10'
+ - value <= 10
docstring: |
Total size of computational domain given as 2**(``world_size``). Always cubic
and symmetric around the origin. Negative value means it will be computed
from the molecular geometry.
- name: world_unit
type: str
- default: "bohr"
+ default: bohr
predicates:
- - 'value.lower() in ["bohr", "angstrom"]'
+ - value.lower() in ["bohr", "angstrom"]
docstring: |
Length unit for *all* coordinates given in user input. Everything will be
converted to atomic units (bohr) before the main executable is launched,
@@ -27,7 +27,7 @@ keywords:
type: List[float]
default: [0.0, 0.0, 0.0]
predicates:
- - 'len(value) == 3'
+ - len(value) == 3
docstring: |
Global gauge origin of the calculation.
sections:
@@ -37,16 +37,16 @@ sections:
keywords:
- name: nuclear_prec
type: float
- default: "user['world_prec']"
+ default: user['world_prec']
predicates:
- - '1.0e-10 < value < 1.0'
+ - 1.0e-10 < value < 1.0
docstring: |
Precision parameter used in smoothing and projection of nuclear potential.
- name: poisson_prec
type: float
- default: "user['world_prec']"
+ default: user['world_prec']
predicates:
- - '1.0e-10 < value < 1.0'
+ - 1.0e-10 < value < 1.0
docstring: |
Precision parameter used in construction of Poisson operators.
- name: exchange_prec
@@ -80,7 +80,7 @@ sections:
type: int
default: 6
predicates:
- - '0 < value < 10'
+ - 0 < value < 10
docstring: |
Number of digits in property output (energies will get twice this
number of digits).
@@ -88,9 +88,14 @@ sections:
type: int
default: 75
predicates:
- - '50 < value < 100'
+ - 50 < value < 100
docstring: |
Line width of printed output (in number of characters).
+ - name: print_constants
+ type: bool
+ default: false
+ docstring: |
+ Print table of physical constants used by MRChem.
- name: Plotter
docstring: |
Give details regarding the density and orbital plots. Three types of
@@ -102,54 +107,54 @@ sections:
keywords:
- name: path
type: str
- default: "plots"
+ default: plots
predicates:
- - "value[-1] != '/'"
+ - value[-1] != '/'
docstring: |
File path to plot directory.
- name: type
type: str
- default: "cube"
+ default: cube
predicates:
- - "value.lower() in ['line', 'surf', 'cube']"
+ - value.lower() in ['line', 'surf', 'cube']
docstring: |
Type of plot: line (1D), surface (2D) or cube (3D).
- name: points
type: List[int]
default: [20, 20, 20]
predicates:
- - "all(p > 0 for p in value)"
- - "not (user['Plotter']['type'] == 'line' and len(value) < 1)"
- - "not (user['Plotter']['type'] == 'surf' and len(value) < 2)"
- - "not (user['Plotter']['type'] == 'cube' and len(value) < 3)"
+ - all(p > 0 for p in value)
+ - not (user['Plotter']['type'] == 'line' and len(value) < 1)
+ - not (user['Plotter']['type'] == 'surf' and len(value) < 2)
+ - not (user['Plotter']['type'] == 'cube' and len(value) < 3)
docstring: |
Number of points in each direction on the cube grid.
- name: O
type: List[float]
default: [0.0, 0.0, 0.0]
predicates:
- - 'len(value) == 3'
+ - len(value) == 3
docstring: |
Origin of plotting ranges.
- name: A
type: List[float]
default: [1.0, 0.0, 0.0]
predicates:
- - 'len(value) == 3'
+ - len(value) == 3
docstring: |
First boundary vector for plot.
- name: B
type: List[float]
default: [0.0, 1.0, 0.0]
predicates:
- - 'len(value) == 3'
+ - len(value) == 3
docstring: |
Second boundary vector for plot.
- name: C
type: List[float]
default: [0.0, 0.0, 1.0]
predicates:
- - 'len(value) == 3'
+ - len(value) == 3
docstring: |
Third boundary vector for plot.
- name: MPI
@@ -199,9 +204,9 @@ sections:
be set automatically based on the world precision.
- name: type
type: str
- default: 'interpolating'
+ default: interpolating
predicates:
- - "value.lower() in ['interpolating', 'legendre']"
+ - value.lower() in ['interpolating', 'legendre']
docstring: |
Polynomial type of multiwavelet basis.
- name: Derivatives
@@ -210,22 +215,22 @@ sections:
keywords:
- name: kinetic
type: str
- default: 'abgv_55'
+ default: abgv_55
docstring: |
Derivative used in kinetic operator.
- name: h_b_dip
type: str
- default: 'abgv_00'
+ default: abgv_00
docstring: |
Derivative used in magnetic dipole operator.
- name: h_m_pso
type: str
- default: 'abgv_00'
+ default: abgv_00
docstring: |
Derivative used in paramagnetic spin-orbit operator.
- name: zora
type: str
- default: 'abgv_00'
+ default: abgv_00
docstring: |
Derivative used ZORA potential.
- name: Molecule
@@ -240,6 +245,8 @@ sections:
- name: multiplicity
type: int
default: 1
+ predicates:
+ - value > 0
docstring: |
Spin multiplicity of molecule.
- name: translate
@@ -260,29 +267,9 @@ sections:
- name: method
type: str
predicates:
- - "value.lower() in
- ['core',
- 'hartree',
- 'hf',
- 'hartreefock',
- 'hartree-fock',
- 'dft',
- 'lda',
- 'svwn3',
- 'svwn5',
- 'pbe',
- 'pbe0',
- 'bpw91',
- 'bp86',
- 'b3p86',
- 'b3p86-g',
- 'blyp',
- 'b3lyp',
- 'b3lyp-g',
- 'olyp',
- 'kt1',
- 'kt2',
- 'kt3']"
+ - value.lower() in ['core', 'hartree', 'hf', 'hartreefock', 'hartree-fock',
+ 'dft', 'lda', 'svwn3', 'svwn5', 'pbe', 'pbe0', 'bpw91', 'bp86', 'b3p86',
+ 'b3p86-g', 'blyp', 'b3lyp', 'b3lyp-g', 'olyp', 'kt1', 'kt2', 'kt3']
docstring: |
Wavefunction method. See predicates for valid methods. ``hf``,
``hartreefock`` and ``hartree-fock`` all mean the same thing, while ``lda``
@@ -296,12 +283,9 @@ sections:
Use spin restricted wavefunction.
- name: relativity
type: str
- default: "none"
+ default: none
predicates:
- - "value.lower() in
- ['none',
- 'zora',
- 'nzora']"
+ - value.lower() in ['none', 'zora', 'nzora']
docstring: |
Set method for relativistic treatment. ``ZORA`` for fully self-consistent ZORA potential, by default
including all potentials (``V_nuc``, ``J``, ``V_xc``) but this can be overwritten in the ``ZORA`` section.
@@ -311,11 +295,6 @@ sections:
docstring: |
Define required parameters for the ZORA Hamiltonian.
keywords:
- - name: light_speed
- type: float
- default: -1.0
- docstring: |
- Adjust speed of light.
- name: include_nuclear
type: bool
default: true
@@ -337,7 +316,7 @@ sections:
keywords:
- name: spin
type: bool
- default: "not(user['WaveFunction']['restricted'])"
+ default: not(user['WaveFunction']['restricted'])
docstring: |
Use spin separated density functionals.
- name: density_cutoff
@@ -407,7 +386,7 @@ sections:
type: List[float]
default: []
predicates:
- - 'len(value) == 0 or len(value) == 3'
+ - len(value) == 0 or len(value) == 3
docstring: |
Strength of external electric field.
- name: Polarizability
@@ -441,102 +420,102 @@ sections:
keywords:
- name: guess_basis
type: str
- default: "initial_guess/mrchem.bas"
+ default: initial_guess/mrchem.bas
docstring: |
File name for GTO basis set, used with ``gto`` guess.
- name: guess_gto_p
type: str
- default: "initial_guess/mrchem.mop"
+ default: initial_guess/mrchem.mop
docstring: |
File name for paired orbitals, used with ``gto`` guess.
- name: guess_gto_a
type: str
- default: "initial_guess/mrchem.moa"
+ default: initial_guess/mrchem.moa
docstring: |
File name for alpha orbitals, used with ``gto`` guess.
- name: guess_gto_b
type: str
- default: "initial_guess/mrchem.mob"
+ default: initial_guess/mrchem.mob
docstring: |
File name for beta orbitals, used with ``gto`` guess.
- name: guess_phi_p
type: str
- default: "initial_guess/phi_p"
+ default: initial_guess/phi_p
docstring: |
File name for paired orbitals, used with ``mw`` guess.
Expected path is ``/phi_p_scf_idx_<0...Np>_.mw
- name: guess_phi_a
type: str
- default: "initial_guess/phi_a"
+ default: initial_guess/phi_a
docstring: |
File name for alpha orbitals, used with ``mw`` guess.
Expected path is ``/phi_a_scf_idx_<0...Na>_.mw
- name: guess_phi_b
type: str
- default: "initial_guess/phi_b"
+ default: initial_guess/phi_b
docstring: |
File name for beta orbitals, used with ``mw`` guess.
Expected path is ``/phi_b_scf_idx_<0...Nb>_.mw
- name: guess_x_p
type: str
- default: "initial_guess/X_p"
+ default: initial_guess/X_p
docstring: |
File name for paired response orbitals, used with ``mw`` guess.
Expected path is ``/x_p_rsp_idx_<0...Np>_.mw
- name: guess_x_a
type: str
- default: "initial_guess/X_a"
+ default: initial_guess/X_a
docstring: |
File name for alpha response orbitals, used with ``mw`` guess.
Expected path is ``/x_a_rsp_idx_<0...Na>_.mw
- name: guess_x_b
type: str
- default: "initial_guess/X_b"
+ default: initial_guess/X_b
docstring: |
File name for beta response orbitals, used with ``mw`` guess.
Expected path is ``/x_b_rsp_idx_<0...Nb>_.mw
- name: guess_y_p
type: str
- default: "initial_guess/Y_p"
+ default: initial_guess/Y_p
docstring: |
File name for paired response orbitals, used with ``mw`` guess.
Expected path is ``/y_p_rsp_idx_<0...Np>_.mw
- name: guess_y_a
type: str
- default: "initial_guess/Y_a"
+ default: initial_guess/Y_a
docstring: |
File name for alpha response orbitals, used with ``mw`` guess.
Expected path is ``/y_a_rsp_idx_<0...Na>_.mw
- name: guess_y_b
type: str
- default: "initial_guess/Y_b"
+ default: initial_guess/Y_b
docstring: |
File name for beta response orbitals, used with ``mw`` guess.
Expected path is ``/y_b_rsp_idx_<0...Nb>_.mw
- name: guess_cube_p
type: str
- default: "initial_guess/phi_p"
+ default: initial_guess/phi_p
docstring: |
File name for paired orbitals, used with ``cube`` guess.
Expected path is ``/phi_p_scf_idx_<0...Np>_.cube
- name: guess_cube_a
type: str
- default: "initial_guess/phi_a"
+ default: initial_guess/phi_a
docstring: |
File name for alpha orbitals, used with ``cube`` guess.
Expected path is ``/phi_a>_scf_idx_<0...Na>_.cube
- name: guess_cube_b
type: str
- default: "initial_guess/phi_b"
+ default: initial_guess/phi_b
docstring: |
File name for beta orbitals, used with ``cube`` guess.
Expected path is ``/phi_b_scf_idx_<0...Nb>_.cube
- name: cube_vectors
type: str
- default: "cube_vectors/"
+ default: cube_vectors/
docstring: |
Directory where cube vectors are stored for mrchem calculation.
-
+
- name: SCF
docstring: |
Includes parameters related to the ground state SCF orbital optimization.
@@ -568,7 +547,7 @@ sections:
Use canonical or localized orbitals.
- name: orbital_thrs
type: float
- default: "10 * user['world_prec']"
+ default: 10 * user['world_prec']
docstring: |
Convergence threshold for orbital residuals.
- name: energy_thrs
@@ -580,7 +559,7 @@ sections:
type: float
default: 1.0e-3
predicates:
- - '1.0e-10 < value < 1.0'
+ - 1.0e-10 < value < 1.0
docstring: |
Precision parameter used in construction of initial guess.
- name: guess_screen
@@ -604,22 +583,10 @@ sections:
Incremental precision in SCF iterations, final value.
- name: guess_type
type: str
- default: 'sad_dz'
+ default: sad_dz
predicates:
- - "value.lower() in
- ['mw',
- 'chk',
- 'gto',
- 'core_sz',
- 'core_dz',
- 'core_tz',
- 'core_qz',
- 'sad_sz',
- 'sad_dz',
- 'sad_tz',
- 'sad_qz',
- 'sad_gto',
- 'cube']"
+ - value.lower() in ['mw', 'chk', 'gto', 'core_sz', 'core_dz', 'core_tz',
+ 'core_qz', 'sad_sz', 'sad_dz', 'sad_tz', 'sad_qz', 'sad_gto', 'cube']
docstring: |
Type of initial guess for ground state orbitals.
``chk`` restarts a previous calculation which was dumped using the
@@ -645,9 +612,9 @@ sections:
there are slashes in the path "path/to/checkpoint".
- name: path_checkpoint
type: str
- default: "checkpoint"
+ default: checkpoint
predicates:
- - "value[-1] != '/'"
+ - value[-1] != '/'
docstring: |
Path to checkpoint files during SCF, used with ``write_checkpoint``
and ``chk`` guess.
@@ -660,9 +627,9 @@ sections:
Can be used as ``mw`` initial guess in subsequent calculations.
- name: path_orbitals
type: str
- default: "orbitals"
+ default: orbitals
predicates:
- - "value[-1] != '/'"
+ - value[-1] != '/'
docstring: |
Path to where converged orbitals will be written in connection with
the ``write_orbitals`` keyword. Note: must be given in quotes if
@@ -688,12 +655,12 @@ sections:
Length of KAIN iterative history.
- name: localize
type: bool
- default: "user['SCF']['localize']"
+ default: user['SCF']['localize']
docstring: |
Use canonical or localized unperturbed orbitals.
- name: orbital_thrs
type: float
- default: "10 * user['world_prec']"
+ default: 10 * user['world_prec']
docstring: |
Convergence threshold for orbital residuals.
- name: property_thrs
@@ -720,17 +687,14 @@ sections:
type: float
default: 1.0e-3
predicates:
- - '1.0e-10 < value < 1.0'
+ - 1.0e-10 < value < 1.0
docstring: |
Precision parameter used in construction of initial guess.
- name: guess_type
type: str
- default: 'none'
+ default: none
predicates:
- - "value.lower() in
- ['none',
- 'chk',
- 'mw']"
+ - value.lower() in ['none', 'chk', 'mw']
docstring: |
Type of initial guess for response.
``none`` will start from a zero guess for the response functions.
@@ -748,9 +712,9 @@ sections:
initial guess in subsequent calculations.
- name: path_checkpoint
type: str
- default: "checkpoint"
+ default: checkpoint
predicates:
- - "value[-1] != '/'"
+ - value[-1] != '/'
docstring: |
Path to checkpoint files during SCF, used with ``write_checkpoint``
and ``chk`` guess.
@@ -763,9 +727,9 @@ sections:
Can be used as ``mw`` initial guess in subsequent calculations.
- name: path_orbitals
type: str
- default: "orbitals"
+ default: orbitals
predicates:
- - "value[-1] != '/'"
+ - value[-1] != '/'
docstring: |
Path to where converged orbitals will be written in connection with
the ``write_orbitals`` keyword.
@@ -787,9 +751,9 @@ sections:
interaction between environment and molecule.
- name: algorithm
type: str
- default: "scrf"
+ default: scrf
predicates:
- - "value.lower() in ['scrf']"
+ - value.lower() in ['scrf']
docstring: |
What algorithm to use for the reaction field
``scrf`` runs a nested algorithm where the generalized Poisson
@@ -797,9 +761,9 @@ sections:
convergence threshold.
- name: convergence_criterion
type: str
- default: "dynamic"
+ default: dynamic
predicates:
- - "value.lower() in ['dynamic', 'static']"
+ - value.lower() in ['dynamic', 'static']
docstring: |
Adjust the convergence threshold for the nested procedure.
``dynamic`` Uses the absolute value of the latest orbital update as
@@ -815,15 +779,15 @@ sections:
charge distribution in the convergence acceleration.
- name: kain
type: int
- default: "user['SCF']['kain']"
+ default: user['SCF']['kain']
docstring: |
Number of previous reaction field iterates kept for convergence
acceleration during the nested precedure.
- name: density_type
type: str
- default: "total"
+ default: total
predicates:
- - "value.lower() in ['total', 'nuclear', 'electronic']"
+ - value.lower() in ['total', 'nuclear', 'electronic']
docstring: |
What part of the total molecular charge density to use in the algorithm.
``total`` uses the total charge density.
@@ -836,7 +800,7 @@ sections:
keywords:
- name: spheres
type: str
- default: ""
+ default: ''
docstring: |
Coordinates and radii of the spheres written as
$spheres
@@ -868,10 +832,65 @@ sections:
solvent used.
- name: formulation
type: str
- default: "exponential"
+ default: exponential
predicates:
- - "value.lower() in ['exponential']"
+ - value.lower() in ['exponential']
docstring: |
Formulation of the Permittivity function. Currently only the
exponential is used.
-
+
+ - name: Constants
+ docstring: Physical and mathematical constants used by MRChem
+ keywords:
+ - name: hartree2simagnetizability
+ default: 78.9451185
+ type: float
+ docstring: '| Conversion factor for magnetizability from atomic units to SI
+ units (unit: J T^-2). Affected code: Printed value of the magnetizability
+ property.'
+ - name: light_speed
+ default: 137.035999084
+ type: float
+ docstring: '| Speed of light in atomic units (unit: au). Affected code: Relativistic
+ Hamiltonians (ZORA, etc.)'
+ - name: angstrom2bohrs
+ default: 1.8897261246257702
+ type: float
+ docstring: '| Conversion factor for Cartesian coordinates from Angstrom to
+ Bohr (unit: Å^-1). Affected code: Parsing of input coordinates, printed
+ coordinates'
+ - name: hartree2kjmol
+ default: 2625.4996394798254
+ type: float
+ docstring: '| Conversion factor from Hartree to kJ/mol (unit: kJ mol^-1).
+ Affected code: Printed value of energies.'
+ - name: hartree2kcalmol
+ default: 627.5094740630558
+ type: float
+ docstring: '| Conversion factor from Hartree to kcal/mol (unit: kcal mol^-1).
+ Affected code: Printed value of energies.'
+ - name: hartree2ev
+ default: 27.211386245988
+ type: float
+ docstring: '| Conversion factor from Hartree to eV (unit: ev). Affected code:
+ Printed value of energies.'
+ - name: hartree2wavenumbers
+ default: 219474.6313632
+ type: float
+ docstring: '| Conversion factor from Hartree to wavenumbers (unit: cm^-1).
+ Affected code: Printed value of frequencies.'
+ - name: fine_structure_constant
+ default: 0.0072973525693
+ type: float
+ docstring: '| Fine-structure constant in atomic units (unit: au). Affected
+ code: Certain magnetic interaction operators.'
+ - name: electron_g_factor
+ default: -2.00231930436256
+ type: float
+ docstring: '| Electron g factor in atomic units (unit: au). Affected code:
+ Certain magnetic interaction operators.'
+ - name: dipmom_au2debye
+ default: 2.5417464739297717
+ type: float
+ docstring: '| Conversion factor for dipoles from atomic units to Debye (unit:
+ ?). Affected code: Printed value of dipole moments.'
diff --git a/src/analyticfunctions/HarmonicOscillatorFunction.h b/src/analyticfunctions/HarmonicOscillatorFunction.h
index a139a5c4d..b31dee6e0 100644
--- a/src/analyticfunctions/HarmonicOscillatorFunction.h
+++ b/src/analyticfunctions/HarmonicOscillatorFunction.h
@@ -28,6 +28,7 @@
#include
#include "MRCPP/MWFunctions"
+#include "chemistry/PhysicalConstants.h"
#include "mrchem.h"
@@ -42,7 +43,7 @@ class HarmonicOscillator1D final {
double operator()(double x) const {
double ax = this->alpha * (x - this->origin);
- double ap = std::sqrt(this->alpha / MATHCONST::sqrt_pi);
+ double ap = std::sqrt(this->alpha / std::sqrt(mrcpp::pi));
double N_nu = std::sqrt(N2(this->nu));
return ap * N_nu * H(this->nu, ax) * exp(-ax * ax / 2.0);
}
diff --git a/src/analyticfunctions/HydrogenFunction.cpp b/src/analyticfunctions/HydrogenFunction.cpp
index a7eaafe2a..464b1d33b 100644
--- a/src/analyticfunctions/HydrogenFunction.cpp
+++ b/src/analyticfunctions/HydrogenFunction.cpp
@@ -26,6 +26,7 @@
#include "MRCPP/Printer"
#include "HydrogenFunction.h"
+#include "chemistry/PhysicalConstants.h"
#include "mrchem.h"
#include "utils/math_utils.h"
@@ -113,6 +114,7 @@ double AngularFunction::evalf(const mrcpp::Coord<3> &r) const {
// clang-format off
double AngularFunction::calcConstant() const {
double c = 0.0;
+ double sqrt_pi = std::sqrt(mrcpp::pi);
if (L == 0 and M == 0) { c = std::sqrt( 1.0/1.0);
} else if (L == 1 and M == 0) { c = std::sqrt( 3.0/1.0);
} else if (L == 1 and M == 1) { c = std::sqrt( 3.0/1.0);
@@ -130,7 +132,7 @@ double AngularFunction::calcConstant() const {
} else if (L == 3 and M == 5) { c = std::sqrt(105.0/4.0);
} else if (L == 3 and M == 6) { c = std::sqrt(105.0/4.0);
} else { NOT_IMPLEMENTED_ABORT; }
- return c/(2.0*MATHCONST::sqrt_pi);
+ return c/(2.0*sqrt_pi);
}
// clang-format on
diff --git a/src/chemistry/CMakeLists.txt b/src/chemistry/CMakeLists.txt
index fc01af372..0dcd10265 100644
--- a/src/chemistry/CMakeLists.txt
+++ b/src/chemistry/CMakeLists.txt
@@ -2,4 +2,5 @@ target_sources(mrchem PRIVATE
${CMAKE_CURRENT_SOURCE_DIR}/chemistry_utils.cpp
${CMAKE_CURRENT_SOURCE_DIR}/PeriodicTable.cpp
${CMAKE_CURRENT_SOURCE_DIR}/Molecule.cpp
+ ${CMAKE_CURRENT_SOURCE_DIR}/PhysicalConstants.cpp
)
diff --git a/src/chemistry/PhysicalConstants.cpp b/src/chemistry/PhysicalConstants.cpp
new file mode 100644
index 000000000..91ef82c7e
--- /dev/null
+++ b/src/chemistry/PhysicalConstants.cpp
@@ -0,0 +1,66 @@
+/*
+ * MRChem, a numerical real-space code for molecular electronic structure
+ * calculations within the self-consistent field (SCF) approximations of quantum
+ * chemistry (Hartree-Fock and Density Functional Theory).
+ * Copyright (C) 2022 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors.
+ *
+ * This file is part of MRChem.
+ *
+ * MRChem is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * MRChem is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with MRChem. If not, see .
+ *
+ * For information on the complete list of contributors to MRChem, see:
+ *
+ */
+
+#include "PhysicalConstants.h"
+#include "MRCPP/Printer"
+#include "utils/print_utils.h"
+#include
+
+using json = nlohmann::json;
+
+namespace mrchem {
+
+bool PhysicalConstants::initialized = false;
+json PhysicalConstants::constants_ = json();
+
+// clang-format off
+json PhysicalConstants::testConstants = {
+ {"angstrom2bohrs", 1.8897261246257702},
+ {"dipmom_au2debye", 2.5417464739297717},
+ {"electron_g_factor", -2.00231930436256},
+ {"fine_structure_constant", 0.0072973525693},
+ {"hartree2ev", 27.211386245988},
+ {"hartree2kcalmol", 627.5094740630558},
+ {"hartree2kjmol", 2625.4996394798254},
+ {"hartree2simagnetizability", 78.9451185},
+ {"hartree2wavenumbers", 219474.6313632},
+ {"light_speed", 137.035999084}
+};
+// clang-format on
+
+PhysicalConstants &PhysicalConstants::Initialize(const json &constants) {
+ initialized = true;
+ static PhysicalConstants obj(constants);
+ return obj;
+}
+
+/** @brief Pretty print physical constants */
+void PhysicalConstants::Print() {
+ mrcpp::print::header(0, "Physical Constants");
+ print_utils::json(0, constants_, true);
+ mrcpp::print::separator(0, '=', 2);
+}
+
+} // namespace mrchem
diff --git a/src/chemistry/PhysicalConstants.h b/src/chemistry/PhysicalConstants.h
new file mode 100644
index 000000000..430e4902f
--- /dev/null
+++ b/src/chemistry/PhysicalConstants.h
@@ -0,0 +1,65 @@
+/*
+ * MRChem, a numerical real-space code for molecular electronic structure
+ * calculations within the self-consistent field (SCF) approximations of quantum
+ * chemistry (Hartree-Fock and Density Functional Theory).
+ * Copyright (C) 2022 Stig Rune Jensen, Luca Frediani, Peter Wind and contributors.
+ *
+ * This file is part of MRChem.
+ *
+ * MRChem is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * MRChem is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with MRChem. If not, see .
+ *
+ * For information on the complete list of contributors to MRChem, see:
+ *
+ */
+
+#pragma once
+#include "MRCPP/Printer"
+#include
+
+using json = nlohmann::json;
+
+namespace mrchem {
+
+class PhysicalConstants {
+public:
+ static PhysicalConstants &Initialize(const json &constants);
+ static double get(const std::string &key) {
+ try {
+ if (initialized) {
+ return constants_[key];
+ } else {
+ return testConstants[key];
+ }
+ } catch (...) { MSG_ABORT("Error getting constant with name: " + key); }
+ }
+
+ static void Print();
+
+ PhysicalConstants() = default;
+ ~PhysicalConstants() = default;
+
+ PhysicalConstants(const PhysicalConstants &) = delete;
+ PhysicalConstants &operator=(const PhysicalConstants &) = delete;
+ PhysicalConstants &operator=(const PhysicalConstants &&) = delete;
+ PhysicalConstants(PhysicalConstants &&) = delete;
+
+ static bool initialized;
+
+private:
+ PhysicalConstants(const json &constants) { constants_ = constants; }
+ static json constants_;
+ static json testConstants;
+};
+
+} // namespace mrchem
diff --git a/src/chemistry/chemistry_utils.cpp b/src/chemistry/chemistry_utils.cpp
index 22b985088..70913ad41 100644
--- a/src/chemistry/chemistry_utils.cpp
+++ b/src/chemistry/chemistry_utils.cpp
@@ -24,10 +24,12 @@
*/
#include "chemistry_utils.h"
#include "Nucleus.h"
+#include "PhysicalConstants.h"
#include "qmfunctions/Density.h"
#include "qmfunctions/density_utils.h"
#include "utils/math_utils.h"
#include
+
namespace mrchem {
/** @brief computes the repulsion self energy of a set of nuclei
@@ -61,7 +63,7 @@ double chemistry::get_total_charge(const Nuclei &nucs) {
}
Density chemistry::compute_nuclear_density(double prec, const Nuclei &nucs, double alpha) {
- auto beta = std::pow(alpha / MATHCONST::pi, 3.0 / 2.0);
+ auto beta = std::pow(alpha / mrcpp::pi, 3.0 / 2.0);
int nNucs = nucs.size();
auto gauss = mrcpp::GaussExp<3>();
diff --git a/src/driver.cpp b/src/driver.cpp
index 65b166f8f..d766b677b 100644
--- a/src/driver.cpp
+++ b/src/driver.cpp
@@ -31,6 +31,7 @@
#include "chemistry/Molecule.h"
#include "chemistry/Nucleus.h"
+#include "chemistry/PhysicalConstants.h"
#include "initial_guess/chk.h"
#include "initial_guess/core.h"
@@ -140,15 +141,19 @@ void driver::init_molecule(const json &json_mol, Molecule &mol) {
auto xyz = coord["xyz"];
nuclei.push_back(atom, xyz);
}
- std::vector radii;
- std::vector> spheres;
- for (const auto &coord : json_mol["cavity_coords"].get()) {
- radii.push_back(coord["radius"].get());
- spheres.push_back(coord["center"].get>());
- }
- auto cavity_width = json_mol["cavity_width"].get();
- mol.initCavity(spheres, radii, cavity_width);
+ if (json_mol.contains("cavity")) {
+ auto json_cavity = json_mol["cavity"];
+ std::vector radii;
+ std::vector> coords;
+ for (const auto &sphere : json_cavity["spheres"]) {
+ radii.push_back(sphere["radius"]);
+ coords.push_back(sphere["center"]);
+ }
+ auto width = json_cavity["width"];
+
+ mol.initCavity(coords, radii, width);
+ }
mol.printGeometry();
}
@@ -981,8 +986,7 @@ void driver::build_fock_operator(const json &json_fock, Molecule &mol, FockBuild
////////////////////// Zora Operator //////////////////
///////////////////////////////////////////////////////////
if (json_fock.contains("zora_operator")) {
- auto c = json_fock["zora_operator"]["light_speed"];
- if (c <= 0.0) c = PHYSCONST::alpha_inv;
+ auto c = PhysicalConstants::get("light_speed");
F.setLightSpeed(c);
auto include_nuclear = json_fock["zora_operator"]["include_nuclear"];
diff --git a/src/environment/Cavity.cpp b/src/environment/Cavity.cpp
index 1af51fced..1bd5b5225 100644
--- a/src/environment/Cavity.cpp
+++ b/src/environment/Cavity.cpp
@@ -23,6 +23,7 @@
*
*/
#include "Cavity.h"
+#include "chemistry/PhysicalConstants.h"
#include "utils/math_utils.h"
namespace mrchem {
@@ -60,6 +61,7 @@ Cavity::Cavity(std::vector> ¢ers, std::vector &radii
auto gradCavity(const mrcpp::Coord<3> &r, int index, const std::vector> ¢ers, std::vector &radii, double width) -> double {
double C = 1.0;
double DC = 0.0;
+ double sqrt_pi = std::sqrt(mrcpp::pi);
for (int i = 0; i < centers.size(); i++) {
double s = math_utils::calc_distance(centers[i], r) - radii[i];
double ds = (r[index] - centers[i][index]) / (math_utils::calc_distance(centers[i], r));
@@ -67,7 +69,7 @@ auto gradCavity(const mrcpp::Coord<3> &r, int index, const std::vector
+#include "chemistry/PhysicalConstants.h"
#include "chemistry/chemistry_utils.h"
#include "qmfunctions/density_utils.h"
#include "qmfunctions/qmfunction_utils.h"
@@ -42,6 +43,7 @@ using DerivativeOperator_p = std::shared_ptr>;
using OrbitalVector_p = std::shared_ptr;
namespace mrchem {
+
SCRF::SCRF(Permittivity e,
const Nuclei &N,
PoissonOperator_p P,
@@ -113,7 +115,7 @@ void SCRF::computeGamma(QMFunction &potential, QMFunction &out_gamma) {
auto d_V = mrcpp::gradient(*derivative, potential.real());
resetQMFunction(out_gamma);
mrcpp::dot(this->apply_prec, out_gamma.real(), d_V, this->d_cavity);
- out_gamma.rescale(std::log((epsilon.getEpsIn() / epsilon.getEpsOut())) * (1.0 / (4.0 * MATHCONST::pi)));
+ out_gamma.rescale(std::log((epsilon.getEpsIn() / epsilon.getEpsOut())) * (1.0 / (4.0 * mrcpp::pi)));
mrcpp::clear(d_V, true);
}
diff --git a/src/mrchem.cpp b/src/mrchem.cpp
index b7adf662c..6a8698405 100644
--- a/src/mrchem.cpp
+++ b/src/mrchem.cpp
@@ -35,6 +35,7 @@
#include "version.h"
#include "chemistry/Molecule.h"
+#include "chemistry/PhysicalConstants.h"
// Initializing global variables
mrcpp::MultiResolutionAnalysis<3> *mrchem::MRA;
@@ -55,6 +56,11 @@ int main(int argc, char **argv) {
const auto &mol_inp = json_inp["molecule"];
const auto &scf_inp = json_inp["scf_calculation"];
const auto &rsp_inp = json_inp["rsp_calculations"];
+ const auto &con_inp = json_inp["constants"];
+
+ // Instantiate the physical constants singleton
+ PhysicalConstants::Initialize(con_inp);
+ if (json_inp["printer"]["print_constants"]) PhysicalConstants::Print();
Timer timer;
Molecule mol;
diff --git a/src/properties/DipoleMoment.h b/src/properties/DipoleMoment.h
index e099c8518..c30ad0f2f 100644
--- a/src/properties/DipoleMoment.h
+++ b/src/properties/DipoleMoment.h
@@ -29,6 +29,8 @@
#include "mrchem.h"
+#include "chemistry/PhysicalConstants.h"
+
#include "utils/math_utils.h"
#include "utils/print_utils.h"
@@ -53,9 +55,9 @@ class DipoleMoment final {
auto nuc_au = getNuclear().norm();
auto tot_au = getTensor().norm();
- auto el_db = el_au * PHYSCONST::Debye;
- auto nuc_db = nuc_au * PHYSCONST::Debye;
- auto tot_db = tot_au * PHYSCONST::Debye;
+ auto el_db = el_au * PhysicalConstants::get("dipmom_au2debye");
+ auto nuc_db = nuc_au * PhysicalConstants::get("dipmom_au2debye");
+ auto tot_db = tot_au * PhysicalConstants::get("dipmom_au2debye");
mrcpp::print::header(0, "Dipole Moment (" + id + ")");
print_utils::coord(0, "r_O", getOrigin());
diff --git a/src/properties/Magnetizability.h b/src/properties/Magnetizability.h
index 1acc61b94..f9c87d92c 100644
--- a/src/properties/Magnetizability.h
+++ b/src/properties/Magnetizability.h
@@ -29,6 +29,8 @@
#include "mrchem.h"
+#include "chemistry/PhysicalConstants.h"
+
#include "utils/math_utils.h"
#include "utils/print_utils.h"
@@ -53,7 +55,7 @@ class Magnetizability final {
void print(const std::string &id) const {
auto w_au = getFrequency();
- auto w_cm = PHYSCONST::cm_m1 * w_au;
+ auto w_cm = PhysicalConstants::get("hartree2wavenumbers") * w_au;
auto dynamic = (w_au > mrcpp::MachineZero);
auto l_nm = (dynamic) ? (1.0e7 / w_cm) : 0.0;
@@ -62,9 +64,9 @@ class Magnetizability final {
auto iso_au_t = iso_au_d + iso_au_p;
// SI units (J/T^2 10^{-30})
- auto iso_si_t = iso_au_t * PHYSCONST::JT_m2;
- auto iso_si_d = iso_au_d * PHYSCONST::JT_m2;
- auto iso_si_p = iso_au_p * PHYSCONST::JT_m2;
+ auto iso_si_t = iso_au_t * PhysicalConstants::get("hartree2simagnetizability");
+ auto iso_si_d = iso_au_d * PhysicalConstants::get("hartree2simagnetizability");
+ auto iso_si_p = iso_au_p * PhysicalConstants::get("hartree2simagnetizability");
mrcpp::print::header(0, "Magnetizability (" + id + ")");
if (dynamic) print_utils::scalar(0, "Wavelength", l_nm, "(nm)");
diff --git a/src/properties/Polarizability.h b/src/properties/Polarizability.h
index 10a6877bf..85ce66608 100644
--- a/src/properties/Polarizability.h
+++ b/src/properties/Polarizability.h
@@ -27,6 +27,7 @@
#include
+#include "chemistry/PhysicalConstants.h"
#include "mrchem.h"
#include "utils/math_utils.h"
@@ -50,7 +51,7 @@ class Polarizability final {
void print(const std::string &id) const {
auto w_au = getFrequency();
- auto w_cm = PHYSCONST::cm_m1 * w_au;
+ auto w_cm = PhysicalConstants::get("hartree2wavenumbers") * w_au;
auto dynamic = (w_au > mrcpp::MachineZero);
auto l_nm = (dynamic) ? (1.0e7 / w_cm) : 0.0;
auto iso_au = getTensor().trace() / 3.0;
diff --git a/src/properties/SCFEnergy.h b/src/properties/SCFEnergy.h
index cd8b41288..c82dd6b78 100644
--- a/src/properties/SCFEnergy.h
+++ b/src/properties/SCFEnergy.h
@@ -27,7 +27,9 @@
#include
+#include "chemistry/PhysicalConstants.h"
#include "mrchem.h"
+
#include "utils/print_utils.h"
/** @class SCFEnergy
@@ -68,9 +70,9 @@ class SCFEnergy final {
void print(const std::string &id) const {
auto E_au = E_nuc + E_el;
- auto E_eV = E_au * PHYSCONST::eV;
- auto E_kJ = E_au * PHYSCONST::kJ;
- auto E_kcal = E_au * PHYSCONST::kcal;
+ auto E_eV = E_au * PhysicalConstants::get("hartree2ev");
+ auto E_kJ = E_au * PhysicalConstants::get("hartree2kjmol");
+ auto E_kcal = E_au * PhysicalConstants::get("hartree2kcalmol");
auto pprec = 2 * mrcpp::Printer::getPrecision();
mrcpp::print::header(0, "Molecular Energy (" + id + ")");
diff --git a/src/qmoperators/one_electron/DeltaOperator.h b/src/qmoperators/one_electron/DeltaOperator.h
index a3bfe2156..066e3e149 100644
--- a/src/qmoperators/one_electron/DeltaOperator.h
+++ b/src/qmoperators/one_electron/DeltaOperator.h
@@ -29,6 +29,8 @@
#include "tensor/RankZeroOperator.h"
+#include "chemistry/PhysicalConstants.h"
+
#include "qmfunctions/qmfunction_utils.h"
#include "qmoperators/QMPotential.h"
@@ -47,7 +49,7 @@ class DeltaOperator final : public RankZeroOperator {
// Define analytic potential
double beta = 1.0 / smooth_prec;
- double alpha = std::pow(beta / MATHCONST::pi, 3.0 / 2.0);
+ double alpha = std::pow(beta / mrcpp::pi, 3.0 / 2.0);
mrcpp::GaussFunc<3> f(beta, alpha, o);
// Project analytic potential
diff --git a/src/qmoperators/one_electron/H_BM_dia.h b/src/qmoperators/one_electron/H_BM_dia.h
index 275a5bc14..a6b80d136 100644
--- a/src/qmoperators/one_electron/H_BM_dia.h
+++ b/src/qmoperators/one_electron/H_BM_dia.h
@@ -25,6 +25,7 @@
#pragma once
+#include "chemistry/PhysicalConstants.h"
#include "tensor/RankTwoOperator.h"
#include "NuclearGradientOperator.h"
@@ -53,7 +54,7 @@ class H_BM_dia final : public RankTwoOperator<3, 3> {
: H_BM_dia(PositionOperator(o), NuclearGradientOperator(1.0, k, proj_prec, smooth_prec)) {}
H_BM_dia(PositionOperator r_o, NuclearGradientOperator r_rm3) {
- const double alpha_2 = PHYSCONST::alpha * PHYSCONST::alpha;
+ const double alpha_2 = PhysicalConstants::get("fine_structure_constant") * PhysicalConstants::get("fine_structure_constant") * 1000000.0;
RankZeroOperator &o_x = r_o[0];
RankZeroOperator &o_y = r_o[1];
diff --git a/src/qmoperators/one_electron/H_B_spin.h b/src/qmoperators/one_electron/H_B_spin.h
index d43a29caa..bbdc1b2df 100644
--- a/src/qmoperators/one_electron/H_B_spin.h
+++ b/src/qmoperators/one_electron/H_B_spin.h
@@ -25,6 +25,7 @@
#pragma once
+#include "chemistry/PhysicalConstants.h"
#include "tensor/RankOneOperator.h"
#include "SpinOperator.h"
@@ -51,13 +52,13 @@ class H_B_spin final : public RankOneOperator<3> {
: H_B_spin(SpinOperator()) {}
explicit H_B_spin(SpinOperator s) {
- const double g_e = PHYSCONST::g_e;
+ const double g_e = PhysicalConstants::get("electron_g_factor");
// Invoke operator= to assign *this operator
RankOneOperator<3> &h = (*this);
- h[0] = (-g_e / 2.0) * s[0];
- h[1] = (-g_e / 2.0) * s[1];
- h[2] = (-g_e / 2.0) * s[2];
+ h[0] = (g_e / 2.0) * s[0];
+ h[1] = (g_e / 2.0) * s[1];
+ h[2] = (g_e / 2.0) * s[2];
h[0].name() = "h_B_spin[x]";
h[1].name() = "h_B_spin[y]";
h[2].name() = "h_B_spin[z]";
diff --git a/src/qmoperators/one_electron/H_MB_dia.h b/src/qmoperators/one_electron/H_MB_dia.h
index 342915af5..dbe2a52ba 100644
--- a/src/qmoperators/one_electron/H_MB_dia.h
+++ b/src/qmoperators/one_electron/H_MB_dia.h
@@ -25,6 +25,7 @@
#pragma once
+#include "chemistry/PhysicalConstants.h"
#include "tensor/RankTwoOperator.h"
#include "NuclearGradientOperator.h"
@@ -53,7 +54,7 @@ class H_MB_dia final : public RankTwoOperator<3, 3> {
: H_MB_dia(PositionOperator(o), NuclearGradientOperator(1.0, k, proj_prec, smooth_prec)) {}
H_MB_dia(PositionOperator r_o, NuclearGradientOperator r_rm3) {
- const double alpha_2 = PHYSCONST::alpha * PHYSCONST::alpha;
+ const double alpha_2 = PhysicalConstants::get("fine_structure_constant") * PhysicalConstants::get("fine_structure_constant") * 1000000.0;
RankZeroOperator &o_x = r_o[0];
RankZeroOperator &o_y = r_o[1];
diff --git a/src/qmoperators/one_electron/H_M_fc.h b/src/qmoperators/one_electron/H_M_fc.h
index 83af6eff9..a2a933398 100644
--- a/src/qmoperators/one_electron/H_M_fc.h
+++ b/src/qmoperators/one_electron/H_M_fc.h
@@ -27,6 +27,8 @@
#include "tensor/RankOneOperator.h"
+#include "chemistry/PhysicalConstants.h"
+
#include "DeltaOperator.h"
#include "H_B_spin.h"
@@ -52,8 +54,8 @@ class H_M_fc final : public RankOneOperator<3> {
: H_M_fc(H_B_spin(), DeltaOperator(o, proj_prec, smooth_prec)) {}
H_M_fc(H_B_spin s, DeltaOperator delta) {
- const double coef = -(8.0 / 3.0) * MATHCONST::pi;
- const double alpha_2 = PHYSCONST::alpha * PHYSCONST::alpha;
+ const double coef = -(8.0 / 3.0) * mrcpp::pi;
+ const double alpha_2 = PhysicalConstants::get("fine_structure_constant") * PhysicalConstants::get("fine_structure_constant") * 1000000.0;
// Invoke operator= to assign *this operator
RankOneOperator<3> &h = (*this);
diff --git a/src/qmoperators/one_electron/H_M_pso.h b/src/qmoperators/one_electron/H_M_pso.h
index 0d4be6149..3f17d7c71 100644
--- a/src/qmoperators/one_electron/H_M_pso.h
+++ b/src/qmoperators/one_electron/H_M_pso.h
@@ -25,6 +25,7 @@
#pragma once
+#include "chemistry/PhysicalConstants.h"
#include "tensor/RankOneOperator.h"
#include "MomentumOperator.h"
@@ -52,7 +53,7 @@ class H_M_pso final : public RankOneOperator<3> {
: H_M_pso(MomentumOperator(D), NuclearGradientOperator(1.0, k, proj_prec, smooth_prec)) {}
H_M_pso(MomentumOperator p, NuclearGradientOperator r_rm3) {
- const double alpha_2 = PHYSCONST::alpha * PHYSCONST::alpha;
+ const double alpha_2 = PhysicalConstants::get("fine_structure_constant") * PhysicalConstants::get("fine_structure_constant") * 1000000.0;
RankZeroOperator &p_x = p[0];
RankZeroOperator &p_y = p[1];
diff --git a/src/scf_solver/HelmholtzVector.cpp b/src/scf_solver/HelmholtzVector.cpp
index cec8a815d..6559b3c0d 100644
--- a/src/scf_solver/HelmholtzVector.cpp
+++ b/src/scf_solver/HelmholtzVector.cpp
@@ -30,6 +30,7 @@
#include "parallel.h"
#include "HelmholtzVector.h"
+#include "chemistry/PhysicalConstants.h"
#include "qmfunctions/Orbital.h"
#include "qmfunctions/orbital_utils.h"
#include "qmfunctions/qmfunction_utils.h"
@@ -40,6 +41,7 @@ using mrcpp::Printer;
using mrcpp::Timer;
namespace mrchem {
+
extern mrcpp::MultiResolutionAnalysis<3> *MRA; // Global MRA
/** @brief HelmholtzVector constructor
@@ -148,13 +150,13 @@ Orbital HelmholtzVector::apply(int i, Orbital &phi) const {
if (phi.hasReal()) {
out.alloc(NUMBER::Real);
mrcpp::apply(this->prec, out.real(), H, phi.real(), -1, true); // Absolute prec
- out.real().rescale(-1.0 / (2.0 * MATHCONST::pi));
+ out.real().rescale(-1.0 / (2.0 * mrcpp::pi));
}
if (phi.hasImag()) {
out.alloc(NUMBER::Imag);
mrcpp::apply(this->prec, out.imag(), H, phi.imag(), -1, true); // Absolute prec
double sign = (phi.conjugate()) ? -1.0 : 1.0;
- out.imag().rescale(sign / (2.0 * MATHCONST::pi));
+ out.imag().rescale(sign / (2.0 * mrcpp::pi));
}
return out;
}
diff --git a/src/utils/print_utils.cpp b/src/utils/print_utils.cpp
index aeb03e37b..d2bb985fb 100644
--- a/src/utils/print_utils.cpp
+++ b/src/utils/print_utils.cpp
@@ -102,6 +102,51 @@ void print_utils::text(int level, const std::string &txt, const std::string &val
println(level, o.str());
}
+void print_utils::json(int level, const nlohmann::json &j, bool ralign) {
+ // Determine longest name
+ // This will be used for the spacing left of :
+ // if the name is too long to be aligned with
+ // the other sections
+ int lshift = 0;
+ for (const auto &item : j.items()) {
+ if (item.key().size() > lshift) lshift = item.key().size();
+ }
+
+ // Loop over json items
+ for (const auto &item : j.items()) {
+ // Extract key and value from json
+ std::string key = item.key();
+ std::stringstream o_val;
+ o_val << item.value();
+ std::string val = o_val.str();
+
+ // Remove quotes from val
+ val.erase(std::remove(val.begin(), val.end(), '\"'), val.end());
+
+ // Standard shift to align all colons
+ int w0 = Printer::getWidth() - 2; // Defines the printable width
+ int w1 = w0 * 2 / 9; // Space dedicated to the json key
+ int w2 = w0 - 3 * w1;
+ int w3 = w2 - (val.size() + 1);
+
+ // Some paddings for book keeping
+ int frontEndPadding = 2; // Empty space at beginning and end
+ int colonPadding = 3; // Two empty spaces around a single colon
+
+ // Use standard spacing if longest name fits
+ if (w2 > lshift) lshift = w2 - frontEndPadding;
+
+ // Calculate the shift needed for right-aligning
+ int rshift = (ralign) ? Printer::getWidth() - lshift - val.size() - frontEndPadding - colonPadding : 0;
+
+ // Check that rshift is not negative (caused by very long names)
+ if (rshift < 0) rshift = 0;
+
+ // Print line
+ std::printf(" %-*s%-s%-s%-s \n", lshift, key.c_str(), " : ", std::string(rshift, ' ').c_str(), val.c_str());
+ }
+}
+
void print_utils::coord(int level, const std::string &txt, const mrcpp::Coord<3> &val, int p, bool s) {
if (p < 0) p = Printer::getPrecision();
int w0 = Printer::getWidth() - 2;
diff --git a/src/utils/print_utils.h b/src/utils/print_utils.h
index 60d405320..25127a510 100644
--- a/src/utils/print_utils.h
+++ b/src/utils/print_utils.h
@@ -28,11 +28,14 @@
#include "mrchem.h"
#include "qmfunctions/qmfunction_fwd.h"
+#include
+
namespace mrchem {
namespace print_utils {
void headline(int level, const std::string &txt);
void text(int level, const std::string &txt, const std::string &val);
+void json(int level, const nlohmann::json &j, bool ralign = false);
void coord(int level, const std::string &txt, const mrcpp::Coord<3> &val, int p = -1, bool s = false);
void scalar(int level, const std::string &txt, double val, const std::string &unit = "", int p = -1, bool s = false);
void vector(int level, const std::string &txt, const DoubleVector &val, int p = -1, bool s = false);
diff --git a/tests/he_zora_scf_lda/he.inp b/tests/he_zora_scf_lda/he.inp
index 564cf5ee3..379cf0d3a 100644
--- a/tests/he_zora_scf_lda/he.inp
+++ b/tests/he_zora_scf_lda/he.inp
@@ -4,6 +4,6 @@
"MPI": {"numerically_exact": true},
"Molecule": {"coords": "He 0.0 0.0 0.0"},
"WaveFunction": {"method": "lda", "relativity": "zora"},
- "ZORA": {"light_speed": -1.0, "include_nuclear": true, "include_coulomb": true, "include_xc": true},
+ "ZORA": {"include_nuclear": true, "include_coulomb": true, "include_xc": true},
"SCF": {"kain": 3, "max_iter": 10, "guess_type": "sad_gto"}
}
diff --git a/tests/li_solv/li.inp b/tests/li_solv/li.inp
index 4b46f72c3..e7b37af5f 100644
--- a/tests/li_solv/li.inp
+++ b/tests/li_solv/li.inp
@@ -1,47 +1,37 @@
-# vim:syntax=sh:
-world_prec = 1.0e-3 # Overall relative precision
-world_size = 5 # Size of simulation box 2^n
-
-MPI {
- numerically_exact = true # Guarantee identical results in MPI
-}
-
-Molecule {
- charge = 1
-$coords
-Li 0.0 0.0 0.0
-$end
-}
-
-WaveFunction {
- method = pbe0 # Wave function method (HF or DFT)
-}
-
-Environment {
- run_environment = true
- algorithm = scrf
- kain = 5
- max_iter = 100
- convergence_criterion = dynamic
- extrapolate_Vr = true
-
- Cavity {
- cavity_width = 0.5
-$spheres
-0.0 0.0 0.0 4.0
-$end
- }
-
- Permittivity {
- epsilon_in = 1.0
- epsilon_out = 2.0
- formulation = exponential
+{
+"world_prec": 1.0e-3,
+"world_size": 5,
+"MPI": {
+ "numerically_exact": true
+},
+"Molecule": {
+ "charge": 1,
+ "coords": "Li 0.0 0.0 0.0"
+},
+"WaveFunction": {
+ "method": "pbe0"
+},
+"Environment": {
+ "run_environment": true,
+ "algorithm": "scrf",
+ "kain": 5,
+ "max_iter": 100,
+ "convergence_criterion": "dynamic",
+ "extrapolate_Vr": true,
+ "Cavity": {
+ "cavity_width": 0.5,
+ "spheres": "0.0 0.0 0.0 4.0"
+ },
+ "Permittivity": {
+ "epsilon_in": 1.0,
+ "epsilon_out": 2.0,
+ "formulation": "exponential"
}
+},
+"SCF": {
+ "run": true,
+ "kain": 5,
+ "orbital_thrs": 1.0e-1,
+ "guess_type": "sad_dz"
}
-
-SCF {
- run = true
- kain = 5
- orbital_thrs = 1.0e-1
- guess_type = sad_dz
}
diff --git a/tests/li_solv/reference/li.json b/tests/li_solv/reference/li.json
index f02f5d70b..8b6ff47a9 100644
--- a/tests/li_solv/reference/li.json
+++ b/tests/li_solv/reference/li.json
@@ -2,14 +2,6 @@
"input": {
"molecule": {
"cavity_coords": [
- {
- "center": [
- 0.0,
- 0.0,
- 0.0
- ],
- "radius": 4.0
- },
{
"center": [
0.0,
@@ -121,6 +113,7 @@
"localize": false,
"method": "DFT (PBE0)",
"prec": 0.001,
+ "relativity": "None",
"restricted": true,
"screen": 12.0,
"type": "sad",
@@ -153,6 +146,7 @@
"method": "DFT (PBE0)",
"orbital_thrs": 0.1,
"proj_prec": 0.001,
+ "relativity": "None",
"rotation": 0,
"shared_memory": false,
"smooth_prec": 0.001,
@@ -172,7 +166,7 @@
"charge": 1,
"dipole_moment": {
"dip-1": {
- "magnitude": 7.124597251131804e-14,
+ "magnitude": 2.5407415880983347e-13,
"r_O": [
0.0,
0.0,
@@ -208,7 +202,7 @@
"multiplicity": 1,
"orbital_energies": {
"energy": [
- -2.21888856419452
+ -2.210659505993009
],
"occupation": [
2.0
@@ -216,23 +210,23 @@
"spin": [
"p"
],
- "sum_occupied": -4.43777712838904
+ "sum_occupied": -4.421319011986018
},
"scf_energy": {
- "E_ee": 3.354625825769642,
+ "E_ee": 3.3576788742677035,
"E_eext": 0.0,
- "E_el": -7.135135523826116,
- "E_en": -16.47750955461854,
- "E_kin": 7.571525854087563,
+ "E_el": -7.127946930407876,
+ "E_en": -16.4914760516857,
+ "E_kin": 7.584212118791978,
"E_next": 0.0,
"E_nn": 0.0,
- "E_nuc": -0.18173935348374062,
- "E_tot": -7.316874877309856,
- "E_x": -0.4192040376909349,
- "E_xc": -1.2857400054483297,
- "Er_el": 0.12116639407448325,
- "Er_nuc": -0.18173935348374062,
- "Er_tot": -0.06057295940925737
+ "E_nuc": -0.19212873389394816,
+ "E_tot": -7.320075664301824,
+ "E_x": -0.4195941455371656,
+ "E_xc": -1.2868490616228483,
+ "Er_el": 0.12808133537815772,
+ "Er_nuc": -0.19212873389394816,
+ "Er_tot": -0.06404739851579039
}
},
"provenance": {
@@ -246,17 +240,17 @@
"rsp_calculations": null,
"scf_calculation": {
"initial_energy": {
- "E_ee": 3.332799383894304,
+ "E_ee": 3.3327993838943355,
"E_eext": 0.0,
- "E_el": -7.222649094979309,
- "E_en": -16.530755287695275,
- "E_kin": 7.675564301900197,
+ "E_el": -7.222649094979323,
+ "E_en": -16.530755287695392,
+ "E_kin": 7.675564301900276,
"E_next": 0.0,
"E_nn": 0.0,
"E_nuc": 0.0,
- "E_tot": -7.222649094979309,
- "E_x": -0.4164844953175967,
- "E_xc": -1.2837729977609398,
+ "E_tot": -7.222649094979323,
+ "E_x": -0.4164844953176007,
+ "E_xc": -1.2837729977609427,
"Er_el": 0.0,
"Er_nuc": 0.0,
"Er_tot": 0.0
@@ -266,28 +260,28 @@
"cycles": [
{
"energy_terms": {
- "E_ee": 3.354625825769642,
+ "E_ee": 3.3576788742677035,
"E_eext": 0.0,
- "E_el": -7.135135523826116,
- "E_en": -16.47750955461854,
- "E_kin": 7.571525854087563,
+ "E_el": -7.127946930407876,
+ "E_en": -16.4914760516857,
+ "E_kin": 7.584212118791978,
"E_next": 0.0,
"E_nn": 0.0,
- "E_nuc": -0.18173935348374062,
- "E_tot": -7.316874877309856,
- "E_x": -0.4192040376909349,
- "E_xc": -1.2857400054483297,
- "Er_el": 0.12116639407448325,
- "Er_nuc": -0.18173935348374062,
- "Er_tot": -0.06057295940925737
+ "E_nuc": -0.19212873389394816,
+ "E_tot": -7.320075664301824,
+ "E_x": -0.4195941455371656,
+ "E_xc": -1.2868490616228483,
+ "Er_el": 0.12808133537815772,
+ "Er_nuc": -0.19212873389394816,
+ "Er_tot": -0.06404739851579039
},
- "energy_total": -7.316874877309856,
- "energy_update": 0.09422578233054679,
- "mo_residual": 0.060097601021086314,
- "wall_time": 17.953042328
+ "energy_total": -7.320075664301824,
+ "energy_update": 0.0974265693225016,
+ "mo_residual": 0.06037084934246604,
+ "wall_time": 11.150312992
}
],
- "wall_time": 17.953267517
+ "wall_time": 11.155897119
},
"success": true
},
diff --git a/tests/li_solv/reference/li.out b/tests/li_solv/reference/li.out
index 7bc82fbaa..bacb33154 100644
--- a/tests/li_solv/reference/li.out
+++ b/tests/li_solv/reference/li.out
@@ -11,10 +11,10 @@
*** ***
*** VERSION 1.1.0-alpha ***
*** ***
-*** Git branch rebase-zora ***
-*** Git commit hash faecb83984c1d010ec90-dirty ***
-*** Git commit author Stig Rune Jensen ***
-*** Git commit date Thu Mar 10 17:22:51 2022 +0100 ***
+*** Git branch molecule-class-cleanup-input-parser ***
+*** Git commit hash d073311b4510a57f986e ***
+*** Git commit author Anders Brakestad ***
+*** Git commit date Mon Mar 28 13:46:33 2022 +0200 ***
*** ***
*** Contact: luca.frediani@uit.no ***
*** ***
@@ -46,14 +46,14 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
---------------------------------------------------------------------------
- MRCPP version : 1.4.1
+ MRCPP version : 1.4.0-alpha
Git branch : HEAD
- Git commit hash : 75d41879b1908a94a452
+ Git commit hash : 2797a7a18efb69942a6e
Git commit author : Stig Rune Jensen
- Git commit date : Thu Jan 6 11:38:53 2022 +0100
+ Git commit date : Fri Mar 12 12:48:07 2021 +0100
- Linear algebra : EIGEN v3.4.0
- Parallelization : OpenMP (1 threads)
+ Linear algebra : EIGEN v3.3.7
+ Parallelization : NONE
---------------------------------------------------------------------------
@@ -133,6 +133,7 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Calculation : Compute initial energy
Method : DFT (PBE0)
+ Relativity : None
Precision : 1.00000e-03
Localization : Off
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -178,6 +179,7 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Calculation : Optimize ground state orbitals
Method : DFT (PBE0)
+ Relativity : None
Checkpointing : Off
Max iterations : 100
KAIN solver : 5
@@ -188,7 +190,6 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
Helmholtz precision : Dynamic
Energy threshold : Off
Orbital threshold : 1.00000e-01
- ZORA potential : None
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -196,7 +197,7 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
Iter MO residual Total energy Update
---------------------------------------------------------------------------
0 1.000000e+00 -7.222649094979 -7.222649e+00
- 1 6.009760e-02 -7.316874877310 -9.422578e-02
+ 1 6.037085e-02 -7.320075664302 -9.742657e-02
---------------------------------------------------------------------------
SCF converged in 1 iterations!
===========================================================================
@@ -227,26 +228,26 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
===========================================================================
Molecular Energy (final)
---------------------------------------------------------------------------
- Kinetic energy : (au) 7.571525854088
- E-N energy : (au) -16.477509554619
- Coulomb energy : (au) 3.354625825770
- Exchange energy : (au) -0.419204037691
- X-C energy : (au) -1.285740005448
+ Kinetic energy : (au) 7.584212118792
+ E-N energy : (au) -16.491476051686
+ Coulomb energy : (au) 3.357678874268
+ Exchange energy : (au) -0.419594145537
+ X-C energy : (au) -1.286849061623
Ext. field (el) : (au) 0.000000000000
---------------------------------------------------------------------------
N-N energy : (au) 0.000000000000
Ext. field (nuc) : (au) 0.000000000000
- Tot. Reac. Energy : (au) -0.060572959409
- El. Reac. Energy : (au) 0.121166394074
- Nuc. Reac. Energy : (au) -0.181739353484
+ Tot. Reac. Energy : (au) -0.064047398516
+ El. Reac. Energy : (au) 0.128081335378
+ Nuc. Reac. Energy : (au) -0.192128733894
---------------------------------------------------------------------------
- Electronic energy : (au) -7.135135523826
- Nuclear energy : (au) -0.181739353484
+ Electronic energy : (au) -7.127946930408
+ Nuclear energy : (au) -0.192128733894
---------------------------------------------------------------------------
- Total energy : (au) -7.316874877310e+00
- : (kcal/mol) -4.591408269000e+03
- : (kJ/mol) -1.921045220996e+04
- : (eV) -1.991022996492e+02
+ Total energy : (au) -7.320075664302e+00
+ : (kcal/mol) -4.593416793146e+03
+ : (kJ/mol) -1.921885587500e+04
+ : (eV) -1.991893974965e+02
===========================================================================
@@ -255,9 +256,9 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
---------------------------------------------------------------------------
n Occ Spin : Epsilon
---------------------------------------------------------------------------
- 0 2 p : (au) -2.218888564195
+ 0 2 p : (au) -2.210659505993
---------------------------------------------------------------------------
- Sum occupied : (au) -4.437777128389
+ Sum occupied : (au) -4.421319011986
===========================================================================
@@ -266,7 +267,7 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
---------------------------------------------------------------------------
r_O : 0.000000 0.000000 0.000000
---------------------------------------------------------------------------
- Electronic vector : -0.000000 -0.000000 -0.000000
+ Electronic vector : -0.000000 0.000000 -0.000000
Magnitude : (au) 0.000000
: (Debye) 0.000000
---------------------------------------------------------------------------
@@ -286,7 +287,7 @@ J.Chem.Theor.Comp. 2010, DOI: 10.1021/ct100117s
*** ***
*** Exiting MRChem ***
*** ***
-*** Wall time : 0h 1m 5s ***
+*** Wall time : 0h 0m 48s ***
*** ***
***************************************************************************
diff --git a/tests/li_solv/test b/tests/li_solv/test
index de8fd98c5..81a10475d 100755
--- a/tests/li_solv/test
+++ b/tests/li_solv/test
@@ -24,6 +24,6 @@ filters = {
ER_NUC: rel_tolerance(1.0e-6),
}
-ierr = run(options, input_file="li", filters=filters)
+ierr = run(options, input_file="li", filters=filters, extra_args=['--json'])
sys.exit(ierr)