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

• 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 NIST Chemistry WebbookExternal 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.

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:

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:

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:

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

AJMPublic/camcasp/ionization-potentials (last edited 2021-04-28 11:56:06 by apw185)