CBS optimization didn't converge and we got imaginary frequencies

izxle · May 30, 2024, 10:03pm

Hi all,

We’re interested in running a CBS calculation with CCSD(T) DT coupled-cluster correction atop a TQ MP2 extrapolation geometry optimization. Followed by a vibrational analysis.

We tested it with glycine’s most stable conformer, to compare with CCSD(T) results [ref], with this input file:

import psi4
psi4.set_memory('100 GB')
psi4.set_num_threads(32)

mol = psi4.geometry("""
  C   -0.53563417153848      0.11199726133568     -0.00218270919096
  N   1.96831088448357      0.01407984608721     -0.00463200595537
  H   2.00194456014429      0.60779655219457     -0.82701948797671
  H   1.99678824661351      0.64509640477476      0.78972730649050
  C   0.72799222755498     -0.72883454326320      0.00883469941762
  O   -1.64638748186916     -0.66599224513834      0.00123189653112
  O   -0.57631328419839      1.31773503135563     -0.01129121640155
  H   0.68464763461558     -1.40639931931541     -0.85043325500567
  H   0.68643497927198     -1.37685378837988      0.89095623842695
  H   -2.40923559507788     -0.06493519965101     -0.00379646633592
""")

psi4.optimize("mp2/aug-cc-pv[tq]z + d:ccsd(t)/aug-cc-pvdz")

psi4.frequency("mp2/aug-cc-pv[tq]z + d:ccsd(t)/aug-cc-pvdz")

The problem is that the optimization stoped at 5 iterations where the energy didn’t converge. We also get a couple of warnings from the vibrational analysis because of this.

$ grep '~' gly.out
             1    -284.36311131   -2.84e+02      5.53e-03      1.80e-03 o    1.04e-02      5.32e-03 o  ~
             2    -284.36325211   -1.41e-04      7.32e-04      2.17e-04 o    1.42e-02      5.01e-03 o  ~
             3    -284.36323893    1.32e-05      5.26e-04      2.40e-04 o    1.47e-02      5.54e-03 o  ~
             4    -284.36325114   -1.22e-05      8.90e-04      2.82e-04 o    4.56e-03      1.53e-03 o  ~
             5    -284.36325689   -5.75e-06      1.04e-04 *    2.38e-05 o    1.09e-03 *    5.09e-04 o  ~
$ tail -3 gly.out
Optimizer: Optimization complete!
Warning: thermodynamics relations excluded imaginary frequencies: ['47.7412i']
 Warning: used thermodynamics relations inappropriate for low-frequency modes: ['197.1585' '250.1161' '460.5099' '515.6141' '613.9087']

Is it that we’re using some default input parameter that we should be changing or maybe we’re missing an option?

We are new to psi4 so we would greatly appreciate the community’s help in resolving this.

jmisiewicz · May 31, 2024, 11:49am

My first thought is to try tightening the geometry convergence. gau_tight and gau_verytight are my go-tos.

If frequencies are at all suspicious, it’s usually good to write the molden file for visualization.

izxle · May 31, 2024, 11:22pm

Thank you for your suggestion. I want to try it but unfortunately, I don’t understand how to use those options from the python API and I didn’t find a way to do it in the documentation. Could you point me in the right direction?

Also, why would it be that the optimization finished without meeting the convergence criteria? As I showed, it only did 5 optimization steps.

jmisiewicz · June 1, 2024, 4:04am

Use set_options, as described in the PythonAPI tutorial.
As described in the documentation on convergence criteria, " Convergence achieved when Max Force and one of Max Energy or Max Disp are fulfilled. … For ultimate control, specifying a value for any of the five monitored options activates that criterium and overwrites/appends it to the criteria set by G_CONVERGENCE. Note that this revokes the special convergence arrangements detailed in notes [5] and [6] and instead requires all active criteria to be fulfilled to achieve convergence."

izxle · June 17, 2024, 11:26pm

Thank you for your answer.

I was able to change the convergence criteria to GAU_VERYTIGHT, but unfortunately the geometry didn’t converge in 50 iterations so I changed it to GAU_TIGHT. This did converge but with one imaginary frequency. So I think I might need an intermediate convergence criteria or changing some other parameter.

I have a couple of questions that I wasn’t able to find how to change from the documentation, could you help me pointing to the keywords to modify this?

How to tighten the tolerance for the energy calculation during the SCF
How to increase the number of optimization iterations
How to restart a calculation, or get the last geometry because I see several on the output and don’t know which one is the correct one.

Thank you in advance for your help.

jmisiewicz · June 18, 2024, 1:10pm

Changing d_convergence is your best bet. You can see the list of keywords associated with any module at the top of its page in the documentation.
geom_maxiter. This is in the geometry optimization page of the manual.
The last geometry in the output is probably the last geometry? I’m not sure where the confusion is coming from.

