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)$.