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. psi4.energy(‘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)
    geometry(’’’
    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
maths

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 psi4.energy(‘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!!
BR
maths

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 input.py; 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') --> psi4.energy('b3lyp')

Many Thanks for the answers!

I think now I almost know how to implement it!

I have some additional “understanding Questions”.

  1. I am not familiar with the orbital nomenclature used in psi 4. For XPS we use the atomic core level (e.g. C - C1s, O1s, O2s…). Psi4 uses of course the MOs. I do not know how to attribute each MO to the atomic core level (e.g. the C1s energies of benzol are obviously given as: 1Ag, 1B3u, 1B2u, 1B1g, 2Ag, 2B2u ) and even more to the correct atom (how they are sorted). In benzol they are equal, but how is it e.g. for tnt featuring three NO2 groups and a -CH3 group. From lith. I know that it has to do with point symmetry and spin, but I could not find a good explanation on it.
    Beside benzol, I already tried to optimize the structure of tnt, but it gives strange results:
    input:
    tnt = psi4.geometry(’’’
    C 4.8526 1.6386 0.7866
    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
    units angstrom
    ‘’’)
    psi4.optimize(‘scf/6-31g**’, molecule=tnt)

this gives only A orbitals
comparison:
For TNT
==> Pre-Iterations <==


Irrep   Nso     Nmo     Nalpha   Nbeta   Ndocc  Nsocc

 A        265     265       0       0       0       0

Total     265     265      58      58      58       0

   1A    -20.639122     2A    -20.638605     3A    -20.636646  
   4A    -20.636587     5A    -20.633679     6A    -20.633510  
   7A    -15.899000     8A    -15.898586     9A    -15.896068  
  10A    -11.363120    11A    -11.361702    12A    -11.361546  
  13A    -11.355752    14A    -11.352386    15A    -11.352260  
  16A    -11.234695    17A     -1.668072    18A     -1.667601  
  19A     -1.664384    20A     -1.450530    21A     -1.449943  
  22A     -1.445561    23A     -1.274092    24A     -1.172317  
  25A     -1.168165    26A     -1.053417    27A     -1.024789  
  28A     -1.012802    29A     -0.946826    30A     -0.873982  
  31A     -0.863392    32A     -0.826574    33A     -0.796022  
  34A     -0.783772    35A     -0.781998    36A     -0.773287  
  37A     -0.767190    38A     -0.764606    39A     -0.757666  
  40A     -0.754640    41A     -0.721514    42A     -0.711405  
  43A     -0.622044    44A     -0.620294    45A     -0.606262  
  46A     -0.580878    47A     -0.573568    48A     -0.542807  
  49A     -0.529508    50A     -0.524357    51A     -0.518670  
  52A     -0.518076    53A     -0.487131    54A     -0.481305  
  55A     -0.481205    56A     -0.464211    57A     -0.437337  
  58A     -0.425390

For benzol
using
benz = psi4.geometry("""
C 0.000 1.396 0.000
C 1.209 0.698 0.000
C 1.209 -0.698 0.000
C 0.000 -1.396 0.000
C -1.209 -0.698 0.000
C -1.209 0.698 0.000
H 0.000 2.479 0.000
H 2.147 1.240 0.000
H 2.147 -1.240 0.000
H 0.000 -2.479 0.000
H -2.147 -1.240 0.000
H -2.147 1.240 0.000
units angstrom
“”")
I got:
==> Pre-Iterations <==


Irrep   Nso     Nmo     Nalpha   Nbeta   Ndocc  Nsocc

 Ag        31      31       0       0       0       0
 B1g       23      23       0       0       0       0
 B2g        7       7       0       0       0       0
 B3g       11      11       0       0       0       0
 Au         7       7       0       0       0       0
 B1u       11      11       0       0       0       0
 B2u       31      31       0       0       0       0
 B3u       23      23       0       0       0       0

Total     144     144      21      21      21       0

   1Ag   -11.235152     1B3u  -11.234600     1B2u  -11.234572  
   1B1g  -11.233377     2Ag   -11.233348     2B2u  -11.232763  
   3Ag    -1.147365     3B2u   -1.012503     2B3u   -1.012470  
   2B1g   -0.821931     4Ag    -0.821929     5Ag    -0.707274  
   4B2u   -0.642661     3B3u   -0.616421     5B2u   -0.586390  
   4B3u   -0.586350     1B1u   -0.499085     6Ag    -0.493154  
   3B1g   -0.493044     1B2g   -0.334290     1B3g   -0.334281

Does it have to do with the choice of the origin?.. Does it influence the energies (does not look like)

  1. Actually I do not understand why I have to swap the orbitals using MOM. Why is it not only removing the core electron of interest. In the paper I cited, they swap the HOMO with the core level. But the I do not understand the drawing in the paper. It looks like that the HOMO has now the energy of the core hole but it is drawn on top of the energy MO-plot (How is this possible, the only axis in a MO plot is Energy. If i swap only this two i am getting something like: -20 … -11 …-0.5…-0.4 --> -0.4 … -11 …-0.5…-20 ??).

How I would implement it:

a) def struc.: stuc = psi4.geometry("""…""")
b) E0, WF = psi4.optimize(‘scf/6-311G(d,p)’, …)
This gives the energy of the relaxed molecule E0 (C1s, O1s,…)
c) def MOM and docc,socc: psi4.set_options({‘MOM_start’: (e.g. 10), ‘docc’: ?, ‘socc’:?})
The problem is at this step docc and socc are not energy dependent (e.g. docc= [Ag, Bg1,…] but there are several Ag with different energies). What I whant is removing an electron from a double occ. orbital (docc - 1), but with a specific energy (e.g. from 1Ag with energy -11.235152 …corresponding to a specific atom in the molecule) and therefore i get a single occ. orbital socc + 1.
d) def. new structure with one electron less:
stuc_ion = psi4.geometry("""
1 1
struc_after_relax
“”")
e) EI, WF = psi4.energy(‘scf/6-311G(d,p)’, molecule=stuc_ion,…)

For different carbon atoms C1…Cn:
E(C1s(C1))=EI_c1(C1s(C1)) - E0(C1s(C1))

E(C1s(Cn))=EI_cn(C1s(Cn)) - E0(C1s(Cn))

again, many thanks again!
BR
maths

If anything is unclear in my question plz. tell me and I will try to improve it.
many thx!
maths

The post is long and complicated. Maybe we go step by step?

The A or Ag, B1g, etc. labels are for the point group of the molecule. Having only A means C1.

I dont think psi4 has a AO decomposition of the MOs implemented. That might be a hurdle for what you want to do.

Thx for your answer.

  1. Ok the topic of point groups is quit complicated. But what is the reason i getting only A in the case of tnt and different for benzol?

  2. Ok… but it should be possible to relate the MO to the corresponding position in the molecule? I think the rest could may be done by comparison with measurements (at least for most cases).

to 1) Benzene has D6h symmetry, but psi4 can only do D2h. If you look here: http://symmetry.jacobs-university.de/cgi-bin/group.cgi?group=602&option=4 you will find the orbital labels and their symmetry elements. TNT has point groups C1 in your case, which means only A.

to 2) MOs are delocalized and often have multiple centers (=atoms) with large coefficients. Just try to visualize the MOs for TNT. Core-orbitals might be localized enough.

