Simulating core photoelectron spectroscopy using the equation-of-motion coupled-cluster method for ionization potentials with the frozen-core and core-valence-separation approximations ((fc)-CVS-EOM-IP-CC)

X-ray lasers can trigger core ionization in a molecular system. The leaving electrons feel the electronic interaction in the molecular environment and thus carry the information of the electronic structure. Thanks to energy separation and locality of core ionizations, core photoelectron spectroscopy with X-ray has excellent element selectivity and spatial resolution, which makes it a powerful experimental tool in modern chemistry, physics, materials science and biology.

As correlation plays an important role in core ionization, high-level correlated methods such as coupled cluster (CC) theory are usually used to simulate core photoelectron spectroscopy (PES). CC and its equation-of-motion (EOM) variant are among the most accurate and versatile tools in studying core level spectroscopy. To consider ionizations, the EOM-IP-CC approach is adopted, in which the ground state is treated at the CC level and the ionized states are accessed by applying an ionization operator to the ground-state wavefunction. In order to reduce the computational cost, the CC calculations are done under the frozen core (fc) approximation, which means the ground state amplitudes and multipliers with occupied indices referring to a core orbital are set to zero and hence the core-core and core-valence correlation in the ground state is neglected. Since we only care about core ionizations here, the EOM states are represented in a restricted Fock space spanned by the many-electron determinants with at least one occupied core orbital index (the core-valence-separation approximation). The above sketches the main profile of the fc-CVS-EOM-CC approach.

Both ground states and excited state core PES can be studied by fc-CVS-EOM-CCSD. In ground state studies, the ionized states are obtained by fc-CVS-EOM-IP-CCSD, and the corresponding left- and right-Dyson orbitals are calculated, as the transformed Hamiltonian in the CC equations are not hermitian. The ionization intensities can be estimated using the product of norms of the left- and right-Dyson orbitals. The case of excited state core PES is different from that of ground states, since the final states of core ionization have double excitation character (2-hole-1-particle), which are not well described by at the level of EOM-CCSD. A scheme to circumvent this difficulty is to first obtain an excited triplet state with a valence hole through MOM, and then use this excited triplet state as a reference to calculate the core-ionized states and the corresponding left- and right-Dyson orbitals. Core ionizations from both alpha and beta channels show similar results.

In the following examples we show how to calculate the ground state PES of adenine and the excited state PES of uracil through fc-CVS-EOM-IP-CCSD.

Example. Simulating the C1s PES of the ground state of adenine. Geometry optimized at the MP2/cc-pVTZ level of theory.

Q-Chem input and comments (after “!”. Comments outside of the $rem section should be removed in order to run the qchem calculation.):

$molecule
0 1
N       1.663948    -1.127531    0.000000
N      -0.643862    -1.804088    0.000000
N      -0.636159     1.773662    0.000000
N      -2.146840     0.119893    0.000000
N       2.353671     1.084678    0.000000
C       0.668787    -2.022720    0.000000
C      -0.919981    -0.499328    0.000000
C       0.000000     0.546991    0.000000
C       1.356490     0.175959    0.000000
C      -1.907009     1.474843    0.000000
H      -3.039931    -0.341646    0.000000
H      -2.712908     2.189994    0.000000
H       0.978346    -3.060960    0.000000
H       2.144452     2.065199    0.000000
H       3.305020     0.766647    0.000000
$end

$rem
print_general_basis = true
jobtype = sp
method = eom-ccsd
cvs_ip_states    [5,0]    ! calculate 5 ionized states in the first irrep, none for the second.
basis = general    ! a general 6-311(2+,+)G** basis set is used (uncontracted core on C)
n_frozen_core = fc    ! the frozen core approximation
cc_trans_prop = true    ! calculate the transition properties of CC states
cc_do_dyson = true    ! calculate Dyson orbitals
print_general_basis = true
molden_format = true    ! output the orbitals in the MOLDEN format
$end

$basis
H    0
6-311(2+,+)G**
****
O    0
6-311(2+,+)G**
****
N    0
6-311(2+,+)G**
****
C    0
S    1    1.000000
   4.56324000E+03    1.00000000E+00    ! uncontracted core on C
S    1    1.000000
   6.82024000E+02    1.00000000E+00
S    1    1.000000
   1.54973000E+02    1.00000000E+00
S    1    1.000000
   4.44553000E+01    1.00000000E+00
