Projecting out rotors in Arkane leads to zero frequency

[sevyharris]

Bug Description

I am trying to compute the thermodynamics for a species (smiles='CCCCOO') by running an Arkane statmech + thermo job, but there is a very large difference between my entropy and the reference data I'm using from the AramcoMech3.0, leading to a 27 kcal/mol error in Gibbs free energy at 1000K.

This led me to look at the rotors. I ran 1-D hindered rotor scans (at 20 degree increments) for the molecule and they looked reasonable in Gaussview, but @rwest noticed that after Arkane projects out the rotors, one of the vibrational frequencies is approximately zero (1e-5).

HarmonicOscillator(
        frequencies = ([4.79104e-05, 424.798, 521.468, 584.421, 778.857, 820.924, 923.676, 952.077, 1015.85, 1044.41, 1084.27, 1107.11, 1156.83, 1241.45, 1259.22, 1294.35, 1298.53, 1335.14, 1346.66, 1375.05, 1392.76, 1419.19, 1430.83, 1464.16, 1477.48, 1501.56, 1553.74, 2611.14, 2672.06, 2713.94, 2729.57, 2761.52, 2804.99, 2828.13, 2855.04, 2875.21, 2981.48], 'cm^-1'),
)

You can see all the statmech modes in the full output.py

I'm trying to figure out whether 1. I've messed up my input to Arkane, 2. Arkane is projecting the rotors out incorrectly, or 3. the rotation and vibrational modes are too tightly coupled for this kind of calculation to work. This is known to be a problem for OOQOOH radicals but my species is a little simpler. I have successfully calculated dozens of species thermodynamics with my current workflow, but I also see this rotors/zero-frequency issue with a handful of other species: (smiles='CC(CC=O)OO', 'CC(CCOO)O[O]', 'CC(CCO[O])OO', and 'O=CCCCOO').

I would appreciate any and all hints for how to tackle debugging this issue.

How To Reproduce

These are the files I'm using (.py extensions changed to .txt). For simplicity, I have excluded my single point energy calculation (used to generate the plots above), and energy, frequencies and geometry are all taken from the conformer_0002.log file.

Arkane Input Files:

input.py
conformer_0002.py

Gaussian Log Files

conformer_0002.log
rotor_0000.log
rotor_0001.log
rotor_0002.log
rotor_0003.log
rotor_0004.log

Expected Behavior

I would expect the lowest frequency to be about 100 cm^-1 or larger after projecting out the hindered rotors, instead of zero.

Installation Information

I'm on my experimental my_bacs branch of RMG-database and the latest RMG-Py.

[mjohnson541]

1. From statmech: lim v_i -> 0 q_vib_i(v_i) = Inf and lim q_vib_i -> Inf S(q_vib_i) = Inf.

2. Near zero frequency implies near zero eigenvalue and that along the associated eigenvector there is little or no barrier to moving from the well in that direction.

3. A Hessian at a well is positive definite which implies that all of the eigenvalues are positive and it cannot have a zero eigenvalue motion like that

4. Projecting out motions from a positive definite Hessian should never make it not positive definite so it shouldn't be possible to create a "flat" motion if the original Hessian doesn't have it

Some possibilities:

a) The "flat" motion exists in the original Hessian => lowest frequency in gaussian file is 62 cm^-1 so no

b) The "flat" motion observed here is actually a projected out motion that hasn't been full removed yet in the code (either a bug or a result of the debugging point) => 42 freq before - 5 rotors = 37 freq after so no

c) Something in the input rotor information is messing up rotor projection...I don't see anything obvious to me, but I would double check your pivots and top indexing against an example you know is right, before digging into the code

d) Bug in rotor projection...I feel like there have been a couple times we thought there was a bug here before and it turned out we were wrong...I would very carefully check the inputs to the project_rotors function before trying to dig into the details...

[sevyharris]

Thanks @mjohnson541 for these possibilities!

I don't think the atom numbering is a problem.

I tried to break it down into a smaller example, including just a single free rotor for the OH at the end. The numbering for the tops/pivots is about as simple an example as I could come up with (you can see the atom numbering in the screenshot) and it's not even reading the Gaussian rotor scan file. This leads to three zero frequencies after projecting out the single free rotor, which makes me wonder if the project_rotors function is capable of doing only a subset of the full number of rotors.

My experiments so far have shown that including any one of the rotors for this species leads to at least one zero frequency after projecting the rotors out.

spinMultiplicity = 1

energy = {'M06-2X/cc-pVTZ': Log('conformer_0002.log')}
geometry = Log('conformer_0002.log')

frequencies = Log('conformer_0002.log')

rotors = [
     FreeRotor(pivots=[1, 2], top=[2, 16], symmetry=1),
]

Output:

