Estimation of XPS core level electron binding energies


Hello all!
I would like to estimate the XPS core level binding energies (BEs) by means of ab initio calculation methods (e.g. using Koopmans theorem: BE=-E and delta-SFC method: BE=-E+Rel with Rel=relaxation of the orbitals). Publications I found so far used throughout the commercial Gaussian-package. However, from what I read so far, I am almost sure it is possible to do this with PSI4.
How I would start, and where I have some troubles (as an ex. I chose vinyl alcohol):
(PSI4 is installed on windows 10 subsystem linux via anaconda, and I am using Jupyter for scripting the calculations)
Q1(Initial Question): I do not 100% how psi4 is implemented in python. If I would have implemented the code the most intuitive way for me would to make a class instance for each molecule e.g. benzen = instance_psi4(), and than performing all the operations on the class instance. benzen.geomentry/energy/… So if i perform several calc. on different molecules i do not want to overwrite the data every time. So what is the best way to work?

  1. Def. the structure: vinyl_alc = psi4.geometry(’’’… ‘’’)
  2. Optimize the structure self consistently (i.e. find a minimum local or global): psi4.optimize( ‘scf/6-311G(d,p)’,molecule=vinyl_alc) (e…g using sfc and a gaussian basis set)
    This gives me a total energy of the system (Etot). And in vinyl_alc the optimized structure is now stored.
    To obtain electron core levels binding energies, I need now to calculate also the energy of the ionized system. And for delta SFC it is also necessary to swap the specific core hole (e.g C1s for C or O1s for O in vinyl alc) with the HOMO (Source: DOI: 10.1002/sia.6319)
    Q2: how can i perform this swap in psi4? (I think this step stands for the Rel. part in the equation above and Koopmans theory does not need this step)
  3.‘scf/6-311G(d,p)’,molecule=vinyl_alc) Now I need to perform this on the system with a hole in the C1s at the alpha pos. ( E1) and one more time with a hole in O1s (E2).
    The Binding energies (BE) should result from BE(C1s)=Etot-E1 and BE(O1s)=Etot-E2.
    Q3: However, I am not jet sure how to correctly ionize a specific atom in the molecule.
    I think somehow like this:
    3.1. load relaxed (and orbital swapt structure)
    1 2 #first no is the charge second no. the spin multipl.
    C … …
    Q3 ad: If i want this for a specific atom in the geometry, where/how do I have to place it??

Thx allot!
Many thx for your time


The problem with core-excited calculations is that the wave function will want to collapse the excitation by moving electrons into the vacant core orbital. To avoid this from happening, one needs to do special tricks.

It appears the paper you sent using TCTSNBM uses an obsolete density matrix minimization technique to avoid the collapse from taking place. (If you use DIIS, the core hole will collapse in the very first iteration.)

Other alternatives include freezing the core orbital occupation number to zero (done by setting the diagonal Fock matrix elements to zero), or not by updating the core orbital occupation number. The last alternative should be feasible in Psi4: this is the maximum overlap method (MOM), which has been known since something like the 1970s.

Looks like there is already a sample for how to run MOM


Many thanks for your answer!
I cold not find the abbreviation TCTSNBM in this context. Could you give me a different name or more infos about that method.
I didn’t understand the idea of core level swap (HOMO <-> core hole).
I already read a paper about the MOM method but I also did not fully understand the idea behind it.
Maybe you could give me also a short really easy explanation about the MOM. Do I still have to swap the orbitals using MOM.
As mentioned I installed psi4 as python module and the manual is not that detailed about the functions. So maybe you could help me a bit with the implementation.
How I would implement it (obviously wrong :slight_smile: ) --example toluol:

  1.   Def. the structure: toluol= psi4.geometry(''' coordinates ''') 
  2.   Optimize the structure self consistently psi4.optimize() (e.g. args: ‘scf/6311G(d,p)’,molecule=toluol)

(As I understand optimize uses

  1.  Create a new molecule with a hole in the specific core level and specific atom under investigation 
    toluol_ionized = psi.geometry(‘’’
    C 4.8526 1.6386 0.7866
    1 2 # ionization
    C 4.9336 2.9467 0.3081
    C 3.9063 3.4662 -0.4792
    C 2.7979 2.6773 -0.7894
    C 2.7172 1.3694 -0.3114
    C 3.7445 0.8501 0.4770
    H 1.9881 3.0871 -1.4103
    H 3.6805 -0.1809 0.8539
    C 6.1567 3.8172 0.6511
    H 5.9778 4.3381 1.5685
    H 7.0197 3.1937 0.7575
    H 6.3225 4.5247 -0.1343
    N 3.9912 4.8447 -0.9828
    O 4.9565 5.4972 -0.6987
    O 3.0916 5.2618 -1.6575
    N 5.9352 1.0915 1.6170
    O 6.8744 1.7947 1.8655
    O 5.8359 -0.0363 2.0128
    N 1.5498 0.5379 -0.6380
    O 1.5103 -0.5823 -0.2115
    O 0.6835 1.0130 -1.3178
    (here I don’t know how to create a (hole in the core level of a specific atom)/(position dependent ion).
  2. Set MOM condition using psi4.set_options({‘reference’: ‘uhf’, ‘docc’: [1,0,2,1,0,…],….}, ‘MOM_START’: ?,molecule=toluol_ionized) (docc should swap the orbitals but i don’t how (what docc is doing?) and if it is even necessary)
    (I think MOM_start needs an int number for the iteration?)
  3. Calculate the energy of the ionized molecule‘scf/6-311G(d,p)’,molecule=benzol) and
  4. EC1s(specific Ion)-EC1s(neutral) = BE=Ionization energy

So as I understood the delta SCF method described in the paper: In order to get BE’s with delta-SCF for a molecule e.g H3C-CH2-OH (with interesting XPS spectra C1s and O1s) one has to perform 4 calculation to obtain the spectra.
One Calc. to get all the core hole energies of the neutral molecule (E(neutral), which I get from optimize. One for C1s(1) –( C -C-O) one for C1s(2) –(C- C -O) and one for O1s (C-C- O ), where each time the core hole is in the specific orbital of the atom of interest in the molecule.
However I don’t know if this is also the approach using MOM.
many thanks again for your help!!


I guess it should actually be The Code That Shall Not Be Named, or TCTSNBN for short, which they used in the paper you linked.

Actually, what you want to do is swap the core orbital with the LUMO. This way the occupation shifts from the core orbital to the LUMO, and you get a core excited state. However, the core hole and the LUMO still need to be relaxed, without the core hole collapsing.

The point about MOM is just that instead of filling in the orbitals in increasing energy (i.e. the Aufbau principle), you occupy the orbitals that have the biggest overlap to the previous density. The Aufbau principle would typically fill in the core hole as the first thing it does, whereas the MOM is less likely to place an electron in the core hole.

As for the Python stuff, I have no idea…


Any Psi4 installation can act as a python module (import psi4; python; PsiAPI) or as an executable (psi4 input.dat; PSIthon), so to get started from the MOM input file, you may as well run in the latter mode. In general to convert PSIthon to PsiAPI:

  • molecule {...} --> psi4.geometry("""...""")
  • set {option1 5\noption2 sto-3g} --> psi4.set_options({'option1': 5, 'option2': 'sto-3g'})
  • energy('b3lyp') -->'b3lyp')