TDDFT For Bulk Tutorial

The program mbpt_lcao is capable of performing a linear response TDDFT calculation for extended systems, i.e. for solids, chains, surfaces etc. The output of the program is macroscopic dielectric function $\epsilon_{\mathrm{M}}(\mathbf{q}, \omega)$ for a given set of frequencies $\omega$ and for a given momentum transfer $\mathbf{q}$.

The iterative method of computing of macroscopic dielectric function is described elsewhere NIMB paper 2015. The purpose of this tutorial is mainly to describe relevant input parameters of the input file.

In the first step, we will need to perform a DFT calculation. The SIESTA calculation should deliver Hamiltonian and Overlap matrices without so-called ''folding'' feature. The feature affects only Gamma-only calculations. As a result of ''folding'', the matrix elements get squizzed into rectangular matrices, i.e. Hamiltonian of the system seems similar to a finite-system calculation. A simple way to avoid ''folding'' is to demand band structure plotting in the SIESTA input file (.fdf file). For instance, one may insert into the .fdf file a block BandLines

%block BandLines
1   0.0    0.0    0.0      Gamma
10  0.5    0.0    0.0      X ?
%endblock BandLines


In mbpt_lcao a simple Monkhorst-Pack (MP) sampling of Brilloine zone (BZ) will be used. Therefore, if you want to be totally consistent with DFT—TDDFT calculation, you should use MP sampling also in SIESTA. For instance, if you plan to use 3x3x3 MP sampling in tddft_bulk, then a block kgrid_Monkhorst_Pack can be used in .fdf file

%block kgrid_Monkhorst_Pack
3 0 0 0
0 3 0 0
0 0 3 0
%endblock kgrid_Monkhorst_Pack


However, mbpt_lcao is using generally its own BZ sampling and one only should worry about convergence of density with BZ sampling in SIESTA calculation. Only real-space data of Hamiltonian matrix and overlap matrix will be used on TDDFT side. The export of data must be caused by two options in .fdf file

COOP.Write     .true.  # Write .WFSX and .HSX files
WriteDenchar   .true.  # Write .DIM and .PLD files


Below you will find an example of .fdf file for calculation of Graphene

SystemName          graphene

DM.UseSaveDM        true     # to use continuation files
NumberOfAtoms       2
NumberOfSpecies     1

%block ChemicalSpeciesLabel
1  6  C      # Species index, atomic number, species label
%endblock ChemicalSpeciesLabel

LatticeConstant  2.46 Ang  ## This is only valid for graphene case

%block LatticeVectors
0.86600   0.500000  0.00000
0.86600  -0.500000  0.00000
0.00000   0.000000  8.08081
%endblock LatticeVectors

AtomicCoordinatesFormat  Ang

%block AtomicCoordinatesAndAtomicSpecies
0.000000     0.00000  0.00000 1
1.420282     0.00000  0.00000 1
%endblock AtomicCoordinatesAndAtomicSpecies

%block kgrid_Monkhorst_Pack
64 0  0  0.0
0 64  0  0.0
0  0  1  0.0
%endblock kgrid_Monkhorst_Pack

ElectronicTemperature  10 K

MeshCutoff       200.0 Ry

BandLinesScale      pi/a               # Default Value
%block BandLines
1  0.577350269   1.00000    0.0000   \M
30  0.0000        0.0000     0.0000  \Gamma
30  0.0000    2.309401077    0.0000  \K
30  0.577350269   1.00000    0.0000  \M
%endblock BandLines

PAO.EnergyShift       50 meV
MaxSCFIterations      500            # Maximum number of SCF iter
DM.Tolerance          1e-5          # Tolerance in maximum difference
DM.MixingWeight       0.12           # New DM amount for next SCF cycle
DM.NumberPulay        4
COOP.Write               .true. ### Exports .WFSX, .HSX
WriteDenchar             .true. ### Exports .DIM and .PLD files


Please, copy this text into siesta.fdf file and run

siesta <siesta.fdf |tee siesta.out


After a successful SIESTA run, a set of output files will appear in the current directory. On TDDFT side, the files siesta.EIG, siesta.HSX, siesta.WFSX (or siesta.fullBZ.WFSX), siesta.DIM, siesta.PLD and C.ion will be used. You may not bother with copying of these files into other place and just run mbpt_lcao standing in the current directory. mbpt_lcao is controlled through another, own, input file. The input file must be tddft_lr.inp.

Besides the other parameters, one may define the vector of momentum transfer $\mathbf{q}$, set of frequencies $\omega$, broadening constant $\varepsilon$, BZ sampling for TDDFT and other more technical parameters which affect TDDFT calculation.

