# 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()
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 epsilon_a[] lists are not available and I have not figured out what they are called at present. This version of the input will work but it will not compute the AC-shift automatically. To complete the calculation you will need to find the energy of the HOMO orbital (the highest occupied molecular orbital) from the Psi4 output, and then compute the AC-shift using:

$$\text{AC-Shift} = \text{IP} + \epsilon_{\rm HOMO}.$$

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

<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()
# The epsilon_a list is not accessible like this in Psi4 1.3.2 and later versions:
# epsilon_a = wfn.epsilon_a()
# e_homo    = epsilon_a[nelec_a-1]
# So we set e_homo = 0. This means that the AC-shift must be computed manually
e_homo = 0.0

# 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.000000   a.u.
AC-shift              0.577954   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.}$

• You can also use the experimental IP if you have it. You can find these at the NIST Chemistry WebBookhttps://webbook.nist.gov/chemistry/

AJMPublic/camcasp/psi4-ip-ac-shift (last edited 2021-04-14 12:51:46 by apw109)