# CP2K Open Source Molecular Dynamics

### Sidebar

#### For Developers

exercises:2019_conexs_newcastle:ex2

# Linear-response TDDFT for XAS in CP2K: the XAS_TDP method

In this exercise, you will get familiar with the XAS_TDP method in CP2K. We will simulate the Aluminium K-edge XANES of trimethylaluminium with chemical formula Al$_2$(CH$_3$)$_6$:

The molecule is made of 26 atoms and we will use medium quality double-zeta basis sets. This way, most of the exercises can be done with little to no parallelization of the calculations. All the concepts discussed below are still valid for larger system though.

You should create a new directory for each of the 3 parts of the exercise

You can find documentation for the XAS_TDP method here and refer to the slides.

Remember that the LR-TDDFT code is not yet included in the main CP2K code, hence you have to modify the module load line in the cp2k.sh file in order to load the correct version:

module load CP2K/xas_tdp_libint2

To make sure you are using the proper version, type module unload CP2K/6.1-foss-2017b before launching any calculation.

## Part 1: the standard workflow

The first step to all atomistic simulations is to get hold of the molecule/material structure. The latter can come from different sources: a careful geometry optimization, frames of a molecular dynamics run or experiments. In order to keep the focus on XANES, a geometry optimization was run ahead of time and the structure stored in the following geometry.xyz file:

geometry.xyz
      26
Trimethylaluminium structure
Al         4.2675935634        5.5407111338        4.7089610588
Al         6.4608798572        4.2997143768        5.3967873091
C         5.1793535458        3.8631556857        3.7238051743
H         5.6698412406        4.2214449576        2.8092452108
H         4.1254815388        3.6491894229        3.4661776259
C         5.5489970130        5.9771338270        6.3820518382
H         5.1739704174        6.9992805163        6.1899010533
H         5.0579905439        5.6189992181        7.2963828304
H         5.5549255131        2.8411776275        3.9158135330
H         6.6028509125        6.1906122100        6.6401744970
C         8.1577533120        4.8462637633        4.5673348260
H         8.0626764249        5.7568375257        3.9554443238
H         8.5548879421        4.0566420550        3.9092516157
H         8.9349946502        5.0418279813        5.3241717501
C         6.3493279327        2.8490785738        6.7195500278
H         5.3140757011        2.6112574534        7.0097068947
H         6.8924547937        3.1044799443        7.6434085417
H         6.8018718123        1.9194670696        6.3386170476
C         4.3795530867        6.9915820838        3.3864761343
H         5.4149333409        7.2291739602        3.0966185981
H         3.9272669851        7.9212012690        3.7677052718
H         3.8365057581        6.7367174672        2.4624343308
C         2.5706506209        4.9941975061        5.5382826840
H         2.6657281632        4.0834902522        6.1499880352
H         1.7933806464        4.7988174956        4.7814409011
H         2.1735548094        5.7836967331        6.1965426058


Download/copy the above file in your working directory. This is a XYZ file, which structure is quite simple: the first line states the number of atoms in the structure, the second line is reserved for a description and the rest contains the Cartesian coordinates in Angstroms.

### Making a XAS_TDP compatible ground state calculation

To accurately describe the core states from which electrons are excited in XAS, all electron calculations are required. In CP2K, this means using the GAPW method. The following part1_gs.inp input file contains all required keywords for such a calculation. Pay attention to the comments (starting with exclamation marks) explaining what is going on.

part1_gs.inp
!This input file is meant for a ground state calculation of Al2(CH3)6 before XAS_TDP
&GLOBAL
PROJECT  part1             !project name, useful to keep track of produced files
RUN_TYPE ENERGY
&END GLOBAL

&FORCE_EVAL
METHOD QS
&DFT
!where to find all-electron basis sets and potentials
BASIS_SET_FILE_NAME  EMSL_BASIS_SETS
POTENTIAL_FILE_NAME  POTENTIAL

&QS
METHOD GAPW          !GAPW for all-electron
&END QS

