diff --git a/docs/methods/properties/x-ray/delta-kick.md b/docs/methods/properties/x-ray/delta-kick.md new file mode 100644 index 0000000000..3246c92ccf --- /dev/null +++ b/docs/methods/properties/x-ray/delta-kick.md @@ -0,0 +1,295 @@ +# X-Ray Absorption from δ-Kick + +This tutorial shows how to run a Real-Time Time-Dependent DFT calculation using the so-called δ-kick +approach to compute absorption electronic spectra. + +We warmly recommend having a look first at the Linear-Response approach, for instance at +[the X-Ray frequencies](./tddft). See also this article and the supplementary information: LINK + +This tutorial is structured as follows: + +- Brief theory overview +- CP2K input file overview + +Then, the link between the CP2K output (time-dependent dipole moment) and the quantity targetted +(the absorption spectra in the frequency domain) is performed in this jupyter notebook: +[download this file](https://www.cp2k.org/_media/howto:delta_kick_analyse.zip). Of course, one can +use its own script and method. Yet, as these analyses present some subtilities, we propose our which +shows quantitative results with the Linear Response one. + +## Theory Overview + +### Response theory reminders + +In the perturbative regime and at the linear order, the time-dependent response of an electronic +cloud can be decomposed in the Fourier space: the response at a given frequency is proportional to +the perturbation acting on the system at the same frequency. In particular, the induced dipole +moment (or equivalently the induced current) at a given frequency, $\mu^\omega$, is given by the +product of the polarizability matrix, $\alpha(\omega, \omega)$, with the perturbative electric field +at the same frequency, $F^\omega$, within the dipolar approximation. + +$$ +\mu^\omega = \alpha(\omega, \omega) \cdot F^\omega +$$ + +This quantity involved a dot product between the field vector and the polarizability matrix. For the +next paragraphs, we will drop the vectorial behavior and discuss it later. + +In Linear-Response-TDDFT approaches, one computes the excited state promoted by an electric field +oscillating at a specific frequency $\omega$ and then computes the transition dipole moment +associated with such electronic transition. This transition dipole moment is then used to compute +the polarizability for resonant or non-resonant frequency. + +In the Real-Time approach presented here, we will excite **all** the possible electronic transitions +at the beginning of the simulation and then propagate this excited electronic cloud to observe the +induced dipole moment. We will deduce the polarizability tensor for any frequency from this induced +dipole moment. + +### The δ-kick technic + +In the δ-kick technic, the electronic structure is first obtained at the ground state, for instance +using DFT, and then perturbed with an instantaneous electric field named a δ-kick: + +$$ +F(t) = F^0 \delta(t) +$$ + +This δ-kick perturbes the ground state wave-function at $t=0^-$ using either the dipole moment +operator (length gauge) or using the momentum one (velocity gauge). Then, this excited wave-function +is propagated in real time by numerically integrating the DFT equivalent of the Schrodinger +equation. For more information, have a look for instance at the recent work in CP2K by +[](#Mattiat2022). + +The electric δ-kick can be written in the frequency domain: + +$$ +F^\omega = \frac{F^0}{2 \pi} \int_{-\infty}^{+ \infty} \delta(t) e^{i \omega t} dt = \frac{F^0}{2 \pi} e^{i \omega \times 0} = \frac{F^0}{2 \pi} +$$ + +The field amplitude in the Fourier space is $F^0 / 2 \pi$ for all frequencies: this instantaneous +perturbation does indeed contain all the frequencies. The time-dependent wave-function can be thus +described to be the one at the ground state plus all the possible excited states. During the time +propagation, the dipole moment of the molecule will be given by the superposition of all possible +oscillations related to all the excited states. This complex behavior in the time domain became +simple in the frequency one since we know the amplitude of the perturbation applied at each +frequency: + +$$ +\mu^\omega = \frac{1}{2 \pi} \alpha(\omega, \omega) F^0 +$$ + +Therefore, in order to get the polarizability $\alpha(\omega, \omega)$ one prepares the perturbed +state at $t=0^-$ for a given field amplitude, then propagates the wave-function in real time and +compute the time-dependent dipole moment. The dipole moment is Fourier transformed and the +polarizability is extracted using: + +$$ +\text{Re} \left[ \alpha(\omega, \omega) \right] = 2 \pi \frac{\text{Re} \left[ \mu^\omega \right] }{F^0} \\ +\text{Im} \left[ \alpha(\omega, \omega) \right] = 2 \pi \frac{\text{Im} \left[ \mu^\omega \right] }{F^0} +$$ + +The amplitude of the dipole moment in the Fourier space may be very small if the frequencies are far +from resonances. But, in principle, one can extract the whole spectra from one Real-Time +calculation. Of course, numerically speaking this is not the case since the time scale involved for +UV frequencies is quite different from the one involved for X-Rays: a Real-Time propagation +calculation is thus run for a specific range of frequency. + +### Absorption spectrum + +Using the time-dependent dipole moment, we got the frequency-dependent polarizability $\alpha$ +assuming a perturbative and linear response regime. The real and imaginary part of the +polarizability defines the nature of the response. If the frequency is at an electronic resonance, +then the polarizability has an imaginary part. The polarizability is pure real off resonances. + +From one Real-Time propagation, one can obtain 3 components of the polarizability tensor. For +instance, applying a δ-kick along the $x$ direction will provides the components $\alpha_{xx}$, +$\alpha_{yx}$ and $\alpha_{zx}$ by Fourier Transform the dipole moment along the $x$, $y$ and $z$ +direction respectively. Therefore, to get the full polarizability tensor, one should run 3 Real-Time +propagations using a perturbation along 3 orthogonal directions, for instance, $x$, $y$, and $z$. + +To compare with experiments, one often assumes that the system is averaged over all possible +orientations in space. Then, the absorption spectrum $I(\omega)$ is proportional to: + +$$ +I(\omega) \propto \sum_{i=x,y,z} \text{Im} \left[ \alpha_{ii}(\omega, \omega) \right] +$$ + +## CP2K Input + +In this tutorial, we will study the response of a carbon-monoxide in the gas phase upon a δ-kick +along its perpendicular direction. We will focus on the X-Ray range, but in principle, other ranges +of frequency can be sampled using this technic. + +This has been done for isolated carbon-monoxide at the DFT/PBEh level using a PCSEG-2 basis set. One +can look at the input file contained in RTP.inp, in the following we will emphasize a few important +points related to the δ-kick technic. + +```none +&GLOBAL + PROJECT RTP + RUN_TYPE RT_PROPAGATION +&END GLOBAL + +&MOTION + &MD + ENSEMBLE NVE + STEPS 50000 + TIMESTEP [fs] 0.00078 + TEMPERATURE [K] 0.0 + &END MD +&END MOTION + +&FORCE_EVAL + METHOD QS + &DFT + &REAL_TIME_PROPAGATION + APPLY_DELTA_PULSE .TRUE. + DELTA_PULSE_DIRECTION 1 0 0 + DELTA_PULSE_SCALE 0.001 + MAX_ITER 100 + MAT_EXP ARNOLDI + EPS_ITER 1.0E-11 + INITIAL_WFN SCF_WFN + PERIODIC .FALSE. + &END REAL_TIME_PROPAGATION + BASIS_SET_FILE_NAME BASIS_PCSEG2 + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + CUTOFF 1000 + NGRIDS 5 + REL_CUTOFF 60 + &END MGRID + &QS + METHOD GAPW + EPS_FIT 1.0E-6 + &END QS + &SCF + MAX_SCF 500 + SCF_GUESS RESTART + EPS_SCF 1.0E-8 + &END SCF + &POISSON + POISSON_SOLVER WAVELET + PERIODIC NONE + &END POISSON + &XC + &XC_FUNCTIONAL PBE + &PBE + SCALE_X 0.55 + &END + &END XC_FUNCTIONAL + &HF + FRACTION 0.45 + &INTERACTION_POTENTIAL + POTENTIAL_TYPE TRUNCATED + CUTOFF_RADIUS 7.0 + &END INTERACTION_POTENTIAL + &END HF + &END XC + &PRINT + &MULLIKEN OFF + &END MULLIKEN + &HIRSHFELD OFF + &END HIRSHFELD + &MOMENTS + PERIODIC .FALSE. + FILENAME =dipole + COMMON_ITERATION_LEVELS 100000 + &EACH + MD 1 + &END EACH + &END MOMENTS + &END PRINT + &END DFT + + &SUBSYS + &CELL + ABC 10 10 10 + ALPHA_BETA_GAMMA 90 90 90 + PERIODIC NONE + &END CELL + &TOPOLOGY + COORD_FILE_NAME carbon-monoxide_opt.xyz + COORD_FILE_FORMAT XYZ + &END TOPOLOGY + &KIND C + BASIS_SET pcseg-2 + POTENTIAL ALL + &END KIND + &KIND O + BASIS_SET pcseg-2 + POTENTIAL ALL + &END KIND + &END SUBSYS +&END FORCE_EVAL +``` + +### Time step + +The time step used to propagate the wave-function is related to the frequency of interest. +Typically, one should do a trade-off between the maximal frequency that can be sampled, this gives +the time step value, and the accuracy of the sampling itself, related to the total simulation time. +Typically, the time step should be about one order of magnitude lower than the frequency of +interest. Here, we will tackle the Oxygen K-edge: the electronic transition between the Oxygen 1s +and the first available excited state. + +According to Linear Response TDDFT calculation (using the XAS_TDP module, see +[this tutorial](./tddft)), the first transition occurs at 529 eV with an oscillator strength of +0.044 a.u. In our δ-kick approach, it means that the induced dipole moment should have an +oscillatory component around the frequency corresponding to 529 eV. Hence, we should have a time +step smaller than this frequency. For this tutorial, we have set the +[TIMESTEP](#CP2K_INPUT.MOTION.MD.TIMESTEP) to `[fs] 0.00078`, which leads to about 10 time steps per +period for an oscillation at 529 eV. + +The total time of the simulation is then given by the number of time steps required, here 500000, +and is related to the precision of the calculation: the more total time, the better the accuracy. +This value is less clear to determine: it depends on the strength of the perturbation, the error +thresholds used for the Real Time simulation, the amplitude of the polarizability value itself, and +the overall noise of the calculation. A good way to check that the total simulation time is enough +is to compute the polarizability for the targetted frequency using either 80% of the simulation time +or 100%. + +### Field properties + +The field perturbation (the δ-kick) is defined by its amplitude and polarization. + +The amplitude depends on the system of interest, typically $10^{-3}$ would do, and a good way to +check it is to run 2 simulations: one with once the amplitude and another with twice the amplitude. +Within the linear regime, the response of the system should have doubled as well. If this is not the +case, lower the field amplitude. For instance, by one or two orders of magnitude. Note that applying +a too-low field may lead to numerical noise. For isolated carbon monoxide, we have checked that +$10^{-3}$ is within the linear regime. + +**Please note that the actual perturbation applied in CP2K is not +[DELTA_PULSE_SCALE](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.DELTA_PULSE_SCALE) and depends +on the cell size. Look at the output file just before the Real Time propagation part to know the +exact perturbation amplitude used** + +The polarization of the field determines which excited state is triggered: the electric field and +the transition dipole moment should be non-perpendicular. One way to know how to select the +polarization is to run a Linear-Response TDDFT calculation and to look at the transition dipole +moment decomposition along the laboratory axis (or the oscillator strength decomposition). For +instance, have a look at the first part of [this tutorial](./resonant_excitation) on Real Time TDDFT +with a time-dependent field for isolated carbon monoxide. Relying on this calculation, the +electronic transition at 529 eV is perpendicular to the CO bond. For the geometry we are using, it +means that the field polarization should be along x (or y) by setting +[DELTA_PULSE_DIRECTION](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.DELTA_PULSE_DIRECTION) to +`1 0 0` and [DELTA_PULSE_SCALE](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.DELTA_PULSE_SCALE) +to `0.001`. + +### Note about the choice of the Gauge + +For an isolated system, use the length implementation because the perturbation is applied up to all +perturbative orders by setting [PERIODIC](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PERIODIC) +to `.FALSE.`. + +For condensed phase systems, use the velocity gauge instead by setting +[PERIODIC](#CP2K_INPUT.FORCE_EVAL.DFT.REAL_TIME_PROPAGATION.PERIODIC) to `.TRUE.`. In this case, the +perturbation will be applied only within the first order. We have used the length one for this +tutorial. + +## Analyzing the results + +To analyze the Real Time propagation simulation, we propose to go through a jupyter notebook: +[download this file](https://www.cp2k.org/_media/howto:delta_kick_analyse.zip). Note that this zip +file also contains the data needed to perform the analyses (the output of the CP2K run). diff --git a/docs/methods/properties/x-ray/index.md b/docs/methods/properties/x-ray/index.md index 799b53c261..620f5d1c50 100644 --- a/docs/methods/properties/x-ray/index.md +++ b/docs/methods/properties/x-ray/index.md @@ -6,5 +6,6 @@ titlesonly: maxdepth: 1 --- tddft +delta-kick correction_scheme ```