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)
```