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.

IL150-IP-calculation.psi4

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

IL150 snippet

<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

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

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

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