Simulating Photoelectron Spectroscopy using Q-Chem with ezDyson and ezSpectrum


In an approximate photoelectron spectroscopy simulation, the norm square of Dyson orbitals are usually taken as a quick estimation of the ionization intensities. However, for quantitative calculations of total and differential ionization cross sections, the photoelectron matrix element should be calculated and evaluating both Dyson orbitals and the wave functions of the leaving photoelectrons is necessary. Q-Chem can calculate Dyson orbitals at the EOM-CC level of theory, which makes calculating absolute photodetachment/photoionization cross sections, photoelectron angular distributions (PADs), and anisotropy parameters (β) using the ezDyson program package possible. In addition, the Franck-Condon factors between the vibrational states of the ionized and neutral species can be calculated by another program package ezSpectrum, therefore the vibrational structures of photoelectron spectroscopy can be studied.

In the following we will illustrate how to simulate photoelectron spectroscopy using Q-Chem with ezDyson and ezSpectrum step-by-step.

Quantum Chemistry Calculations with Q-Chem

Before we run ezDyson and ezSpectrum, Q-Chem calculations need to be carried out in order to produce the necessary Dyson orbitals and other relevant information. The following is an example qchem input file:

Example 1. An EOM-IP-CCSD calculation of a methanol molecule. Geometry optimized at the 𝜔B97X-D/aug-cc-pVTZ level of theory (comments after the symbol “!”).

0 1
C         0.660492       -0.036257    0.000000
O        -0.742798     0.118380    0.000000
H         1.019207      -0.559731      -0.891530
H         1.019207       -0.559731    0.891530
H         1.087065     0.964343    0.000000
H        -1.152635       -0.745352    0.000000
jobtype = SP
correlation = CCSD
BASIS = aug-cc-pVTZ
purecart = 111        ! Required for ezDyson to know if pure or cartesian functions are used.
scf_convergence = 7
ip_states = [1,1]    ! Requests Dyson orbital for states with specific symmetries
ccman2 = true        ! Invoke the ccman2 code
cc_trans_prop = true    ! Request the transitions between all EOM states and the reference state
cc_do_dyson = true    ! Calculate Dyson orbitals
cc_t_conv = 7
eom_davidson_convergence = 6
print_general_basis = true    ! Prints basis set information in output, needed for ezDyson
molden_format = 1    ! Print out Dyson orbitals in the MOLDEN format, needed for visualization

After the calculation is done, the Dyson orbitals can be visualized by IQmol, the graphical frontend software of Q-Chem. Just open the qchem output file with IQmol, and select “Dyson orbital” from ’Surfaces’ on the left side panel; then select the transition for which the Dyson orbital will be plotted. An example calculated Dyson orbital of methanol in IQmol is shown in the following figure:

Fig. 1: Visualization of Dyson orbitals with IQmol.

Photoionization Cross Section Calculation with ezDyson

  1. Obtain ezDyson

Download the current version of ezDyson from the link:

and run (under Linux)

tar -zvxf ezdyson.v4.tar.gz

to extract the files. Usually the pre-compiled executable (exedys) can run under Linux. If not, please copy the makefiles for different platforms provided in the main directory (copy the “mk” file for the proper platform into makefile), and run:

make clean



to compile the code and generate the executable. For more details about the package, please refer to the manual from the link:

  1. Generate ezDyson input file from qchem output file.

Suppose the qchem output file in the last section is “methanol-qchem.out”, and the path to the extracted ezDyson directory is PATHDYSON, one can run the command:

PATHDYSON/ezDyson/xmlsamples/makexml_script/ methanol-ezdy.xml methanol-qchem.out

Here the python script generates the ezDyson input file methanol-ezdy.xml from the qchem output file.

  1. Modify the ezDyson input file if necessary.

The head part of the ezDyson input file methanol-ezdy.xml is as follows:

Example 2. The head of ezDyson input file of methanol (notes after the symbol “!”) .

