lcaotddft.rst 13.9 KB
 askhl committed Mar 24, 2015 1 2 .. _lcaotddft:  kuismam committed Jun 25, 2015 3 ===================================================  Tuomas Rossi committed Jan 16, 2018 4 Time-propagation TDDFT with LCAO: Theory and usage  kuismam committed Jun 25, 2015 5 ===================================================  tprossi committed Jun 18, 2015 6   askhl committed Mar 24, 2015 7 This page documents the use of time-propagation TDDFT in :ref:LCAO  askhl committed Jun 29, 2015 8 mode . The implementation is described in Ref. [#Kuisma2015]_.  tprossi committed Jun 18, 2015 9   Tuomas Rossi committed Jan 16, 2018 10 Real-time propagation of LCAO functions  kuismam committed Jun 25, 2015 11 12 =======================================  askhl committed Jun 29, 2015 13 14 15 In the real-time LCAO-TDDFT approach, the time-dependent wave functions are represented using the localized basis functions \tilde{\phi}_{\mu}(\mathbf r) as  kuismam committed Jun 25, 2015 16 17 18  .. math::  askhl committed Jun 29, 2015 19  \tilde{\psi}_n(\mathbf{r},t) = \sum_{\mu} \tilde{\phi}_{\mu}(\mathbf{r}-\mathbf{R}^\mu) c_{\mu n}(t) .  kuismam committed Jun 25, 2015 20   askhl committed Jun 29, 2015 21 The time-dependent Kohn--Sham equation in the PAW formalism can be written as  kuismam committed Jun 25, 2015 22 23 24 25 26  .. math:: \left[ \widehat T^\dagger \left( -i \frac{{\rm d}}{{\rm d}t} + \hat H_{\rm KS}(t) \right) \widehat T \right] \tilde{\Psi(\mathbf{r},t)} = 0.  askhl committed Jun 29, 2015 27 28 From this, the following matrix equation can be derived for the LCAO wave function coefficients:  kuismam committed Jun 25, 2015 29 30 31 32  .. math:: {\rm i}\mathbf{S} \frac{{\rm d}\mathbf{C}(t)}{{\rm d}t} = \mathbf{H}(t) \mathbf{C}(t).  askhl committed Jun 29, 2015 33 34 35 36 37 38 In the current implementation, \mathbf C, \mathbf S and \mathbf H are dense matrices which are distributed using ScaLAPACK/BLACS. Currently, a semi-implicit Crank--Nicolson method (SICN) is used to propagate the wave functions. For wave functions at time t, one propagates the system forward using \mathbf H(t) and solving the linear equation  kuismam committed Jun 25, 2015 39 40 41  .. math::  askhl committed Jun 29, 2015 42  \left( \mathbf{S} + {\rm i} \mathbf{H}(t) {\rm d}t / 2 \right) \mathbf{C}'(t+{\rm d}t) = \left( \mathbf{S} - {\rm i} \mathbf{H}(t) {\rm d}t / 2 \right) \mathbf{C}(t).  kuismam committed Jun 25, 2015 43   askhl committed Jun 29, 2015 44 45 46 Using the predicted wave functions C'(t+\mathrm dt), the Hamiltonian H'(t+\mathrm dt) is calculated and the Hamiltonian at middle of the time step is estimated as  kuismam committed Jun 25, 2015 47 48 49  .. math::  kuismam committed Jun 26, 2015 50  \mathbf{H}(t+{\rm d}t/2) = (\mathbf{H}(t) + \mathbf{H}'(t+{\rm d}t)) / 2  kuismam committed Jun 25, 2015 51   askhl committed Jun 29, 2015 52 53 With the improved Hamiltonian, the wave functions are again propagated from t to t+\mathrm dt by solving  kuismam committed Jun 09, 2015 54   kuismam committed Jun 25, 2015 55 56 .. math::  askhl committed Jun 29, 2015 57  \left( \mathbf{S} + {\rm i} \mathbf{H}(t+{\rm d}t/2) {\rm d}t / 2 \right) \mathbf{C}(t+{\rm d}t) = \left( \mathbf{S} - {\rm i} \mathbf{H}(t+{\rm d}t/2) {\rm d}t / 2 \right) \mathbf{C}(t).  tprossi committed Jun 11, 2015 58   askhl committed Jun 29, 2015 59 60 This procedure is repeated using 500--2000 time steps of 5--40 as to calculate the time evolution of the electrons.  tprossi committed Jun 11, 2015 61   Tuomas Rossi committed Jan 18, 2018 62 63 64 65 .. _example: Example usage =============  kuismam committed Jun 09, 2015 66   Tuomas Rossi committed Jan 16, 2018 67 First do a standard ground-state calculation with the GPAW calculator:  kuismam committed Jun 09, 2015 68   Tuomas Rossi committed Jan 16, 2018 69 70 .. literalinclude:: lcaotddft.py :lines: 3-21  kuismam committed Jun 09, 2015 71   askhl committed Jun 29, 2015 72 Some important points are:  kuismam committed Jun 09, 2015 73   Tuomas Rossi committed Jan 16, 2018 74 75 76 77 78 79 80 81 82 83 84 85 86 87 * The grid spacing is only used to calculate the Hamiltonian matrix and therefore a coarser grid than usual can be used. * Completely unoccupied bands should be left out of the calculation, since they are not needed. * The density convergence criterion should be a few orders of magnitude more accurate than in usual ground-state calculations. * The convergence tolerance of the Poisson solver should be at least 1e-16, but 1e-20 does not hurt (note that this is the **quadratic** error). * One should use multipole-corrected Poisson solvers or other advanced Poisson solvers in any TDDFT run. See the documentation on :ref:advancedpoisson. Next the calculation proceeds as in the grid mode with TDDFT object. We kick the system in the z direction and propagate 3000 steps of 10 as:  kuismam committed Jun 09, 2015 88   Tuomas Rossi committed Jan 16, 2018 89 90 .. literalinclude:: lcaotddft.py :lines: 23-35  kuismam committed Jun 09, 2015 91   Tuomas Rossi committed Jan 16, 2018 92 After the time propagation, the spectrum can be calculated:  kuismam committed Jun 09, 2015 93   Tuomas Rossi committed Jan 16, 2018 94 95 .. literalinclude:: lcaotddft.py :lines: 37-39  kuismam committed Jun 09, 2015 96   Tuomas Rossi committed Jan 18, 2018 97 98 The previous example as a complete script can be downloaded here: :download:lcaotddft.py.  kuismam committed Jun 09, 2015 99 100 101 102 103  General notes about basis sets ==============================  askhl committed Jun 29, 2015 104 In time-propagation LCAO-TDDFT, it is much more important to think  Tuomas Rossi committed Jan 16, 2018 105 about the basis sets compared to ground-state LCAO calculations. It  askhl committed Jun 29, 2015 106 is required that the basis set can represent both the occupied  kuismam committed Jul 30, 2015 107 (holes) and relevant unoccupied states (electrons) adequately. Custom  askhl committed Jun 29, 2015 108 109 basis sets for the time propagation should be generated according to one's need, and then benchmarked.  kuismam committed Jun 09, 2015 110   jensj committed Jun 26, 2015 111 **Irrespective of the basis sets you choose, ALWAYS, ALWAYS, benchmark LCAO  askhl committed Jun 29, 2015 112 113 114 115 116 117 results with respect to the FD time-propagation code** on the largest system possible. For example, one can create a prototype system which consists of similar atomic species with similar roles as in the parent system, but small enough to calculate with grid propagation mode. See Figs. 4 and 5 of Ref. [#Kuisma2015]_ for example benchmarks.  Tuomas Rossi committed Jan 16, 2018 118 119 120 After these remarks, we describe two packages of basis sets that can be used as a starting point for choosing a suitable basis set for your needs. Namely, :ref:pvalence basis sets and :ref:coopt basis sets.  tprossi committed Jun 18, 2015 121   Jens Jørgen Mortensen committed Mar 10, 2016 122   tprossi committed Jun 18, 2015 123 124 125 126 127 .. _pvalence basis sets: p-valence basis sets --------------------  Jens Jørgen Mortensen committed Mar 10, 2016 128 129 130 131 132 133 134 135 136 137 The so-called p-valence basis sets are constructed for transition metals by replacing the p-type polarization function of the default basis sets with a bound unoccupied p-type orbital and its split-valence complement. Such basis sets correspond to the ones used in Ref. [#Kuisma2015]_. These basis sets significantly improve density of states of unoccupied states. The p-valence basis sets can be easily obtained for the appropriate elements with the :command:gpaw install-data tool using the following options:: $gpaw install-data {} --basis --version=pvalence  138   Jens Jørgen Mortensen committed Mar 15, 2016 139 See :ref:installation of paw datasets for more information on basis set  Jens Jørgen Mortensen committed Mar 10, 2016 140 141 142 installation. It is again reminded that these basis sets are not thoroughly tested and **it is essential to benchmark the performance of the basis sets for your application**.  tprossi committed Jun 18, 2015 143 144 145 146 147 148 149  .. _coopt basis sets: Completeness-optimized basis sets ---------------------------------  jensj committed Jun 26, 2015 150 151 152 A systematic approach for improving the basis sets can be obtained with the so-called completeness-optimization approach. This approach is used in Ref. [#Rossi2015]_ to generate basis set series for TDDFT calculations of copper,  kuismam committed Jun 25, 2015 153 silver, and gold clusters.  tprossi committed Jun 18, 2015 154   jensj committed Jun 26, 2015 155 156 For further details of the basis sets, as well as their construction and performance, see [#Rossi2015]_. For convenience, these basis sets can be easily  Jens Jørgen Mortensen committed Mar 10, 2016 157 obtained with::  158   Jens Jørgen Mortensen committed Mar 10, 2016 159 $ gpaw install-data {} --basis --version=coopt  160   Jens Jørgen Mortensen committed Mar 15, 2016 161 See :ref:installation of paw datasets for basis set installation. Finally,  Jens Jørgen Mortensen committed Mar 10, 2016 162 163 it is again emphasized that when using the basis sets, **it is essential to benchmark their suitability for your application**.  tprossi committed Jun 18, 2015 164   kuismam committed Jun 09, 2015 165 166 167 168  Parallelization ===============  Tuomas Rossi committed Jan 16, 2018 169 170 LCAO-TDDFT is parallelized using ScaLAPACK. It runs without ScaLAPACK, but in this case only a single core is used for linear alrebra.  kuismam committed Jun 25, 2015 171   Tuomas Rossi committed Jan 16, 2018 172 173 174 175 176 177 178 * Use parallel={'sl_default':(N, M, 64)}; See :ref:manual_parallel. * It is necessary that N*M equals the total number of cores used by the calculator, and max(N,M)*64 < nbands, where 64 is the used block size. The block size can be changed to, e.g., 16 if necessary. * Apart from parallelization of linear algrebra, normal domain and band parallelizations can be used. As in ground-state LCAO calculations, use band parallelization to reduce memory consumption.  kuismam committed Jun 09, 2015 179   tprossi committed Jun 11, 2015 180   Tuomas Rossi committed Jan 16, 2018 181 182 183 184 185 186 .. TODO Timing ====== Add ParallelTimer example  tprossi committed Jun 18, 2015 187   tprossi committed Jun 11, 2015 188   Tuomas Rossi committed Jan 18, 2018 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 Advanced analysis tools ======================= In :ref:example it was demonstrated how to calculate photoabsorption spectrum from the time-dependent dipole moment data collected with DipoleMomentWriter() observer. The code is not limited to this analysis but any (also user-written) analysis tools can be embedded in the general time-propagation framework. Here we describe some analysis tools available in GPAW. We use a finite sodium atom chain as an example system. First, let's do the ground-state calculation: .. literalinclude:: lcaotddft_Na8/gs.py Note the recommended use of :ref:advancedpoisson and :ref:pvalence basis sets. Recording the wave functions and replaying the time propagation --------------------------------------------------------------- We can record the time-dependent wave functions during the propagation with WaveFunctionWriter() observer: .. literalinclude:: lcaotddft_Na8/td.py .. tip:: The time propagation can be in the same manner from the restart file: .. literalinclude:: lcaotddft_Na8/tdc.py The created wfw.ulm file contains the time-dependent wave functions C_{\mu n}(t) that define the state of the system at each time instance. We can use that file to replay the time propagation: .. literalinclude:: lcaotddft_Na8/td_replay.py The update keyword in propagator has following options: ============== =============================== update variables updated during replay ============== =============================== 'all' Hamiltonian and density 'density' density 'none' nothing ============== ===============================  Tuomas Rossi committed Jan 16, 2018 237 238   Tuomas Rossi committed Jan 18, 2018 239 240 241 242 243 244 245 246 247 248 249 Kohn--Sham decomposition of the density matrix ---------------------------------------------- Kohn--Sham decomposition is an illustrative way of analyzing electronic excitations in Kohn--Sham electron-hole basis. For demonstration and description of the implementation, see Ref. [#Rossi2017]_. Here we demonstrate how to construct the photoabsorption decomposition at a specific frequency in Kohn--Sham electon-hole basis.  Tuomas Rossi committed Jan 18, 2018 250 251 First, let's calculate and :download:plot  the spectrum:  Tuomas Rossi committed Jan 18, 2018 252 253 254  .. literalinclude:: lcaotddft_Na8/spectrum.py  Tuomas Rossi committed Jan 18, 2018 255 256 257 258 .. image:: lcaotddft_Na8/spec.png :scale: 70% The two main resonances are analyzed in the following.  Tuomas Rossi committed Jan 18, 2018 259 260 261 262 263 264 265 266 267 268  Frequency-space density matrix """""""""""""""""""""""""""""" We generate the density matrix for the frequencies of interest: .. literalinclude:: lcaotddft_Na8/td_fdm_replay.py .. tip::  Tuomas Rossi committed Jan 18, 2018 269 270 271  Instead of replaying the propagation, one can do the same analysis on-the-fly by attaching the analysis tools to the usual time-propagation calculation.  Tuomas Rossi committed Jan 18, 2018 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292  Transform the density matrix to Kohn--Sham electron-hole basis """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" First, we construct the Kohn--Sham electron-hole basis. For that we need to calculate unoccupied Kohn--Sham states, which is conveniently done by restarting from the earlier ground-state file: .. literalinclude:: lcaotddft_Na8/ksd_init.py Next, we can use the created objects to transform the LCAO density matrix to the Kohn--Sham electron-hole basis and visualize the photoabsorption decomposition as a transition contribution map (TCM): .. literalinclude:: lcaotddft_Na8/tcm_plot.py Note that the sum over the decomposition (the printed total absorption) equals to the value of the photoabsorption spectrum at the particular frequency in question.  Tuomas Rossi committed Jan 18, 2018 293 294 We obtain the following transition contributions for the resonances (presented both as tables and TCMs):  Tuomas Rossi committed Jan 18, 2018 295 296 297 298  .. literalinclude:: lcaotddft_Na8/table_1.12.txt :language: none  Tuomas Rossi committed Jan 18, 2018 299 300 301 .. image:: lcaotddft_Na8/tcm_1.12.png :scale: 70%  Tuomas Rossi committed Jan 18, 2018 302 303 304 .. literalinclude:: lcaotddft_Na8/table_2.48.txt :language: none  Tuomas Rossi committed Jan 18, 2018 305 306 307 .. image:: lcaotddft_Na8/tcm_2.48.png :scale: 70%  Tuomas Rossi committed Jan 16, 2018 308   kuismam committed Jun 25, 2015 309 Induced density  Tuomas Rossi committed Jan 16, 2018 310 ---------------  kuismam committed Jun 25, 2015 311   Tuomas Rossi committed Jan 18, 2018 312 313 The density matrix gives access to any other quantities. For instance, the induced density can be conveniently obtained:  kuismam committed Jun 25, 2015 314   Tuomas Rossi committed Jan 18, 2018 315 .. literalinclude:: lcaotddft_Na8/ksd_ind.py  kuismam committed Jun 25, 2015 316   Tuomas Rossi committed Jan 18, 2018 317 318 319 The resulting cube files can be visualized, for example, with :download:this script  producing the figures:  kuismam committed Jun 25, 2015 320   Tuomas Rossi committed Jan 18, 2018 321 322 .. image:: lcaotddft_Na8/ind_1.12.png :scale: 70%  kuismam committed Jun 25, 2015 323   Tuomas Rossi committed Jan 18, 2018 324 325 .. image:: lcaotddft_Na8/ind_2.48.png :scale: 70%  kuismam committed Jun 25, 2015 326 327   Tuomas Rossi committed Jan 18, 2018 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 Advanced tutorials ================== Plasmon resonance of silver cluster ----------------------------------- One should think about what type of transitions of interest are present, and make sure that the basis set can represent such Kohn-Sham electron and hole wave functions. The first transitions in silver clusters will be 5s \rightarrow 5p like. We require 5p orbitals in the basis set, and thus, we must generate a custom basis set. Here is how to generate a double-zeta basis set with 5p orbital in valence for silver for GLLB-SC potential. Note that the extra 5p valence state effectively improves on the ordinary polarization function, so this basis set is **better** than the default double-zeta polarized one. We will use the 11-electron Ag setup, since the semi-core p states included in the default setup are not relevant here. .. literalinclude:: lcaotddft_basis.py We calculate the icosahedral Ag55 cluster: :download:ag55.xyz This code uses ScaLAPACK parallelization with 48 cores. .. literalinclude:: lcaotddft_ag55.py Code runs for approximately two hours wall-clock time. The resulting spectrum shows already emerging plasmonic excitation around 4 eV. For more details, see [#Kuisma2015]_. .. image:: fig1.png .. TODO Large organic molecule ---------------------- On large organic molecules, on large conjugated systems, there will \pi \rightarrow \pi^*, \sigma \rightarrow \sigma^*. These states consist of only the valence orbitals of carbon, and they are likely by quite similar few eV's below and above the fermi lavel. These are thus a reason to believe that these states are well described with hydrogen 1s and carbon 2s and 2p valence orbitals around the fermi level. Here, we will calculate a small and a large organic molecule with lcao-tddft.  tprossi committed Jun 11, 2015 377 378 379 380 References ========== .. [#Kuisma2015]  Tuomas Rossi committed Jan 18, 2018 381 382 383 384 385  M. Kuisma, A. Sakko, T. P. Rossi, A. H. Larsen, J. Enkovaara, L. Lehtovaara, and T. T. Rantala, Localized surface plasmon resonance in silver nanoparticles: Atomistic first-principles time-dependent density functional theory calculations,  tprossi committed Jun 18, 2015 386  *Phys. Rev. B* **69**, 245419 (2004).  Tuomas Rossi committed Jan 18, 2018 387  doi:10.1103/PhysRevB.91.115431 _  tprossi committed Jun 18, 2015 388 389 390  .. [#Rossi2015] T. P. Rossi, S. Lehtola, A. Sakko, M. J. Puska, and R. M. Nieminen,  Tuomas Rossi committed Jan 18, 2018 391 392  Nanoplasmonics simulations at the basis set limit through completeness-optimized, local numerical basis sets,  tprossi committed Jun 18, 2015 393  *J. Chem. Phys.* **142**, 094114 (2015).  Tuomas Rossi committed Jan 18, 2018 394 395 396 397 398 399 400 401 402  doi:10.1063/1.4913739 _ .. [#Rossi2017] T. P. Rossi, M. Kuisma, M. J. Puska, R. M. Nieminen, and P. Erhart, Kohn--Sham Decomposition in Real-Time Time-Dependent Density-Functional Theory: An Efficient Tool for Analyzing Plasmonic Excitations, *J. Chem. Theory Comput.* **13**, 4779 (2017). doi:10.1021/acs.jctc.7b00589 _