Source code for qcelemental.physical_constants.context

"""
Contains relevant physical constants
"""

import collections
from decimal import Decimal
from functools import lru_cache
from typing import TYPE_CHECKING, Union

from ..datum import Datum, print_variables
from .ureg import build_units_registry

if TYPE_CHECKING:
    from pint import Quantity as _Quantity  # lgtm: [py/unused-import]
    from pint import UnitRegistry  # lgtm: [py/unused-import]


[docs]class PhysicalConstantsContext: r"""CODATA physical constants set from NIST. Parameters ---------- context : str {'CODATA2014', 'CODATA2018'} Origin of loaded data. Attributes ---------- doi : str The DOI of the current context. name : str The name of the context ('CODATA2014') pc : Dict[str, Datum] Each physical constant is an entry in `pc`, where key is the lowercased string of the NIST name (or any alias) and the value is a Datum object with `lbl` the exact NIST name string, `units`, `data` value as Decimal object, and any uncertainty in the `comment` field. raw_codata : Dict[str, Any] A dictionary representation of the raw context data. year : int The year the context was created. """ _transtable = str.maketrans(" -/{", "__p_", ".,()}") # Alias Typing h: float c: float kb: float R: float bohr2angstroms: float bohr2cm: float amu2g: float amu2kg: float hartree2J: float hartree2aJ: float cal2J: float dipmom_au2si: float dipmom_au2debye: float c_au: float hartree2ev: float hartree2wavenumbers: float hartree2kcalmol: float hartree2kJmol: float hartree2MHz: float kcalmol2wavenumbers: float e0: float na: float me: float def __init__(self, context="CODATA2014"): self.pc = collections.OrderedDict() if context == "CODATA2014": from ..data import nist_2014_codata self.doi = nist_2014_codata["doi"] self.raw_codata = nist_2014_codata["constants"] elif context == "CODATA2018": from ..data import nist_2018_codata self.doi = nist_2018_codata["doi"] self.raw_codata = nist_2018_codata["constants"] else: raise KeyError("Context set as '{}', only contexts {'CODATA2014', 'CODATA2018', } are currently supported") # physical constant loop for k, v in self.raw_codata.items(): self.pc[k] = Datum( v["quantity"], v["unit"], Decimal(v["value"]), comment="uncertainty={}".format(v["uncertainty"]), doi=self.doi, ) self.name = context self.year = int(context.replace("CODATA", "")) self._ureg = None # Extra relationships self.pc["calorie-joule relationship"] = Datum( "calorie-joule relationship", "J", Decimal("4.184"), comment="uncertainty=(exact)" ) rename_2018_from_2014 = { "atomic unit of momentum": "atomic unit of mom.um", "reduced Planck constant": "Planck constant over 2 pi", "reduced Planck constant in eV s": "Planck constant over 2 pi in eV s", "reduced Planck constant times c in MeV fm": "Planck constant over 2 pi times c in MeV fm", "natural unit of momentum": "natural unit of mom.um", "natural unit of momentum in MeV/c": "natural unit of mom.um in MeV/c", "electron gyromag. ratio in MHz/T": "electron gyromag. ratio over 2 pi", "vacuum mag. permeability": "mag. constant", "lattice spacing of ideal Si (220)": "{220} lattice spacing of silicon", "Planck constant in eV/Hz": "Planck constant in eV s", "Bohr magneton in inverse meter per tesla": "Bohr magneton in inverse meters per tesla", "Boltzmann constant in inverse meter per kelvin": "Boltzmann constant in inverse meters per kelvin", "Copper x unit": "Cu x unit", "Molybdenum x unit": "Mo x unit", "proton gyromag. ratio in MHz/T": "proton gyromag. ratio over 2 pi", "shielded proton gyromag. ratio in MHz/T": "shielded proton gyromag. ratio over 2 pi", "reduced proton Compton wavelength": "proton Compton wavelength over 2 pi", "reduced tau Compton wavelength": "tau Compton wavelength over 2 pi", "tau energy equivalent": "tau mass energy equivalent in mev", "reduced neutron Compton wavelength": "neutron Compton wavelength over 2 pi", "neutron gyromag. ratio in MHz/T": "neutron gyromag. ratio over 2 pi", "nuclear magneton in inverse meter per tesla": "nuclear magneton in inverse meters per tesla", "shielded helion gyromag. ratio in MHz/T": "shielded helion gyromag. ratio over 2 pi", "reduced Compton wavelength": "Compton wavelength over 2 pi", "vacuum electric permittivity": "electric constant", "reduced muon Compton wavelength": "muon Compton wavelength over 2 pi", } # fmt: off if context == "CODATA2014": aliases = [] elif context == "CODATA2018": for new_name, old_name in rename_2018_from_2014.items(): dm = self.pc[new_name.lower()] self.pc[old_name.lower()] = Datum(old_name, dm.units, dm.data, comment=dm.comment, doi=dm.doi) aliases = [ ("molar Planck constant times c", "J m mol^{-1}", self.pc["molar planck constant"].data * self.pc["speed of light in vacuum"].data, ""), ("Faraday constant for conventional electric current", "C_{90} mol^{-1}", self.pc["faraday constant"].data / self.pc["conventional value of coulomb-90"].data, ""), ("elementary charge over h", "A J^{-1}", self.pc["elementary charge over h-bar"].data / (2 * _get_pi(from_scratch=False)), ""), ] aliases.extend([ ('h', 'J s', self.pc['hertz-joule relationship'].data, 'The Planck constant (Js)'), ('hbar', 'J s', self.pc['planck constant over 2 pi'].data, 'Reduced Planck constant (Js)'), ('c', 'm s^-1', self.pc['inverse meter-hertz relationship'].data, 'Speed of light (ms$^{-1}$)'), ('kb', 'J K^-1', self.pc['kelvin-joule relationship'].data, 'The Boltzmann constant (JK$^{-1}$)'), ('R', 'J mol^-1 K^-1', self.pc['molar gas constant'].data, 'Universal gas constant (JK$^{-1}$mol$^{-1}$)'), ('bohr2angstroms', 'AA', self.pc['bohr radius'].data * Decimal('1.E10'), 'Bohr to Angstroms conversion factor'), ('bohr2m', 'm', self.pc['bohr radius'].data, 'Bohr to meters conversion factor'), ('bohr2cm', 'cm', self.pc['bohr radius'].data * Decimal('100'), 'Bohr to centimeters conversion factor'), ('amu2g', 'g', self.pc['atomic mass constant'].data * Decimal('1000'), 'Atomic mass units to grams conversion factor'), ('amu2kg', 'kg', self.pc['atomic mass constant'].data, 'Atomic mass units to kg conversion factor' ), ('au2amu', 'u', self.pc['electron mass in u'].data, 'Atomic units (m$@@e$) to atomic mass units conversion factor'), ('hartree2J', 'J', self.pc['hartree energy'].data, 'Hartree to joule conversion factor'), ('hartree2aJ', 'aJ', self.pc['hartree energy'].data * Decimal('1.E18'), 'Hartree to attojoule (10$^{-18}$J) conversion factor'), ('cal2J', 'J', self.pc['calorie-joule relationship'].data, 'Calorie to joule conversion factor'), ('dipmom_au2si', 'C m', self.pc['atomic unit of electric dipole mom.'].data, 'Atomic units to SI units (Cm) conversion factor for dipoles'), ('dipmom_au2debye', '???', self.pc['atomic unit of electric dipole mom.'].data * Decimal('1.E21') / self.pc['hertz-inverse meter relationship'].data, 'Atomic units to Debye conversion factor for dipoles'), ('dipmom_debye2si', 'C m', self.pc['hertz-inverse meter relationship'].data * Decimal('1.E-21'), 'Debye to SI units (Cm) conversion factor for dipoles'), ('c_au', '', self.pc['inverse fine-structure constant'].data, 'Speed of light in atomic units'), ('hartree2ev', 'eV', self.pc['hartree energy in ev'].data, 'Hartree to eV conversion factor'), ('hartree2wavenumbers', 'cm^-1', self.pc['hartree-inverse meter relationship'].data * Decimal('0.01'), 'Hartree to cm$^{-1}$ conversion factor'), ('hartree2kcalmol', 'kcal mol^-1', self.pc['hartree energy'].data * self.pc['avogadro constant'].data * Decimal('0.001') / self.pc['calorie-joule relationship'].data, 'Hartree to kcal mol$^{-1}$ conversion factor'), ('hartree2kJmol', 'kJ mol^-1', self.pc['hartree energy'].data * self.pc['avogadro constant'].data * Decimal('0.001'), 'Hartree to kilojoule mol$^{-1}$ conversion factor'), ('hartree2MHz', 'MHz', self.pc['hartree-hertz relationship'].data * Decimal('1.E-6'), 'Hartree to MHz conversion factor'), ('na', 'mol^-1', self.pc['avogadro constant'].data, "Avogadro's number"), ('me', 'kg', self.pc['electron mass'].data, 'Electron rest mass (in kg)'), ]) if context == "CODATA2014": aliases.extend([ ('kcalmol2wavenumbers', 'kcal cm mol^-1', Decimal('10') * self.pc['calorie-joule relationship'].data / self.pc['molar planck constant times c'].data, 'kcal mol$^{-1}$ to cm$^{-1}$ conversion factor',), ('e0', 'F m^-1', self.pc['electric constant'].data, 'Vacuum permittivity (Fm$^{-1}$)'), ]) elif context == "CODATA2018": aliases.extend([ ("kcalmol2wavenumbers", "kcal cm mol^-1", Decimal("10") * self.pc["calorie-joule relationship"].data / (self.pc["molar planck constant"].data * self.pc["speed of light in vacuum"].data), "kcal mol$^{-1}$ to cm$^{-1}$ conversion factor",), ("e0", "F m^-1", self.pc["vacuum electric permittivity"].data, "Vacuum permittivity (Fm$^{-1}$)"), ]) # fmt: on # add alternate names for constants or derived values to help QC programs for alias in aliases: ident, units, value, comment = alias self.pc[ident.lower()] = Datum(ident, units, value, comment=comment) # add constants as directly callable member data for qca in self.pc.values(): callname = qca.label.translate(self._transtable) setattr(self, callname, float(qca.data)) def __str__(self) -> str: return "PhysicalConstantsContext(context='{}')".format(self.name) @property def ureg(self) -> "UnitRegistry": r"""Returns the internal Pint units registry. Returns ------- UnitRegistry The pint context """ if self._ureg is None: self._ureg = build_units_registry(self) return self._ureg
[docs] def get(self, physical_constant: str, return_tuple: bool = False) -> Union[float, Datum]: r"""Access a physical constant, `physical_constant`. Parameters ---------- physical_constant Case-insensitive string of physical constant with NIST name. return_tuple See below. Returns ------- Union[float, Datum] When ``return_tuple=False``, value of physical constant. When ``return_tuple=True``, Datum with units, description, uncertainty, and value of physical constant as Decimal. """ qca = self.pc[physical_constant.lower()] if return_tuple: return qca else: return float(qca.data)
# h 'hertz-joule relationship' = 6.62606896E-34 # The Planck constant (Js) # c 'inverse meter-hertz relationship' = 2.99792458E8 # Speed of light (ms$^{-1}$) # kb 'kelvin-joule relationship' = 1.3806504E-23 # The Boltzmann constant (JK$^{-1}$) # R 'molar gas constant' = 8.314472 # Universal gas constant (JK$^{-1}$mol$^{-1}$) # bohr2angstroms 'Bohr radius' * 1.E10 = 0.52917720859 # Bohr to Angstroms conversion factor # bohr2m 'Bohr radius' = 0.52917720859E-10 # Bohr to meters conversion factor # bohr2cm 'Bohr radius' * 100 = 0.52917720859E-8 # Bohr to centimeters conversion factor # amu2kg 'atomic mass constant' = 1.660538782E-27 # Atomic mass units to kg conversion factor # au2amu 'electron mass in u' = 5.485799097E-4 # Atomic units (m$@@e$) to atomic mass units conversion factor # hartree2J 'Hartree energy' = 4.359744E-18 # Hartree to joule conversion factor # hartree2aJ 'Hartree energy' * 1.E18 = 4.359744 # Hartree to attojoule (10$^{-18}$J) conversion factor # cal2J = 4.184 # Calorie to joule conversion factor # dipmom_au2si 'atomic unit of electric dipole mom.' = 8.47835281E-30 # Atomic units to SI units (Cm) conversion factor for dipoles # dipmom_au2debye 'atomic unit of electric dipole mom.' / ('hertz-inverse meter relationship' * 1.E-21) = 2.54174623 # Atomic units to Debye conversion factor for dipoles # dipmom_debye2si 'hertz-inverse meter relationship' * 1.E-21 = 3.335640952E-30 # Debye to SI units (Cm) conversion factor for dipoles # c_au 'inverse fine-structure constant' = 137.035999679 # Speed of light in atomic units # hartree2ev 'Hartree energy in eV' = 27.21138 # Hartree to eV conversion factor # hartree2wavenumbers 'hartree-inverse meter relationship' * 0.01 = 219474.6 # Hartree to cm$^{-1}$ conversion factor # hartree2kcalmol hartree2kJmol / cal2J = 627.5095 # Hartree to kcal mol$^{-1}$ conversion factor # hartree2kJmol 'Hartree energy'*'Avogadro constant'*0.001 = 2625.500 # Hartree to kilojoule mol$^{-1}$ conversion factor # hartree2MHz 'hartree-hertz relationship' = 6.579684E9 # Hartree to MHz conversion factor # kcalmol2wavenumbers 10. / 'molar Planck constant times c'*4.184 = 349.7551 # kcal mol$^{-1}$ to cm$^{-1}$ conversion factor # e0 'electric constant' = 8.854187817E-12 # Vacuum permittivity (Fm$^{-1}$) # na 'Avogadro constant' = 6.02214179E23 # Avogadro's number # me 'electron mass' = 9.10938215E-31 # Electron rest mass (in kg)
[docs] def Quantity(self, data: str) -> "_Quantity": """Returns a Pint Quantity.""" return self.ureg.Quantity(data)
[docs] @lru_cache() def conversion_factor(self, base_unit: Union[str, "_Quantity"], conv_unit: Union[str, "_Quantity"]) -> float: r"""Provides the conversion factor from one unit to another. The conversion factor is based on the current contexts CODATA. Parameters ---------- base_unit The original units conv_unit The units to convert to Examples -------- >>> conversion_factor("meter", "picometer") 1e-12 >>> conversion_factor("feet", "meter") 0.30479999999999996 >>> conversion_factor(10 * ureg.feet, "meter") 3.0479999999999996 Returns ------- float The requested conversion factor """ # Add a little magic in case the incoming values have scalars from pint import Quantity as _Quantity factor = 1.0 # First make sure we either have Units or Quantities if isinstance(base_unit, str): base_unit = self.ureg.parse_expression(base_unit) if isinstance(conv_unit, str): conv_unit = self.ureg.parse_expression(conv_unit) # Need to play with prevector if we have Quantities if isinstance(base_unit, _Quantity): factor *= base_unit.magnitude base_unit = base_unit.units if isinstance(conv_unit, _Quantity): factor /= conv_unit.magnitude conv_unit = conv_unit.units return self.ureg.convert(factor, base_unit, conv_unit)
[docs] def string_representation(self) -> str: """Print name, value, and units of all physical constants.""" return print_variables(self.pc)
def run_comparison(context: str) -> None: """Compare the existing physical constant information for Psi4 (in checkup_data folder) to object instantiation from `context`. Specialized use.""" self = PhysicalConstantsContext(context) try: from .. import checkup_data except ImportError: # pragma: no cover print("Info for comparison (directory checkup_data) not installed. Run from source.") raise class bcolors: HEADER = "\033[95m" OKBLUE = "\033[94m" OKGREEN = "\033[92m" WARNING = "\033[93m" FAIL = "\033[91m" ENDC = "\033[0m" BOLD = "\033[1m" UNDERLINE = "\033[4m" tol = 1.0e-8 print(bcolors.OKBLUE + "\nChecking ({}) physconst vs. Psi4 ...".format(tol) + bcolors.ENDC) for pc in dir(checkup_data.physconst): if not pc.startswith("__"): ref = self.get(pc) val = getattr(checkup_data.physconst, pc) assert isinstance(ref, (int, float)) rat = abs(1.0 - float(ref) / val) if rat > 1.0e-4: print( bcolors.FAIL + "Physical Constant {} ratio differs by {:12.8f}: {} (this) vs {} (psi)".format(pc, rat, ref, val) + bcolors.ENDC ) if rat > tol: print("Physical Constant {} ratio differs by {:12.8f}: {} (this) vs {} (psi)".format(pc, rat, ref, val)) def run_internal_comparison(old_context: str, new_context: str) -> None: """Compare two physical constant versions. Useful when adding new version.""" self = PhysicalConstantsContext(old_context) other = PhysicalConstantsContext(new_context) class bcolors: HEADER = "\033[95m" OKBLUE = "\033[94m" OKGREEN = "\033[92m" WARNING = "\033[93m" FAIL = "\033[91m" ENDC = "\033[0m" BOLD = "\033[1m" UNDERLINE = "\033[4m" both = set(self.pc.keys()).intersection(other.pc.keys()) self_only = set(self.pc.keys()).difference(both) print(bcolors.WARNING + f"\nPhysConst in old {self.name} missing in new {other.name} ..." + bcolors.ENDC) for pc in sorted(self_only): print(pc, self.get(pc)) other_only = set(other.pc.keys()).difference(both) print(bcolors.WARNING + f"\nPhysConst in new {other.name} missing in old {self.name} ..." + bcolors.ENDC) for pc in sorted(other_only): print(pc, other.get(pc)) tol = 1.0e-8 print(bcolors.OKBLUE + f"\nChecking ({tol}) physconst {self.name} vs. {other.name} ..." + bcolors.ENDC) for pc in sorted(both): if not pc.startswith("__"): ref = self.get(pc) val = other.get(pc) assert isinstance(ref, (int, float)) rat = abs(1.0 - float(ref) / val) if rat > 1.0e-4: print( bcolors.FAIL + f"Physical Constant {pc} ratio differs by {rat:12.8f}: {ref} (this, {self.name}) vs {val} ({other.name})" + bcolors.ENDC ) if rat > tol: print( f"Physical Constant {pc} ratio differs by {rat:12.8f}: {ref} (this, {self.name}) vs {val} ({other.name})" ) refu = self.pc[pc].units valu = other.pc[pc].units if refu != valu: print( bcolors.BOLD + f"Physical Constant {pc} units differs: {refu} (this, {self.name}) vs {valu} ({other.name})" + bcolors.ENDC ) def _get_pi(from_scratch: bool = False) -> "Decimal": r"""Get pi to 36 digits (or more with mpmath). Parameters ---------- from_scratch If True, recomputes Pi from mpmath. Returns ------- Decimal A representation of Pi """ if from_scratch: # pragma: no cover from mpmath import mp mp.dps = 36 return mp.pi else: return Decimal("3.14159265358979323846264338327950288") def write_c_header(context, filename="physconst.h", prefix="pc_"): """Write C header file defining physical constants and pi, constants names are prefixed.""" self = PhysicalConstantsContext(context) pi = _get_pi(from_scratch=False) tau = 2 * pi text = [ "#ifndef _qcelemental_physconst_h_", "#define _qcelemental_physconst_h_", "", f"/* This file is autogenerated from the QCElemental python module from {self.name} */", "", "/* clang-format off */", f"#define {prefix}pi {pi}", f"#define {prefix}twopi {tau}", ] for pc, qca in sorted(self.pc.items()): callname = qca.label.translate(self._transtable) noncomment = "#define {}{} {}".format(prefix, callname, qca.data) text.append("{:80} /*- {} [{}] {} -*/".format(noncomment, qca.label, qca.units, qca.comment)) text.append("/* clang-format on */") text.append("") text.append("/* For Cray X1 compilers */") text.append("#ifndef M_PI") text.append("#define M_PI 3.14159265358979323846") text.append("#endif") text.append("") text.append("#endif /* header guard */") text.append("") with open(filename, "w") as handle: handle.write("\n".join(text)) print("File written ({}). Remember to add license and clang-format it.".format(filename)) def write_fortran_header(context, filename="physconst.fh", prefix="pc_", kind=None): """Write Fortran include file defining physical constants and pi, constants names are prefixed. The name of the kind parameter can be passed as well.""" self = PhysicalConstantsContext(context) text = [ f"! This file is autogenerated from the QCElemental python module from {self.name}", ] pi = _get_pi(from_scratch=False) tau = 2 * pi if kind is None: declare = f"real, parameter ::" text.append(f"{declare} {prefix}pi = {pi}") text.append(f"{declare} {prefix}twopi = {tau}") else: declare = f"real({kind}), parameter ::" text.append(f"{declare} {prefix}pi = {pi}_{kind}") text.append(f"{declare} {prefix}twopi = {tau}_{kind}") for pc, qca in sorted(self.pc.items()): callname = qca.label.translate(self._transtable) if kind is None: noncomment = "{} {}{} = {}".format(declare, prefix, callname, qca.data) else: noncomment = "{} {}{} = {}_{}".format(declare, prefix, callname, qca.data, kind) text.append("{:80} !! {} [{}] {}".format(noncomment, qca.label, qca.units, qca.comment)) with open(filename, "w") as handle: handle.write("\n".join(text)) print("File written ({}). Remember to add license header.".format(filename)) # singleton constants = PhysicalConstantsContext("CODATA2014")