S    1    1.000000
   1.30290000E+01    1.00000000E+00
S    1    1.000000
   1.82773000E+00    1.00000000E+00
SP   3    1.000000
   2.09642000E+01    1.14660000E-01   4.02487000E-02
   4.80331000E+00    9.19999000E-01   2.37594000E-01
   1.45933000E+00   -3.03068000E-03   8.15854000E-01
SP   1    1.000000
   4.83456000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   1.45585000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   4.38000000E-02    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   1.31927711E-02    1.00000000E+00   1.00000000E+00
D    1    1.000000
   6.26000000E-01    1.00000000E+00
****
$end

The simulation results are shown in the following figure:

Example. Simulating the O1s PES of the excited state (S1) of uracil. Geometry optimized at the CCSD/aug-cc-pVDZ level of theory.

Q-Chem input and comments (after “!” in. Comments outside of the $rem section should be removed in order to run the qchem calculation.):

$comment
First step: a HF single point calculation of the neutral singlet state
$end

$molecule
0 1
C        -0.0042112405          -1.7549930536        0.0000000000
H        -2.0949460216          -1.4184138214        0.0000000000
H        -0.1300971688          -2.8358204634        0.0000000000
H          0.0374086733    2.0243048927        0.0000000000
C        -1.2210197500    0.4010705948        0.0000000000
O        -2.2535781024    1.0567627444        0.0000000000
C          1.2399442445    0.2839850904        0.0000000000
O          2.3361887279    1.0691923037        0.0000000000
N          0.0349480223    1.0125648826        0.0000000000
N        -1.1853832587          -0.9770870742        0.0000000000
C          1.2251302704          -1.1165150648        0.0000000000
H          2.1596039787          -1.6767971209        0.0000000000
$end

$rem
jobtype sp
basis general
method hf
n_frozen_core = 2    ! freeze 2 O1s electrons
$end

$basis
H    0
6-311(2+,+)G**
****
O    0
S    1    1.000000
   8.58850000E+03    1.00000000E+00
S    1    1.000000
   1.29723000E+03    1.00000000E+00
S    1    1.000000
   2.99296000E+02    1.00000000E+00
S    1    1.000000
   8.73771000E+01    1.00000000E+00
S    1    1.000000
   2.56789000E+01    1.00000000E+00
S    1    1.000000
   3.74004000E+00    1.00000000E+00
SP   3    1.000000
   4.21175000E+01    1.13889000E-01   3.65114000E-02
   9.62837000E+00    9.20811000E-01   2.37153000E-01
   2.85332000E+00   -3.27447000E-03   8.19702000E-01
SP   1    1.000000
   9.05661000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   2.55611000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   8.45000000E-02    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   2.54518072E-02    1.00000000E+00   1.00000000E+00
D    1    1.000000
   1.29200000E+00    1.00000000E+00
****
N    0
S    1    1.000000
   6.29348000E+03    1.00000000E+00
S    1    1.000000
   9.49044000E+02    1.00000000E+00
S    1    1.000000
   2.18776000E+02    1.00000000E+00
S    1    1.000000
   6.36916000E+01    1.00000000E+00
S    1    1.000000
   1.88282000E+01    1.00000000E+00
S    1    1.000000
   2.72023000E+00    1.00000000E+00
SP   3    1.000000
   3.06331000E+01    1.11906000E-01   3.83119000E-02
   7.02614000E+00    9.21666000E-01   2.37403000E-01
   2.11205000E+00   -2.56919000E-03   8.17592000E-01
SP   1    1.000000
   6.84009000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   2.00878000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   6.39000000E-02    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   1.92469880E-02    1.00000000E+00   1.00000000E+00
D    1    1.000000
   9.13000000E-01    1.00000000E+00
****
C    0
S    1    1.000000
   4.56324000E+03    1.00000000E+00
S    1    1.000000
   6.82024000E+02    1.00000000E+00
S    1    1.000000
   1.54973000E+02    1.00000000E+00
S    1    1.000000
   4.44553000E+01    1.00000000E+00
S    1    1.000000
   1.30290000E+01    1.00000000E+00
S    1    1.000000
   1.82773000E+00    1.00000000E+00
SP   3    1.000000
   2.09642000E+01    1.14660000E-01   4.02487000E-02
   4.80331000E+00    9.19999000E-01   2.37594000E-01
   1.45933000E+00   -3.03068000E-03   8.15854000E-01
