Dielectric function from Bloch-states dynamics
Step 0: Theoretical background
The dynamics of Bloch-states is found by integrating a set of Time Dependent Effective Schrödinger Equations (TD-ESEs). Effective means that an effective Hamiltonian is considered. We consider the many-body effective Hamiltonians seen in the Linear response from real time simulations. In addition, one can also consider effective Hamiltonians based on density-functional theory. The latter choice gives time-dependent density-functional theory.
At difference with what seen in Linear response from real time simulations, the coupling between electrons and the external field is described by means of the Modern Theory of Polarization:
Where [math]\displaystyle{ | v_{i \mathbf{k}} \rangle }[/math] are the time-dependent Bloch states (corresponding to the filled Bloch-states at zero-field), [math]\displaystyle{ \mathcal E(t) }[/math] is the external field and the term [math]\displaystyle{ i e \partial k }[/math] is the dipole operator consistent with periodic boundary condition. For more details on this last term, see Ref. [1] and [2].
[math]\displaystyle{ h_{rt} }[/math] is the Hamiltonian,
.
This contains the equilibrium Hamiltonian, that is the zero-field eigenvalues evaluated through DFT; the quasiparticle corrections, evaluated from GW or estimated from experiment, and the variation of the self-energy. The latter is a functional of the density matrix [math]\displaystyle{ \rho }[/math], which is obtained from the Bloch states.
From the solution of the equation of motion for the Bloch states, the polarization is calculated by means of Berry's phase
.
One can decide the level of theory by selecting the terms to include in the Hamiltonian:
- by including only the first term, one selects the independent-particle approximation
- by including the first and second terms, one selects the independent-quasiparticle approximation
- by including the first, second terms and the Hartree part of the self-energy, one selects the time-dependent Hartree or random-phase approximation
- by including all terms, one selects the time-dependent Screened-exchange or Bethe-Salpeter equation level of theory.
The type of "experiment" to perform, or of the response one would like to extract, depend on the time-dependent field in input and the post-processing of the Berry's phase polarization. The induced polarization is the sum of the responses of the system to all orders in the external field: [math]\displaystyle{ \bf{P}(t) = \sum_{n=-\infty}^{+\infty} \bf{P}_n (t) }[/math]
In the case we are interested in linear response only, we use a delta-like electric field that excites all frequencies of the system. We choose the intensity of the field [math]\displaystyle{ E }[/math] to be weak enough so that the response is dominated by the first order term. The response at all frequencies is obtained by taking the Fourier-transform of the polarization as
[math]\displaystyle{ \chi(\omega) = \frac{P(\omega)}{E} }[/math].
The response [math]\displaystyle{ \chi }[/math] is related to the dielectric function through the relation [math]\displaystyle{ \epsilon(\omega) = 1 + 4 \pi \chi(\omega) }[/math]. See also Ref. [2].
Notes:
- Length vs Velocity gauge: the equations of motion of Yambo are in the length gauge. Other codes choose instead to work in the velocity gauge. If you want to know more about the advantages and disadvantages of the two gauges read section 2.7 of Ref. [3].
- Advantages and disadvantages with respect to the density-matrix based formalism: Yambo implements two distinct formalisms for studying the real-time response, the one based on density-matrix and the one based on Boch states. Which formalism is best to use depends on what one wants to simulate. The dynamics formulation in terms of density matrix or non-equilibrium Green's functions allows a systematic treatment of correlation effects and electron-phonon coupling[4], but the price to pay is that the coupling with the external field is correct only to the linear order. This formalism is important in all those phenomena where electronic relaxation plays a major role. On the other hand, the formulation in terms of the Bloch states treats coupling with the external field in an exact way even beyond the linear regime, and for this reason, it allows the description of phenomena coherent with the external fields and to access, for example, non-linear responses.
Step 1: Prerequisites
In this example, we will consider a single layer of hexagonal boron nitride (hBN). Before running real time simulations in yambo you need to:
- download input files and Yambo databases for this tutorial here: hBN-2D-RT.tar.gz.
- perform the steps described in Prerequisites for Real Time propagation with Yambo.
Step 2: Independent-particle approximation
Similar to what you have seen in the tutorial for the density matrix formalism (Real time approach to linear response) we start from the lowest level of theory. Here it is assumed you ran the tutorial on the density-matrix formalism and so you are familiar with the optical spectrum of 2D h-BN at this level of theory. You will recognise that the input and output for the Bloch-state formalism are very similar to those for the density-matrix formalism.
Run the simulations
First of all, create a directory for the inputs to be run by yambo_nl:
mkdir Inputs_nl
Then, use the command:
yambo_nl -nl n -F Inputs_nl/01_td_ip.in
to generate the input:
nloptics # [R] Non-linear spectroscopy % NLBands 3 | 6 | # [NL] Bands range % NLverbosity= "high" # [NL] Verbosity level (low | high) NLtime= 55.000000 fs # [NL] Simulation Time NLintegrator= "INVINT" # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC") NLCorrelation= "IPA" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX") NLLrcAlpha= 0.000000 # [NL] Long Range Correction % NLEnRange -1.000000 |-1.000000 | eV # [NL] Energy range % NLEnSteps= 1 # [NL] Energy steps NLDamping= 0.000000 eV # [NL] Damping (or dephasing) RADLifeTime=-1.000000 fs # [RT] Radiative life-time (if negative Yambo sets it equal to Phase_LifeTime in NL) #EvalCurrent # [NL] Evaluate the current #FrPolPerdic # [DIP] Force periodicity of polarization respect to the external field HARRLvcs= 18475 RL # [HA] Hartree RL components EXXRLvcs= 18475 RL # [XX] Exchange RL components % Field1_Freq 0.100000 | 0.100000 | eV # [RT Field1] Frequency % Field1_NFreqs= 1 # [RT Field1] Frequency Field1_Int= 1000.00 kWLm2 # [RT Field1] Intensity Field1_Width= 0.000000 fs # [RT Field1] Width Field1_kind= "DELTA" # [RT Field1] Kind(SIN|COS|RES|ANTIRES|GAUSS|DELTA|QSSIN) Field1_pol= "linear" # [RT Field1] Pol(linear|circular) % Field1_Dir 0.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor % Field1_Tstart= 0.010000 fs # [RT Field1] Initial Time
Note that the standard input of yambo_nl has default settings for nonlinear response (e.g. the SOFTSIN for the Field1_kind (that set the time-dependence 'kind' of the field).
Set the field direction Field1_Dir along y, the field kind to DELTA, the length of the simulation NLtime to 55 fs, select the bands from 3 to 6 (NLBands) and set the NLDamping to zero. The fields you need to change with respect to the default settings are shown above in red.
Now run
yambo_nl -F Inputs_nl/01_td_ip.in -J TD-IP_nl -C TD-IP_nl
The code produces different files (in the TD-IP_nl directory, all with the suffix TD-IP_nl: o.polarization_F1 that contains the polarization, o.external_potential_F1 the external field we used, and finally r_optics_nloptics a report with all information about the simulation.
The simulation takes longer than the density-matrix counterpart. This is the price to pay for computing the polarization as a Berry Phase. While there is no gain for calculating the linear response, this allows computing nonlinear coefficients, which cannot be obtained within the density based approach. Have a look at the time profile of the run towards the end of the r_optics_nloptics (absolute time may differ):
NL Berry Pol NEQ : 2.2409s CPU ( 5502 calls, 0.000 sec avg) io_WF : 3.8788s CPU ( 757 calls, 0.005 sec avg) DIPOLE_overlaps : 4.1834s CPU DIPOLE_buil_covariants : 5.5611s CPU Dipoles : 7.3264s CPU
Most of the time is spent computing the covariant dipoles and the non-equilibrium Berry Polarization, which involves the inversion of small matrices via the Lapack inversion routine.
Some time is also required by the integrator NL_Integrator of the equation of motion inversion, which is coded via the solution of a linear system of equations as implemented in the Lapack library. This integrator is more accurate than the Runge-Kutta 2nd order employed for the corresponding calculation within the density matrix formalism.
Plot the third column of o.polarization_F1 versus the first one (time-variable) to get the time-dependent polarization along the y-direction:
For comparison, in the figure above, the polarization is obtained either through the Bloch-states or the density-matrix. There are few differences which could be due to the following reasons:
- Coupling with the external field: the
yambo_rtevaluates the coupling with the external field through the dipoles (by default computing the commutator [x,Hnl]; theyambo_nlevaluates the coupling using an expression based on the Berry Phase. For the latter, a numerical differentiation in reciprocal space is performed. The 10x10 2D k-mesh may not be fine enough. For denser k-mesh, the Berry phase approach converges to the one based on the dipoles. To avoid this one can useDipApproach="Covariant"in the scheme used by yambo_rt. - Integrator: the integrator used for the
yambo_nlis more accurate. The fact that the differences are more pronounced for long times may hint to inaccuracies in the integration of the equation of motions inyambo_rt. One can use a similar integrator with yambo_rt, by setting in inputIntegrator= "INV RK2"
Output postprocessing: the dielectric function
Once we obtained the polarization, we Fourier transform the polarization to obtain the dielectric function.
Use the command:
ypp_nl -nl -F Inputs_nl/ypp_abs.in
to generate the input file:
nonlinear # [R] NonLinear Optics Post-Processing Xorder= 1 # Max order of the response functions % TimeRange -1.000000 |-1.000000 | fs # Time-window where processing is done % ETStpsRt= 1001 # Total Energy steps % EnRngeRt 0.00000 | 10.00000 | eV # Energy range % DampMode= "LORENTZIAN" # Damping type ( NONE | LORENTZIAN | GAUSSIAN ) DampFactor= 0.10000 eV # Damping parameter
where we set a Lorentzian smearing corresponding to 0.1 eV. Note that due to the finite time of our simulation, we need to introduce a finite smearing to Fourier transform the result.
To run the post-processing, use the command:
ypp_nl -F Inputs_nl/ypp_abs.in -J TD-IP_nl -C TD-IP_nl
and obtain the files for the dielectric constant along with the field direction, the EELS along with the same direction, and the damped polarization.
o-TD-IP_nl.YPP-eps_along_E o-TD-IP_nl.YPP-eels_along_E o-TD-IP_nl.YPP-damped_polarization
inside the folder TD-IP_nl
Plot the dielectric constant and compare it with the linear response and the density-matrix formalism (available from the Real time approach to linear response tutorial):
gnuplot plot "CHI_IP/o-CHI_IP.eps_q1_ip" u 1:2 w l rep "TD-IP_rt/o-TD-IP_rt.YPP-eps_along_E" u 1:2 w l rep "TD-IP_nl/o-TD-IP_nl.YPP-eps_along_E" u 1:2 w l
The differences observed with the spectra obtained from either the density-matrix and linear-response formalisms are due to how the dipoles are evaluated. Within the density-matrix and linear-response formalisms, dipoles are computed in the same way. The Bloch-states formalism uses instead covariant dipoles evaluated numerically on the k-mesh. The differences are due to numerical errors in evaluating the covariant dipoles. These differences disappear with a denser k-mesh.
It is a good practice to converge the dielectric function obtained from the Bloch-state dynamics with respect to the k-points using the linear-response spectrum as a reference. This gives the minimum number of k-points to be used to obtain accurate results.
Step 3: Hartree+Screened exchange approximation
Evaluating the collisions
The evaluation of the collisions. (see Introduction to Real Time propagation in Yambo to remember what are the collisions)
This part is common in between the yambo_nl and the yambo_rt. If you have computed the collisions for the corresponding yambo_rt, you can re-used those results. If not, follow the steps in Real time Bethe-Salpeter Equation (density matrix only).
Run the simulations
To generate the input for the real-time simulation you can run
yambo_nl -nl n -V qp -F Inputs_nl/04_td-sex.in
The input file is
nloptics # [R NL] Non-linear optics % NLBands 4 | 5 | # [NL] Bands % NLverbosity= "high" # [NL] Verbosity level (low | high) NLtime= 20.00000 fs # [NL] Simulation Time NLintegrator= "CRANKNIC" # [NL] Integrator ("EULEREXP/RK2/RK4/RK2EXP/HEUN/INVINT/CRANKNIC") NLCorrelation= "SEX" # [NL] Correlation ("IPA/HARTREE/TDDFT/LRC/LRW/JGM/SEX") NLLrcAlpha= 0.000000 # [NL] Long Range Correction % NLEnRange -1.000000 | -1.000000 | eV # [NL] Energy range % NLEnSteps= 1 # [NL] Energy steps NLDamping= "0.10000 eV # [NL] Damping #EvalCurrent # [NL] Evaluate the current #FrPolPerdic # [DIP] Force periodicity of polarization respect to the external field HARRLvcs= 1000 mHa # [HA] Hartree RL components EXXRLvcs= 1000 mHa # [XX] Exchange RL components
The next section of the input is about external quasiparticle corrections (EXTQP G). There change only
% GfnQP_E
3.000000 | 1.000000 | 1.000000 | # [EXTQP G] E parameters (c/v) eV|adim|adim
%
The latter line introduces a scissor operator (a rigid shift of the conduction bands) of 3.0 eV. In principle, it is possible to perform a G0W0 calculation with Yambo and use the Quasi-particle band structure instead of the rigid shift.
The last section regards the applied field (RT Field1). There change these two lines:
% Field1_Dir 0.000000 | 1.000000 | 0.000000 | # [RT Field1] Versor % Field1_kind= "DELTA" # [RT Field1] Kind(SIN|SOFTSIN|RES|ANTIRES|GAUSS|DELTA|QSSIN)
Run the calculation:
nohup yambo_nl -F Inputs_nl/04_td_hsex.in -J "TD-SEX_nl,COLL_SEX" -C TD-SEX_nl
Output postprocessing: dielectric function
The post-processing of the output does not depend on the level of approximation. The same steps are taken as in the independent-particle approximation.
Use the command:
ypp_nl -F Inputs_nl/ypp_abs.in -J TD-SEX_nl -C TD-SEX_nl
Plot the results (compare with results obtained from the density-matrix formalism):
gnuplot> plot "TD-HSEX_rt/o-TD-HSEX_rt.YPP-eps_along_E" u 1:2 w l gnuplot> rep "TD-HSEX_nl/o-TD-HSEX_nl.YPP-eps_along_E" u 1:2 w l gnuplot> set xrange [3:10] gnuplot> set yrange [0:12] gnuplot> rep "TD-IP_rt/o-TD-IP_rt.YPP-eps_along_E" u ($1+3):2 w l
Similarly to what we have seen in the independent-particle case, the differences between the two approaches are due to the accuracy in evaluating numerically the covariant dipoles in the Bloch-state dynamics approach.
For comparison, the linear response results at the BSE level can be obtained following the BSE tutorial.
References
- ↑ I. Souza, J. Íñiguez, and D. Vanderbilt, PRB 69, 085106 (2004)
- ↑ 2.0 2.1 C. Attaccalite and M. Gruning Rev. B, 88, 235113 (2013)
- ↑ C. Attaccalite, arXiv 1609.09639 (2017)
- ↑ D. Sangalli, and A. Marini EPL 110, 47004 (2015)