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.

Note Asymptotic Correction 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:

Ionization potential

First see if there is an experimental IP available from the NIST Chemistry WebbookExternal Link. Convert it to atomic units.

If you need to compute the IP use the $\Delta$DFT approach: $$

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

Important

For notes on using NWChem see my tutorials at:

In particular, see the following:

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 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:

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.

  $ rm *.db
  $ rm *.movecs

Note down the following information from the output file:

...
...
 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):

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:

Here's what this should look like:

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.

  $ rm *.db
  $ rm *.movecs

Note down the following information from the output file:

Analysis

Energies are all in a.u.

Now we can calculate the IP: $$

$$

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:

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