mbpt_lcao supports several calculations which serve as cross checks. The most ''advanced'' would be so-called iterative calculation. This calculation must be demanded in the input explicitly. This must be done by setting parameter do_eels_iter4 to a positive value

do_eels_iter4  1


The calculation invoked by do_eels_iter4 is mainly realized with single-precision arithmetics. Single precision is most advanced option in terms of computational speed and required memory.

Momentum transfer $\mathbf{q}$ is defined with a 3-vector parameter qvec (three numbers must follow). For example, in order to define momentum transfer 0.046 $\mathrm{Ang}^{-1}$ along $\Gamma-\mathrm{M}$ direction, in graphene plane (this has been used in the publication we reproducing PhysRevB.83.245122), one could specify in the input file

qvec_type     CARTESIAN
qvec 0.0243422     0.0     0.0


Currently, the units of momentum transfer must be $\mathrm{Bohr}^{-1}$ only.

The set of frequencies $\omega$ is defined through number of frequencies nff and frequency spacing d_omega_win1. Optionally, one may also define minimum and maximum limits for the set of frequencies with parameters omega_min_tddft and omega_max_tddft. Default units for these parameters electron volts (eV). The broadening constant $\varepsilon$ is given through parameter freq_eps_win1 (by default also in eV). For instance, to reproduce the electron-energy loss spectrum (EELS) and dielectric function of graphene from PhysRevB.83.245122, we found was necessary a broadening of $0.24$eV, frequency spacing $0.12$eV, number of frequencies 292

nff                      292
d_omega_win1    0.12
freq_eps_win1     0.24


The sensible values of parameters omega_min_tddft, omega_max_tddft must be in the range $0 \ldots$nff$\times$d_omega_win1. The frequency broadening freq_eps_win1 must be twice or more larger than frequency spacing d_omega_win1.

TDDFT has its own sampling of BZ. This is controlled through a 3-vector parameter BZ_sampling. For instance, for graphene example it would be sensible to set

BZ_sampling 64 64 1


The plane waves (PW) are used in the computation of (periodic) Hartree kernel $K^{\mu\nu}(\mathbf{q})$. The set of PW is defined by BZ and by an energy cutoff. Similarly to GPAW package, we use a parameter Ecut to define PW cutoff. Sensible values for PW cutoff in SIESTA, LCAO calculations would be from 150 eV to 300 eV. Larger cutoff is not necessary because pseudo-atomic orbitals used in SIESTA do not cause steep changes in space. Typically, the calculation is well converged with Ecut equal 200 eV

Ecut 200


The product basis which will be used in the iterative calculation is a local basis, not PW basis. There are several options to generate this basis. We recommend to use an atom-centered basis set, which can be demanded via parameter prod_basis_type

prod_basis_type MIXED


The size of atom-centered basis is controlled through an eigenvalue cutoff for local products of atomic orbitals

eigmin_local    1E-5   !! Threshold for local products.


During generation of atom-centered basis, we also use so-called bilocal dominant functions. These functions can be more numerous. The number of bilocal functions is controlled with another eigenvalue cutoff

eigmin_bilocal  1E-8   !! Threshold for bilocal products.


The whole input file tddft_lr.inp may look as following

iv           1    # Verbosity

do_eels_iter4  1

eigmin_local    1E-5   !! Threshold for local products.
eigmin_bilocal  1E-8   !! Threshold for bilocal products.

prod_basis_type  MIXED

qvec_type     CARTESIAN
qvec 0.0243422     0.0     0.0  ! must be Thygesen's 0.046 A^{-1}

BZ_sampling 64 64 1

Ecut          200 eV

nff             292
d_omega_win1    0.12
freq_eps_win1    0.24


Please, copy this text into the file tddft_lr.inp and start mbpt_lcao.
After a while (can be an hour in this graphene example), you will obtain files eels_nlfe_iter4.txt, eels_lfe_iter4.txt, epsm_nlfe_iter4.txt, epsm_lfe_iter4.txt. The files contain EELS without local field effects (nlfe), EELS with local field effects (lfe), macroscopic dielectric function without local field effects and macroscopic dielectric function with local field effects, correspondingly. The first two files contain two columns: set of frequencies $\omega$ and EELS signal $-\mathrm{Im}\frac{1}{\epsilon_{\mathrm{M}}(\mathbf{q}, \omega)}$. The files with epsm_ prefix contain three columns: for set of frequencies, real and imaginary parts of $\epsilon_{\mathrm{M}}(\mathbf{q}, \omega)$.

Below is a comparison of EELS computed with GPAW with the calculation we did in this tutorial