CBS optimization didn't converge and we got imaginary frequencies

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.

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.

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.

  1. Use set_options, as described in the PythonAPI tutorial.
  2. 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."

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.

  1. 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.
  2. geom_maxiter. This is in the geometry optimization page of the manual.
  3. 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.