conformer(
    label = 'C4H10O2',
    E0 = (
        -226.736,
        'kJ/mol',
    ),
    modes = [
        IdealGasTranslation(mass=(90.0681, 'amu')),
        NonlinearRotor(
            inertia = ([83.9394, 242.509, 246.662], 'amu*angstrom^2'),
            symmetry = 1,
        ),
        HarmonicOscillator(
            frequencies = ([1.72094e-05, 1.73751e-05, 3.94123e-05, 492.042, 528.128, 557.903, 588.439, 733.637, 819.461, 877.679, 933.604, 964.708, 1037.86, 1065.64, 1092.01, 1115.61, 1159.76, 1251.67, 1269.09, 1298.71, 1302.19, 1340.69, 1355.65, 1382.61, 1398.4, 1421.2, 1432.78, 1466.02, 1478.26, 1508.59, 1557.52, 2711.73, 2716.45, 2743.62, 2762.49, 2796.45, 2817.41, 2830.6, 2860.77, 2876.46, 2995.24], 'cm^-1'),
        ),
        FreeRotor(inertia=(0.355768, 'amu*angstrom^2'), symmetry=1),
    ],
    spin_multiplicity = 1,
    optical_isomers = 2,
)

[sevyharris]

Here are the plots of the rotor fits. I had to specify "best" fit to generate the plots, but I'm using the Fourier fit for my calculation.

[mjohnson541]

If the issue is with every rotor, I would check the Hessian, particularly what Hessian is being fed into project_rotors. In particular, what its eigenvalues and/or frequencies are.

[sevyharris]

Not sure if this is useful, but I made a set of diagnostic plots. These are the frequencies before and after the projection out:

The frequencies overlap much better on other example species (title is the species smiles):

And I can find other OO species that have bad before/after matches, but don't have zero frequencies:

[rwest]

I've spent some time in the statmech code....
main...rwest:RMG-Py:arkanebug

I haven't solved anything yet, but did find that when I replaced the previous Gram-Schmidt orthonormalization with a QR decomposition from numpy, I found column 47 of the QR decomposition is nearly zero, so we could lose a basis.
(A small diagonal element of R during a QR decomposition signals that the corresponding column in the original matrix is almost entirely explained by the previous columns, implying near-linear dependence and a potential loss of full rank in the matrix.)
But I don't know the cause or consequence.

This is in the project_rotors function, so only gets called when doing rotors, but happens before we actually reach the rotors: at that point it is just projecting out the translation and external rotation. Could those be wrong somehow? Or something weird about the hessian in the gaussian log file?

[mjohnson541]

Yeah, that should imply a near zero frequency/eigenvalue. Gram-Schmidt is an odd algorithm choice for orthonomalization since its numerically unstable, so it would probably not be a bad thing to change it.

Could it be related to this old issue? or something similar? #981 :
"""
The fundamental issue is still that sometimes the molecule geometry and Hessian are input in different Cartesian coordinates. The Hessian that Cantherm reads (and that is output by the iop(7/33=1) Gaussian option) has the header "Force constants in Cartesian coordinates: " in the Gaussian log file. However, this description does not say which Cartesian coordinate system. After talking with Gaussian tech support, they confirmed that with the iop(7/33=1) keyword, that Hessian is always defined in the Input Orientation, i.e., in whatever coordinate system the molecule was in the input, not the Standard orientation. This is true whether the user specifies "nosymm" or not. In fact, if you don't specify nosymm, at a certain point in the Gaussian log it will say: "Axes restored to original set". Above this line, everything is defined in Standard orientation, but below this line everything (including the Hessian) gets rotated back into the Input orientation. So clearly, the molecule geometry that Cantherm reads in should be consistent with the Input Orientation coordinate system of the freq job. However, currently Cantherm does not check if this is the case and so it is entirely up to the user to input a geometry and Hessian in consistent coordinates.
"""

[ttpekkan]

If the inputted Cartesians and Hessians are inconsistent, I think the projected frequencies should disagree with the ones found in the Gaussian output file even when no internal rotors are included. Is this the case here?

As an aside, MESMER has a similar projection method implemented as Arkane, and to ensure I'm using consistent Cartesians and Hessians (obtained with Gaussian), I take them from an .fchk file.

[sevyharris]

If the inputted Cartesians and Hessians are inconsistent, I think the projected frequencies should disagree with the ones found in the Gaussian output file even when no internal rotors are included. Is this the case here?

As an aside, MESMER has a similar projection method implemented as Arkane, and to ensure I'm using consistent Cartesians and Hessians (obtained with Gaussian), I take them from an .fchk file.

That might be it!

Here's the result when I call the rotor projection function with no rotors:

And here's a different species after calling the rotor projection function with no rotors:

[sevyharris]

Per one of @rwest 's suggestions, I reran the optimization with IOP(2/9=2000),

  | Do print the Cartesian coordinates in the input orientation. -- | --

This produces the correct matchup in frequencies after projecting out with no rotors

And here it is after projecting out with rotors included!

[sevyharris]

Don't want to celebrate too early, but it looks like it's the cartesians being inconsistent with the Hessian. Thank you @mjohnson541 Timo and @rwest for all your help debugging this!

[sevyharris]

The default in Gaussian is to print the Cartesian coordinates in the input orientation only for small cases, which explains why it was working for many species and then failing on bigger molecules. I will add a warning to Arkane's statmech.py either as a comment or in code to check for iop(2/9=2000) for larger molecules (or in general).