It looks like this system has a symmetry plane. What is the symmetry of the vibrational mode? It should be either A or A'. If it’s A', tell me immediately.

If you’re going to do more frequency computations, I recommend you use normal_modes_write.

izxle · June 19, 2024, 11:20pm

I tried using the GAU_VERYTIGHT option for G_CONVERGENCE but it didn’t converge. Looking at the values, it seemed like it was oscillating.

            37    -284.36325702    1.05e-09 o    1.38e-06 *    3.97e-07 *    1.86e-04      9.35e-05    ~
            38    -284.36325702   -1.49e-09 o    9.31e-07 *    3.72e-07 *    2.20e-04      1.03e-04    ~
            39    -284.36325702    5.07e-10 o    1.20e-06 *    4.35e-07 *    8.78e-05      4.64e-05    ~
            40    -284.36325702    7.19e-10 o    1.21e-06 *    3.17e-07 *    8.57e-05      4.27e-05    ~
            41    -284.36325702   -8.24e-10 o    1.24e-06 *    3.04e-07 *    9.51e-05      4.35e-05    ~
            42    -284.36325702    2.98e-10 o    7.91e-07 *    2.56e-07 *    5.88e-05      2.64e-05    ~
            43    -284.36325702   -9.41e-10 o    8.38e-07 *    2.73e-07 *    4.24e-05      2.17e-05    ~
            44    -284.36325702    5.76e-10 o    1.10e-06 *    3.07e-07 *    7.41e-05      3.76e-05    ~
            45    -284.36325702    8.75e-10 o    9.99e-07 *    3.49e-07 *    3.90e-05      1.39e-05    ~
            46    -284.36325702   -1.15e-09 o    1.61e-06 *    5.09e-07 *    1.34e-04      5.66e-05    ~
            47    -284.36325702    7.86e-10 o    1.28e-06 *    3.73e-07 *    8.87e-05      4.17e-05    ~
            48    -284.36325702   -1.25e-09 o    9.29e-07 *    2.75e-07 *    8.35e-05      4.41e-05    ~
            49    -284.36325702    7.38e-10 o    5.29e-07 *    2.11e-07 *    7.26e-05      2.99e-05    ~
            50    -284.36325702    7.92e-10 o    7.00e-07 *    2.60e-07 *    5.18e-05      1.79e-05    ~

So I took the last geometry in the file and did a vibrational analysis. Here’s the python script for the calculation:

# gly.py
import psi4
 psi4.set_output_file('gly.out')
 psi4.set_memory('100 GB')
 psi4.set_num_threads(32)

 mol = psi4.geometry("""
     C           -0.496720819429     0.085973520038    -0.002581303779
     N            1.999463902179    -0.005407262029    -0.002584180948
     H            2.026257392296     0.592951952784    -0.817352809296
     H            2.027102166863     0.618729964820     0.792582526584
     C            0.763407304381    -0.755356018637     0.010077751184
     O           -1.605159173664    -0.684347933019     0.010203570329
     O           -0.535179581955     1.289631891684    -0.021612796934
     H            0.718858385973    -1.424410405149    -0.849370445104
     H            0.719870290390    -1.396689509161     0.890458724804
     H           -2.359843677532    -0.082215176293     0.001099823752
 """)

 # psi4.optimize("mp2/aug-cc-pv[tq]z + d:ccsd(t)/aug-cc-pvdz", molecule=mol,
 #                         optimizer_keywords={'G_CONVERGENCE': 'GAU_TIGHT'})
 psi4.set_options({'NORMAL_MODES_WRITE': True})
 psi4.frequency("mp2/aug-cc-pv[tq]z + d:ccsd(t)/aug-cc-pvdz", molecule=mol)

I still get an imaginary frequency. I will try changing the D_CONVERGENCE as you suggested, I’ll set it to 10 and will report with the results, I’ll use the same geometry as above but setting G_CONVERGENCE to GAU_VERYTIGHT.

Where can I check the symmetry of the vibrational mode? I think it’s A but I paste the beginning of the Harmonic Vibrational Analysis for you to confirm:

