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.
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 https://doi.org/10.1103/PhysRevB.5.844 and 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 10.1007/s002140050523, 10.1039/B615522G.
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 WATER_AIMD-pos-1.xyz
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
../LIB_TOOLS/get_rand_snapshot.sh
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:
&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
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.
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:
../LIB_TOOLS/get_average_spectrum.sh
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
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 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.
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