H$_2$ : the potential energy landscape

In this practical session we will perform various types of calculations on H$_2$. We will not bother excessively about basis sets, but more on the methods used. From our discussion in the lectures, we have seen how different RHF, UHF and FCI are for this system. The RHF wavefunction dissociates into the wrong states (H and 1/2 H$^-$), while the UHF and FCI wavefunctions dissociate correctly into two H atoms, though the UHF wavefunction doesn't result in correct energies in the intermediate region.

To see all of this we will not only calculate energies at the H$_2$ equilibrium separation of $R_e=1.4$ a.u. (nearly), but at other separations so as to get the complete potential energy shape. I suggest you choose the separations (in multiples of $R_e$): 0.85, 0.90, 0.95, 1.0, 1.1, 1.2, 1.5, 2.0, 4.0, 8.0, 16.0.

RHF

Here is the restricted Hartree-Fock input file for the STO-3G basis set:

Memory 500 mb

charge 0
Geometry units bohr
  H     0.0    0.0    0.0
  H     0.0    0.0    1.4
End

Basis "ao basis" spherical
  H  library  STO-3G
End

Title "H2 STO-3G "

scf
  Singlet
  RHF
  Direct
end

task scf energy

Save this file in

h2-rhf-sto3g.nw

This file is run as follows:

  mpirun.mpich -np 1 nwchem h2-rhf-sto3g.nw >& h2-rhf-sto3g-1Re.out

Notice the name of the output file includes a 1Re. I've done this to indicate that this particular output is for 1 times $R_e$, the equilibrium bond length of H$_2$. It is good practice to include such information in the file names.

Batch job

Save this file as H2_batch.bash and make it executable using

• $ chmod +x H2_batch.bash

You can now run it using

• $ ./H2_batch.bash

There will be many input and output file created. And the final results will be in a file called data.csv. You may delete all files ending in

• .c
• .b
• .db
• .movecs
• .b^-1

These are temporary files created by NWChem and are not needed here. }}}

UHF

You should have seen that the RHF dissociation energy was not very good. So we now turn to the unrestricted Hartree-Fock (UHF) method. Here we allow the spin up and down orbitals to have different spatial parts. How do we get NWChem to do this? Well, here is the input file:

Memory 500 mb

charge 0
Geometry units bohr
  H     0.0    0.0    0.0
  H     0.0    0.0    1.4
End

Basis "ao basis" spherical
  H  library  STO-3G
End

Title "H2 STO-3G "

scf
  Singlet
  Nopen 2
  UHF
  Sym   Off
  Adapt Off
  Direct
end

task scf energy

The only difference is the scf block. Let's look at it in detail:

• Singlet: We need this as we need a spin $S=0$ state.

• Nopen 2: Tells NWChem we have two open-shell orbitals. This is correct in the dissociation limit, but it isn't in the small-$R$ limit where we would like both electrons in the same orbital (i.e. the RHF solution). So we should expect the UHF energy to be dodgy at short H-H bond lengths. There should be a way of getting NWChem to handle both the short and long-$R$ cases seamlessly, but as of yet I have not found it. So let me know if you figure it out.

• UHF : Should be obvious!

• Sym Off and Adapt Off: This sets symmetry and symmetry-adaptation off. The NWChem wiki suggests this as the dissociated state - as we know from class work - places electrons in non-symmetrized spatial orbitals. But I'm not sure if we need it.

• Direct : As before, don't use the disk for I/O.

...basis geometry etc as before...

Title "H2: RHF energy"

scf
  Singlet
  RHF
  Sym   On
  Adapt On
  Direct
end

task scf energy

Title "H2: UHF energy"

scf
  Singlet
  Nopen 2
  UHF
  Sym   Off
  Adapt Off
  Direct
end

task scf energy

CI : Configuration Interaction

The only way to get the dissociation energy right at all bond lengths is to use full configuration interaction (FCI). But NWChem will not let us do this. Instead we use something called CISD - we will see more of this in subsequent lectures, but for now, CISD = configuration interaction with single and double excitations only. The CISD wavefunctions however, do not dissociate correctly either and for a sufficiently large bond length, will break down. However, let's see just how well it fares.

The STO-3G basis is a minimal basis set. It contains just one $s$-function: that is, one spherical Gaussian orbital. This means the CI space contains just two orbitals of gerade symmetry (exactly as we had in our lectures). While this should be OK, NWChem 6.1.1 seems to have a problem with such a small CI space. So we use a larger basis: cc-pVDZ.

Here is a CISD NWChem command file:

Memory 1000 mb

charge 0
Geometry units bohr
  Symmetry group c2v
  H     0.0    0.0    0.0
  H     0.0    0.0    1.4
End

Basis "ao basis" spherical
  H  library  cc-pVDZ
End

Title "H2 cc-pVDZ "

scf
  Singlet
  RHF
  Direct
end
TCE
  CISD
END

task TCE energy

There are a few differences from what we have used before:

• The geometry block inlcudes the line Symmetry group c2v. NWChem's CI module cannot deal with the infinite point groups $D_{\infty h}$ and $C_{\infty v}$. So we use a sub-group: in this case $C_{2v}$.

• The CI wavefunction is calculated using the '''TCE''' module. This is the Tensor Contraction Engine.

• I have used a RHF wavefunction in the above. For longer bond lengths use the UHF wavefunctions.

How do your results compare with those from Helgaker et al.:

and here are the energies from different methods: