Skip to content

2. Core Hole Localisation: Azulene

In this part, you will learn how to perform core-hole calculations for highly symmetrical systems. In such cases, the core-state wave functions are not localised on a single atom but rather form linear combinations of the single core-orbitals. However, The experimental evidence shows that the core-hole is localised on one atom. Therefore, we must also force the core-hole of our excited wave function to localise on one single atom.

As an example, to show the principal approach to this problem, we chose the azulene molecule. This molecule only contains carbon and hydrogen atoms, but due to its non-alternant topology, the chemical shift between the carbon atoms is sufficiently large to see differences. Azulene possesses C2v point symmetry, and several pairs of atoms are symmetry equivalent. There is experimental data for the multilayer of azulene, to which we can compare our results later.1


Scheme 3. Structure of the azulene molecule and numbering scheme used in this tutorial.

The following geometry.in file can be used for setting up the azulene calculations:

atom    0.00000000      0.00000000      2.15610643      C
atom    1.15055599      0.00000000      1.34809461      C
atom    -1.15055599     0.00000000      1.34809461      C
atom    0.74965762      0.00000000      -0.00057279     C
atom    -0.74965762     0.00000000      -0.00057279     C
atom    1.59583767      0.00000000      -1.10633748     C
atom    -1.59583767     0.00000000      -1.10633748     C
atom    1.26699680      0.00000000      -2.46635052     C
atom    -1.26699680     0.00000000      -2.46635052     C
atom    0.00000000      0.00000000      -3.05928669     C
atom    0.00000000      0.00000000      3.24440762      H
atom    2.18076818      0.00000000      1.69493655      H
atom    -2.18076818     0.00000000      1.69493655      H
atom    2.66541089      0.00000000      -0.87579776     H
atom    -2.66541089     0.00000000      -0.87579776     H
atom    2.11174769      0.00000000      -3.15756315     H
atom    -2.11174769     0.00000000      -3.15756315     H
atom    0.00000000      0.00000000      -4.15206631     H

In general, we follow the same steps as in Part 1.

Ground State Calculation

First, we perform the single point ground state calculation:

  1. Use the functional PBE and add a tight-tier-2 basis set for carbon and hydrogen.
  2. Add the basis set augmentation to the basis set definition of carbon
  3. Take care to use the keyword spin collinear
  4. Once again, add output cube eigenstate n with n = 1 to 10 to the input file to check for the localisation of the wave functions
  5. Run the calculation

By looking at the exported cube files, we can see that the wave functions of the core states are now mixed and not localised properly at single atoms anymore.

Using deltascf_projector and Breaking the Symmetry

We can again use the deltascf_projector command, but due to the symmetry equivalent atoms in azulene, we must now break the symmetry of the wave function first to force a proper localisation.

To break the symmetry, we need to use a trick: We differentiate the atom, where the wavefunction is to be localised, by increasing its nuclear charge in an intermittent calculation and carrying the localised wave function over to a later calculation with the physically correct core charge. This later calculation then keeps to the localised wave function but yields the sought-after physically meaningful energy.

It turns out that the best results in the localisation are obtained if we perform two so-called initialisation calculations before the final hole calculation, yielding a three-step procedure:

  1. The init_1 calculation with the nuclear charge increased by +0.1 and the system charge set to +0.1. A full SCF cycle is performed and a restart_file is written)

  2. The init_2 calculation with the nuclear charge increased by +0.1 and the system charge set to +1.1. deltascf_projector is used to initialise the core hole at the correct atom. It is enough to calculate 1 SCF step, the restart_file from the init-a calculation is read and overwritten)

  3. The hole calculation with the nuclear charge back to the original value and the system charge set to 1.0. deltascf_projector is used to de-occupy the core-hole, a full SCF cycle is performed and the restart_file from the init_2 calculation is read. This calculation produces the core-hole energy used to determine the binding energy according to equation (3)

To determine which atom the core-hole is to be localised at, it has to be differentiated in the geometry.in file, e.g. by calling it C1. As described above, the same basis set is used as for all the other atoms, but the control.in file must now contain a second definition for the carbon basis set, called C1. We will use the C1 basis set definition to change the nuclear charge only for this specific atom.

Note on Boys Localisation: We are currently investigating using Boys' Localisation here to converge the core hole rather than using an additional partial nuclear charge. Until this is thoroughly tested, the partial core hole should suffice as a method in most circumstances.

