Counterpoise correction and projected Hessian

1

Hello eveyone,

Summary: I found a curious result when running a Counterpoise (CP) corrected Hessian for a known global minimum structure. The projected Hessian differ significantly from the unprojected results, something that did not happen without CP.

Context: I made a second-derivative run at a global minimum structure (C_2h point group symmetry) which was optimized at the Frozen-Core CCSD(T)/aug-cc-pVDZ level (using other quantum package).

  1. This is supposed to be a minimum (thought it may not be at this level of theory).
  2. The imaginary frequency is very large, even if this stationary point, at this level of theory, is a first-order TS.
  3. Since the forces should be zero, before and after the projection the frequencies should be the same, which was the case when no counterpoise correction was included in the job. What surprised me the most is that the nuclear repulsion energy and the electronic energy computed by PSI4 are “exactly” the same as the ones I obtained with the other quantum package.

Input:
#counterpoise vibfreq
psi4.set_memory(‘50 GB’)
psi4.set_num_threads(10)

dimer = psi4.geometry(“”"
0 1
C 0.0000000000 0.6890807464 0.9171748908
O 0.0000000000 0.0000000000 0.0000000000

0 1
C 0.0000000000 -0.6890807457 4.3410689806
O 0.0000000000 0.0000000000 5.2582438714
units angstrom
“”")

psi4.set_options({
“basis”: “AUG-CC-PVDZ”,
“freeze_core”: True,
“scf_type”: “pk”
})

psi4.frequency(‘ccsd(t)’, bsse_type = ‘cp’, molecule=dimer, return_wfn=True)

Output:
==> N-Body: Counterpoise Corrected (CP) energies <==

    n-Body     Total Energy            Interaction Energy                          N-body Contribution to Interaction Energy
               [Eh]                    [Eh]                  [kcal/mol]            [Eh]                  [kcal/mol]
         1     -226.148105592093        0.000000000000        0.000000000000        0.000000000000        0.000000000000

FULL/RTN 2 -226.148590186506 -0.000484594413 -0.304087585344 -0.000484594413 -0.304087585344

==> 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.0000 [cm^-1]
pre-proj low-frequency mode: 0.0000 [cm^-1]
pre-proj low-frequency mode: 2.2855 [cm^-1]
pre-proj low-frequency mode: 2.3337 [cm^-1]
pre-proj low-frequency mode: 11.9213 [cm^-1]
pre-proj low-frequency mode: 23.9077 [cm^-1]
pre-proj low-frequency mode: 25.4804 [cm^-1]
pre-proj low-frequency mode: 69.4484 [cm^-1]
pre-proj all modes:[‘94.3053i’ ‘0.0000i’ ‘0.0000’ ‘0.0000’ ‘2.2855’ ‘2.3337’ ‘11.9213’
‘23.9077’ ‘25.4804’ ‘69.4484’ ‘2104.1167’ ‘2105.6230’]
projected mass-weighted Hessian: Symmetric? True Hermitian? True Lin Dep Dim? 6 (6)
post-proj low-frequency mode: 94.3053i [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.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: 23.8892 [cm^-1] (V)
post-proj low-frequency mode: 25.4769 [cm^-1] (V)
post-proj low-frequency mode: 69.4474 [cm^-1] (V)
post-proj all modes:[‘94.3053i’ ‘0.0000i’ ‘0.0000i’ ‘0.0000i’ ‘0.0000i’ ‘0.0000’ ‘0.0000’
‘23.8892’ ‘25.4769’ ‘69.4474’ ‘2104.1167’ ‘2105.6230’]

Vibration 1 8 9
Freq [cm^-1] 94.3053i 23.8892 25.4769
Irrep Au
Reduced mass [u] 13.4388 13.4388 14.6891
Force const [mDyne/A] -0.0704 0.0045 0.0056
Turning point v=0 [a0] 0.0000 0.6124 0.5672
RMS dev v=0 [a0 u^1/2] 0.0000 1.5874 1.5372
Char temp [K] 0.0000 34.3712 36.6557

  1   C                0.45 -0.34  0.00    0.00 -0.00  0.57    0.10  0.39  0.00
  2   O               -0.34  0.26  0.00    0.00 -0.00 -0.42   -0.13  0.57 -0.00
  3   C                0.45 -0.34 -0.00   -0.00  0.00  0.57   -0.11 -0.39  0.00
  4   O               -0.34  0.25  0.00   -0.00  0.00 -0.42    0.13 -0.57 -0.00
2

Thanks for the report. Can you send me full output files? This reminds me of finite difference displacements landing on a bad state. You normally only see that with open-shell structures or very pathological geometries, though…

3

output.txt (3.2 MB)

Hello,
I noticed my previous upload did not work, so I am uploading it again.
Please let me know what do I need to improve my simulation.

Thank you very much for your support and patience.

4

I recently talked to the original poster and ran some tests with her input plus slight modifications. It does seem like there might be a bug in the CP-corrected frequencies pathway while the CP-corrected gradient seems all fine.

For the CP-corrected calculation for a complex AB, 5 different energy/gradient/Hessian contributions are needed: AB, A in AB basis, A in A basis, B in AB basis, and B in B basis.

  1. In the original calculation in this thread (with Psi4 1.9.1), the displacements selected to evaluate numerical derivatives for those 5 contributions were not all the same. In my calculations (with Psi4 1.10), these displacements are all the same. This is definitely a good thing, but it hasn’t solved the problem - I still observe the same imaginary frequency.

  2. I can confirm that the geometry is a CP-corrected minimum - the coordinates of the CP-corrected geometric gradient don’t exceed 6e-5. I ran a geometry optimization to converge to the minimum even tighter (using G_CONVERGENCE = "GAU_TIGHT"). Now my CP-corrected gradient coordinates don’t exceed 9e-6 and, instead of one frequency of 94.3i cm^-1, I get two smaller ones, 42.7i and 16.4i cm^-1. This lack of stability worries me too. I wasn’t able to converge a CP-corrected geometry optimization with G_CONVERGENCE = "GAU_VERYTIGHT".

  3. The above calculations used point-group symmetry. I ran everything without symmetry as well, and got three imaginary frequencies instead of one. However, those frequencies were stable - they changed minimally when I converged the geometry optimization tighter.

That’s all I was able to find out - let me know if I can be of any further help!