Many Thanks for the answers!
I think now I almost know how to implement it!
I have some additional “understanding Questions”.
- 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)
- 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