<?xml version="1.0" encoding="ISO-8859-1"?>
  job = "dyson"

 text = "
      C       0.6610636725        -0.0210744946    -0.0000000000
      O      -0.7435042327        0.1214892640      0.0000000000
      H       1.0242669977        -0.5414443867    -0.8915300000
      H       1.0242669977        -0.5414443867      0.8915300000
      H       1.0790162890         0.9831568008    -0.0000000000
      H      -1.1458984577       -0.7457351718      0.0000000000

   l_max = "5" ! l_max=5 is usually enough for low energy ionization. See ezDyson manual.
   charge_of_ionized_core = "0" > ! 0 is plane wave, 1 is coulomb wave, you can also use <1 as effective charge.
     <k_grid n_points="11" min="0.01" max="1.01" />     ! Kinetic energies for which to compute cross sections, in eV. Should be changed if necessary.

   method= "avg" ! How to average over molecular orientations. Avg is faster as long as l_max is converged.
   method_possible_values="noavg, avg, num" >

   ionization_energy = "10.84" > ! Ionization energy of the molecule. Needs to be specified. Here the number “10.84” is for methanol.
   <laser_polarization x="0.0" y="0.0" z="1.0" /> ! For “avg”, this does nothing. Only for “noavg”.

<lab_xyz_grid> ! Numerical integration grid. This is ok if the integrate norm ~1.
 <axis n_points="201" min="-10.0" max="10.0" />
 <axis n_points="201" min="-10.0" max="10.0" />
 <axis n_points="201" min="-10.0" max="10.0" />

The other parts of the .xml file just copies information from the qchem output file such as the basis functions and Dyson orbitals, which can be left untouched.

  1. Run ezDyson.

Once the .xml file is well prepared, you can run ezDyson with the command: PATHDYSON/ezDyson/exedys methanol-ezdy.xml > methanol-ezdy.out

  1. Result analysis.

In the output file “methanol-ezdy.out” before the content copied from the .xml input file, one can see a summary of the calculated ionization cross sections and anisotropy:

Example 3. Summary of the calculated ionization cross sections and anisotropy of methanol.

E=IE+E_k,eV   Sigma_par     Sigma_perp    Sigma_tot     beta
10.850000      0.000012      0.000009      0.000125      0.190750
10.950000      0.000396      0.000306      0.004223      0.177615
11.050000      0.000990      0.000789      0.010759      0.156109
11.150000      0.001683      0.001385      0.018654      0.134058
11.250000      0.002427      0.002061      0.027432      0.111715
11.350000      0.003189      0.002797      0.036794      0.089172
11.450000      0.003949      0.003580      0.046529      0.066489
11.550000      0.004691      0.004396      0.056480      0.043714
11.650000      0.005405      0.005239      0.066529      0.020886
11.750000      0.006082      0.006100      0.076580       -0.001958
11.850000      0.006718      0.006974      0.086561       -0.024785

Plotting the E (first column) vs. Sigma_tot, you get the “photoionization spectrum”. Plotting E vs. beta you get the information on how the anisotropy parameter varies with energy. PADs are printed in a separate file named pad.dat.

Franck-Condon Factor (FCF) Calculation with ezSpectrum

  1. Run Q-Chem to obtain the optimized geometries and vibrational frequencies and normal mode vectors of both the ground and ionized states.

The example Q-Chem input files for the geometry optimization and frequency calculations of the ionized state of methanol are as the following:

Example 4. Q-Chem input files for the geometry optimization and frequency calculations of the ionized state of methanol.

Geometry optimization

1 2
 C         0.660492       -0.036257    0.000000
 O       -0.742798     0.118380    0.000000
 H         1.019207       -0.559731      -0.891530
 H         1.019207       -0.559731    0.891530
 H         1.087065     0.964343    0.000000
 H       -1.152635       -0.745352    0.000000

   BASIS  =  aug-cc-pVTZ
   JOB_TYPE  =  Optimization
   exchange  =  omegaB97X-D
   sym_ignore = true
   symmetry = false

Geometry optimization