==> Harmonic Vibrational Analysis <==

   non-mass-weighted Hessian:       Symmetric? True   Hermitian? True   Lin Dep Dim?  3 (0)
   projection of translations (True) and rotations (True) removed 6 degrees of freedom (6)
   total projector:                 Symmetric? True   Hermitian? True   Lin Dep Dim?  6 (6)
   mass-weighted Hessian:           Symmetric? True   Hermitian? True   Lin Dep Dim?  3 (0)
   pre-proj  low-frequency mode:    0.0000i [cm^-1]
   pre-proj  low-frequency mode:    0.0000i [cm^-1]
   pre-proj  low-frequency mode:    0.0000i [cm^-1]
   pre-proj  low-frequency mode:    0.0000  [cm^-1]
   pre-proj  low-frequency mode:    0.0000  [cm^-1]
   pre-proj  low-frequency mode:   66.6287  [cm^-1]
   pre-proj  all modes:['130.2130i' '38.0162i' '0.0000i' '0.0000' '0.0000' '66.6287' '105.6617'
  '211.6142' '264.1203' '463.5174' '519.8997' '619.8253' '658.9138'
  '835.6966' '913.0734' '940.6887' '1140.9267' '1175.5874' '1191.3013'
  '1305.8826' '1389.7999' '1414.1836' '1430.7359' '1684.3767' '1821.8944'
  '3089.4105' '3127.3956' '3518.9315' '3602.9946' '3772.9338']
   projected mass-weighted Hessian: Symmetric? True   Hermitian? True   Lin Dep Dim?  6 (6)
   post-proj low-frequency mode:   41.4351i [cm^-1] (V)
   post-proj low-frequency mode:    0.0000i [cm^-1] (TR)
   post-proj low-frequency mode:    0.0000i [cm^-1] (TR)
   post-proj low-frequency mode:    0.0000i [cm^-1] (TR)
   post-proj low-frequency mode:    0.0000  [cm^-1] (TR)
   post-proj low-frequency mode:    0.0000  [cm^-1] (TR)
   post-proj low-frequency mode:    0.0000  [cm^-1] (TR)
   post-proj  all modes:['41.4351i' '0.0000i' '0.0000i' '0.0000i' '0.0000' '0.0000' '0.0000'
  '209.2949' '263.5180' '463.0805' '519.8472' '619.3245' '658.6918'
  '835.6828' '913.0530' '940.6859' '1140.9125' '1175.5551' '1191.3008'
  '1305.8220' '1389.7970' '1414.1778' '1430.5999' '1684.3759' '1821.8905'
  '3089.3939' '3127.3893' '3518.9314' '3602.9946' '3772.9316']

   Vibration                       1                   8                   9
   Freq [cm^-1]                 41.4351i            209.2949            263.5180
   Irrep                           A                   A                   A
   Reduced mass [u]              2.0091              1.1814              4.4036
   Force const [mDyne/A]        -0.0020              0.0305              0.1802
   Turning point v=0 [a0]        0.0000              0.6978              0.3221
   RMS dev v=0 [a0 u^1/2]        0.0000              0.5363              0.4780
   IR activ [km/mol]            16.6662             192.5582            48.5721
   Char temp [K]                 0.0000             301.1287            379.1438
   ----------------------------------------------------------------------------------
       1   C                0.00 -0.02  0.03    0.00 -0.01  0.03    0.03 -0.14 -0.01
       2   N               -0.01  0.02 -0.12   -0.01  0.02 -0.08   -0.18  0.24  0.04
       3   H               -0.03 -0.29 -0.35   -0.42  0.50  0.26   -0.33  0.21  0.02
       4   H               -0.06  0.33 -0.36    0.37 -0.45  0.28   -0.55  0.34 -0.02
       5   C                0.01 -0.01  0.16    0.01 -0.01  0.07    0.02 -0.11 -0.05
       6   O               -0.01  0.01 -0.14   -0.01  0.01 -0.03   -0.12  0.08  0.04
       7   O                0.02 -0.01  0.11    0.02 -0.01 -0.02    0.28 -0.14 -0.02
       8   H               -0.03 -0.29  0.39   -0.05 -0.09  0.13    0.18 -0.03 -0.12
       9   H                0.07  0.26  0.36    0.09  0.07  0.12    0.13 -0.20 -0.11
      10   H                0.03  0.07 -0.17    0.01  0.04 -0.12    0.01  0.25 -0.04

My confusion from the geometries came from trying to visualize the changes in geometry from the optimization. I didn’t know know which of all the geometries reported I should look at. By checking the output file I noticed that only one pattern (Next Geometry) repeated the same number of times as the number of optimization steps so I got that one. I’m visualizing with ASE, here’s my code for that:

import re

from ase import Atoms
from ase.visualize import view

with open('gly.out', 'r') as f:
    text = f.read()

atoms_list = list()
m_list = re.findall(r'Next Geometry.+?Geometry \(in Angstrom\).+?\s+((?:[A-Z][a-z]?(?: +-?\d+\.\d+){3}\s+)+)', text, re.DOTALL)
for m in m_list:
    symbols = list()
    positions = list()
    for line in m.strip().split('\n'):
        element, x, y, z = line.strip().split()
        symbols.append(element)
        positions.append([float(x), float(y), float(z)])
    atoms = Atoms(symbols=symbols, positions=positions)
    atoms_list.append(atoms)
view(atoms_list)