Initialisation and Hole Calculations for deltascf_projector

Following the procedure laid out in the last section, now we perform the initialisation and hole calculations:

  1. Create 3 folders for each atom: init_1 , init_2, and hole
  2. Create a special geometry.in file for each atom where that atom is labelled as C1

Running the init_1 calculation requires the following settings:

  1. In the C1 basis set information, the nuclear charge is increased by 0.1 by setting nucleus from 6 to 6.1
  2. To compensate this increased nuclear charge, the valence electron count of the C1 basis set is also increased by 0.1 by setting valence 2 p from 2 to 2.1.
  3. In the general part of the control.in file, the keyword charge = 0.1 is added to ensure charge neutrality is maintained
  4. add restart_write_only restart_file and force_single_restartfile .true. to write out a restart file
  5. Run all ten init_1 calculations

Now we move on to the init_2 calculations

  1. Copy the restart file from each init_1 calculation into the corresponding init_2 folder
  2. the changes in the C1 basis set definition remain the same
  3. we change the system charge to charge = 1.1 to remove one complete electron in addition to the compensation charge
  4. we add deltascf_projector 1 1 0.0 1 10 to introduce the core hole. Note that we do not need to change the KS_state value, because for this system, the added nuclear charge will always force the corresponding 1s state to the lowest energy.
  5. set sc_iter_limit to 1, because we only need to initialise the correct wave function
  6. add restart restart_file and force_single_restartfile .true. to both load the restart file from the previous init_1 calculation and overwrite it for use in the next calculation
  7. Run all ten init_2 calculations

Finally, we can perform the hole calculations

  1. Copy the restart_file from the init_2 calculations into the corresponding hole folder
  2. Remove the changes in the C1 basis set definition. nucleus goes back to 6 and valence 2 p goes back to 2
  3. Change the system charge to charge = 1.0 to have a regularly charged system
  4. deltascf_projector 1 1 0.0 1 10 to introduce the core hole remains the same.
  5. Set sc_iter_limit back to a normal value, e.g. 500, to perform a normal SCF cycle
  6. Add restart_read_only restart_file and force_single_restartfile .true. to load the restart file from the previous init_2 calculation
  7. To check the localisation of the core hole after calculation, we can again use the output cube spin_density keyword
  8. Run all ten hole calculations

With the total energies of these hole calculations and the total energy of the reference ground state single point calculation, we again calculate the core electron binding energies using equation (3).

deltascf_basis Command

As an alternative to the deltascf_projector procedure described above, there is an additional method implemented in FHI-aims to force localisation of the core-hole: the deltascf_basis command. This command also uses the maximum overlap method (MOM) methodology to track the state with the occupation constraint. However, deltascf_basis tries to force the occupation to a state similar to a specified atomic basis function in the system rather than a previously determined KS state.

The deltascf_basis keyword follows the following syntax:

deltascf_basis atom spin basis_type basis_n basis_l basis_m occ max_KS_state

  • atom indicates the number of the specified atom (in the geometry.in file), on which the basis function will be centred
  • spin is again arbitrarily set to 1
  • basis_type is chosen to be atomic, because we want to constrain the hole to an atomic orbital like hole
  • basis_n, basis_l, and basis_m now define the basis function via quantum numbers, to specify our 1s state we choose accordingly basis_n = 1, basis_l = 0 and basis_m = 0
  • occ is again set to 0.0 to introduce a full core hole
  • max_KS_state is again set to the highest number of 1s states present, which is 10 for the azulene molecule

The deltascf_basis command is often more convenient than the deltascf_projector keyword because it needs less calculations and attention to detail when setting up the input files. However, deltascf_basis is only available for the aperiodic case and can be problematic when used with hybrid functionals.

Hole Calculations using deltascf_basis

Performing the calculations is much simpler with the deltascf_basis keyword:

  1. Generate one hole folder for each carbon atom
  2. It is not required to specify a restart file
  3. The control.in file should contain:
  4. charge = 1.0 to remove one electron
  5. deltascf_basis N 1 atomic 1 0 0 0.0 10 where N is chosen from 1 to 10 for each hole calculation
  6. output cube spin_density to check for the correct localisation of the core hole
  7. Run all core-hole calculations