&POISSON                !Necessary to define POISSON section in non-perdiodic
PERIODIC NONE        !boudary conditions as the default assumes PBCs
POISSON_SOLVER MT
&END POISSON

&SCF
MAX_SCF    50
EPS_SCF    1.0E-08  !high quality ground state required for XAS_TDP
SCF_GUESS  RESTART  !to avoid recomputing when we can
&END SCF

&MGRID
CUTOFF 600
&END

&XC
&XC_FUNCTIONAL                !the stendard PBE functional, the LIBXC way
&LIBXC
FUNCTIONAL GGA_C_PBE
&END LIBXC
&LIBXC
FUNCTIONAL GGA_X_PBE
&END LIBXC
&END XC_FUNCTIONAL
&END XC

&END DFT

&SUBSYS
&CELL
ABC 10.0 10.0 10.0           !Big enough cell, we don't want density to spread out of it
PERIODIC NONE                !Need to specify non-periodicity here too
&END CELL
&TOPOLOGY
!define the structure externally => reuse the same file all the way
COORD_FILE_NAME geometry.xyz
COORD_FILE_FORMAT xyz
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&KIND H
BASIS_SET Ahlrichs-pVDZ         !Kinds are described by all-electron potentials and
POTENTIAL ALL                   !double-zeta quality all-electron basis sets
&END KIND
&KIND C
BASIS_SET Ahlrichs-pVDZ
POTENTIAL ALL
&END
&KIND Al
BASIS_SET Ahlrichs-pVDZ
POTENTIAL ALL
&END
&END SUBSYS
&END FORCE_EVAL


Download this input and place it into your working directory. Make sure that the geometry.xyz file is also there. Because it is a relatively small system, we do not need to run CP2K in parallel. To launch the calculation, write in the terminal (after modifying the paths in the cp2k.sh file):

sbatch cp2k.sh

Where cp2k.sopt is the serial version of CP2K, the -i flag indicates the input file and the -o flag the output. You can follow the progress of your calculation by typing the following in the terminal (press ctrl+C to stop):

tail -f part1_gs.out

Once the calculation finishes, you will notice that a new file named part1-RESTART.wfn was created. This file contains information about the converged ground state wave function, such that further calculation sharing the same project name will be able to skip the SCF convergence step and be much faster.

Note that in principle, you should make sure your ground state calculation is converged. For this exercise, we assume that it has been done ahead of time.

### Checking the core donor states

The next step in the workflow consists in verifying that we can find the proper donor core states, which are 2 Aluminium 1s in our case. To do that, we need to trigger a XAS_TDP calculation on top of our previously computed ground state with the CHECK_ONLY keyword. The latter makes sure that the programs only analyzes the potential donor molecular orbitals (MOs) and stops before doing any expensive operation. The next input file does just that:

part1_check.inp
&GLOBAL
PROJECT  part1  !This is important to keep the same project name to be able to RESTART
RUN_TYPE ENERGY
&END GLOBAL

&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME  EMSL_BASIS_SETS
POTENTIAL_FILE_NAME  POTENTIAL

&QS
METHOD GAPW
&END QS

&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON

&SCF
MAX_SCF    50
EPS_SCF    1.0E-08
SCF_GUESS  RESTART
&END SCF

&MGRID
CUTOFF 600
&END

&XC
&XC_FUNCTIONAL
&LIBXC
FUNCTIONAL GGA_C_PBE
&END LIBXC
&LIBXC
FUNCTIONAL GGA_X_PBE
&END LIBXC
&END XC_FUNCTIONAL
&END XC

!The section controlling the XAS calculations
&XAS_TDP

CHECK_ONLY   !in order to verify our donor MOs

&DONOR_STATES
DEFINE_EXCITED BY_KIND
KIND_LIST Al
STATE_TYPES 1s
N_SEARCH 2       !We know that Al 1s are the 2 MOs with the lowest energy
#LOCALIZE        !rmove the # to uncomment LOCALIZE
&END DONOR_STATES

&END XAS_TDP

&END DFT

