<<TableOfContents()>>

= Ionization potentials and the asymptotic correction =


== Ionization energy and AC-SHIFT ==
Eventually we will use LC-PBE/LC-PBE0 for SAPT(DFT) calculations, but till then, we will continue to use an asymptotically corrected hybrid functional like PBE0. 

{{{#!wiki  note
Note
Asymptotic Correction
[[AJMPublic/teaching/intermolecular-interactions/dft_in_saptdft|Click here for an introduction to the asymptotic correction]], but bear in mind that the Fermi-Amaldi correction that I talk about on that page is not the one we now recommend. We now think that the GRLB correction is best. More on this some other time. 
}}}

For the asymptotic correction you need to know the shift needed to correct the XC potential. This shift is defined as
$$  \Delta = {\rm IP} + \epsilon_{\rm HOMO}.  $$
So you need two quantities to define the shift:
  * The vertical ionization energy, and
  * the energy of the HOMO eigenvalue.
  * NOTE: all quantities must be in atomic units!


=== Ionization potential ===
First see if there is an experimental IP available from the [[http://webbook.nist.gov/chemistry/|NIST Chemistry Webbook]][[http://example.com|External Link]]. Convert it to atomic units.

If you need to compute the IP use the $\Delta$DFT approach:
$$  {\rm IP} = E[N-1] - E[N],  $$
that is, calculate the energies of the $N$ and $N-1$ electron systems and determine the IP by subtraction.
Use PBE0 or LC-PBE0 to determine these energies. Basis sets could be Sadlej-pVTZ, cc-pVDZ or cc-pVTZ. These basis sets may be augmented, but augmentation is typically needed for anions.   

NWChem and DALTON2016 also need the AC-SHIFT. This is just the $\Delta$ in the equation defined above. DALTON2006 uses only the IPs and determines the shift on-the-fly. This is better but is more time-consuming. To compute the shift you need to calculate the HOMO eigenvalue //without// an asymptotic correction.

{{{#!wiki  important
Important

For notes on using NWChem see my tutorials at:
  * [[AJMPublic/teaching/electronic-structure|Electronic Structure Methods]]
In particular, see the following:
  * [[AJMPublic/teaching/electronic-structure/practical1|Linux : Using Comanche]]
  * [[AJMPublic/teaching/electronic-structure/practical-ar2-dft|Ar$_2$ : DFT calculations of the PES]]
Of course, the one on Using Comanche only applies if you have an account there.

It may be easiest to use the Psi4 code for this calculation. Psi4 has nice tutorials on all sorts of calculations [[http://www.psicode.org/|that you can find on their wedsite.]] 
}}}

== Calculations with NWChem ==
You will need two NWChem command files. We will use Cluster to generate one and then get the other using a simple edit.

Here is the Cluster file we will use for the NWChem calculation:

[[attachment:HCN.clt]]

{{{
Title HCN 
 
Global
  Units Bohr Degrees
  Overwrite Yes
End
 
Molecule HCN
   ! HCN from 99_BenzeneHCNcomplex.xyz
   Ion-Charges Auto
   Units Angstrom
   N   -0.0034118   3.5353926   0.0000000
   C    0.0751963   2.3707040   0.0000000
   H    0.1476295   1.3052847   0.0000000
End

Show HCN in XYZ format  ( you could also use PDB as the format )
 
Run-Type
  Properties
 
  Molecule    HCN
 
  Main-Basis   aVTZ   Type  MC
  Func         PBE0
  AC           NONE            ( Must be NONE )
  SCF-code     NWChem          ( DALTON2015, NWCHEM, etc )
  File-Prefix  HCN_aTZ         ( change this to something reasonable )
End

Finish
}}}

You may want to perform these calculations in a special directory. I tend to call mine ''IP'':

{{{
  $ mkdir IP
  $ cd IP  
}}}

Create the Cluster file here.

Run this using either
{{{
  $ cluster --job HCN_aTZ < HCN.clt  
}}}

or, if you have the newer Cluster code

{{{
  $ cluster < HCN.clt  
}}}
  
Now you should have the following files:
{{{
  $ ls
  HCN_aTZ_A.nw          HCN_aTZ.cks           HCN_aTZ.ornt  HCN_aTZ.sites  HCN.xyz
  HCN_aTZ_casimir.prss  HCN_aTZ_display.ornt  HCN_aTZ.prss  HCN.clt  
}}}

You need only the ''HCN_aTZ_A.nw'' file. 

Make two copies of this file:

{{{
  $ cp HCN_aTZ_A.nw HCN_neutral.nw
  $ cp HCN_aTZ_A.nw HCN_cation.nw  
}}}

  
=== Neutral ===

The first, neutral, file needs only one edit:
{{{
  * Change the ''Scratch_dir <SCRATCHDIR>'' line to contain your actual scratch directory path: ''Scratch_dir /scratch/alston/''.  
}}}

This file should look like:

{{{
Title "Properties:  HCN Method: PBE0 "
Scratch_dir /scratch/alston/
Memory     4096 mb
Start HCN_aTZ_A

Charge     0
Geometry Units Bohr nocenter noautoz noautosym
   N            -0.00644737       6.68092328       0.00000000
   C             0.14210040       4.47998096       0.00000000
   H             0.27897930       2.46663042       0.00000000
End
Basis "ao basis" SPHERICAL
      N          library   aug-cc-pVTZ
      C          library   aug-cc-pVTZ
      H          library   aug-cc-pVTZ
End
set lindep:tol 1e-5
dft
   direct
   grid fine
   iterations 60
   convergence energy 1e-6
   xc PBE0
   Mult 1
end
task dft energy
}}}

Run this using the command:

{{{
  $ mpirun.mpich -np 2 nwchem HCN_neutral.nw >& HCN_neutral.out &  
}}}

For larger systems you may use more processors using ''-np 4'' for 4 cores, etc. Just don't use more cores than you have. 

  * The output will be in the file ''HCN_neutral.out''. 
  * Remove the *.db and *.movecs files as these can use a problem for the next calculation.

{{{
  $ rm *.db
  $ rm *.movecs
}}}

Note down the following information from the output file:
  * ''Total DFT energy'': $E[N] = -93.342497219381$ (units are Hartree, or atomic units). This is the energy of the $N$-electron system.
  * We also need the energy of the highest occupied molecular orbital (HOMO). This can be found after the Total DFT energy is printed. You will see text like the following:
{{{
...
...
 Vector    6  Occ=2.000000D+00  E=-3.794473D-01
              MO Center=  3.3D-02,  3.0D+00,  3.0D-05, r^2= 1.1D+00
   Bfn.  Coefficient  Atom+Function         Bfn.  Coefficient  Atom+Function
  ----- ------------  ---------------      ----- ------------  ---------------
    11      0.301797  1 N  pz                57      0.288181  2 C  pz
     8      0.245963  1 N  pz                54      0.225239  2 C  pz
    14      0.210149  1 N  pz

 Vector    7  Occ=2.000000D+00  E=-3.794473D-01
              MO Center=  3.3D-02,  3.0D+00, -3.0D-05, r^2= 1.1D+00
   Bfn.  Coefficient  Atom+Function         Bfn.  Coefficient  Atom+Function
  ----- ------------  ---------------      ----- ------------  ---------------
     9      0.301122  1 N  px                55      0.287522  2 C  px
     6      0.245413  1 N  px                52      0.224720  2 C  px
    12      0.209677  1 N  px

 Vector    8  Occ=0.000000D+00  E=-5.255275D-03
              MO Center=  7.3D-02,  2.4D+00, -1.4D-05, r^2= 3.5D+00
   Bfn.  Coefficient  Atom+Function         Bfn.  Coefficient  Atom+Function
  ----- ------------  ---------------      ----- ------------  ---------------
    61      0.598636  2 C  px                15     -0.496325  1 N  px
    12     -0.341957  1 N  px                58      0.335482  2 C  px
     9     -0.288065  1 N  px                55      0.272180  2 C  px
     6     -0.227074  1 N  px                52      0.191189  2 C  px
   103      0.175193  3 H  px
...
}}}
Vector 7 happens to be the HOMO level here. Vector 8 has no electrons in it (occupancy=0.0):
  Vector    7  Occ=2.000000D+00  E=-3.794473D-01
The energy of this vector is $\epsilon_{\rm HOMO} = -3.794473\times 10^{-01}$ a.u. Note this down too.

=== Cation ===

The second, cation, file needs to be edited. You need to change the charge and multplicity:
  * Change the ''Scratch_dir <SCRATCHDIR>'' line to contain your actual scratch directory path: ''Scratch_dir /scratch/alston/''.
  * charge should be 1
  * multplicity should be 2
Here's what this should look like:

[[attachment:HCN_cation.nw]]

{{{
Title "Properties:  HCN Method: PBE0 "
Scratch_dir /scratch/alston/
Memory     4096 mb
Start HCN_aTZ_A

Charge     1
Geometry Units Bohr nocenter noautoz noautosym
   N            -0.00644737       6.68092328       0.00000000
   C             0.14210040       4.47998096       0.00000000
   H             0.27897930       2.46663042       0.00000000
End
Basis "ao basis" SPHERICAL
      N          library   aug-cc-pVTZ
      C          library   aug-cc-pVTZ
      H          library   aug-cc-pVTZ
End
set lindep:tol 1e-5
dft
   direct
   grid fine
   iterations 60
   convergence energy 1e-6
   xc PBE0
   Mult 2
end
task dft energy
}}}

Run this using the command:
{{{
  $ mpirun.mpich -np 2 nwchem HCN_cation.nw >& HCN_cation.out &
}}}

For larger systems you may use more processors using ''-np 4'' for 4 cores, etc. Just don't use more cores than you have. 
  * The output will be in the file ''HCN_cation.out''. 
  * Remove the *.db and *.movecs files.

{{{
  $ rm *.db
  $ rm *.movecs
}}}
  
Note down the following information from the output file:
  * ''Total DFT energy'': $E[N-1] = -92.846349549335$ (units are Hartree, or atomic units). This is the energy of the $N-1$-electron system.

=== Analysis ===
  * $E[N-1] = -92.846349549335$
  * $E[N] = -93.342497219381$
  * $\epsilon_{\rm HOMO} = -3.794473\times 10^{-01}$
Energies are all in a.u.

Now we can calculate the IP:
$$  {\rm IP} = E[N-1] - E[N] = -92.846349549335 - (-93.342497219381) = 0.4962 {\rm a.u.}            $$

For the asymptotic correction you need to know the shift needed to correct the XC potential. This shift is defined as
$$  \Delta = {\rm IP} + \epsilon_{\rm HOMO} = 0.4962 + (-0.3794) = 0.1168 {\rm a.u.}  $$

Insert these bits of information in the ''MOLECULE'' block in the Cluster file:

[[attachment:HCN.clt]]

{{{
Title HCN 
 
Global
  Units Bohr Degrees
  Overwrite Yes
End
 
Molecule HCN
   IP       0.4962 a.u.
   AC-Shift 0.1168 a.u.
   ! HCN from 99_BenzeneHCNcomplex.xyz
   Ion-Charges Auto
   Units Angstrom
   N   -0.0034118   3.5353926   0.0000000
   C    0.0751963   2.3707040   0.0000000
   H    0.1476295   1.3052847   0.0000000
End
...
...
}}}