xtb
xtb copied to clipboard
Published
20 hours ago •
grimme-lab
grimme-lab
Fix imaginary frequencies
Fix #693
Instead of trying to detect how many frequencies should be rot/trans, this patch counts number of low freqs and sort them.
For linear H2O molecule:
3
H 0.0 0.0 0.0
O 0.0 0.0 1.0
H 0.0 0.0 2.0
the diff between old behaviour and new (--hess):
vibrational frequencies (cm⁻¹)
-eigval : 0.00 0.00 0.00 0.00 0.00 935.71
-eigval : 935.71 2867.70 3045.49
+eigval : -1840.78 -1840.78 935.71 935.71 2867.70 3045.49
+eigval : 0.00 0.00 0.00
...
vibrational frequencies (cm⁻¹)
-eigval : 0.00 0.00 0.00 0.00 0.00 935.71
-eigval : 935.71 2867.70 3045.49
+eigval : -1840.78 -1840.78 935.71 935.71 2867.70 3045.49
+eigval : 0.00 0.00 0.00
reduced masses (amu)
- 1: 2.69 2: 2.69 3: 14.32 4: 14.32 5: 14.32 6: 1.01 7: 1.01 8: 1.01
- 9: 2.69
+ 1: 2.69 2: 2.69 3: 1.01 4: 1.01 5: 1.01 6: 2.69 7: 14.32 8: 14.32
+ 9: 14.32
IR intensities (km·mol⁻¹)
- 1:726.33 2:726.33 3: 0.00 4: 0.00 5: 0.00 6: 0.00 7: 0.00 8: 0.00
- 9:217.07
+ 1:726.33 2:726.33 3: 0.00 4: 0.00 5: 0.00 6:217.07 7: 0.00 8: 0.00
+ 9: 0.00
Raman intensities (Ä⁴*amu⁻¹)
1: 0.00 2: 0.00 3: 0.00 4: 0.00 5: 0.00 6: 0.00 7: 0.00 8: 0.00
9: 0.00
...
: SETUP :
:.................................................:
: # frequencies 4 :
- : # imaginary freq. 0 :
+ : # imaginary freq. 2 :
: linear? true :
: only rotor calc. false :
: symmetry din :
...
mode ω/cm⁻¹ T·S(HO)/kcal·mol⁻¹ T·S(FR)/kcal·mol⁻¹ T·S(vib)
------------------------------------------------------------------------
+ 1 935.71 -0.03610 (100.00%) -0.18483 ( 0.00%) -0.03611
+ 2 935.71 -0.03610 (100.00%) -0.18483 ( 0.00%) -0.03611
+ 3 2867.70 -0.00001 (100.00%) 0.14430 ( 0.00%) -0.00001
+ 4 3045.49 -0.00000 (100.00%) 0.16204 ( 0.00%) -0.00000
------------------------------------------------------------------------
...
:::::::::::::::::::::::::::::::::::::::::::::::::::::
:: THERMODYNAMIC ::
:::::::::::::::::::::::::::::::::::::::::::::::::::::
- :: total free energy -5.009429056662 Eh ::
+ :: total free energy -5.009429056672 Eh ::
::.................................................::
:: total energy -5.010690547588 Eh ::
- :: zero point energy 0.017734665620 Eh ::
+ :: zero point energy 0.017734665611 Eh ::
:: G(RRHO) w/o ZPVE -0.016473174695 Eh ::
- :: G(RRHO) contrib. 0.001261490926 Eh ::
+ :: G(RRHO) contrib. 0.001261490916 Eh ::
:::::::::::::::::::::::::::::::::::::::::::::::::::::
+imag cut-off (cm-1) : 5.00
+ found 2 significant imaginary frequencies
+ writing imag mode distorted coords to xtbhess.xyz
+ for further optimization.
CO.xyz:
2
O 1.128 0.0 0.0
C 0.0 0.0 0.0
The diff (running with --ohess):
@@ -402,7 +402,7 @@ RMS gradient : 0.00032
writing file <hessian>, containing the non-mass-weighted Hessian matrix in atomic units (Eₕ/Bohr²).
vibrational frequencies (cm⁻¹)
-eigval : 0.00 0.00 0.00 0.00 0.00 2223.01
+eigval : -5.74 -5.74 2223.01 0.00 0.00 0.00
-------------------------------------------------
| Property Printout |
-------------------------------------------------
@@ -480,11 +480,11 @@ It seems to be the Cinfv point group
| Frequency Printout |
-------------------------------------------------
vibrational frequencies (cm⁻¹)
-eigval : 0.00 0.00 0.00 0.00 0.00 2223.01
+eigval : -5.74 -5.74 2223.01 0.00 0.00 0.00
reduced masses (amu)
- 1: 13.72 2: 13.72 3: 14.29 4: 14.29 5: 14.29 6: 13.72
+ 1: 13.72 2: 13.72 3: 13.72 4: 14.29 5: 14.29 6: 14.29
IR intensities (km·mol⁻¹)
- 1: 1.84 2: 1.84 3: 0.00 4: 0.00 5: 0.00 6: 26.98
+ 1: 1.84 2: 1.84 3: 26.98 4: 0.00 5: 0.00 6: 0.00
Raman intensities (Ä⁴*amu⁻¹)
1: 0.00 2: 0.00 3: 0.00 4: 0.00 5: 0.00 6: 0.00
output can be read by thermo (or use thermo option).
@@ -495,11 +495,13 @@ eigval : 0.00 0.00 0.00
-------------------------------------------------
cin symmetry found (for desy threshold: 0.10E+00) used in thermo
+ inverting freq 1 5.7385921955661043
+ inverting freq 2 5.7384824331576896
...................................................
: SETUP :
:.................................................:
- : # frequencies 1 :
+ : # frequencies 3 :
: # imaginary freq. 0 :
: linear? true :
: only rotor calc. false :
@@ -512,39 +514,46 @@ cin symmetry found (for desy threshold:
mode ω/cm⁻¹ T·S(HO)/kcal·mol⁻¹ T·S(FR)/kcal·mol⁻¹ T·S(vib)
------------------------------------------------------------------------
+ 1 5.74 -2.71750 ( 0.02%) -1.57666 ( 99.98%) -1.57686
+ 2 5.74 -2.71751 ( 0.02%) -1.57667 ( 99.98%) -1.57686
+ 3 2223.01 -0.00015 (100.00%) 0.06796 ( 0.00%) -0.00015
------------------------------------------------------------------------
temp. (K) partition function enthalpy heat capacity entropy
cal/mol cal/K/mol cal/K/mol J/K/mol
- 298.15 VIB 1.00 0.139 0.005 0.001
+ 298.15 VIB 0.134E+04 1168.810 1.993 10.578
ROT 107. 592.502 1.987 11.276
- INT 107. 592.641 1.992 11.277
+ INT 0.144E+06 1761.312 3.980 21.855
TR 0.143E+27 1481.254 4.968 35.908
- TOT 2073.8949 6.9604 47.1852 197.4228
+ TOT 3242.5660 8.9480 57.7631 241.6810
T/K H(0)-H(T)+PV H(T)/Eh T*S/Eh G(T)/Eh
------------------------------------------------------------------------
- 298.15 0.330496E-02 0.836936E-02 0.224192E-01 -0.140498E-01
+ 298.15 0.516736E-02 0.102579E-01 0.274451E-01 -0.171872E-01
------------------------------------------------------------------------
:::::::::::::::::::::::::::::::::::::::::::::::::::::
:: THERMODYNAMIC ::
:::::::::::::::::::::::::::::::::::::::::::::::::::::
- :: total free energy -6.135900073382 Eh ::
+ :: total free energy -6.139037456355 Eh ::
::.................................................::
:: total energy -6.121850225907 Eh ::
- :: zero point energy 0.005064394619 Eh ::
- :: G(RRHO) w/o ZPVE -0.019114242094 Eh ::
- :: G(RRHO) contrib. -0.014049847475 Eh ::
+ :: zero point energy 0.005090541319 Eh ::
+ :: G(RRHO) w/o ZPVE -0.022277771767 Eh ::
+ :: G(RRHO) contrib. -0.017187230448 Eh ::
:::::::::::::::::::::::::::::::::::::::::::::::::::::
+imag cut-off (cm-1) : 5.00
+ found 2 significant imaginary frequencies
+ writing imag mode distorted coords to xtbhess.xyz
+ for further optimization.
For bigger system:
-------------------------------------------------
| Frequency Printout |
-------------------------------------------------
projected vibrational frequencies (cm⁻¹)
-eigval : -0.00 -0.00 -0.00 -0.00 0.00 0.00
-eigval : 53.81 134.56 182.61 284.67 371.44 655.69
-eigval : 709.53 796.89 888.75 941.11 994.10 1024.29
+eigval : 53.81 134.55 182.61 284.67 371.44 655.69
+eigval : 709.53 796.89 888.74 941.11 994.10 1024.29
eigval : 1143.84 1191.69 1270.75 1337.11 1346.31 1371.68
-eigval : 1785.08 1814.93 2775.63 2855.93 2911.01 2978.60
+eigval : 1785.08 1814.93 2775.64 2855.93 2911.01 2978.60
+eigval : 0.00 0.00 0.00 0.00 0.00 0.00
reduced masses (amu)
- 1: 12.58 2: 13.67 3: 13.65 4: 13.11 5: 14.05 6: 13.75 7: 14.58 8: 9.29
- 9: 12.26 10: 10.39 11: 12.03 12: 8.22 13: 9.23 14: 11.13 15: 10.44 16: 5.72
- 17: 8.68 18: 8.53 19: 11.85 20: 5.45 21: 4.21 22: 2.75 23: 2.82 24: 2.00
- 25: 13.42 26: 13.48 27: 1.70 28: 1.75 29: 1.73 30: 1.73
+ 1: 14.58 2: 9.29 3: 12.26 4: 10.39 5: 12.03 6: 8.22 7: 9.23 8: 11.13
+ 9: 10.44 10: 5.72 11: 8.68 12: 8.53 13: 11.85 14: 5.45 15: 4.21 16: 2.75
+ 17: 2.82 18: 2.00 19: 13.42 20: 13.48 21: 1.70 22: 1.75 23: 1.73 24: 1.73
+ 25: 12.25 26: 12.82 27: 13.75 28: 13.97 29: 14.59 30: 13.43
IR intensities (km·mol⁻¹)
- 1: 2.37 2: 1.24 3: 0.01 4: 3.95 5: 1.04 6: 1.92 7: 0.87 8: 8.59
- 9: 11.13 10: 20.62 11: 14.21 12: 7.87 13: 28.63 14: 31.72 15: 19.04 16: 13.81
- 17: 35.02 18: 4.17 19:170.86 20: 47.05 21: 14.71 22: 32.04 23: 15.06 24: 25.25
- 25:427.51 26:280.50 27:108.09 28: 74.70 29: 7.14 30: 2.52
+ 1: 0.87 2: 8.59 3: 11.13 4: 20.62 5: 14.21 6: 7.87 7: 28.63 8: 31.72
+ 9: 19.04 10: 13.81 11: 35.02 12: 4.17 13:170.86 14: 47.04 15: 14.71 16: 32.04
+ 17: 15.06 18: 25.25 19:427.51 20:280.50 21:108.09 22: 74.70 23: 7.14 24: 2.52
+ 25: 2.93 26: 1.59 27: 1.70 28: 0.23 29: 2.54 30: 1.52
Raman intensities (Ä⁴*amu⁻¹)
1: 0.00 2: 0.00 3: 0.00 4: 0.00 5: 0.00 6: 0.00 7: 0.00 8: 0.00
9: 0.00 10: 0.00 11: 0.00 12: 0.00 13: 0.00 14: 0.00 15: 0.00 16: 0.00
@@ -989,69 +989,84 @@ c1 symmetry found (for desy threshold:
mode ω/cm⁻¹ T·S(HO)/kcal·mol⁻¹ T·S(FR)/kcal·mol⁻¹ T·S(vib)
------------------------------------------------------------------------
- 1 53.81 -1.39307 ( 57.28%) -1.03419 ( 42.72%) -1.23977
- 2 134.56 -0.85863 ( 98.13%) -0.76303 ( 1.87%) -0.85684
+ 1 53.81 -1.39306 ( 57.28%) -1.03419 ( 42.72%) -1.23977
+ 2 134.55 -0.85864 ( 98.13%) -0.76303 ( 1.87%) -0.85685
3 182.61 -0.68621 ( 99.44%) -0.67264 ( 0.56%) -0.68614
- 4 284.67 -0.44885 ( 99.90%) -0.54118 ( 0.10%) -0.44894
+ 4 284.67 -0.44886 ( 99.90%) -0.54118 ( 0.10%) -0.44894
+ 5 371.44 -0.32017 ( 99.97%) -0.46239 ( 0.03%) -0.32021
+ 6 655.69 -0.10828 (100.00%) -0.29407 ( 0.00%) -0.10828
+ 7 709.53 -0.08795 (100.00%) -0.27070 ( 0.00%) -0.08795
+ 8 796.89 -0.06257 (100.00%) -0.23630 ( 0.00%) -0.06257
+ 9 888.74 -0.04354 (100.00%) -0.20399 ( 0.00%) -0.04354
+ 10 941.11 -0.03533 (100.00%) -0.18703 ( 0.00%) -0.03533
+ 11 994.10 -0.02856 (100.00%) -0.17081 ( 0.00%) -0.02856
+ 12 1024.29 -0.02528 (100.00%) -0.16195 ( 0.00%) -0.02529
+ 13 1143.84 -0.01553 (100.00%) -0.12925 ( 0.00%) -0.01554
+ 14 1191.69 -0.01276 (100.00%) -0.11711 ( 0.00%) -0.01276
+ 15 1270.75 -0.00920 (100.00%) -0.09808 ( 0.00%) -0.00920
+ 16 1337.11 -0.00697 (100.00%) -0.08300 ( 0.00%) -0.00697
+ 17 1346.31 -0.00671 (100.00%) -0.08097 ( 0.00%) -0.00671
+ 18 1371.68 -0.00603 (100.00%) -0.07544 ( 0.00%) -0.00603
+ 19 1785.08 -0.00103 (100.00%) 0.00259 ( 0.00%) -0.00103
+ 20 1814.93 -0.00091 (100.00%) 0.00750 ( 0.00%) -0.00091
+ 21 2775.64 -0.00001 (100.00%) 0.13335 ( 0.00%) -0.00001
+ 22 2855.93 -0.00001 (100.00%) 0.14180 ( 0.00%) -0.00001
+ 23 2911.01 -0.00001 (100.00%) 0.14746 ( 0.00%) -0.00001
+ 24 2978.60 -0.00001 (100.00%) 0.15426 ( 0.00%) -0.00001
------------------------------------------------------------------------
temp. (K) partition function enthalpy heat capacity entropy
cal/mol cal/K/mol cal/K/mol J/K/mol
- 298.15 VIB 29.2 2157.784 13.674 13.425
+ 298.15 VIB 29.2 2157.786 13.674 13.425
ROT 0.123E+06 888.752 2.981 26.273
- INT 0.360E+07 3046.537 16.655 39.698
+ INT 0.360E+07 3046.538 16.655 39.698
TR 0.800E+27 1481.254 4.968 39.321
- TOT 4527.7906 21.6227 79.0193 330.6168
+ TOT 4527.7919 21.6227 79.0193 330.6169
T/K H(0)-H(T)+PV H(T)/Eh T*S/Eh G(T)/Eh
------------------------------------------------------------------------
- 298.15 0.721549E-02 0.751504E-01 0.375446E-01 0.376058E-01
+ 298.15 0.721550E-02 0.751504E-01 0.375446E-01 0.376058E-01
------------------------------------------------------------------------
:::::::::::::::::::::::::::::::::::::::::::::::::::::
:: THERMODYNAMIC ::
:::::::::::::::::::::::::::::::::::::::::::::::::::::
- :: total free energy -20.585323405221 Eh ::
+ :: total free energy -20.585323411084 Eh ::
::.................................................::
- :: total energy -20.622929187851 Eh ::
- :: zero point energy 0.067934918791 Eh ::
- :: G(RRHO) w/o ZPVE -0.030329136161 Eh ::
- :: G(RRHO) contrib. 0.037605782630 Eh ::
+ :: total energy -20.622929187960 Eh ::
+ :: zero point energy 0.067934920100 Eh ::
+ :: G(RRHO) w/o ZPVE -0.030329143225 Eh ::
+ :: G(RRHO) contrib. 0.037605776876 Eh ::
:::::::::::::::::::::::::::::::::::::::::::::::::::::
I see the problem for linear water, where you want to have two large imaginary modes ( -1840.78 cm-1) and the smaller modes as the rotation (935.71 cm-1). Here, xtb assigns them the wrong way around. But I don't think that not assigning translation/rotation is better here. For example, the relaxed CO2 geometry
3
energy: -10.308452263543 gnorm: 0.000250996772 xtb: 6.7.1 (0a1c9e1)
O -0.00000000000000 0.00000000000000 -0.14372775733893
C 0.00000000000000 -0.00000000000000 0.99999999982051
O -0.00000000000000 0.00000000000000 2.14372775751843
yields the following frequencies:
6.7.1:
< eigval : 0.00 0.00 0.00 0.00 0.00 600.60
< eigval : 600.60 1424.42 2592.42
This PR:
> eigval : -0.12 -0.12 -0.09 15.64 15.64 600.60
> eigval : 600.60 1424.40 2592.37
However, the first five frequencies really should be zero, since they are translation and rotation where the Hessian calculation applies a (wrong) harmonic approximation to the degree of freedom. The same holds for the linear water molecule, where five frequencies must again be zero. For the user, it is highly unintuitive if xtb provides the wrong number of non-zero vibrational frequencies. Therefore, we must assign linear and bent molecules to find how many frequencies should be zero.
How we decide, which frequencies are translations/rotations is a different question. Maybe, we should always take the five/six frequencies with the lowest magnitude, but I am unsure whether this is always correct (not sure how other programs do this).
As a side note, I think we should retain the ordering of the frequencies, where the imaginary modes come first followed by the translation/rotation and the positive frequencies. There are likely many scripts based on this ordering, which would immediately break :D