&SUBSYS
&CELL
ABC 10.0 10.0 10.0
PERIODIC NONE
&END CELL
&TOPOLOGY
COORD_FILE_NAME geometry.xyz
COORD_FILE_FORMAT xyz
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&KIND H
BASIS_SET Ahlrichs-pVDZ
POTENTIAL ALL
&END KIND
&KIND C
BASIS_SET Ahlrichs-pVDZ
POTENTIAL ALL
&END
&KIND Al
BASIS_SET Ahlrichs-pVDZ
POTENTIAL ALL
&END
&END SUBSYS
&END FORCE_EVAL


You can download/copy the above file in your working directory. Note that the project name is the same as in part1_gs.inp; this ensure that we can restart and skip the ground state calculation. Only a minimal &XAS_TDP section was added, with the CHECK_ONLY keyword and the &DONOR_STATES subsection. Note that the LOCALIZE keyword is commented for now. Adapt the paths in cp2k.sh and run this input.

Once the calculation is over, open part1_check.out and look for XAS_TDP. This will lead you to the relevant part of the output file, in which you can check the quality of the donor MOs. Remember: for the theory to apply correctly, the donor states need to be local in space and of 1s character. Look at the overlap and charge indicators, are we good to proceed ?

No, we are not, because the Aluminium atoms are equivalent under symmetry and the 2 lowest energy MOs will be arbitrary linear combinations of the 1s states. To get the needed donor state quality, uncomment the LOCALIZE keyword in the input file and rerun.

### Computing the XANES spectrum

Now that we have a proper ground state and we know that the donor MOs are of good quality, we can do an actual XANES calculation on top. To do so, we, have to complete the &XAS_TDP section to include information about the integration grid, the dipole form and most importantly the kernel:

part1_xas.inp
&GLOBAL
PROJECT  part1
RUN_TYPE ENERGY
&END GLOBAL

&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME  EMSL_BASIS_SETS
POTENTIAL_FILE_NAME  POTENTIAL
AUTO_BASIS RI_XAS MEDIUM            !automatically generates an RI basis

&QS
METHOD GAPW
&END QS

&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON

&SCF
MAX_SCF    50
EPS_SCF    1.0E-08
SCF_GUESS  RESTART
&END SCF

&MGRID
CUTOFF 600
&END

&XC
&XC_FUNCTIONAL
&LIBXC
FUNCTIONAL GGA_C_PBE
&END LIBXC
&LIBXC
FUNCTIONAL GGA_X_PBE
&END LIBXC
&END XC_FUNCTIONAL
&END XC

!The section controlling the XAS calculations
&XAS_TDP

CHECK_ONLY  F !set CHECK_ONLY to false for actual computation

&DONOR_STATES
DEFINE_EXCITED BY_KIND
KIND_LIST Al
STATE_TYPES 1s
N_SEARCH 2
LOCALIZE
&END DONOR_STATES

TAMM_DANCOFF         !For a full-TDDFT calculation, set TAMM_DANCOFF F
DIPOLE_FORM LENGTH   !LENGTH or VELOCITY, it does not really matter
GRID Al 100 100      !100x100 grid rather small but enough here

&KERNEL
&XC_FUNCTIONAL                    !In principle, one could use a
&LIBXC                         !different functional for the
FUNCTIONAL GGA_C_PBE        !kernel, but to be consistent
&END LIBXC                     !one usually uses the same as
&LIBXC                         !for the ground state. Only
FUNCTIONAL GGA_X_PBE        !LIBXC functionals are available
&END LIBXC
&END XC_FUNCTIONAL
&END KERNEL

&END XAS_TDP

&END DFT

&SUBSYS
&CELL
ABC 10.0 10.0 10.0
PERIODIC NONE
&END CELL
&TOPOLOGY
COORD_FILE_NAME geometry.xyz
COORD_FILE_FORMAT xyz
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&KIND H
BASIS_SET Ahlrichs-pVDZ
POTENTIAL ALL
&END KIND
&KIND C
BASIS_SET Ahlrichs-pVDZ
POTENTIAL ALL
&END
&KIND Al
BASIS_SET Ahlrichs-pVDZ
POTENTIAL ALL
&END
&END SUBSYS
&END FORCE_EVAL


