accidentally I observed that in some cases GAMESS calculates lower energies for ROHF wavefunctions than MOLPRO or ORCA. Below are the energies obtained with SCFTYP=ROHF, MULT=2, ICHARG=1 (i.e. the radical cation) of the monomer, dimer, and trimer of glycine:
gly1(+) gly2(+) gly3(+)
E_rohf_orca -282.53982397 -489.36557909 -696.16963524
E_rohf_molpro -282.53982408 -489.38065114 -696.16963031
E_rohf_gamess -282.53982408 -489.38313835 -696.22387017
In all cases the cc-pVDZ basis was used, with D5=.T. set in $CONTRL in Firefly.
Whereas all three programs agree for the monomer, the results differ for the dimer and the trimer.
Since all three programs agree for the monomer, I conclude that the problem is not a difference in the basis set (e.g. spherical harmonics vs. cartesian).
Since in all three cases the geometry is fixed to planar, the unpaired electron could be in an a' or an a'' orbital. GAMESS always puts it in a'. It seems that one can not influence this choice in ROHF, but one can set up a MCSCF run that mimicks a ROHF calculation. I.e., if I read in the ROHF orbitals and then make a CAS(1,1) calculation by setting NDOC=0 NALP=1 NVAL=0, I get exactly the same result as ROHF. If I now switch the singly occupied orbital of a'-symmetry with either an occupied or a virtual a''-orbital, the MCSCF yields energies of -696.190903226 and -695.90391231 H for the trimer. The second state is apparently an excited state, the first is lower than the ORCA and MOLPRO result.
So I suspect that the algorithm in GAMESS is different from ORCA and MOLPRO, but even MOLPRO and ORCA dont always agree. So is ROHF defined differently by different groups? Or could it be a result of different initial guesses? In the latter case, ORCA and MOLPRO could run into local minima, and a different guess would go to the same energy as GAMESS did.
Does anybody know about such a difference between GAMESS and other programs with respect to ROHF?