CP2K Open Source Molecular Dynamics

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== 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.

<note>
For more information see [[doi>10.1007/s002140050523]], [[doi>10.1039/B615522G]],

[[https://www.cp2k.org/_media/events:2015_cecam_tutorial:iannuzzi_gpw_gapw_b.pdf|iannuzzi_gpw_gapw_b.pdf]] and
[[https://www.cp2k.org/_media/events:2015_cecam_tutorial:iannuzzi_spectroscopy_b.pdf|iannuzzi_spectroscopy_b.pdf]]
</note>



===== 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 [[doi>https://doi.org/10.1103/PhysRevB.5.844]] and [[doi>https://doi.org/10.1103/PhysRevB.64.115107]]) 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 [[doi>10.1007/s002140050523]], [[doi>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 {{ :exercises:2017_uzh_cp2k-tutorial:tutorial_gapw.tar.gz|}}
where you will be running the calculations. Also, load the CP2K module 
as explained in [[exercises:2017_uzh_cp2k-tutorial:login|]].

Go to the directory ''gapw_water_xas'', where there is
a trajectory of liquid water ''WATER_AIMD-pos-1.xyz'' 
at room temperature that was obtained using
the GPW method (see the previous tutorial for more info [[exercises:2016_summer_school:aimd|Ab initio molecular dynamics]]).
Now run the following script to extract a random snapshot from this trajectory
<code>
../LIB_TOOLS/get_rand_snapshot.sh
</code>

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:

<code - water_xas_gapw.inp>
&FORCE_EVAL
  METHOD QS
  &DFT
    ! the EMSL basis set file contains all electron basis sets
    BASIS_SET_FILE_NAME  EMSL_BASIS_SETS
    POTENTIAL_FILE_NAME  POTENTIAL
    ! spin polarization required for the XAS calculation
    LSD
    &MGRID
      NGRIDS 5
      CUTOFF 400
      REL_CUTOFF 60
    &END MGRID
    &QS
      ! Task: insert METHOD keyword to use gaussian and augmented plane wave method
      !METHOD GAPW
      EXTRAPOLATION ASPC
      EXTRAPOLATION_ORDER 3
      MAP_CONSISTENT
      EPS_DEFAULT 1.0E-10
      ! algorithm to construct the atomic radial grid for GAPW
      QUADRATURE   GC_LOG
      ! 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
    &SCF
      MAX_SCF 20
      EPS_SCF 1.0E-5
      SCF_GUESS RESTART
     
      &MIXING
         METHOD BROYDEN_MIXING
         ALPHA 0.5
      &END MIXING

      &OUTER_SCF
        EPS_SCF 1.0E-5
        MAX_SCF 50
      &END OUTER_SCF
    &END SCF

    &XC
      &XC_FUNCTIONAL
        &PBE
        &END PBE
      &END XC_FUNCTIONAL
      &XC_GRID
         XC_SMOOTH_RHO NN50
         XC_DERIV NN50_SMOOTH
      &END XC_GRID
         &VDW_POTENTIAL
            POTENTIAL_TYPE PAIR_POTENTIAL
            &PAIR_POTENTIAL
               PARAMETER_FILE_NAME dftd3.dat
               TYPE DFTD3
               REFERENCE_FUNCTIONAL PBE
               R_CUTOFF [angstrom] 16
            &END
         &END VDW_POTENTIAL
    &END XC

    &XAS
      RESTART F
      ! 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
      
      DIPOLE_FORM   VELOCITY

      ! 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     
      &SCF 
         EPS_SCF 1.0E-5
         MAX_SCF 100
         &MIXING
            METHOD BROYDEN_MIXING
            ALPHA 0.6
         &END MIXING

         &OUTER_SCF
           EPS_SCF 1.0E-5
           MAX_SCF 50
         &END OUTER_SCF
      &END SCF
      &LOCALIZE
      &END
      &PRINT
         &PROGRAM_RUN_INFO
         &END
         &RESTART
             FILENAME ./water
             &EACH
               XAS_SCF 15
             &END
             ADD_LAST NUMERIC
         &END
         &XAS_SPECTRUM
           FILENAME ./water
         &END
         &XES_SPECTRUM
           FILENAME ./water
         &END
      &END
    &END
  &END DFT
  &SUBSYS
    &TOPOLOGY
      COORD_FILE_NAME  water.xyz
      COORDINATE XYZ 
    &END TOPOLOGY
    &CELL
             ABC 9.8528 9.8528 9.8528 
    &END CELL
    ! 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*
      ! POTENTIAL ALL
      ! number of points for the angular part of the grid, needed for GAPW
      LEBEDEV_GRID 80
      ! number of points for the radial part of the grid, needed for GAPW
      RADIAL_GRID 200
    &END KIND
    &KIND O
      ! BASIS_SET 6-31G*
      ! POTENTIAL ALL
      LEBEDEV_GRID 80
      RADIAL_GRID 200
    &END KIND
  &END SUBSYS
&END FORCE_EVAL

&GLOBAL
  PROJECT_NAME WATER_XAS
  RUN_TYPE ENERGY
  PRINT_LEVEL LOW
  FLUSH_SHOULD_FLUSH T
&END GLOBAL
</code>

After having edited the input perform the XAS simulation by running
<code>
mpirun -n 4 cp2k.popt  -i water_xas_gapw.inp -o out.out &
</code>
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:

<code - ice1h-xas_at1.spectrum>
 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
</code>

Now run the following command to convolute the spectrum so we can compare it to
previously published calculations and NEXAFS experiments for ice 1h:
<code>
../LIB_TOOLS/get_average_spectrum.sh
</code>

Open ''gnuplot'' and plot, respectively, the bare and the convoluted spectrum by typing
<code>
plot "./spectrum.inp" using 2:6 with l,  "./spectrum.out" using 1:5 with l
</code>

**Question:**

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 [[doi> 10.1063/1.2928842]]
or Fig 2 of [[doi> 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
<code>
plot "./spectrum.out" using 1:5 with l,  "../gapw_ice1h_xas/spectrum.out" using 1:5 with l
</code>

  * Compare the spectra between each other and with the paper [[doi>10.1063/1.2928842]] where the spectra for both ice 1h and single water snapshots have been calculated.
  * Fig.5 of this review [[doi>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?


<note> 
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.
</note>

<note>
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.
</note>