You can run this XAS calculation the usual way.

Note that again, the program reuses the previous ground state calculation. The spectral data can be found in the part1.spectrum file. The excitation energies are given in eV and the oscillator strength are within the dipole approximation (arbitrary units). There are two blocks of data, one for each Al 1s. Note that the values are practically the same since the Al atoms are indistinguishable.

If you wish, you can plot the XANES spectrum using this python script. Make sure the script is located in the same directory as your part1.spectrum file and type: Note: before using this script, you have to load a proper python installation by typing: module load Anaconda3/

python plot_spectrum.py part1.spectrum Al 1s 1.0 1500 1525 1000

Where Al and 1s tell the script which data to read in part1.spectrum, 1.0 is the width for the gaussian broadening, 1500 and 1525 are the x-axis (energy) bounds and 1000 is for the resolution (1000 points). You can read the first few lines of the script for documentaion. Remember that TDDFT results usually need to be shifted to match experiments (either by an empirical amount or by a ΔSCF calculation)

If you have time, feel free to:

• Try different broadening when plotting. If you use a small enough value, you will see the individual contributions
• Try different basis sets (*e.g. 6-31G*, Ahlrichs-def2-TZVP, etc.)
• Change the precision parameters such as GRID, AUTO_BASIS or add a RI_RADIUS. You will notice that the calculations we did were already well converged.

## Part 2: Higher quality hybrid functional calculation

It is generally accepted that standard GGA level TDDFT is not enough for accurately describing X-ray absorption spectrosocpy. Instead, most calculations in the literature are done at the hybdrid level. In this case, a fraction of the exchange part of the XC functional is replaced by the so-called exact exchange or Hartree-Fock exchange:

$$E_{xc}[n] = E_c^{DFT}[n] + (1-s_x)E_x^{DFT}[n] + s_x E_x^{HF}[n]$$

where $s_x$ is the fraction of Hartree-Fock exchange. This exact exchange notably allows to catch non-local effects, which cannot be described by a semi-local description such as GGA.

In this part of the exercise, we will simulate the Aluminium K-edge of Al$_2$(CH$_3$)$_6$ at the hybrid level, using the BHandHLYP functional. Note that such calculations are significantly more expensive and it is recommended that you run CP2K in parallel.

hybrid.inp
&GLOBAL
PROJECT  hybrid
RUN_TYPE ENERGY
&END GLOBAL

&FORCE_EVAL
METHOD QS
&DFT
BASIS_SET_FILE_NAME  EMSL_BASIS_SETS
POTENTIAL_FILE_NAME  POTENTIAL
AUTO_BASIS RI_XAS LARGE                   !want high precision

&QS
METHOD GAPW
&END QS

&POISSON
PERIODIC NONE
POISSON_SOLVER MT
&END POISSON

&SCF
MAX_SCF    50
EPS_SCF    1.0E-08
SCF_GUESS  RESTART
&END SCF

&MGRID
CUTOFF 600
&END

&XC
&XC_FUNCTIONAL                !The BHandHLYP functional, popular with XAS TDDFT
&LIBXC
FUNCTIONAL HYB_GGA_XC_BHandHLYP
&END LIBXC
&END XC_FUNCTIONAL
&HF
FRACTION 0.5               !the functional requires 50% of exact exchange
&END HF
&END XC

!The section controlling the XAS calculations
&XAS_TDP

&DONOR_STATES
DEFINE_EXCITED BY_KIND
KIND_LIST Al
STATE_TYPES 1s
N_SEARCH 2
LOCALIZE
&END DONOR_STATES

TAMM_DANCOFF
DIPOLE_FORM LENGTH
GRID Al 300 300              !Denser grid for integration