SP   1    1.000000
   4.83456000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   1.45585000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   4.38000000E-02    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   1.31927711E-02    1.00000000E+00   1.00000000E+00
D    1    1.000000
   6.26000000E-01    1.00000000E+00
****
$end

@@@

$comment
Second step: read in the MOs from the previous step, and obtain a neutral triplet excited state through MOM, then run CVS-EOM-IP-CCSD to calculate the ionized states
$end

$molecule
0 3    ! triplet  configuration is considered as reference for core ionization
C        -0.0042112405       -1.7549930536    0.0000000000
H        -2.0949460216       -1.4184138214    0.0000000000
H        -0.1300971688       -2.8358204634    0.0000000000
H         0.0374086733        2.0243048927    0.0000000000
C        -1.2210197500        0.4010705948    0.0000000000
O        -2.2535781024       1.0567627444     0.0000000000
C         1.2399442445        0.2839850904    0.0000000000
O         2.3361887279       1.0691923037     0.0000000000
N         0.0349480223       1.0125648826     0.0000000000
N       -1.1853832587       -0.9770870742     0.0000000000
C        1.2251302704       -1.1165150648    0.0000000000
H        2.1596039787       -1.6767971209    0.0000000000
$end

$rem
job_type  =  sp
basis  = general
print_general_basis = true
method  =  eom-ccsd
Scf_guess = read    ! read in the MOs from the previous step 
mom_start = 1        ! invoke MOM to obtain an excited state
cvs_ip_states = [4,0]    ! calculate 4 ionized states in the first irrep, none for the second.
cvs_eom_shift = 19500    ! specifies energy shift in CVS-EOM calculations.
n_frozen_core = 2    ! freeze 2 O1s electrons
cc_trans_prop = true    ! calculate the transition properties of CC states
cc_do_dyson = true    ! calculate Dyson orbitals
molden_format = true    ! output the orbitals in the MOLDEN format
$end

$basis
H    0
6-311(2+,+)G**
****
O    0
S    1    1.000000
   8.58850000E+03    1.00000000E+00
S    1    1.000000
   1.29723000E+03    1.00000000E+00
S    1    1.000000
   2.99296000E+02    1.00000000E+00
S    1    1.000000
   8.73771000E+01    1.00000000E+00
S    1    1.000000
   2.56789000E+01    1.00000000E+00
S    1    1.000000
   3.74004000E+00    1.00000000E+00
SP   3    1.000000
   4.21175000E+01    1.13889000E-01   3.65114000E-02
   9.62837000E+00    9.20811000E-01   2.37153000E-01
   2.85332000E+00   -3.27447000E-03   8.19702000E-01
SP   1    1.000000
   9.05661000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   2.55611000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   8.45000000E-02    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   2.54518072E-02    1.00000000E+00   1.00000000E+00
D    1    1.000000
   1.29200000E+00    1.00000000E+00
****
N    0
S    1    1.000000
   6.29348000E+03    1.00000000E+00
S    1    1.000000
   9.49044000E+02    1.00000000E+00
S    1    1.000000
   2.18776000E+02    1.00000000E+00
S    1    1.000000
   6.36916000E+01    1.00000000E+00
S    1    1.000000
   1.88282000E+01    1.00000000E+00
S    1    1.000000
   2.72023000E+00    1.00000000E+00
SP   3    1.000000
   3.06331000E+01    1.11906000E-01   3.83119000E-02
   7.02614000E+00    9.21666000E-01   2.37403000E-01
   2.11205000E+00   -2.56919000E-03   8.17592000E-01
SP   1    1.000000
   6.84009000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   2.00878000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   6.39000000E-02    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   1.92469880E-02    1.00000000E+00   1.00000000E+00
D    1    1.000000
   9.13000000E-01    1.00000000E+00
****
C    0
S    1    1.000000
   4.56324000E+03    1.00000000E+00
S    1    1.000000
   6.82024000E+02    1.00000000E+00
S    1    1.000000
   1.54973000E+02    1.00000000E+00
S    1    1.000000
   4.44553000E+01    1.00000000E+00
S    1    1.000000
   1.30290000E+01    1.00000000E+00
S    1    1.000000
   1.82773000E+00    1.00000000E+00
SP   3    1.000000
   2.09642000E+01    1.14660000E-01   4.02487000E-02
   4.80331000E+00    9.19999000E-01   2.37594000E-01
   1.45933000E+00   -3.03068000E-03   8.15854000E-01