Your grammar doesn’t help. I have a very difficult time reading your posts. Please use proper English to the best of your ability.

As for question (1), spin is irrelevant here. It sounds like you’ve been reading about term symbols for electronic states. That’s related to the symbols you’re asking about, but term symbols are not the place to start.

The condensed version is:
Molecules have certain “symmetry operations,” e.g., you can rotate a benzene molecule by 60 degrees and get the same molecule. In fact, if you rotate it by 60 degrees 6 times, you get every atom in the exact same place as it was before. Wavefunctions obey similar rules, and we can enforce this computationally by making the orbitals obey those rules. Different sets of symmetry operations have different possibilities for what those rules are. 1B1g means you are talking about the first orbital that obeys the B1g rule, and 2B2u means you are talking about the second orbital that obeys the B2u rule.

For (2), try setting print_MOs true in your input file that should print out the decomposition of MOs in terms of AOs.

Hi there,

if you want to calculate core electron binding energies of organic molecules, you may want to check my PSI4 plugin https://github.com/Masterluke87/psixas. Though it is more written to calculate NEXAFS spectra (using TP-DFT), you can use it to calculate XPS positions as well. For example, the 1s CEBE of water can be obtained by using:

import psixas

molecule {
O 0.27681793323501 0.00000014791107 0.00000000000000
H 0.86159097690242 0.76505117501585 0.00000000000000
H 0.86159108986257 -0.76505132292693 0.00000000000000
symmetry c1
}

set {
basis def2-TZVP
}

set scf {
reference uks
scf_type MEM_DF
}

set psixas {
print 1
prefix WATER
MODE GS+EX
ORBS [0 ]
OCCS [0.0]
SPIN [b ]
DAMP 0.2
OVL [T]
FREEZE [T]
}
energy(‘psixas’,functional=‘PBE’)

Then:

grep “FINAL” output.dat
FINAL GS SCF ENERGY: -76.37667773 [Ha]
FINAL EX SCF ENERGY: -56.56555074 [Ha]

And you can calculate the XPS position, which is in this case 539.08 eV. Hope that helps!

1 Like

Interesting. What is the best way to install your plugin – I am using the conda install of psi4. Thanks.