&KERNEL
NPROCS_GRID 4             !Use 4 procs per grid, ignored if serial run
RI_RADIUS 3.0             !Use also neighbor's RI basis functions for density projection
&XC_FUNCTIONAL
&LIBXC
FUNCTIONAL HYB_GGA_XC_BHandHLYP
&END LIBXC
&END XC_FUNCTIONAL
&EXACT_EXCHANGE
SCALE 0.5
&END EXACT_EXCHANGE
&END KERNEL

&END XAS_TDP

&END DFT

&SUBSYS
&CELL
ABC 10.0 10.0 10.0
PERIODIC NONE
&END CELL
&TOPOLOGY
COORD_FILE_NAME geometry.xyz
COORD_FILE_FORMAT xyz
&CENTER_COORDINATES
&END CENTER_COORDINATES
&END TOPOLOGY
&KIND H
BASIS_SET Ahlrichs-def2-TZVP   !moved to higher quality triple zeta basis sets
POTENTIAL ALL
&END KIND
&KIND C
BASIS_SET Ahlrichs-def2-TZVP
POTENTIAL ALL
&END
&KIND Al
BASIS_SET Ahlrichs-def2-TZVP
POTENTIAL ALL
&END
&END SUBSYS
&END FORCE_EVAL


Create a new directory for this part of the exercise and download the above input file into it. Make also sure that you copy the gemoetry.xyz file from the previous part in there. Have a look at the changes made in the input and their corresponding comments.

In particular, notice how the new keyword RI_RADIUS was added in the &KERNEL subsection. This defines a sphere around the excited atoms in which all RI basis functions contribute to the projection of the density. Remeber:

$$n(\mathbf{r}) \approx \sum_{pq}\sum_{\mu\nu} P_{pq} \ \big(\varphi_p \varphi_q \chi_\mu \big) \ S_{\mu\nu}^{-1} \ \chi_\nu(\mathbf{r})$$

where $P_{pq}$ is the density matrix, $\varphi_p,\varphi_q$ are atomic orbitals and $\chi_\mu, \chi_\nu$ are RI basis functions. Extending RI_RADIUS allows for a greater number of RI basis and a higher quality projection of the density, to which hybrid kernels have shown to be very sensitive.

With the keyword NPROCS_GRID set to 4, you can run the calculation on 8 cores and have simultaneous integration of the XC kernel for both Al atoms using all available resources. cp2k.popt -i hybrid.inp -o hybrid.out & </code>

If you copy the plot_spectrum.py script from part one in your working directory, you can use it the same way to plot the XANES spetrum:

python plot_spectrum.py hybrid.spectrum Al 1s 1.0 1545 1570 1000

If you have time, feel free to:

• Use smaller values for the GRID, AUTO_BASIS and especially RI_RADIUS keywords and survey the effects on the excitation energies
• Use other hybrid functionals (B3LYP, PBE0, etc.) and/or other basis sets (6-31G*, Ahlrichs-def2-QZVP, etc.) and how this affects the spectrum

## Part 3: L-edge spectroscopy (extra)

If you have enough time or a particular interest in L-edge spectroscopy, you can try and setup a calculation for it.

You should first create a new directory for this part. Then, starting from either part1_xas.inp or hybrid.inp, you can tweak the input file to look for 2s and/or 2p donor states by modifying the &DONOR_STATES subsection:

&DONOR_STATES
DEFINE_EXCITED BY_KIND
KIND_LIST Al
STATE_TYPES 2s          !Note that the STATE_TYPES keyword can be used multiple times such
STATE_TYPES 2p          !that both 2s and 2p excitations are done in a single calculation
N_SEARCH -1             !We do not know exactly where are the 2s/2p states in terms of
LOCALIZE                !energy => N_SEARCH -1 will look among all occupied MOs
&END DONOR_STATES

Once you have copied the above in your input file, you can launch your calculation as usual. Make sure that the gemoetry.xyz file is present in the same directory too. The plot_spectrum.py script can be used again.

Note that if you chose to do the hybrid calculation, you might have to increase the integration grid density and the RI_RADIUS even further to obtain converged results.