Fix imaginary frequencies #1210
Conversation
There was a problem hiding this comment.
Thank you for your contribution!
This looks good to me, I played around with a few systems and the results are good.
I am unsure, however, if this could have further impacts that I am possibly unaware of, @thfroitzheim?
Signed-off-by: Igor S. Gerasimov <[email protected]>
Signed-off-by: Igor S. Gerasimov <[email protected]>
Signed-off-by: Igor S. Gerasimov <[email protected]>
Signed-off-by: Igor S. Gerasimov <[email protected]>
Signed-off-by: Igor S. Gerasimov <[email protected]>
Signed-off-by: Igor S. Gerasimov <[email protected]>
Signed-off-by: Igor S. Gerasimov <[email protected]>
|
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 yields the following frequencies: 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 |
|
I think I can use detrotra8 procedure to set rot/tra freqs properly. |
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:
the diff between old behaviour and new (
--hess):CO.xyz:
The diff (running with
--ohess):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 :: :::::::::::::::::::::::::::::::::::::::::::::::::::::