Freq scaling factors from ORCA do not transfer

[ffmulks]

Hello again!

I tried using PBEh-3c as optimisation method and that seems to work great so far. They suggest, however, to use a scaling factor of 0.95 for computing thermodynamics. I added the following bit to the hess-keywords which works fine.

%freq
scalfreq 0.95
end

ORCA 5.0.3 produces the following after autodE called for the calculation.

-----------------------
VIBRATIONAL FREQUENCIES
-----------------------

Scaling factor for frequencies =  0.950000000  (already applied!)

   0:         0.00 cm**-1
   1:         0.00 cm**-1
   2:         0.00 cm**-1
   3:         0.00 cm**-1
   4:         0.00 cm**-1
   5:         0.00 cm**-1
   6:      1620.46 cm**-1
   7:      3708.47 cm**-1
   8:      3815.22 cm**-1

Checking the frequencies in autodE, however, gives:

if molecule.imaginary_frequencies:
    frequencies = molecule.imaginary_frequencies + molecule.vib_frequencies
else:
    frequencies = molecule.vib_frequencies
print(*[f'{(freq).to("cm-1"):.1f}' for freq in frequencies],sep=', ')

Frequencies in cm-1
1705.8, 3903.8, 4016.2

This is weird, because I cannot find the numbers at all in ORCA's output files. Does autodE re-compute the frequencies straight from the found normal modes? Both in the .hess and the .out, ORCA only stores the scaled frequencies.

Do I assume correctly that also the thermodynamics did not get the scaling factor applied? My test computing water several times seems to support that. The deviations appear to be in our margin of error (full cycles with conformer search and opt were run each time).

For H2O:

No scaling
E_opt, G_cont, H_cont, E_sp
-76.267818102248,0.0073088101511278696,0.0257081115702916,-76.385195142966

Strong scaling (0.5)
E_opt, G_cont, H_cont, E_sp
-76.26781810866,0.007308818432288296,0.025708124156216056,-76.385195140162
-76.26781810866,0.0073088378537382,0.02570814342764412,-76.385195140162
-76.267818106119,0.007308850889597047,0.025708153656751487,-76.385195141268
-76.267818102248,0.007308814161775101,0.02570811555037913,-76.385195142966
-76.26781810866,0.007308829535297789,0.025708135175725795,-76.385195140162

Is there a good way to transfer the scaling factor then (and propagate it to the stored H_cont and G_cont)? I could surely overwrite the frequencies with the scaled ones but how would I go best about re-computing H_cont and G_cont without triggering new frequency calculations?

Or better: can you add a snippet to ade.Config to define a scaling factor that is given before the thermodynamics are computed?

[ffmulks]

I fixed it for my specific use case by adding the *0.95 down there in the return line (in hessians.py) That is obviously not a useful solution because it required re-building from source :D

def _eigenvalues_to_freqs(lambdas) -> List[Frequency]:
    """
    Convert eigenvalues of the Hessian matrix (SI units) to
    frequencies in wavenumber units

    -----------------------------------------------------------------------
    Arguments:
        lambdas (np.ndarray):

    Returns:
        (list(autode.values.Frequency)):
    """

    nus = (np.sqrt(np.complex_(lambdas))
           / (2.0 * np.pi * Constants.ang_to_m * Constants.c_in_cm))

    # Cast the purely complex eigenvalues to negative real numbers, as is
    # usual in quantum chemistry codes
    nus[np.iscomplex(nus)] = -np.abs(nus[np.iscomplex(nus)])

    return [Frequency(np.real(nu)*0.95, units=wavenumber) for nu in nus]

With the desired 0.95 scaling this gave -76.267818109491,0.006212212664853274,0.024612544739443177,-76.385195140062

[t-young31]

Hi!

As you found out autodE does its own hessian diagonalisation to get the frequencies from the Hessian matrix and unfortunately doesn't have any way of applying a frequency scaling factor. It will however be pretty simple to implement! Is there a table somewhere of scaling factors for different functionals?

---

Notes for me:

• Config attr
• ORCA warning for unused scalfreq kwd
• Freq scale factor in functional defn (override by Config)

[ffmulks]

Hey,

this is the broadest one that I am aware of. They use the weird Gaussian notation but I suppose they cover your standard PBE0 method.

https://cccbdb.nist.gov/vibscalejust.asp

I am unsure if this is a good idea, but you may improve the calculation times by a couple of seconds with switching ORCA's internal qRRHO off and things like that when you do it yourself anyways. But honestly, it is probably better to leave it in so that people have actually useful output files in case they want to do their own math with the output files.

It is fair to say that the basis set dependence is low, so you may consider for unknown functional/BS combinations to set a functional standard. I would personally round up in order to not overcorrect unknown errors, e.g. 0.97 would be a decent guess for any B3LYP calculation. 0.96 for PBE0. However, it is easier to stick to known combinations because this generalisation surely requires some solid foundation in the methods description and a bunch of references. I would implement the generalisation only when time and resources to write and reference a solid background are available.

[ffmulks]

Tested and works well! Looking forward to the 1.3.0 release. Want this closed or after the release?

[t-young31]

Awesome – thanks @ffmulks! Let's leave this open until the release is out (hopefully soon 🤞🏼)