SP   1    1.000000
   4.83456000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   1.45585000E-01    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   4.38000000E-02    1.00000000E+00   1.00000000E+00
SP   1    1.000000
   1.31927711E-02    1.00000000E+00   1.00000000E+00
D    1    1.000000
   6.26000000E-01    1.00000000E+00
****
$end

$occupied    ! specify the orbital occupation
1:29 41    ! alpha orbitals no. 1-29 and 41 are occupied, which is specific to state S1
1:26 28 29    ! beta orbital no. 1-29 are occupied except no. 27, which is specific to state S1
$end

The simulation results are shown in the following figure:

2 Likes

This howto is prepared with the help from Marta Lopez Vidal at Technical University of Denmark.

Dear yzhang,

thank you for this very interesting tutorial.
May I ask you two questions:

  1. What is the difference to the fc-CVS-EOM-EE-CC approach presented in the manual and in Webinar 36?
  2. I understood that the number of frozen core is very important and edge sensitive. How does one account for this if there are (as it is in your case) several degenerate core orbitals. Are by choosing n_frozen_core = fc all C core orbitals included ?

Best
Michael

Dear Michael,
My answers to your questions are as the following:

  1. I think the fc-CVS-EOM-EE/IP-CC approach were covered in webinar 36. As I presented in the introduction section of the tutorial above, the frozen core (fc) approximation means the ground state amplitudes and multipliers with occupied indices referring to a core orbital are set to zero and hence the core-core and core-valence correlation in the ground state is neglected. And in the core-valence-separation approximation the EOM states are represented in a restricted Fock space spanned by the many-electron determinants with at least one occupied core orbital index.
  2. Yes, all C1s core orbitals are frozen is you set n_frozen_core = fc.

Dear yzhang,

thank you for your prompt reply.
I am sorry I think I did not express myself very clear. I meant what is the difference in calculation ground state XAS via the EOM-EE or EOM-IP approach. Why should one choose one or the other. In the webinar they only showed the calculation of ground state XAS via EOM-EE and so does the manual 5.1. Thus, I am very interested in this new way of calculating XAS with the EOM-IP version you outlined above. Do you expect differences with the two approaches? Or are there computational cost benefits?
Regarding my second question, I understand that all C1s orbitals are frozen and also all N1s orbitals I assume. So overal 6 x 1s orbitals. My question is now which 1s - valence transition are included if I set n_frozen_core = fc. According to the manual it is edge sensitive and only the highest lying frozen orbital is incorporated. So in case I wanted the nitrogen edge, I would need to unfreeze the C1s orbitals. In your example we have several Carbon 1s orbitals. They will not all have the same orbital energy, since their chemical environment is different. So by freezing all 1s orbitals, and since only the 1s orbital highest in energy is regarded for the transition, I was just wondering if you still include XAS transitions from all 4 Carbon atoms or only from the one with the 1s highest in energy?

Thank you very much,
Michael

hi, Michael,
I think in webinar 36 it was stated that CVS-EOM-IP is for XPS and XES, while CVS-EOM-EE is for XAS and RIXS. In the tutorial above I was talking about XPS, so I used CVS-EOM-IP.
For your second question, I will talk to the original author of the code and get back to you.
Thanks
Yu

I have consulted the author of the code and get the following answer to your second question: you will get transitions from all 4 carbon atoms if you do not ask for too less ionized states in your calculation (for example, only 1 or 2 ionized states).

The difference between using EOM-IP and EOM-EE is the type of the target states you are accessing. For example, if you want to compute core and valence ionization energies and XES transitions in benzene molecule, you will need to use EOM-IP and CVS-EOM-IP with a closed-shell neutral reference. If you would like to compute XAS and RIXS of the neutral benzene, then you need to use EOM-EE and fc-EOM-EE, again, with the neutral reference. In all cases, you should use FC=6 to freeze all 6 core orbitals (which is the default setting for both calculations).

You may find this review paper useful:
http://arjournals.annualreviews.org/eprint/ukNyyVKWir2agwdJ5RPx/full/10.1146/annurev.physchem.59.032607.093602

It explains the EOM-CC approach approach and gives examples.

And as a general introduction to the EOM-CC methods in Q-Chem, you can watch this webinar:
https://www.q-chem.com/webinars/3/