# Open SourceMolecular Dynamics

### Sidebar

#### For Developers

exercises:2017_uzh_cp2k-tutorial:gapw

# 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 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.

## 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 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:

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
! 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

! 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
&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
&END KIND
&KIND O
! BASIS_SET 6-31G*
! POTENTIAL ALL
LEBEDEV_GRID 80
&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.

## 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:

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

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.

## 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.