Proper handling of physical constants

[stigrj]

E.g. through https://github.com/MolSSI/QCElemental, but this requires new C++ interface

[robertodr]

What about the mrchem Python script handling conversions to/from AU at the beginning/end of a calculation? mrchem.x could dump the computed data to JSON, re-read it in Python and convert it to whatever units.

[Brakjen]

From @stigrj in #404 (comment)

I think one nice way to handle the constants is to let Python define them and feed them into the C++ program through the JSON input. Then we also have complete control and can tweak the parameters without recompiling, by intercepting the input parser and edit the JSON file before execution.

Is something like this what you had in mind?

from qcelemental.physical_constants.context import PhysicalConstantsContext
from qcelemental.datum import Datum


# Subclass it here so that the same (non-default) context can easily be used
# when importing the class
# We could also define custom aliases here
class MRChemConstants(PhysicalConstantsContext):
    def __init__(self, context="CODATA2018"):
        super().__init__(context)

        customs = [
            ('pi',      '', 3.1415926535897932384, 'Pi'),
            ('pi_sqrt', '', 1.7724538509055160273, 'Square root of pi')
        ]

        # Add aliases to internally stored constants and make them callable
        for ident, units, value, comment in customs:
            self.pc[ident.lower()] = Datum(ident, units, value, comment=comment)
            self.__setattr__(ident, value)

    def to_json(self):
        """Generate a JSON file storing the constants, which can be passed to the c++ program."""
        return {
            qca.label: float(qca.data) for qca in self.pc.values()
        }

And then in api.py

import math

from .helpers import *
from .CUBEparser import *
from .periodictable import PeriodicTable as PT, PeriodicTableByZ as PT_Z
from .physical_constants import MRChemConstants

pc = MRChemConstants()


def translate_input(user_dict):
    # get the origin in the desired units of measure
    origin = user_dict["world_origin"]
    if user_dict["world_unit"] == "angstrom":
        origin = [r / pc.bohr2angstroms for r in origin]

    # prepare bits and pieces
    mol_dict = write_molecule(user_dict, origin)
    mpi_dict = write_mpi(user_dict)
    mra_dict = write_mra(user_dict, mol_dict)
    scf_dict = write_scf_calculation(user_dict, origin)
    rsp_dict = write_rsp_calculations(user_dict, mol_dict, origin)
    constants_dict = pc.to_json()

    # piece everything together
    program_dict = {
        "schema_name": "mrchem_input",
        "schema_version": 1,
        "mpi": mpi_dict,
        "mra": mra_dict,
        "printer": user_dict["Printer"],
        "molecule": mol_dict,
        "scf_calculation": scf_dict,
        "rsp_calculations": rsp_dict,
        "constants": constants_dict
    }
    return program_dict

[stigrj]

Yes, this looks excellent! Then we'll need a structure in the C++ code to house these, perhaps a Singleton class? We should probably do the same for the periodic table data.