Update Physical Constants and Corresponding Unit Tests

[ttpekkan]

Motivation or Problem

RMG currently uses outdated values for many physical constants. Many of these constants were given exact values when they were last updated, so now we have an opportunity to fix the values of these once and for all. This likely will have little to no difference in practice but is the correct thing to do.

Description of Changes

I changed the physical constants to the values reported here: https://arxiv.org/pdf/2409.03787 For atomic weights, I used the most recent values I found from IUPAC: https://iupac.qmul.ac.uk/AtWt/ Note that for some important elements like H and O it's not possible to give an unambiguous isotope-corrected atomic weight as it depends on the environment. I just went with the average value given in the table. Also, converting the dimensionless atomic weights to mol/kg requires multiplication by the molar mass constant (1.0000000010531e-3), which is no longer exactly 1.0e-3 due to how the physical constant were updated. I did this conversion, and although the correct thing to do, does not make a difference in practice. This also slightly affects the atomic mass unit as it is defined as the molar mass constant divided by Avogadro's constant. I did not include unstable elements. These can be added as needed, but I don't think RMG currently supports using them anyway.

Testing

I ran make–test, and as expected, it threw quite a few errors because in the unit tests some RMG output values are compared with pre-defined values. The RMG output values changed a tiny amount because of the physical-constant and atomic-weight updates, but enough for the unit tests to fail. I simply updated the pre-set values in the unit tests to those that RMG outputs with the updated physical constant and atomic weights. In a few cases I adjusted the tolerances for approximate equality comparisons to make the unit tests pass. Most of the errors appear to come from the change to the Hartree --> kJ/mol conversion factor.

Reviewer Tips

I guess what mostly needs to be done is to check that the physical constants do not have typos.

[JacksonBurns]

Hi @ttpekkan welcome to RMG!

For future reference, this is somewhat of a revival of this very old PR: #1627

This PR will resolve the previously abandoned #2435

[JacksonBurns]

This is the current status of the tests:

 =========================== short test summary info ============================
FAILED test/arkane/arkaneOutputTest.py::OutputTest::test_prettify - assert "    E0 = (193.749, 'kJ/mol'),\n" in ['# Coordinates for C7H7 in Input Orientation (angstroms):\n', '#   C    0.0000    0.0000   -1.9093\n', '#   C    0.0000    1.2097   -1.2046\n', '#   C    0.0000    1.2164    0.1782\n', '#   C    0.0000    0.0000    0.9209\n', '#   C    0.0000   -1.2164    0.1782\n', '#   C    0.0000   -1.2097   -1.2046\n', '#   H    0.0000    0.0000   -2.9930\n', '#   H    0.0000    2.1490   -1.7469\n', '#   H    0.0000    2.1576    0.7177\n', '#   H    0.0000   -2.1576    0.7177\n', '#   H    0.0000   -2.1490   -1.7469\n', '#   C    0.0000    0.0000    2.3249\n', '#   H    0.0000    0.9273    2.8840\n', '#   H    0.0000   -0.9273    2.8840\n', 'conformer(\n', "    label = 'C7H7',\n", "    E0 = (193.688, 'kJ/mol'),\n", '    modes = [\n', "        IdealGasTranslation(mass=(91.0548, 'amu')),\n", '        NonlinearRotor(\n', "            inertia = ([91.0567, 186.675, 277.733], 'amu*angstrom^2'),\n", '            symmetry = 2,\n', '        ),\n', '        HarmonicOscillator(\n', "            frequencies = ([199.381, 360.536, 413.795, 480.347, 536.285, 630.723, 687.118, 709.613, 776.662, 830.404, 834.386, 901.841, 973.498, 975.148, 993.349, 998.606, 1040.14, 1120.69, 1179.22, 1189.07, 1292.86, 1332.91, 1357.18, 1479.46, 1495.36, 1507.91, 1583.14, 1604.63, 3156.85, 3170.22, 3172.78, 3185.05, 3189.8, 3203.23, 3253.99], 'cm^-1'),\n", '        ),\n', '        HinderedRotor(\n', "            inertia = (1.70013, 'amu*angstrom^2'),\n", '            symmetry = 2,\n', '            fourier = (\n', '                [\n', '                    [-0.31592, -27.7665, 0.177659, 3.24367, 0.0509783],\n', '                    [-0.00152732, -0.00192155, -0.00354267, -0.000523773, 0.00324541],\n', '                ],\n', "                'kJ/mol',\n", '            ),\n', '            quantum = None,\n', '            semiclassical = None,\n', '        ),\n', '    ],\n', '    spin_multiplicity = 2,\n', '    optical_isomers = 1,\n', ')\n', '\n', '# Thermodynamics for C7H7:\n', '#   Enthalpy of formation (298 K)   =    50.553 kcal/mol\n', '#   Entropy of formation (298 K)    =    75.505 cal/(mol*K)\n', '#    =========== =========== =========== =========== ===========\n', '#    Temperature Heat cap.   Enthalpy    Entropy     Free energy\n', '#    (K)         (cal/mol*K) (kcal/mol)  (cal/mol*K) (kcal/mol)\n', '#    =========== =========== =========== =========== ===========\n', '#            300      25.640      50.604      75.676      27.902\n', '#            400      33.155      53.551      84.101      19.911\n', '#            500      39.702      57.202      92.222      11.092\n', '#            600      45.261      61.458      99.966       1.479\n', '#            800      53.880      71.423     114.239     -19.969\n', '#           1000      59.876      82.836     126.948     -44.112\n', '#           1500      67.906     115.082     152.988    -114.400\n', '#           2000      71.822     150.084     173.100    -196.117\n', '#           2400      74.219     179.310     186.415    -268.086\n', '#    =========== =========== =========== =========== ===========\n', 'thermo(\n', "    label = 'C7H7',\n", '    thermo = NASA(\n', '        polynomials = [\n', '            NASAPolynomial(\n', '                coeffs = [3.95186, 0.0029882, 0.000162149, -2.93429e-07, 1.70812e-10, 23197, 9.643],\n', "                Tmin = (10, 'K'),\n", "                Tmax = (445.411, 'K'),\n", '            ),\n', '            NASAPolynomial(\n', '                coeffs = [-2.82364, 0.0638647, -4.29617e-05, 1.37173e-08, -1.66592e-12, 23800.3, 36.8479],\n', "                Tmin = (445.411, 'K'),\n", "                Tmax = (3000, 'K'),\n", '            ),\n', '        ],\n', "        Tmin = (10, 'K'),\n", "        Tmax = (3000, 'K'),\n", "        E0 = (192.856, 'kJ/mol'),\n", "        Cp0 = (33.2579, 'J/(mol*K)'),\n", "        CpInf = (328.421, 'J/(mol*K)'),\n", '    ),\n', ')\n', '\n']
FAILED test/rmgpy/reactionTest.py::TestChargeTransferReaction::test_get_rate_coeff - assert 0.0078125 < 1e-06
 +  where 0.0078125 = abs((43870307169260.06 - 43870307169260.055))
== 2 failed, 1968 passed, 39 skipped, 3937806 warnings in 2329.92s (0:38:49) ===

These both look like regular old rounding errors. The first one is because E0 = (193.749, 'kJ/mol') is not present in the output, instead we have E0 = (193.688, 'kJ/mol') which is just a tiny difference. The other failure is just an incorrectly set reference value, or a too-tight tolerance (take your pick when issuing a fix).

Thanks again @ttpekkan!