incorrect free rotor entropy due to implicit type conversion #84
Description
Description:
I encountered an issue in GoodVibes 3.2 where the free rotor entropy calculation (calc_freerot_entropy) produces incorrect results due to an type conversion. Specifically, freq_scale_factor is passed as a numpy.float32 and it triggers a type conversion, leading to floating-point underflow. As a result, the computed values for mu are all zero.
Original output:
(.venv) PS > goodvibes .\pyridazine003_conf-1.log --spc link -t 298.15
GoodVibes v3.2 2025/02/05 15:15:32
Citation: Luchini, G.; Alegre-Requena, J. V.; Funes-Ardoiz, I.; Paton, R. S. F1000Research, 2020, 9, 291.
GoodVibes version 3.2 DOI: 10.12688/f1000research.22758.1
o Requested: --spc link -t 298.15
Temperature = 298.15 Kelvin Pressure = 1 atm
All energetic values below shown in Hartree unless otherwise specified.
o Found vibrational scaling factor of 0.985 for B3LYP/def2TZVP level of theory
I. M. Alecu, unpublished (2011).
Entropic quasi-harmonic treatment: frequency cut-off value of 100.0 wavenumbers will be applied.
QS = Grimme: Using a mixture of RRHO and Free-rotor vibrational entropies.
REF: Grimme, S. Chem. Eur. J. 2012, 18, 9955-9964
Combining final single point energy with thermal corrections.\.venv\Lib\site-packages\goodvibes\thermo.py:343: RuntimeWarning: invalid value encountered in scalar divide
factor = [8 * math.pi ** 3 * entry * BOLTZMANN_CONSTANT * temperature / PLANCK_CONSTANT ** 2 for entry in mu_primed]
Structure E_SPC E ZPE H_SPC T.S T.qh-S G(T)_SPC qh-G(T)_SPC
**********************************************************************************************************************************************
o pyridazine003_conf-1 -910.279625 -910.780414 0.226953 -910.036705 0.061099 nan -910.097805 nan
**********************************************************************************************************************************************
So I check the code at thermo.py and all the mu there become np.float32(nan) due to freq_scale_factor itself being numpy.float32 triggle type conversion in
mu = [PLANCK_CONSTANT / (8 * math.pi ** 2 * freq * SPEED_OF_LIGHT * freq_scale_factor) for freq in frequency_wn]
Function changed:
def calc_freerot_entropy(frequency_wn, temperature, freq_scale_factor, fract_modelsys, file, inertia, roconst):
"""
Free rotor entropy evaluation.
Entropic contributions (J/(mol*K)) according to a free-rotor
description for a list of vibrational modes
Sr = R(1/2 + 1/2ln((8pi^3u'kT/h^2))
Parameters:
frequency_wn (list): list of frequencies parsed from file.
temperature (float): temperature for calculations to be performed at.
freq_scale_factor (float): frequency scaling factor based on level of theory and basis set used.
fract_modelsys (list): MM frequency scale factors obtained from ONIOM calculations.
inertia (str): flag for choosing global average moment of inertia for all molecules or computing individually from parsed rotational constants
roconst (list): list of parsed rotational constants for computing the average moment of inertia.
Returns:
float: free rotor entropy of chemical system.
"""
# This is the average moment of inertia used by Grimme
if inertia == "global" or len(roconst) == 0:
bav = 1.00e-44
else:
av_roconst_ghz = sum(roconst)/len(roconst) #GHz
av_roconst_hz = av_roconst_ghz * 1000000000 #Hz
av_roconst_s = 1 / av_roconst_hz #s
av_roconst = av_roconst_s * PLANCK_CONSTANT #kg m^2
bav = av_roconst
if fract_modelsys is not False:
freq_scale_factor = [freq_scale_factor[0] * fract_modelsys[i] + freq_scale_factor[1] * (1.0 - fract_modelsys[i])
for i in range(len(fract_modelsys))]
mu = [PLANCK_CONSTANT / (8 * math.pi ** 2 * frequency_wn[i] * SPEED_OF_LIGHT * freq_scale_factor[i]) for i in
range(len(frequency_wn))]
else:
mu = [PLANCK_CONSTANT / (8 * math.pi ** 2 * freq * SPEED_OF_LIGHT * freq_scale_factor) for freq in frequency_wn]
print(f"----------------------------------DEBUG----------------------------------")
print(f"frequency_wn: {frequency_wn}")
print(f"freq_scale_factor: {freq_scale_factor}")
print(f"temperature: {temperature}")
print(f"fract_modelsys: {fract_modelsys}")
print(f"inertia: {inertia}")
print(f"roconst: {roconst}")
print("typeof freq_scale_factor: ", type(freq_scale_factor))
print(f"----------------------------------")
print(f"mu: {mu}")
mu_primed = [entry * bav / (entry + bav) for entry in mu]
print(f"----------------------------------")
print(f"mu_primed: {mu_primed}")
factor = [8 * math.pi ** 3 * entry * BOLTZMANN_CONSTANT * temperature / PLANCK_CONSTANT ** 2 for entry in mu_primed]
print(f"----------------------------------")
print(f"factor: {factor}")
entropy = [(0.5 + math.log(entry ** 0.5)) * GAS_CONSTANT for entry in factor]
print(f"----------------------------------")
print(f"entropy: {entropy}")
return entropy
Test output:
(.venv) PS > goodvibes .\pyridazine003_conf-1.log --spc link -t 298.15
GoodVibes v3.2 2025/02/05 15:29:02
Citation: Luchini, G.; Alegre-Requena, J. V.; Funes-Ardoiz, I.; Paton, R. S. F1000Research, 2020, 9, 291.
GoodVibes version 3.2 DOI: 10.12688/f1000research.22758.1
o Requested: --spc link -t 298.15
Temperature = 298.15 Kelvin Pressure = 1 atm
All energetic values below shown in Hartree unless otherwise specified.
o Found vibrational scaling factor of 0.985 for B3LYP/def2TZVP level of theory
I. M. Alecu, unpublished (2011).
Entropic quasi-harmonic treatment: frequency cut-off value of 100.0 wavenumbers will be applied.
QS = Grimme: Using a mixture of RRHO and Free-rotor vibrational entropies.
REF: Grimme, S. Chem. Eur. J. 2012, 18, 9955-9964
Combining final single point energy with thermal corrections.----------------------------------DEBUG----------------------------------
frequency_wn: [26.2781, 33.5112, 52.4551, 62.9903, 93.3022, 154.6547, 171.9286, 191.045, 226.4802, 238.6293, 265.5067, 278.9429, 307.11, 368.9971, 398.8422, 411.7189, 413.8857, 458.0241, 504.1258, 527.2091, 599.4305, 621.3473, 645.8964, 670.7584, 732.1868, 735.0064, 814.075, 816.4233, 841.6155, 854.0938, 880.2988, 903.4555, 939.0067, 959.6354, 971.8484, 991.1861, 998.9826, 1027.3866, 1048.1411, 1063.2691, 1084.5765, 1097.6447, 1123.5917, 1137.9219, 1156.5445, 1175.0828, 1192.7165, 1213.5184, 1253.9492, 1295.3426, 1313.5811, 1331.0286, 1346.5684, 1349.9224, 1396.285, 1397.2969, 1417.9172, 1438.4246, 1461.1311, 1486.4614, 1493.7461, 1501.098, 1514.9533, 1523.7389, 1584.4884, 1607.2943, 1622.7974, 1632.8506, 3023.1284, 3033.4975, 3036.7239, 3094.3194, 3096.8561, 3110.9842, 3115.6735, 3176.4438, 3182.9769, 3196.8559, 3200.6226, 3204.9053, 3214.7447]
freq_scale_factor: 0.9850000143051147
temperature: 298.15
fract_modelsys: False
inertia: global
roconst: [2.00498, 0.09379, 0.09153]
typeof freq_scale_factor: <class 'numpy.float32'>
----------------------------------
mu: [np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0)]
----------------------------------
mu_primed: [np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0), np.float32(0.0)]
\.venv\Lib\site-packages\goodvibes\thermo.py:355: RuntimeWarning: invalid value encountered in scalar divide
factor = [8 * math.pi ** 3 * entry * BOLTZMANN_CONSTANT * temperature / PLANCK_CONSTANT ** 2 for entry in mu_primed]
----------------------------------
factor: [np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan), np.float32(nan)]
----------------------------------
entropy: [nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]
Structure E_SPC E ZPE H_SPC T.S T.qh-S G(T)_SPC qh-G(T)_SPC
**********************************************************************************************************************************************
o pyridazine003_conf-1 -910.279625 -910.780414 0.226953 -910.036705 0.061099 nan -910.097805 nan
**********************************************************************************************************************************************
Convert the freq_scale_factor to float fix the issues.