Table of Contents
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.