The total energies of these hole calculations and the total energy of the reference ground state single point calculation are again used to calculate the core electron binding energies via equation (3).

Comparison of deltascf_basis and deltascf_projector

The following table contains the results for the absolute and relative binding energies for the ten atoms of azulene calculated both with the deltascf_basis and deltascf_projector commands.

Table 3. Absolute and relative C 1s core-electron binding energies for azulene were calculated with the keywords deltascf_basis and deltascf_projector (PBE/tight-tier2). All values are given in eV.

C-atom EB (proj.) ΔEB (proj.) EB (basis) ΔEB (basis)
C1 289.03 0.00 289.03 0.00
C2 288.70 -0.33 288.70 -0.33
C3 288.70 -0.33 288.70 -0.33
C4 289.44 0.41 289.44 0.41
C5 289.44 0.41 289.44 0.41
C6 289.64 0.62 289.64 0.62
C7 289.64 0.62 289.64 0.62
C8 289.29 0.26 289.29 0.26
C9 289.29 0.26 289.29 0.26
C10 289.54 0.51 289.54 0.51

We calculated the binding energies of all ten carbon atoms in azulene. As is visible in the resulting energies, there are two unique atoms and four symmetry-equivalent pairs of atoms. Overall, only six hole calculations were necessary after all. However, in this case, it was good practice to calculate the binding energies for all atoms in order to use the symmetry-enforced equivalence of binding energies as an internal sanity check.

We can also clearly see that deltascf_basis and deltascf_projector provide equivalent results when performed correctly.

Comparison with Experimental Data

There exists no experimental XPS data for azulene in the gas phase. However, experimental XPS data for a thin film multilayer of azulene was published by Klein et al. 1 Assuming that the inter-molecular interactions in the multilayer only include weak dispersion interactions, which do not influence the core-hole binding energies much, we can still compare our "gas-phase" DFT data to these experiments.

Figure 2 shows this comparison; again, we used the pseudo-Voigt scheme above to broaden the simulated ΔSCF data and simulate the experimental resolution. We used an FWHM of 0.7 eV and a GL-ratio of 0.4 for the azulene data. Then, we summed up the ten sub-peaks to generate the whole spectrum.


Figure 2. Comparison of the calculated and experimental C 1s XP spectrum of azulene. The experimental data was measured for a thin film multilayer. The calculated data are for a gas-phase molecule. The calculated data was shifted by -4.8 eV to match the experimental energy scale.

In the figure above, the experimental data is taken from Klein et al.1

In Figure 2, we used the absolute energy scale of the experimental data. We needed to rigidly shift the calculated data by -4.8 eV to match this energy scale. Comparing the absolute energy scales of experiment and simulation would necessitate a cautious treatment of the work function for the experimental and simulated data. At the same time, systematic approximations in the DFT calculation would remain. The comparison in Figure 2, however, shows that the peak shape and, therefore, the relative binding energies of carbon atoms in different chemical environments are already in much better agreement than the size of the empirical shift parameter (-4.8 eV) would suggest.

The experimental resolution in these condensed phase measurements allows only a comparison of the overall shape. However, a closer inspection is possible in high-resolution XPS experiments in the gas phase. For example, in the gas phase spectra for naphthalene by Minkov et al., 2 two subpeaks with a difference in binding energy of 0.36 eV are resolved. This relative binding energy shift is in excellent agreement with the value of 0.32 eV obtained by ΔSCF (PBE-tight-tier2).

References


  1. B.P. Klein, N. J. van der Heijden, S. R. Kachel, M. Franke, C. K. Krug, K. K. Greulich, L. Ruppenthal, P. Müller, P. Rosenow, S. Parhizkar, F. C. Bocquet, M. Schmid, W. Hieringer, R. J. Maurer, R. Tonner, C. Kumpf, I. Swart, and J. M. Gottfried, Phys. Rev. X 2019, 9, 011030. https://doi.org/10.1103/PhysRevX.9.011030 

  2. I. Minkov, F. Gel’mukhanov, R. Friedlein, W. Osikowicz, C. Suess, G. Öhrwall, S. L. Sorensen, S. Braun, R. Murdey, W. R. Salaneck, and H. Ågren, J. Chem. Phys., 2004, 121, 5733. https://doi.org/10.1063/1.1784450