1""" 2Contains relevant physical constants 3""" 4 5import collections 6from decimal import Decimal 7from functools import lru_cache 8from typing import TYPE_CHECKING, Union 9 10from ..datum import Datum, print_variables 11from .ureg import build_units_registry 12 13if TYPE_CHECKING: 14 from pint import UnitRegistry, quantity # lgtm: [py/unused-import] 15 16 17class PhysicalConstantsContext: 18 r"""CODATA physical constants set from NIST. 19 20 Parameters 21 ---------- 22 context : str 23 {'CODATA2014', 'CODATA2018'} 24 Origin of loaded data. 25 26 Attributes 27 ---------- 28 doi : str 29 The DOI of the current context. 30 name : str 31 The name of the context ('CODATA2014') 32 pc : Dict[str, Datum] 33 Each physical constant is an entry in `pc`, where key is the 34 lowercased string of the NIST name (or any alias) and the 35 value is a Datum object with `lbl` the exact NIST name string, 36 `units`, `data` value as Decimal object, and any uncertainty 37 in the `comment` field. 38 raw_codata : Dict[str, Any] 39 A dictionary representation of the raw context data. 40 year : int 41 The year the context was created. 42 """ 43 44 _transtable = str.maketrans(" -/{", "__p_", ".,()}") 45 46 # Alias Typing 47 h: float 48 c: float 49 kb: float 50 R: float 51 bohr2angstroms: float 52 bohr2cm: float 53 amu2g: float 54 amu2kg: float 55 hartree2J: float 56 hartree2aJ: float 57 cal2J: float 58 dipmom_au2si: float 59 dipmom_au2debye: float 60 c_au: float 61 hartree2ev: float 62 hartree2wavenumbers: float 63 hartree2kcalmol: float 64 hartree2kJmol: float 65 hartree2MHz: float 66 kcalmol2wavenumbers: float 67 e0: float 68 na: float 69 me: float 70 71 def __init__(self, context="CODATA2014"): 72 self.pc = collections.OrderedDict() 73 74 if context == "CODATA2014": 75 from ..data import nist_2014_codata 76 77 self.doi = nist_2014_codata["doi"] 78 self.raw_codata = nist_2014_codata["constants"] 79 elif context == "CODATA2018": 80 from ..data import nist_2018_codata 81 82 self.doi = nist_2018_codata["doi"] 83 self.raw_codata = nist_2018_codata["constants"] 84 85 else: 86 raise KeyError("Context set as '{}', only contexts {'CODATA2014', 'CODATA2018', } are currently supported") 87 88 # physical constant loop 89 for k, v in self.raw_codata.items(): 90 self.pc[k] = Datum( 91 v["quantity"], 92 v["unit"], 93 Decimal(v["value"]), 94 comment="uncertainty={}".format(v["uncertainty"]), 95 doi=self.doi, 96 ) 97 98 self.name = context 99 self.year = int(context.replace("CODATA", "")) 100 self._ureg = None 101 102 # Extra relationships 103 self.pc["calorie-joule relationship"] = Datum( 104 "calorie-joule relationship", "J", Decimal("4.184"), comment="uncertainty=(exact)" 105 ) 106 107 rename_2018_from_2014 = { 108 "atomic unit of momentum": "atomic unit of mom.um", 109 "reduced Planck constant": "Planck constant over 2 pi", 110 "reduced Planck constant in eV s": "Planck constant over 2 pi in eV s", 111 "reduced Planck constant times c in MeV fm": "Planck constant over 2 pi times c in MeV fm", 112 "natural unit of momentum": "natural unit of mom.um", 113 "natural unit of momentum in MeV/c": "natural unit of mom.um in MeV/c", 114 "electron gyromag. ratio in MHz/T": "electron gyromag. ratio over 2 pi", 115 "vacuum mag. permeability": "mag. constant", 116 "lattice spacing of ideal Si (220)": "{220} lattice spacing of silicon", 117 "Planck constant in eV/Hz": "Planck constant in eV s", 118 "Bohr magneton in inverse meter per tesla": "Bohr magneton in inverse meters per tesla", 119 "Boltzmann constant in inverse meter per kelvin": "Boltzmann constant in inverse meters per kelvin", 120 "Copper x unit": "Cu x unit", 121 "Molybdenum x unit": "Mo x unit", 122 "proton gyromag. ratio in MHz/T": "proton gyromag. ratio over 2 pi", 123 "shielded proton gyromag. ratio in MHz/T": "shielded proton gyromag. ratio over 2 pi", 124 "reduced proton Compton wavelength": "proton Compton wavelength over 2 pi", 125 "reduced tau Compton wavelength": "tau Compton wavelength over 2 pi", 126 "tau energy equivalent": "tau mass energy equivalent in mev", 127 "reduced neutron Compton wavelength": "neutron Compton wavelength over 2 pi", 128 "neutron gyromag. ratio in MHz/T": "neutron gyromag. ratio over 2 pi", 129 "nuclear magneton in inverse meter per tesla": "nuclear magneton in inverse meters per tesla", 130 "shielded helion gyromag. ratio in MHz/T": "shielded helion gyromag. ratio over 2 pi", 131 "reduced Compton wavelength": "Compton wavelength over 2 pi", 132 "vacuum electric permittivity": "electric constant", 133 "reduced muon Compton wavelength": "muon Compton wavelength over 2 pi", 134 } 135 136 # fmt: off 137 if context == "CODATA2014": 138 aliases = [] 139 140 elif context == "CODATA2018": 141 for new_name, old_name in rename_2018_from_2014.items(): 142 dm = self.pc[new_name.lower()] 143 self.pc[old_name.lower()] = Datum(old_name, dm.units, dm.data, comment=dm.comment, doi=dm.doi) 144 145 aliases = [ 146 ("molar Planck constant times c", "J m mol^{-1}", self.pc["molar planck constant"].data * self.pc["speed of light in vacuum"].data, ""), 147 ("Faraday constant for conventional electric current", "C_{90} mol^{-1}", self.pc["faraday constant"].data / self.pc["conventional value of coulomb-90"].data, ""), 148 ("elementary charge over h", "A J^{-1}", self.pc["elementary charge over h-bar"].data / (2 * _get_pi(from_scratch=False)), ""), 149 ] 150 151 aliases.extend([ 152 ('h', 'J s', self.pc['hertz-joule relationship'].data, 'The Planck constant (Js)'), 153 ('hbar', 'J s', self.pc['planck constant over 2 pi'].data, 'Reduced Planck constant (Js)'), 154 ('c', 'm s^-1', self.pc['inverse meter-hertz relationship'].data, 'Speed of light (ms$^{-1}$)'), 155 ('kb', 'J K^-1', self.pc['kelvin-joule relationship'].data, 'The Boltzmann constant (JK$^{-1}$)'), 156 ('R', 'J mol^-1 K^-1', self.pc['molar gas constant'].data, 'Universal gas constant (JK$^{-1}$mol$^{-1}$)'), 157 ('bohr2angstroms', 'AA', self.pc['bohr radius'].data * Decimal('1.E10'), 'Bohr to Angstroms conversion factor'), 158 ('bohr2m', 'm', self.pc['bohr radius'].data, 'Bohr to meters conversion factor'), 159 ('bohr2cm', 'cm', self.pc['bohr radius'].data * Decimal('100'), 'Bohr to centimeters conversion factor'), 160 ('amu2g', 'g', self.pc['atomic mass constant'].data * Decimal('1000'), 'Atomic mass units to grams conversion factor'), 161 ('amu2kg', 'kg', self.pc['atomic mass constant'].data, 'Atomic mass units to kg conversion factor' ), 162 ('au2amu', 'u', self.pc['electron mass in u'].data, 'Atomic units (m$@@e$) to atomic mass units conversion factor'), 163 ('hartree2J', 'J', self.pc['hartree energy'].data, 'Hartree to joule conversion factor'), 164 ('hartree2aJ', 'aJ', self.pc['hartree energy'].data * Decimal('1.E18'), 'Hartree to attojoule (10$^{-18}$J) conversion factor'), 165 ('cal2J', 'J', self.pc['calorie-joule relationship'].data, 'Calorie to joule conversion factor'), 166 ('dipmom_au2si', 'C m', self.pc['atomic unit of electric dipole mom.'].data, 'Atomic units to SI units (Cm) conversion factor for dipoles'), 167 ('dipmom_au2debye', '???', self.pc['atomic unit of electric dipole mom.'].data * Decimal('1.E21') / self.pc['hertz-inverse meter relationship'].data, 168 'Atomic units to Debye conversion factor for dipoles'), 169 ('dipmom_debye2si', 'C m', self.pc['hertz-inverse meter relationship'].data * Decimal('1.E-21'), 'Debye to SI units (Cm) conversion factor for dipoles'), 170 ('c_au', '', self.pc['inverse fine-structure constant'].data, 'Speed of light in atomic units'), 171 ('hartree2ev', 'eV', self.pc['hartree energy in ev'].data, 'Hartree to eV conversion factor'), 172 ('hartree2wavenumbers', 'cm^-1', self.pc['hartree-inverse meter relationship'].data * Decimal('0.01'), 'Hartree to cm$^{-1}$ conversion factor'), 173 ('hartree2kcalmol', 'kcal mol^-1', self.pc['hartree energy'].data * self.pc['avogadro constant'].data * Decimal('0.001') / self.pc['calorie-joule relationship'].data, 174 'Hartree to kcal mol$^{-1}$ conversion factor'), 175 ('hartree2kJmol', 'kJ mol^-1', self.pc['hartree energy'].data * self.pc['avogadro constant'].data * Decimal('0.001'), 'Hartree to kilojoule mol$^{-1}$ conversion factor'), 176 ('hartree2MHz', 'MHz', self.pc['hartree-hertz relationship'].data * Decimal('1.E-6'), 'Hartree to MHz conversion factor'), 177 ('na', 'mol^-1', self.pc['avogadro constant'].data, "Avogadro's number"), 178 ('me', 'kg', self.pc['electron mass'].data, 'Electron rest mass (in kg)'), 179 ]) 180 181 if context == "CODATA2014": 182 aliases.extend([ 183 ('kcalmol2wavenumbers', 'kcal cm mol^-1', Decimal('10') * self.pc['calorie-joule relationship'].data / self.pc['molar planck constant times c'].data, 184 'kcal mol$^{-1}$ to cm$^{-1}$ conversion factor',), 185 ('e0', 'F m^-1', self.pc['electric constant'].data, 'Vacuum permittivity (Fm$^{-1}$)'), 186 ]) 187 elif context == "CODATA2018": 188 aliases.extend([ 189 ("kcalmol2wavenumbers", "kcal cm mol^-1", Decimal("10") * self.pc["calorie-joule relationship"].data / (self.pc["molar planck constant"].data * 190 self.pc["speed of light in vacuum"].data), "kcal mol$^{-1}$ to cm$^{-1}$ conversion factor",), 191 ("e0", "F m^-1", self.pc["vacuum electric permittivity"].data, "Vacuum permittivity (Fm$^{-1}$)"), 192 ]) 193 # fmt: on 194 195 # add alternate names for constants or derived values to help QC programs 196 for alias in aliases: 197 ident, units, value, comment = alias 198 self.pc[ident.lower()] = Datum(ident, units, value, comment=comment) 199 200 # add constants as directly callable member data 201 for qca in self.pc.values(): 202 callname = qca.label.translate(self._transtable) 203 setattr(self, callname, float(qca.data)) 204 205 def __str__(self) -> str: 206 return "PhysicalConstantsContext(context='{}')".format(self.name) 207 208 @property 209 def ureg(self) -> "UnitRegistry": 210 r"""Returns the internal Pint units registry. 211 212 Returns 213 ------- 214 UnitRegistry 215 The pint context 216 """ 217 if self._ureg is None: 218 self._ureg = build_units_registry(self) 219 220 return self._ureg 221 222 def get(self, physical_constant: str, return_tuple: bool = False) -> Union[float, Datum]: 223 r"""Access a physical constant, `physical_constant`. 224 225 Parameters 226 ---------- 227 physical_constant 228 Case-insensitive string of physical constant with NIST name. 229 return_tuple 230 See below. 231 232 Returns 233 ------- 234 Union[float, Datum] 235 When ``return_tuple=False``, value of physical constant. 236 When ``return_tuple=True``, Datum with units, description, uncertainty, and value of physical constant as Decimal. 237 238 """ 239 qca = self.pc[physical_constant.lower()] 240 241 if return_tuple: 242 return qca 243 else: 244 return float(qca.data) 245 246 # h 'hertz-joule relationship' = 6.62606896E-34 # The Planck constant (Js) 247 # c 'inverse meter-hertz relationship' = 2.99792458E8 # Speed of light (ms$^{-1}$) 248 # kb 'kelvin-joule relationship' = 1.3806504E-23 # The Boltzmann constant (JK$^{-1}$) 249 # R 'molar gas constant' = 8.314472 # Universal gas constant (JK$^{-1}$mol$^{-1}$) 250 # bohr2angstroms 'Bohr radius' * 1.E10 = 0.52917720859 # Bohr to Angstroms conversion factor 251 # bohr2m 'Bohr radius' = 0.52917720859E-10 # Bohr to meters conversion factor 252 # bohr2cm 'Bohr radius' * 100 = 0.52917720859E-8 # Bohr to centimeters conversion factor 253 # amu2kg 'atomic mass constant' = 1.660538782E-27 # Atomic mass units to kg conversion factor 254 # au2amu 'electron mass in u' = 5.485799097E-4 # Atomic units (m$@@e$) to atomic mass units conversion factor 255 # hartree2J 'Hartree energy' = 4.359744E-18 # Hartree to joule conversion factor 256 # hartree2aJ 'Hartree energy' * 1.E18 = 4.359744 # Hartree to attojoule (10$^{-18}$J) conversion factor 257 # cal2J = 4.184 # Calorie to joule conversion factor 258 # dipmom_au2si 'atomic unit of electric dipole mom.' = 8.47835281E-30 # Atomic units to SI units (Cm) conversion factor for dipoles 259 # 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 260 # dipmom_debye2si 'hertz-inverse meter relationship' * 1.E-21 = 3.335640952E-30 # Debye to SI units (Cm) conversion factor for dipoles 261 # c_au 'inverse fine-structure constant' = 137.035999679 # Speed of light in atomic units 262 # hartree2ev 'Hartree energy in eV' = 27.21138 # Hartree to eV conversion factor 263 # hartree2wavenumbers 'hartree-inverse meter relationship' * 0.01 = 219474.6 # Hartree to cm$^{-1}$ conversion factor 264 # hartree2kcalmol hartree2kJmol / cal2J = 627.5095 # Hartree to kcal mol$^{-1}$ conversion factor 265 # hartree2kJmol 'Hartree energy'*'Avogadro constant'*0.001 = 2625.500 # Hartree to kilojoule mol$^{-1}$ conversion factor 266 # hartree2MHz 'hartree-hertz relationship' = 6.579684E9 # Hartree to MHz conversion factor 267 # kcalmol2wavenumbers 10. / 'molar Planck constant times c'*4.184 = 349.7551 # kcal mol$^{-1}$ to cm$^{-1}$ conversion factor 268 # e0 'electric constant' = 8.854187817E-12 # Vacuum permittivity (Fm$^{-1}$) 269 # na 'Avogadro constant' = 6.02214179E23 # Avogadro's number 270 # me 'electron mass' = 9.10938215E-31 # Electron rest mass (in kg) 271 272 def Quantity(self, data: str) -> "quantity._Quantity": 273 """Returns a Pint Quantity.""" 274 275 return self.ureg.Quantity(data) 276 277 @lru_cache() 278 def conversion_factor( 279 self, base_unit: Union[str, "quantity._Quantity"], conv_unit: Union[str, "quantity._Quantity"] 280 ) -> float: 281 r"""Provides the conversion factor from one unit to another. 282 283 The conversion factor is based on the current contexts CODATA. 284 285 Parameters 286 ---------- 287 base_unit 288 The original units 289 conv_unit 290 The units to convert to 291 292 Examples 293 -------- 294 295 >>> conversion_factor("meter", "picometer") 296 1e-12 297 298 >>> conversion_factor("feet", "meter") 299 0.30479999999999996 300 301 >>> conversion_factor(10 * ureg.feet, "meter") 302 3.0479999999999996 303 304 Returns 305 ------- 306 float 307 The requested conversion factor 308 """ 309 310 # Add a little magic in case the incoming values have scalars 311 312 from pint import quantity 313 314 factor = 1.0 315 316 # First make sure we either have Units or Quantities 317 if isinstance(base_unit, str): 318 base_unit = self.ureg.parse_expression(base_unit) 319 320 if isinstance(conv_unit, str): 321 conv_unit = self.ureg.parse_expression(conv_unit) 322 323 # Need to play with prevector if we have Quantities 324 if isinstance(base_unit, quantity._Quantity): 325 factor *= base_unit.magnitude 326 base_unit = base_unit.units 327 328 if isinstance(conv_unit, quantity._Quantity): 329 factor /= conv_unit.magnitude 330 conv_unit = conv_unit.units 331 332 return self.ureg.convert(factor, base_unit, conv_unit) 333 334 def string_representation(self) -> str: 335 """Print name, value, and units of all physical constants.""" 336 337 return print_variables(self.pc) 338 339 340def run_comparison(context: str) -> None: 341 """Compare the existing physical constant information for Psi4 (in checkup_data folder) to object instantiation from 342 `context`. Specialized use.""" 343 344 self = PhysicalConstantsContext(context) 345 346 try: 347 from .. import checkup_data 348 except ImportError: # pragma: no cover 349 print("Info for comparison (directory checkup_data) not installed. Run from source.") 350 raise 351 352 class bcolors: 353 HEADER = "\033[95m" 354 OKBLUE = "\033[94m" 355 OKGREEN = "\033[92m" 356 WARNING = "\033[93m" 357 FAIL = "\033[91m" 358 ENDC = "\033[0m" 359 BOLD = "\033[1m" 360 UNDERLINE = "\033[4m" 361 362 tol = 1.0e-8 363 print(bcolors.OKBLUE + "\nChecking ({}) physconst vs. Psi4 ...".format(tol) + bcolors.ENDC) 364 for pc in dir(checkup_data.physconst): 365 if not pc.startswith("__"): 366 ref = self.get(pc) 367 val = getattr(checkup_data.physconst, pc) 368 assert isinstance(ref, (int, float)) 369 rat = abs(1.0 - float(ref) / val) 370 if rat > 1.0e-4: 371 print( 372 bcolors.FAIL 373 + "Physical Constant {} ratio differs by {:12.8f}: {} (this) vs {} (psi)".format(pc, rat, ref, val) 374 + bcolors.ENDC 375 ) 376 if rat > tol: 377 print("Physical Constant {} ratio differs by {:12.8f}: {} (this) vs {} (psi)".format(pc, rat, ref, val)) 378 379 380def run_internal_comparison(old_context: str, new_context: str) -> None: 381 """Compare two physical constant versions. Useful when adding new version.""" 382 383 self = PhysicalConstantsContext(old_context) 384 other = PhysicalConstantsContext(new_context) 385 386 class bcolors: 387 HEADER = "\033[95m" 388 OKBLUE = "\033[94m" 389 OKGREEN = "\033[92m" 390 WARNING = "\033[93m" 391 FAIL = "\033[91m" 392 ENDC = "\033[0m" 393 BOLD = "\033[1m" 394 UNDERLINE = "\033[4m" 395 396 both = set(self.pc.keys()).intersection(other.pc.keys()) 397 398 self_only = set(self.pc.keys()).difference(both) 399 print(bcolors.WARNING + f"\nPhysConst in old {self.name} missing in new {other.name} ..." + bcolors.ENDC) 400 for pc in sorted(self_only): 401 print(pc, self.get(pc)) 402 403 other_only = set(other.pc.keys()).difference(both) 404 print(bcolors.WARNING + f"\nPhysConst in new {other.name} missing in old {self.name} ..." + bcolors.ENDC) 405 for pc in sorted(other_only): 406 print(pc, other.get(pc)) 407 408 tol = 1.0e-8 409 print(bcolors.OKBLUE + f"\nChecking ({tol}) physconst {self.name} vs. {other.name} ..." + bcolors.ENDC) 410 for pc in sorted(both): 411 if not pc.startswith("__"): 412 ref = self.get(pc) 413 val = other.get(pc) 414 assert isinstance(ref, (int, float)) 415 rat = abs(1.0 - float(ref) / val) 416 if rat > 1.0e-4: 417 print( 418 bcolors.FAIL 419 + f"Physical Constant {pc} ratio differs by {rat:12.8f}: {ref} (this, {self.name}) vs {val} ({other.name})" 420 + bcolors.ENDC 421 ) 422 if rat > tol: 423 print( 424 f"Physical Constant {pc} ratio differs by {rat:12.8f}: {ref} (this, {self.name}) vs {val} ({other.name})" 425 ) 426 427 refu = self.pc[pc].units 428 valu = other.pc[pc].units 429 if refu != valu: 430 print( 431 bcolors.BOLD 432 + f"Physical Constant {pc} units differs: {refu} (this, {self.name}) vs {valu} ({other.name})" 433 + bcolors.ENDC 434 ) 435 436 437def _get_pi(from_scratch: bool = False) -> "Decimal": 438 r"""Get pi to 36 digits (or more with mpmath). 439 440 Parameters 441 ---------- 442 from_scratch 443 If True, recomputes Pi from mpmath. 444 445 Returns 446 ------- 447 Decimal 448 A representation of Pi 449 """ 450 451 if from_scratch: # pragma: no cover 452 from mpmath import mp 453 454 mp.dps = 36 455 return mp.pi 456 else: 457 return Decimal("3.14159265358979323846264338327950288") 458 459 460def write_c_header(context, filename="physconst.h", prefix="pc_"): 461 """Write C header file defining physical constants and pi, constants names are prefixed.""" 462 463 self = PhysicalConstantsContext(context) 464 465 pi = _get_pi(from_scratch=False) 466 tau = 2 * pi 467 468 text = [ 469 "#ifndef _qcelemental_physconst_h_", 470 "#define _qcelemental_physconst_h_", 471 "", 472 f"/* This file is autogenerated from the QCElemental python module from {self.name} */", 473 "", 474 "/* clang-format off */", 475 f"#define {prefix}pi {pi}", 476 f"#define {prefix}twopi {tau}", 477 ] 478 479 for pc, qca in sorted(self.pc.items()): 480 callname = qca.label.translate(self._transtable) 481 noncomment = "#define {}{} {}".format(prefix, callname, qca.data) 482 text.append("{:80} /*- {} [{}] {} -*/".format(noncomment, qca.label, qca.units, qca.comment)) 483 text.append("/* clang-format on */") 484 485 text.append("") 486 text.append("/* For Cray X1 compilers */") 487 text.append("#ifndef M_PI") 488 text.append("#define M_PI 3.14159265358979323846") 489 text.append("#endif") 490 text.append("") 491 492 text.append("#endif /* header guard */") 493 text.append("") 494 495 with open(filename, "w") as handle: 496 handle.write("\n".join(text)) 497 print("File written ({}). Remember to add license and clang-format it.".format(filename)) 498 499 500def write_fortran_header(context, filename="physconst.fh", prefix="pc_", kind=None): 501 """Write Fortran include file defining physical constants and pi, constants names are prefixed. 502 The name of the kind parameter can be passed as well.""" 503 504 self = PhysicalConstantsContext(context) 505 506 text = [ 507 f"! This file is autogenerated from the QCElemental python module from {self.name}", 508 ] 509 510 pi = _get_pi(from_scratch=False) 511 tau = 2 * pi 512 if kind is None: 513 declare = f"real, parameter ::" 514 text.append(f"{declare} {prefix}pi = {pi}") 515 text.append(f"{declare} {prefix}twopi = {tau}") 516 else: 517 declare = f"real({kind}), parameter ::" 518 text.append(f"{declare} {prefix}pi = {pi}_{kind}") 519 text.append(f"{declare} {prefix}twopi = {tau}_{kind}") 520 521 for pc, qca in sorted(self.pc.items()): 522 callname = qca.label.translate(self._transtable) 523 if kind is None: 524 noncomment = "{} {}{} = {}".format(declare, prefix, callname, qca.data) 525 else: 526 noncomment = "{} {}{} = {}_{}".format(declare, prefix, callname, qca.data, kind) 527 text.append("{:80} !! {} [{}] {}".format(noncomment, qca.label, qca.units, qca.comment)) 528 529 with open(filename, "w") as handle: 530 handle.write("\n".join(text)) 531 print("File written ({}). Remember to add license header.".format(filename)) 532 533 534# singleton 535constants = PhysicalConstantsContext("CODATA2014") 536