## page was renamed from AJMGrpOnly/teaching/electronic-structure/practical2
## page was renamed from ajm/teaching/electronic-structure/practical2

<<TableOfContents()>>

= 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.

{{{#!wiki  important
Important
NWChem: You will need to get used to using the NWChem wiki that is found at:
  * [[http://www.nwchem-sw.org/index.php/Release62:NWChem_Documentation|NWChem Documentation Wiki]]
Bookmark this page as you will need to use it often. This should be your first port of call should something go wrong. 
}}}

== 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.

{{{#!wiki  warning
Warning
Ideally we should not be running jobs on our home directory. NWChem will create temporary files that are written to disk, but the home directory on '''comanche''' is special: it is a mounted RAID6 disk consisting of 6 or more hard drives. It is a secure disk system, but isn't meant to handle large amounts of disk I/O. So we should ideally run on '''/scratch/HDD_2T/esm/<user>/''' and copy only the results to our home directory. 

To make the disk I/O as little as possible, we use the '''direct''' option in NWChem. This tells the program not to use the disk for storage, but to do everything //directly// in memory.
}}}

{{{#!wiki  tip
Tip
Why have we chosen a '''singlet''' state?
  * This is because we want to study the closed-shell $S=0$ H$_2$ state.
  * To get the open-shell triplet state use '''triplet''' in the '''scf''' block.
}}}

{{{#!wiki  important
ImportantTask: RHF Energy Surface:
  * How will you change the above file to calculate energies at the other configurations?
  * Reminder: configurations needed (in units 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.
Do this. And tabulate the energies.

'''Analysis''': Does the RHF wavefunction dissociate into two H-atoms? How would you know? According to our theory, the RHF dissociation energy should be
\begin{align}
  E({\rm H}_2,{\rm RHF}) \rightarrow E({\rm H}) + \frac{1}{2} E({\rm H}^-)
\end{align}
To verify this we need to calculate the energy of H$^-$. How would you modify the above input file to get this energy? Hint: [[ajm:teaching:electronic-structure:h-ion-nwchem|here's the H$^-$ command file]].

Perform this calculation and see if the theoretical dissociation energy is correct.
}}}
{{{#!wiki  tip
Tip
'''Making Plots'''\\
You will need to make plots energy as a function of separation $R$. The goal here should be to display the data points clearly and make a smooth fit to the data. This can be done using a variety of programs, but I recommend [[http://www.gnuplot.info/|Gnuplot]] as it can do the job very well and without much of a fuss. Gnuplot is available on comanche but you will need to enable X-forwarding (instructions for [[http://www.tldp.org/HOWTO/XDMCP-HOWTO/ssh.html|Putty]], and [[https://wiki.archlinux.org/index.php/Secure_Shell#X11_forwarding|OSX via the Terminal]]) to see the output of Gnuplot on your screen. If this is a problem, consider installing Gnuplot on your laptop. Or else something like [[http://home.gna.org/veusz/|Veusz]]. If all fails, use Excel. 

Plots should be neat and informative. Axes should be clearly labeled, units included, fits should be smooth, lines and points thick/large enough and a legend present.
}}}

{{{#!wiki  tip
Tip
'''Searching for energies in multiple output files'''
The Unix ''grep'' command is great for searching for lines that will always contain a particular string. To see how it can be used try
  man grep
To search for SCF energies in files that are of the form ''JOB_nnn.out'' use
  grep -in 'Total SCF energy' *.out
}}}

[[attachment:Batch job]]

{{{#!wiki  tip
Tip

Here's a BASH script from Matej Veis to do a whole batch of calculations for you and search the output:

#creates input files and executes NWChem, writing energy values for separation n
#version 1.0 by MV 

Rhh=0.0
Re=1.4

for n in {0.85,0.9,0.95,1.0,1.1,1.2,1.5,2.0,4.0,8.0,16.0}; do

Rhh_str="${n} * ${Re}"
Rhh="$(echo "$Rhh_str" | bc -l)"

echo -e "Memory 500 mb

charge 0
Geometry units bohr
  H     0.0    0.0    0.0
  H     0.0    0.0    $Rhh
End

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

Title \"H2 STO-3G\" 

scf
  Singlet
  RHF
  Direct
end

task scf energy" >input$n.nw
mpirun.mpich -np 1 nwchem input$n.nw >& output$n.out

s=`grep -in "Total SCF energy" output$n.out`

echo -e "$s , $n " >>data.csv
done

}}}
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.
}}}

{{{#!wiki  tip
Tip
If you prefer to use Python3 then use this code written by Harry Campion:
{{{
#!/usr/bin/env python3
#
# Created by Harry Campion on 26/01/2019
#
# Runs NWChem batch job and prints output into .dat file

import os
import re
import time
import subprocess
import glob

Re = 1.40
sep = [0.85,0.9,0.95,1.0,1.1,1.2,1.5,2.0,4.0,8.0,16.0]
basis = "sto-3g"
Method = "rhf"

data_filename = "h2-" + Method + "-" + basis + ".dat"
with open(data_filename, "w") as data:
    data.write("""# Created by H2_batch script.
# 
# SCF Energies for H2 computed by NWChem.
#
# SCF Method : %s
# Basis Set  : %s
# Units      : a.u.
#
 R     Energy
""" % (Method, basis))

def maketemplate(Rhh_str, basis):
    template="""Memory 500 mb

charge 0
Geometry units bohr
  H     0.0    0.0    0.0
  H     0.0    0.0    %.2f
End
     
Basis "ao basis" spherical
  H  library  %s
End
        
Title "H2 %s" 
         
scf
  Singlet
  RHF
  Direct
end
                
task scf energy""" % (Rhh_str, basis, basis)
    return template

for n in sep:
    Rhh_str = n*Re
    n = str(n)
    template = maketemplate(Rhh_str, basis)
    dir_name = n + "Re"
    os.mkdir(dir_name)
    os.chdir(dir_name)
    filename_input = "h2-rhf-" + basis + "_" + n + "Re" + ".nw"
    filename_output = "h2-rhf-" + basis + "_" + n + "Re" + ".out"
    with open(filename_input, "w") as file:
        file.write(template)
    runjobs = subprocess.Popen("mpirun.mpich -np 1 nwchem %s > %s &" % (filename_input, fi
lename_output), shell=True)
    
    time.sleep(0.25)

    pattern = re.compile("         Total SCF energy =     (\-\d+\.\d+)")
    with open(filename_output, "r") as output_data:
        for line in output_data:
            match = pattern.match(line)
            if match:
                scf_energy = float(match.group(1))
                os.chdir("../")
                with open(data_filename, "a") as data:
                    data.write("\n%.2f" % (Rhh_str) + "  %.5f" % (scf_energy))

}}}

== 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.

{{{#!wiki  important
Important
Task: UHF energy profile:
  * As with the RHF case, calculate UHF energies for the set of H$_2$ bond lengths. 
  * Compare these with the RHF case. At some bond length you should see the UHF energy get lower than the RHF energy. This is the crossover point.
  * Is the UHF dissociation energy close to that of two H-atoms?
  * It should be possible to do both an RHF and UHF calculation in the same file using a command block like:
}}}
{{{
...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 [[http://www.nwchem-sw.org/index.php/Release62:TCE|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.

{{{#!wiki  important
Important
Tasks:
  * Include CI energies at the bond lengths you have already used.
  * NWChem may not allow convergent energies at all bond lengths. This is a bug in the code. 
}}}

{{{#!wiki  important
Important
Task:
  * Plot the energies. You will not be able to make plots using '''comanche''' due to a problem forwarding the display to Windows machines. 
  * Use '''PhysPlot''' or something else of your choice.
}}}

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

{{attachment:helgaker-fig5p12-h2-pes.png||width=600}}

and here are the energies from different methods:

{{attachment:helgaker-tab5p6-h2-energies-methods.png||width=600}}