1 2
 C         0.6510503166         0.0426377405     0.0000000201
 O       -0.6979023816         0.1143406096       -0.0000000149
 H         0.9951302666       -0.6240943420       -0.8380349100
 H         0.9951302927       -0.6240943394     0.8380349133
 H         1.0882097264         1.0330700697       -0.0000000061
 H       -1.1410802206       -0.7602077383     0.0000000000

   BASIS  =  aug-cc-pVTZ
   JOB_TYPE  =  Freq
   exchange  =  omegaB97X-D
   MAX_SCF_CYCLES  =  150
   sym_ignore = true
   symmetry = false
   threads 4
   mem_static = 2000
   cc_memory = 14000
  1. Obtain ezSpectrum

Download the current version of ezSpectrum from the link:
and run (under Linux)
tar -zvxf ezSpectrum.v3.0.tar.gz
to extract the files. Usually the pre-compiled executable (ezSpectrum.64) can run under Linux.

  1. Prepare the ezSpectrum input file.

Copy the qchem output files of the frequency calculations of both the ground and ionized states of methanol into your working directory. Suppose these files are called “methanol-S0.out” and “methanol-D1.out”. Suppose the path to the ezSpectrum directory is PATHSPEC, copy two file from the ezSpectrum directory to your working directory by running the commands:
cp PATHSPEC/ezSpectrum/atomicMasses.xml .
cp PATHSPEC/ezSpectrum/ .
Then run the python script to generate the input methanol-ezSpec.xml file for ezSpectrum:
./ methanol-ezSpec.xml methanol-S0.out methanol-D1.out
Note: do not change the order of the two qchem output files in the command above.

  1. Modify the .xml input file if necessary.

Open the methanol-ezSpec.xml file and change the temperature on the fifth line from 300 to 0. Then change the ionization potential (ip), which is right under the “<target_state>” section of the xml input file. By default, the xml file has 1; change the 1 to 10.84, the actual experimental IP for the methanol.

  1. Run ezSpectrum.

Run the command:
PATHSPEC/ezSpectrum/ezSpectrum.64 methanol-ezSpec.xml > methanol-ezSpec.out

  1. Use the script xsecFCF to combine the results of ezDyson and ezSpectrum and get the photoelectron data for plotting.

The script xsecFCF can be found in PATHDYSON/xmlsamples/xsecFCF. It requires two files called FCFs.dat and xsec.dat. The first file can be copied from an output file from the ezSpectrum calculation, usually called filename.xml.spectrum_parallel or filename.xml.spectrum_duschinsky (depending on the options used in ezSpectrum). This file needs to be modified by adding only one line in the beginning of the file indicating the total number of FCF transitions included in the file. The second file, xsec.dat, is just a segment copied from the ezDyson output from the section titled “Total cross-section / E (in a.u.) vs. Ek (in eV)”, with again one additional line indicating the number of energies for which cross sections were computed. Note that for ezSpectrum, ionization energies (first column) are reported as negative numbers, so you need to convert them all to positive values first for the script to work correctly. Put FCFs.dat and xsec.dat into the same folder as xsecFCF and run the script:

./xsecFCF > xsecFCF.out

In the end of the output file xsecFCF.out you will see the result summary:

Example 5. Photoionization cross sections with Franck-Condon factors of methanol.

Ionization energy: 10.84000
cross-sections after accounting for Franck-Condon factors:
E, eV     xsec, a.u.
10.85000  0.00003
10.95000  0.00093
11.05000  0.00330
11.15000  0.00722
11.25000  0.01267
11.35000  0.01938
11.45000  0.02703
11.55000  0.03538
11.65000  0.04420
11.75000  0.05330
11.85000  0.06255

The results can be shown in the following picture:

Fig. 1: Absolute cross sections for photoionization of methanol with different screened nuclear charges (Z). Adapted from J. Phys. Chem. Lett., 6, 4532 (2015).


I want to thank Prof. Samer Gozem for his help in the process of preparing this howto.

1 Like

Excellent instructions!

I have two minor comments:

  1. For more details and updates, see ezDyson and ezSpectrum manuals (available from the website). They are very detailed and also include background theory.
  2. In principle, you can use other orbitals in ezDyson calculations, such as Hartree-Fock or Kohn-Sham ones, or Dyson orbitals from ADC calculations.