There are (at least) two types of wave function stability:

(1) ROHF-> ROHF internal stability (usually meaning whether a lower ROHF solution can be found.

(2) ROHF->UHF external stability (usually meaning whether UHF energy is lower than that of ROHF).

The 2nd type of stability can be equivalently checked by computing UHF and comparing the electronic energies. As for the 1st type of stability, it is not supported in Gaussian as far as I know.

But the open source packages PySCF(https://github.com/pyscf/pyscf) and PSI4(https://github.com/psi4/psi4) support that.

You can transfer ROHF MOs from Gaussian into PySCF/PSI4 and check the wave function stability in the latter package. Taking the triplet O2 as an example, the Gaussian input file is

```

%chk=O2_triplet.chk

%mem=4GB

%nprocshared=4

#p ROHF/cc-pVDZ nosymm int=nobasistransform

title

0 3

O 0.0 0.0 0.0

O 0.0 0.0 1.1616

```

The keywords 'nosymm int=nobasistransform' here is for the usage of transferring MOs. After the Gaussian job is accomplished, run the following command in Shell

```

bas_fch2py O2_triplet.chk

```

This will generate the PySCF input file O2_triplet.py. Cartesian coordinates, basis set data and some keywords have been written in this file. Since we want to check the wave function stability, we can open this .py file and append

```

dm = mf.make_rdm1()

mf.max_cycle = 10

mf.kernel(dm)

mf.stability()

```

Then submit the PySCF job

```

python O2_triplet.py >O2_triplet.out

```

You can find in the output that ROHF is converged in 1 cycle and it is stable in the internal stability analysis.

The `bas_fch2py` utility is free and open source, in the MOKIT package(https://gitlab.com/jxzou/mokit). You can also use the `fch2psi` utility to transfer MOs from Gaussian to PSI4, and check ROHF wave function stability.

Jingxiang Zou Note however that when you run UHF with Gaussian, you can reach the RHF stationary point and then you can falsely convinced that everything is OK

Tatiana Korona Hi, what does "you can falsely convinced that everything is OK" mean?

Jingxiang Zou Well, if you run Gaussian for H2 at distance 4 Ang with the input

#N UHF/6-31G(d) SP

you will get

SCF Done: E(UHF) = -0.770202601598 A.U. after 5 cycles NFock= 5 Conv=0.30D-11 -V/T= 1.9607 = 0.0000 = 0.0000 = 0.0000 =-0.0000 S=-0.0000 = 0.00000000000 Annihilation of the first spin contaminant: S**2 before annihilation -0.0000, after 0.0000

while the correct answer for the UHF for this system can be obtained only with guess=mix added:

#N UHF/6-31G(d) SP guess=mix

SCF Done: E(UHF) = -0.996485078493 A.U. after 6 cycles NFock= 6 Conv=0.17D-10 -V/T= 1.9774 = 0.0000 = 0.0000 = 0.0000 = 0.9999 S= 0.6180 = 0.00000000000 Annihilation of the first spin contaminant: S**2 before annihilation 0.9999, after 0.0000

However, in the first version there is no warning in the output that we are actually calculating RHF:

#N RHF/6-31G(d) SP

SCF Done: E(RHF) = -0.770202601598 A.U. after 5 cycles

so somebody could be happy that he/she obtained the UHF results, while in reality these are the RHF results

Tatiana Korona Your detailed explanation is appreciated. I know what you mean, and I know that a singlet UHF calculation may lead to RHF result.

But I'm sorry that the RHF stationary point has nothing to do with this question, and has nothing to do with my answer in April 28. This question is about ROHF wave function stability, not RHF, where we assuming that we are talking about non-singlet calculations. My answer in April 28 focused on the ROHF wave function stability, and a triplet example was shown.