Navigation:
IP and AC-shift calculations with Psi4
The following Psi4 code can be used to easily calculate vertical ionization potentials (IPs) and the shifts needed for the asymptotic correction (AC-Shift). Since we do this sort of calculation often and it can be error-prone, I wrote the following code. It was also a small challenge to me to do everything in one calculation, and use some of the Python constructs in Psi4.
This is an example for a molecule taken from the IL16 database. To enter your molecule just alter the contents of the molecule block. Make sure you set the correct charge and multiplicity. The code should deal with charged systems correctly. You probably want to choose a better basis too. The code defaults to the PBE0 functional. This is a good choice for IPs and is recommended. Though some of the new meta-GGAs (the "//ab initio//" ones) may prove better.
<code | IL150-IP-calculation.psi4>
memory 4000 mb
molecule monomer {
1 1 # charge and multiplicity. In this case the molecule has charge +1
units angstrom
C1 1.25456 0.04517 -0.22319
N2 -0.03032 0.04506 0.31669
C3 0.07512 0.04497 1.66566
N4 1.36986 0.04500 1.99001
C5 2.12628 0.04514 0.83592
C6 -1.29617 0.04507 -0.43477
H7 1.68755 0.04500 2.95634
H8 -0.72290 0.04489 2.39723
H9 3.20705 0.04520 0.85849
H10 1.42695 0.04526 -1.29040
H11 -1.07074 0.04515 -1.50429
H12 -1.87711 -0.85035 -0.18843
H13 -1.87716 0.94042 -0.18831
# Unfortunately, Psi4 1.x does not write out the orbital energies
# correctly if symmetry is used. So we disable it here by forcing C1.
# This will have the consequence that the calculation will proceed
# more slowly, and, for symmetric systems, the calculation may even
# fail to converge. If symmetry is enabled then the e_HOMO and AC-shift
# values may be incorrect and these should be checked and corrected manually.
Symmetry C1
}
# Define the XC functional here
dft_func = 'pbe0'
set basis aug-cc-pVDZ
set freeze_core TRUE
# get the charge and multiplicity of the monomer
charge = monomer.molecular_charge()
mult = monomer.multiplicity()
set reference rks
EN, wfn = energy(dft_func, return_wfn=True)
# We need to find the HOMO eigenvalue of this system
nelec_a = wfn.nalpha()
epsilon_a = wfn.epsilon_a().np # The .np is need to convert the Psi4 Vector object to an np.array
e_homo = epsilon_a[nelec_a-1]
# now remove and electron and increase the multiplicity.
monomer.set_molecular_charge(charge+1)
monomer.set_multiplicity(mult+1)
set reference uks
EN1 = energy(dft_func)
IP = EN1 - EN
AC_shift = IP + e_homo
psi4.print_out("\n")
psi4.print_out("%s energies and IP \n\n" % (dft_func))
psi4.print_out(" Num Electrons Charge E_int [a.u.] \n")
psi4.print_out("-----------------------------------------------------\n")
psi4.print_out(" N %3d %10.6f \n" % (charge,EN))
psi4.print_out(" N-1 %3d %10.6f \n" % (charge+1,EN1))
psi4.print_out("-----------------------------------------------------\n")
psi4.print_out(" IP %10.6f a.u. \n" % (IP))
psi4.print_out(" e_HOMO %10.6f a.u. \n" % (e_homo))
psi4.print_out(" AC-shift %10.6f a.u. \n" % (AC_shift))
psi4.print_out("-----------------------------------------------------\n")The important part of the output will be at the bottom and should look like:
<code | IL150 snippet>
PBE0 energies and IP
Num Electrons Charge E_int [a.u.]
-----------------------------------------------------
N 1 -265.631128
N-1 2 -265.094762
-----------------------------------------------------
IP 0.536366 a.u.
e_HOMO -0.457412 a.u.
AC-shift 0.078954 a.u.
-----------------------------------------------------
Known issues
- For systems with symmetry (atoms) this code is known to pick the wrong energy of the HOMO level. I believe it is a bug in Psi4 1.x. as the orbital energy array misses certain energies associated with symmetry-related degenerate levels.
With newer versions of Psi4 the wfn.epsilon_a() lists are not available directly. The fix is to use wfn.epsilon_a().np to obtain the np.array of up-spin orbital energies.
- The IP could be the computed one (listed at the bottom of the output), or the experimental one. Make sure everything is in atomic units. See the example output after this input file.
Ar_IP-calculation_Psi4_1.3_version.psi4
<code | Ar_IP-calculation_Psi4_1.3_version.psi4>
memory 4 gb
molecule monomer {
0 1 # charge and multiplicity.
units angstrom
Ar 0.0 0.0 0.0
# Unfortunately, Psi4 1.x does not write out the orbital energies
# correctly if symmetry is used. So we disable it here by forcing C1.
# This will have the consequence that the calculation will proceed
# more slowly, and, for symmetric systems, the calculation may even
# fail to converge. If symmetry is enabled then the e_HOMO and AC-shift
# values may be incorrect and these should be checked and corrected manually.
Symmetry C1
}
# Define the XC functional here
dft_func = 'pbe0'
set basis aug-cc-pVTZ
set freeze_core TRUE
# get the charge and multiplicity of the monomer
charge = monomer.molecular_charge()
mult = monomer.multiplicity()
set reference rks
EN, wfn = energy(dft_func, return_wfn=True)
# We need to find the HOMO eigenvalue of this system
nelec_a = wfn.nalpha()
epsilon_a = wfn.epsilon_a().np
e_homo = epsilon_a[nelec_a-1]
# now remove and electron and increase the multiplicity.
monomer.set_molecular_charge(charge+1)
monomer.set_multiplicity(mult+1)
set reference uks
EN1 = energy(dft_func)
IP = EN1 - EN
AC_shift = IP + e_homo
psi4.print_out("\n")
psi4.print_out("%s energies and IP \n\n" % (dft_func))
psi4.print_out(" Num Electrons Charge E_int [a.u.] \n")
psi4.print_out("-----------------------------------------------------\n")
psi4.print_out(" N %3d %10.6f \n" % (charge,EN))
psi4.print_out(" N-1 %3d %10.6f \n" % (charge+1,EN1))
psi4.print_out("-----------------------------------------------------\n")
psi4.print_out(" IP %10.6f a.u. \n" % (IP))
psi4.print_out(" e_HOMO %10.6f a.u. \n" % (e_homo))
psi4.print_out(" AC-shift %10.6f a.u. \n" % (AC_shift))
psi4.print_out("-----------------------------------------------------\n")Here are the relevant parts of the output file. The first part is the energy and orbital energies for the species with $N$ electrons (in this case, the neutral Ar atom). We need to use this and not the data for the species with $N-1$ electrons. This is very important! The last part is the summary table produced by the above code. It will be at the end of the output file.
...
...
==> Geometry <==
Molecular point group: c1
Geometry (in Angstrom), charge = 0, multiplicity = 1:
Center X Y Z Mass
------------ ----------------- ----------------- ----------------- -----------------
AR 0.000000000000 0.000000000000 0.000000000000 39.962383123700
Running in c1 symmetry.
Rotational constants: A = '''''''''''''''''' B = '''''''''''''''''' C = '''''''''''''''''' [cm^-1]
Rotational constants: A = '''''''''''''''''' B = '''''''''''''''''' C = '''''''''''''''''' [MHz]
Nuclear repulsion = 0.000000000000000
Charge = 0
Multiplicity = 1
Electrons = 18
Nalpha = 9
Nbeta = 9
...
...
SCF Guess: Core (One-Electron) Hamiltonian.
==> Iterations <==
Total Energy Delta E RMS |[F,P]|
@DF-RKS iter 1: -511.52636091434454 -5.11526e+02 1.20735e-01 DIIS
@DF-RKS iter 2: -516.43497634350854 -4.90862e+00 8.70805e-02 DIIS
@DF-RKS iter 3: -526.80716133498390 -1.03722e+01 2.38723e-02 DIIS
@DF-RKS iter 4: -527.38114680046090 -5.73985e-01 1.15418e-03 DIIS
@DF-RKS iter 5: -527.38145757657207 -3.10776e-04 4.30197e-04 DIIS
@DF-RKS iter 6: -527.38168212177584 -2.24545e-04 3.76718e-05 DIIS
@DF-RKS iter 7: -527.38168337388311 -1.25211e-06 2.39707e-06 DIIS
@DF-RKS iter 8: -527.38168337947411 -5.59100e-09 5.69820e-08 DIIS
Energy and wave function converged.
==> Post-Iterations <==
Orbital Energies [Eh]
---------------------
Doubly Occupied:
1A -115.306790 2A -11.219611 3A -8.742421
4A -8.742421 5A -8.742421 6A -0.992776
7A -0.440597 8A -0.440597 9A -0.440597
Virtual:
10A 0.063800 11A 0.084352 12A 0.084352
13A 0.084352 14A 0.279236 15A 0.279236
...
...
...
...
pbe0 energies and IP
Num Electrons Charge E_int [a.u.]
-----------------------------------------------------
N 0 -527.381683
N-1 1 -526.803729
-----------------------------------------------------
IP 0.577954 a.u.
e_HOMO -0.440597 a.u.
AC-shift 0.137357 a.u.
-----------------------------------------------------
Psi4 stopped on: Monday, 04 January 2021 11:52AM
Psi4 wall time for execution: 0:00:01.05
'''* Psi4 exiting successfully. Buy a developer a beer!From here we get the following:
$\epsilon_{\rm HOMO} = -0.440597$ a.u. This is the energy of the highest-occupied molecular orbital of the Ar atom (species with $N$ electrons). Note that the output file will also contain the relevant data for the Ar$^+$ cation (species with $N-1$ electrons): we do not need any information from this.
$\text{IP} = 0.577954$ a.u. is the computed ionization potential defined as $\text{IP} = E[N-1] - E[N]$, where $E[N]$ is the energy of the species with $N$ electrons, and $E[N-1]$ is the energy of the ionised species.
In the table above, the e_HOMO and AC-shift values are incorrect. The correct e_HOMO if the one we have extracted above. The AC-shift is computed as follows:$\text{AC-Shift} = \text{IP} + \epsilon_{\rm HOMO}$
- = 0.577954 - 0.440597 = 0.137357
= 0.1374 ~ $\text{a.u. to 4-sig figs.}$
- = 0.577954 - 0.440597 = 0.137357
You can also use the experimental IP if you have it. You can find these at the NIST Chemistry WebBookhttps://webbook.nist.gov/chemistry/