5.3.2 Open-shell Hartree-Fock calculations

In the above example all orbitals are doubly occupied, and therefore the total symmetry of the ground state wavefunction, which is the direct product of the spin-orbital symmetries, is $A_1$. Since all electrons are paired, it is a singlet wavefunction, $^1A_1$.

Even though most stable molecules in the electronic ground states are closed-shell singlet states, this is not always the case. Open-shell treatments are necessary for radicals or ions. As a first simple example we consider the positive ion of formaldehyde, H$_2$CO$^+$. In this case it turns out that the lowest cation state is $^2B_2$, i.e., one electron is removed from the highest occupied orbital in symmetry $b_2$ (denoted 2.3 in MOLPRO, see above). By default MOLPRO assumes that the number of electrons equals the total nuclear charge. Therefore, in order to compute an ion, one has to specify the number of electrons. Alternatively, the total charge of the molecule can be given (see below). Furthermore, one has to specify the symmetry and spin of the wavefunction. This is done using the wf directive (wf stands for wavefunction):

wf,15,3,1

The first entry on a wf card is the number of electrons. The second entry is the total symmetry of the wavefunction. For a doublet this equals the symmetry of the singly occupied molecular orbital. Finally, the third entry specifies the number of singly occupied orbitals, or, more generally the total spin. Zero means singlet, 1 doublet, 2 triplet and so on. Alternatively, the WF card can be written in the form

wf,charge=1,symmetry=3,spin=1

where now the total charge of the molecule instead of the number of electrons is given. In this case the number of electrons is computed automatcially from the nuclear and total charges.

In summary, the input for H$_2$CO$^+$ is

{geometry specification}
{basis specification}

{hf           !invoke spin restricted Hartree-Fock program (rhf can also be used)
occ,5,1,2    !number of occupied orbitals in the irreps a1, b1, b2, respectively
wf,15,3,1}   !define wavefunction: number of electrons, symmetry and spin
---

Note the curley brackets, which are required and enclose the command block hf.

As a second example we consider the ground state of O$_2$, which is $^3\Sigma^-_g$. The geometry specification is simply

geometry     !geometry specification, using z-matrix
o1
o2,o1,r
end

r=2.2 bohr   !bond distance

MOLPRO is unable to use non-abelian point groups, and can therefore only use $D_{2h}$ in the present case. The axis of a linear molecule is placed on the $z$-axis of the coordinate system. Then the symmetries of the $\sigma_g$, $\pi_{u,x}$, $\pi_{u,y}$, $\sigma_u$, $\pi_{g,x}$, $\pi_{g,y}$ orbitals are 1,2,3,5,6,7, respectively. It is easier to remember that the irreducible representations in $D_{2h}$ are carried by the functions (molpro symmetry numbers in parenthesis) $s(1), x(2), y(3), xy(4), z(5), xz(6), yz(7), xyz(8)$. The order in $C_{2v}$ is the same, but then there are only the first four irreducible representations.

The electron configuration of the electronic ground state of O$_2$ is

$1\sigma_g^2 1\sigma_u^2 2\sigma_g^2 2\sigma_u^2 3\sigma_g^2 1\pi_{u,x}^2 1\pi_{u,y}^2 1\pi_{g,x}^1 1\pi_{g,y}^1$.

Thus, the number of occupied orbitals in the 8 different irreducible representations of the $D_{2h}$ point group are specified as

occ,3,1,1,0,2,1,1,0

The product symmetry of the singly occupied orbitals 1.6 and 1.7 is 4 (xy), and therefore the symmetry of the total wavefunction is 4 (please refer to the MOLPRO reference manual for a more complete account of symmetry groups and the numbering of irreducible representations). Thus, the wf card reads

wf,16,4,2      !16 electrons, symmetry 4, triplet (2 singly occupied orbitals)

This is still not unambiguous, since the product symmetry of $\pi_{u,x}$ (2) and $\pi_{u,y}$ (3) is also 4, and therefore the program might not be able to decide if it should singly occupy the $\pi_u$ or $\pi_g$ orbitals. Therefore, the singly occupied orbitals can be specified using the open directive:

open,1.6,1.7

This now defines the wavefunction uniquely. In summary, the input for O$_2$ reads

***,O2
print,basis,orbitals
geometry                   !geometry specification, using z-matrix
o1
o2,o1,r
end

r=2.2 bohr                 !bond distance

basis=vtz                  !triple-zeta basis set

{hf                        !invoke RHF program
wf,16,4,2                  !define wavefunction: 16 electrons, symmetry 4, triplet
occ,3,1,1,,2,1,1           !number of occupied orbitals in each symmetry
open,1.6,1.7}              !define open shell orbitals

In fact, the last 2 lines are not necessary in the present case, since the correct configuration can be automatically determined using the Aufbau principle, but this might not always be true.