Enforcing Symmetry for Dissociated Molecules

[ibosko]

Hello everyone!

When calculating the SCF (STO-3G) for small diatomic molecules (hydrogen molecules or helium dimer), I noticed that as the distance between atoms increases, the MO coefficients change in such a way as if at some point there is a switch from RHF to UHF.

For example, MO coefficients for a hydrogen molecule at 21 rBohr are

[[ 0.70710678 0.70710678]
[ 0.70710678 -0.70710678]]

and for 22 rBohr

[[1. 0.]
[0. 1.]].

Does anyone know exactly what setting needs to be changed to make the program work in RHF even over very long distances? Or maybe something else is going on?

There is a full script:

import psi4
import numpy as np

psi4.set_memory(int(5e8))
numpy_memory = 2

h2 = psi4.geometry("""
H
H 1 21.0
units=bohr
symmetry c1
""")

psi4.set_options({'basis': 'sto-3g',
                  'scf_type': 'pk',
                  'e_convergence': 1e-15
                 },)

e, wfn = psi4.energy('scf', return_wfn=True)
Matrix_C = wfn.Ca()
C = Matrix_C.np

print(C)

Thanks in advance.

[jmisiewicz]

Your question confuses me, and I suspect you are using technical terms incorrectly.

You cannot look at the Ca matrix and say anything about whether spin symmetry (RHF vs UHF) is preserved or not. The test for whether you have switched from “same orbitals for different spin” to “different orbitals for different spin” is whether wfn.Ca() matches wfn.Cb() or not. If the reference keyword is not explicitly invoked, the two will always match.

My best guess is that you mean that the MO coefficients change in a way that breaks spatial symmetry. Your first Ca() matrix has a clear bonding and antibonding orbital, where your second Ca() matrix has one orbital localized on each hydrogen. The simplest fix to this is to remove symmetry c1 from your molecule specification. The role of that keyword is exactly to enable this behavior.

If you are performing a potential energy curve, you should also consider reading in the orbitals from a previous point on the curve as a guess for the next point on the curve.

Please let me know whether this solves your problem. If not, please describe your problem in more detail.

[ibosko]

Thank you for your reply!

You are right, a breaking of spatial symmetry is what I meant.

So, when I remove symmetry c1 from the code above, the program returns

ValidationError: Attempted to call .np on a Psi4 data object with multiple irreps.Please use .nph for objects with irreps.

Therefore I use

Matrix_C = wfn.Ca()
C = Matrix_C.nph

instead of

Matrix_C = wfn.Ca()
C = Matrix_C.np

Here is what it returns:

(array([[1.]]), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([[1.]]), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64))

Seems it does not work correctly without symmetry specification.

[jmisiewicz]

If you don’t want your C matrix to be symmetry-blocked, use wfn.Ca_subset("AO", "ALL") instead of wfn.Ca(). This tells Psi to deliver the C matrix in terms of atomic orbitals, even if a description in terms of symmetry-adapted orbitals is possible. You’ll then use np rather than nph.

[ibosko]

Thank you so much!
Now it all works as needed.