User Tools

Site Tools


Gaussian and Augmented Plane Wave Method

In this tutorial we perform all-electron calculations with GAPW to compute the near-edge X-ray absorption spectra of ice-1h and of liquid water. When performing all-electron calculations GAPW accounts for electron density redistribution of the core electrons, and therefore can be used for the computation of X-ray absorption spectra resulting from the absorption of X-rays exciting the core-shell electrons.

Short Intro

In near edge X-ray absorption fine structure (NEXAFS), the X-rays absorbed by each individual atom excite the core electrons, which therefore leave a hole in the core shell (i.e. a core hole). In the decay process the photo-electron can scatter from electrons of neighboring atoms. This scattering process affects the absorption coefficient that is measured for different energies of the incoming synchrothron radiation. In the case of water and ice we are interested in the spectra with a range of transition energies roughly between 530 and 550 eV.

The absorption coefficient is proportional to the transition probability, which can be computed by applying Fermi's Golden rule, and where the initial unperturbed state is that of the core orbital and the final state is that of the excited photo-electron that is affected by the presence of neighboring atoms.

The transition probability within DFT is computed by modifying the Kohn-Sham equations to include the core-hole potential on the absorbing atom. In practice the transition probability can be computed under the electric dipole approximation and using the full-core hole or half-core hole scheme (see and from the following expression: $$ I_{if} \propto |\big{<} \psi_{f} | \mathbf{\nabla} | \psi_{i} \big{>} |^2, $$

which is referred to as the “velocity” form (i.e. $-i \hbar \mathbf{\nabla}$ is the momentum operator). The excitation energies are obtained as $$ \hbar \omega_{if} = \varepsilon_f -\varepsilon_i. $$

We will compute the transition probabilities and energies in ice1h and liquid water using the GAPW method.

Briefly, in GAPW the total electron density $n(\mathbf{r})$ is written as $$ n(\mathbf{r}) = \tilde{n}(\mathbf{r}) - \sum_A \tilde{n}_A(\mathbf{r}) + \sum_A n_A(\mathbf{r}), $$ where $\tilde{n}(\mathbf{r})$ is the soft and smooth part of the density that is expanded in plane waves, while $\tilde{n}_A(\mathbf{r})$ and $n_A(\mathbf{r})$ are, respectively, the soft and hard densities, which are expanded in gaussian basis functions centred on each atom $A$. In this framework the KS equations can be solved without the use of pseudopotentials, while still accounting for the strongly varying electron density of the core electrons. For more details, see 10.1007/s002140050523, 10.1039/B615522G.

1. Task: compute the XAS spectrum of liquid water

In principle to compute the XAS spectrum of bulk water an average over many uncorrelated structures is required. However, due to limited time and resources in this tutorial we will just compute the XAS spectra of single snapshots for a small (1ps) trajectory of bulk water (32 H2Os). To start, extract the directories tutorial_gapw.tar.gz where you will be running the calculations. Also, load the CP2K module as explained in login.

Go to the directory gapw_water_xas, where there is a trajectory of liquid water at room temperature that was obtained using the GPW method (see the previous tutorial for more info Ab initio molecular dynamics). Now run the following script to extract a random snapshot from this trajectory


Now open the CP2K input file water_xas_gapw.inp with a text editor and add the tags required to perform the XAS simulation using the GAPW method. We have inserted the comment !Task: in parts where edits are required in order to successfully run CP2K.

We include the input file also here for reference:

    ! the EMSL basis set file contains all electron basis sets
    ! spin polarization required for the XAS calculation
      NGRIDS 5
      CUTOFF 400
      REL_CUTOFF 60
      ! Task: insert METHOD keyword to use gaussian and augmented plane wave method
      EPS_DEFAULT 1.0E-10
      ! algorithm to construct the atomic radial grid for GAPW
      ! parameters needed for the GAPW method, look at the manual for more details
      EPSFIT       1.E-4 ! precision to give the extension of a hard gaussian
      EPSISO       1.0E-12
      EPSRHO0      1.E-8
      LMAXN0       4
      LMAXN1       6
      ALPHA0_H     10 ! Exponent for hard compensation charge
    &END QS
      MAX_SCF 20
      EPS_SCF 1.0E-5
         ALPHA 0.5

        EPS_SCF 1.0E-5
        MAX_SCF 50
    &END SCF

        &END PBE
         XC_SMOOTH_RHO NN50
      &END XC_GRID
               PARAMETER_FILE_NAME dftd3.dat
               TYPE DFTD3
               R_CUTOFF [angstrom] 16
    &END XC

      ! Task: specify below the METHOD to use to compute the XAS spectra
      ! half-core hole and the full core-hole are possible methods, choose transition potential half hole
      ! METHOD TP_HH

      ! Task: include the STATE_TYPE keyword to specify the states to compute the spectra
      ! in NEXAFS experiments one looks at the excitation of the innermost-core shell
      ! STATE_TYPE 1s
      ! Task: include the ATOMS_LIST keyword for the calculation of XAS
      ! you can look at the list of atoms to include in the .xyz file for the snapshot
      ! In order to include atoms from X to Y use the syntax X..Y
      ! ATOMS_LIST 1..32
      ! This keyword indicates the number of virtual KS orbitals
      ! to compute the XAS
      ADDED_MOS     60

      ! Below we add another SCF section to perform the calculation of the XAS spectrum
      ! after having optimized the wave-function with GAPW     
         EPS_SCF 1.0E-5
         MAX_SCF 100
            ALPHA 0.6
         &END MIXING

           EPS_SCF 1.0E-5
           MAX_SCF 50
         &END OUTER_SCF
      &END SCF
             FILENAME ./water
               XAS_SCF 15
             ADD_LAST NUMERIC
           FILENAME ./water
           FILENAME ./water
             ABC 9.8528 9.8528 9.8528 
    ! Task: specify the BASIS_SET and POTENTIAL for the gapw calculation
    ! for both O and H we want to use the all-electron 6-31G* basis set
    &KIND H
      ! BASIS_SET 6-31G*
      ! number of points for the angular part of the grid, needed for GAPW
      ! number of points for the radial part of the grid, needed for GAPW
      RADIAL_GRID 200
    &KIND O
      ! BASIS_SET 6-31G*
      RADIAL_GRID 200


After having edited the input perform the XAS simulation by running

mpirun -n 4 cp2k.popt  -i water_xas_gapw.inp -o out.out &

Check the output file as it's being written. You should notice that after the first SCF cycle, where the wavefunction is optimized, there are a number of SCF loops where CP2K is performing the XAS simulation, see for instance in the output XAS_TP_SCF WAVEFUNCTION OPTIMIZATION. It should take about half an hour to run this job so while waiting, let us calculate the XAS spectra for ice. At the end of the tutorial we will try to compare the water and ice spectra.

2. Task: compute the XAS spectrum of ice-1h

Go to the directory gapw_ice1h_xas Run the XAS simulation using only 2 cores. There is no need to edit the input file this time.

The spectrum calculation for ice should be quicker because there is no need to average it over many different molecules. After this calculation is finished you should see several files called ice1h-xas_at*.spectrum, which contain the bare XAS spectra. In the first column there is the index of the virtual KS state, the second column is the transition energy, the third, fourth and fifth column are the transition probabilities projected onto X, Y and Z, and finally the 6th column is the norm of the transition probability.

Below there is an example:

 Absorption spectrum for atom      1, index of excited core MO is     1, # of lines     7
   161    539.43932962     -0.07310887      0.03134486     -0.05316957      0.00915441   0.00000
   162    540.72836460     -0.06587398     -0.04533860     -0.03016382      0.00730483   0.00000
   163    540.80816872      0.02102299     -0.13579881      0.02240266      0.01938516   0.00000
   164    541.59544040     -0.11225534      0.13024845     -0.07569374      0.03529546   0.00000
   165    541.92869553     -0.23353331     -0.06321242     -0.03824229      0.05999609   0.00000
   166    542.18907300      0.03676265      0.04213791      0.19284643      0.04031684   0.00000
   167    542.38178876      0.02306601      0.13887860     -0.07127073      0.02489882   0.00000

Now run the following command to convolute the spectrum so we can compare it to previously published calculations and NEXAFS experiments for ice 1h:


Open gnuplot and plot, respectively, the bare and the convoluted spectrum by typing

plot "./spectrum.inp" using 2:6 with l,  "./spectrum.out" using 1:5 with l


How do your results for the convoluted spectrum compare with previous experiments and simulations? Look for instance at the bottom of Fig.2 of the papers 10.1063/1.2928842 or Fig 2 of 10.1063/1.1879752. There should be several things that do not match with our calculations. Apart from a shift towards larger binding energies compared with the two papers, it looks like not all the states of the conduction band are present.

Now the task is to correct this by editing a tag in the &XAS ..&END XAS section. So remove the *spectrum files, re-run the simulation (you can also input RESTART T in the XAS section) and repeat the tasks just described. How does the spectrum compare now? There should still be clear differences, but the main qualitative features should be observed. The presence of a broad band corresponding to the conduction band of ice should be visible, as well as an apparent shoulder (the main edge) just to the left of the higher and broader (post-edge) peak.

3. Task: compare the XAS spectra of ice-1h and water

Now the calculation of the water spectra should be finished. So from the directory gapw_water_xas run the script to obtain the convoluted spectrum for water averaged over the 32 molecules and plot it in gnuplot together with the ice spectrum

plot "./spectrum.out" using 1:5 with l,  "../gapw_ice1h_xas/spectrum.out" using 1:5 with l
  • Compare the spectra between each other and with the paper 10.1063/1.2928842 where the spectra for both ice 1h and single water snapshots have been calculated.
  • Fig.5 of this review 10.1021/acs.chemrev.5b00672 also shows a comparison between recent ice and liquid water spectra. Can you clearly identify the pre-, main- and post-edge features in the bulk water spectra and in ice?
  • What is the main reason for the different shapes between water and ice?
  • If our results do not match well with previously published experimental or the simulation data, could you think of possible reasons for the discrepancy?
In order to fix the shift of the XAS spectra towards larger binding energies it is possible to use the $\Delta SCF$ approach to compute the first excitation energy and rigidly shift the XAS spectra accordingly. This can be done in CP2K.
Here we have shown simple examples that can be performed on small machines and using a limited time. Therefore, careful tests on basis set size, XC functional, system size, sampling of configuration space, etc. have to be carried out for production runs to get more reliable spectra.
exercises/2017_uzh_cp2k-tutorial/gapw.txt · Last modified: 2017/07/12 22:27 by gtocci