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