lcaotddft.rst 10.9 KB
 askhl committed Mar 24, 2015 1 2 .. _lcaotddft:  kuismam committed Jun 25, 2015 3 4 5 =================================================== Time-propagation TDDFT with LCAO : Theory and usage ===================================================  tprossi committed Jun 18, 2015 6   askhl committed Mar 24, 2015 7 This page documents the use of time-propagation TDDFT in :ref:LCAO  tprossi committed Jun 18, 2015 8 9 mode . The implementation is described in [#Kuisma2015]_.  kuismam committed Jun 25, 2015 10 11 12 Real time propagation of LCAO-functions =======================================  jensj committed Jun 26, 2015 13 In real time LCAO-TDDFT approach, the time-dependent wave functions are  kuismam committed Jun 25, 2015 14 represented using localized basis sets as  kuismam committed Jun 25, 2015 15 16 17  .. math::  kuismam committed Jun 25, 2015 18  \tilde{\Psi}_i(\mathbf{r},t) = \sum_{\mu} C_{\mu i}(t) \tilde{\phi}^{\mu}(\mathbf{r}-\mathbf{R}^\mu).  kuismam committed Jun 25, 2015 19 20 21 22 23 24 25  The TD-Kohn-Sham equation in PAW formalism can be written as .. 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.  jensj committed Jun 26, 2015 26 Using these equations, following matrix equation can be derived for LCAO wave  kuismam committed Jun 25, 2015 27 function coefficients  kuismam committed Jun 25, 2015 28 29 30 31  .. math:: {\rm i}\mathbf{S} \frac{{\rm d}\mathbf{C}(t)}{{\rm d}t} = \mathbf{H}(t) \mathbf{C}(t).  jensj committed Jun 26, 2015 32 33 34 In current implementation in GPAW, C, S and H are full matrices, which are parellelized using ScaLAPACK. Currently semi implicit Crank-Nicholson method (SICN) is used to propagate wave functions. For wave functions at time t, one  kuismam committed Jun 25, 2015 35 propagates the system forward using H(t) and solving a linear equation  kuismam committed Jun 25, 2015 36 37 38  .. math::  kuismam committed Jun 26, 2015 39  \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 40   jensj committed Jun 26, 2015 41 Using the predicted wave functions at C'(t+dt), the Hamiltonian H'(t+dt) is  kuismam committed Jun 25, 2015 42 calculated and the Hamiltonian at middle of the time step is estimated as  kuismam committed Jun 25, 2015 43 44 45  .. math::  kuismam committed Jun 26, 2015 46  \mathbf{H}(t+{\rm d}t/2) = (\mathbf{H}(t) + \mathbf{H}'(t+{\rm d}t)) / 2  kuismam committed Jun 25, 2015 47 48  With the improved Hamiltonian, have functions are again propagated from t to t+dt  kuismam committed Jun 09, 2015 49   kuismam committed Jun 25, 2015 50 51 .. math::  kuismam committed Jun 26, 2015 52  \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 53   jensj committed Jun 26, 2015 54 This procedure is repeated using time step of 5-40as and for 500-2000 times to  kuismam committed Jun 25, 2015 55 obtain time evolution of electrons.  tprossi committed Jun 11, 2015 56   kuismam committed Jun 25, 2015 57 58 59 ===== Usage =====  kuismam committed Jun 09, 2015 60 61 62 63  Create LCAOTDDFT object like a GPAW calculator:: >>> from gpaw.lcaotddft import LCAOTDDFT  askhl committed Jun 29, 2015 64  >>> td_calc = LCAOTDDFT(setups={'Na':'1'}, basis='dzp', xc='LDA', h=0.3, nbands=1,  kuismam committed Jun 09, 2015 65 66 67 68 69 70 71 72 73 74  convergence={'density':1e-7}, poissonsolver=PoissonSolver(eps=1e-20, remove_moment=1+3+5)) Important points are: * Use always compatible setups and basis sets. * Grid spacing is only used to calculate the Hamiltonian matrix and therefore larger grid than usual can be used. * Completely unoccupied bands should be left out of the calculation, since they are not needed. * Convergence of density should be few orders of magnitude more accurate than in ground state calculations. * Convergence of poisson solver should be at least 1e-14, but 1e-20 does not hurt (this is the quadratic error).  kuismam committed Jun 25, 2015 75  * One should use multipole corrected poisson solver in any TDDFT run. For more advanced poisson solvers, see :ref:advancedpoisson.  kuismam committed Jun 09, 2015 76 77 78 79 80 81 82 83 84 85 86 87  Perform a regular ground state calculation, the get the ground state wave functions:: >>> atoms.set_calculator(td_calc) >>> atoms.get_potential_energy() If you wish to save here, write the wave functions also:: >>> td_calc.write('Na2.gpw', mode='all') The calculation proceeds as in grid mode. We kick the system to x-direction and propagate with 10as time steps for 500 steps::  kuismam committed Jun 25, 2015 88  >>> td_calc.absorption_kick([1e-5, 0.0, 0.0])  kuismam committed Jun 09, 2015 89 90 91 92 93 94 95 96 97 98 99 100 101  >>> td_calc.propagate(10, 500, out='Na2.dm') The spectrum is obtained in same manner, as in grid propagation mode. Simple run script ================= .. literalinclude:: lcaotddft.py General notes about basis sets ==============================  jensj committed Jun 26, 2015 102 103 104 105 In time-propagation LCAO-TDDFT, the basis sets are in even more crucial role than in a ground state LCAO calculation. It is required, that basis set can represent both the occupied (electrons) and relevant unoccupied states (holes) adequately. Custom basis sets for the time propagation should be generated according to ones  kuismam committed Jun 25, 2015 106 need, and then benchmarked.  kuismam committed Jun 09, 2015 107   jensj committed Jun 26, 2015 108 109 110 111 **Irrespective of the basis sets you choose, ALWAYS, ALWAYS, benchmark LCAO results with respect to grid time-propagation code** on a largest system possible. For example, one can create a prototype system, which consists of similar atom species with similar roles than in the parent system, but small  kuismam committed Jun 25, 2015 112 enough to calculate with grid propagation mode. For example, in relevance of benchmarking the basis set, see Fig. 4 and 5 of [#Kuisma2015]_.  tprossi committed Jun 18, 2015 113   jensj committed Jun 26, 2015 114 115 After these remarks, we describe two sets of basis sets that can be used as a starting point for choosing suitable basis set for your needs. Namely,  kuismam committed Jun 25, 2015 116 :ref:pvalence basis sets and :ref:coopt basis sets.  tprossi committed Jun 18, 2015 117 118 119 120 121 122  .. _pvalence basis sets: p-valence basis sets --------------------  jensj committed Jun 26, 2015 123 124 125 126 The so-called p-valence basis sets are constructed by replace the p-type polarization function of the default basis sets with 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  kuismam committed Jun 25, 2015 127 states of unoccupied states.  tprossi committed Jun 18, 2015 128   tprossi committed Jun 26, 2015 129 130 131 132 The p-valence basis sets can be easily obtained for appropriate elements with gpaw-install-setups tool with following options: gpaw-install-setups --basis --version=pvalence. See :ref:installationguide_setup_files for basis set installation. It is again reminded that these basis sets  jensj committed Jun 26, 2015 133 are not thoroughly tested and **it is essential to benchmark the performance of  kuismam committed Jun 25, 2015 134 the basis sets for your application**.  tprossi committed Jun 18, 2015 135 136 137 138 139 140  .. _coopt basis sets: Completeness-optimized basis sets ---------------------------------  jensj committed Jun 26, 2015 141 142 143 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 144 silver, and gold clusters.  tprossi committed Jun 18, 2015 145   jensj committed Jun 26, 2015 146 147 For further details of the basis sets, as well as their construction and performance, see [#Rossi2015]_. For convenience, these basis sets can be easily  tprossi committed Jun 26, 2015 148 149 150 151 obtained with gpaw-install-setups tool with following options: gpaw-install-setups --basis --version=coopt. See :ref:installationguide_setup_files for basis set installation. Finally, it is again  jensj committed Jun 26, 2015 152 emphasized that when using the basis sets, **it is essential to benchmark their  kuismam committed Jun 25, 2015 153 suitability for your application**.  tprossi committed Jun 18, 2015 154   kuismam committed Jun 09, 2015 155 156 157 158  Parallelization ===============  kuismam committed Jun 25, 2015 159 160 161 162 163 164 LCAO-TDDFT is parallelized using ScaLAPACK. It runs without scalapack, but in this case only single core is used for linear alrebra. * ScaLAPACK can be enabled by specifying --sl_default=N,M,64 in command line. * Alternatively, use parallization={'sl_default':(N,M,64)}. * It is necessary that N*M equals to total number of cores and max(N,M)*64 < nbands. 64 can be changed to 16 if necessary. * Apart from parallelization of linear algrebra, normal domain and band parallelizations can be used. As in LCAO-mode usually, use band parallelization to reduce memory consumption.  kuismam committed Jun 09, 2015 165   tprossi committed Jun 18, 2015 166 167 Timing ======  tprossi committed Jun 11, 2015 168   tprossi committed Jun 18, 2015 169 170 TODO: add ParallelTimer example  tprossi committed Jun 11, 2015 171   tprossi committed Jun 18, 2015 172 Advanced tutorial - Plasmon resonance of silver cluster  kuismam committed Jun 09, 2015 173 174 =======================================================  jensj committed Jun 26, 2015 175 176 177 178 One should think what type of transitions is there of interest, and make sure that the basis set can represent such Kohn-Sham electron and hole wave functions. The first transitions in silver cluster will be 5s \rightarrow 5p like. We require 5p orbitals in the basis set, and thus, we must generate a custom basis  kuismam committed Jun 25, 2015 179 180 set.  jensj committed Jun 26, 2015 181 Here is how to generate a double-zeta basis set with 5p orbital in valence for  kuismam committed Jun 25, 2015 182 Silver for GLLB-SC potential. We will use GPAW 0.8 setup definition, since semi-core p states are not relevant here.  kuismam committed Jun 25, 2015 183 184  .. literalinclude:: lcaotddft_basis.py  kuismam committed Jun 09, 2015 185   kuismam committed Jun 25, 2015 186 187 188 189 190 We calculate the icosahedral Ag55 cluster: :download:ag55.xyz This code uses ScaLAPACK parallelization with 64 cores. .. literalinclude:: lcaotddft_ag55.py  kuismam committed Jun 09, 2015 191   kuismam committed Jun 25, 2015 192 193 194 195 196 Code runs for approximately two wall hours. The resulting spectrum shows already emerging plasmonic excitation around 4 eV. For more details, see [#Kuisma2015]_. .. image:: fig1.png  kuismam committed Jun 25, 2015 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 Induced density =============== Plotting the induced density is especially interesting in case of plasmon resonances. As an example, we calculate a dummy Na8 wire and write density to a file on every iteration. There is certain advantages in writing the density on every iteration instead of using the predefined frequencies and on-fly Fourier transformation: only one TDDFT run is required as any frequency can be analysed as a post processing operation. Hard disk requirements are large, but tolerable (1-100GB) in usual cases. .. literalinclude:: lcaotddft_induced.py Files with extensions .sG and .asp are created, where .sG contains the density in coarse grid and .asp contains the atomic density matrix. With these, it is possible to reconstruct the full density. This can now be fourier transformed at desired frequency. Here, we look from the produced spectrum file that plasmonic peak, and perform Fourier transform at that frequency. .. literalinclude:: lcaotddft_analyse.py Two cube files are created, one for sin (imag) and cos (real) transform at the frequency. Usually, one is interested in the absorbing part i.e. imaginary part. Below is visualized the plasmon resonance  kuismam committed Jun 25, 2015 212 in Na8 wire. In the current form, these cube files contain just the pseudo part of density.  kuismam committed Jun 25, 2015 213 214 215 216  .. image:: Na8_imag.png  kuismam committed Jun 09, 2015 217 218 219 220 Advanced tutorial - large organic molecule ========================================== General notes  tprossi committed Jun 18, 2015 221 -------------  kuismam committed Jun 09, 2015 222   jensj committed Jun 26, 2015 223 224 225 226 227 On large organic molecules, on large conjugated systems, there will \pi \rightarrow \pi^*, \sigma \rightarrow \sigma^*. These states consists of only the valence orbitals of carbon, and they are likely by quite similar few eV's below and above the fermi lavel. These is thus a reason to believe that these states are well described with hydrogen 1s and carbon 2s and 2p valence orbitals  kuismam committed Jun 25, 2015 228 around the fermi level.  kuismam committed Jun 09, 2015 229 230 231  Here, we will calculate a small and a large organic molecule with lcao-tddft.  kuismam committed Jun 25, 2015 232 **TODO**  kuismam committed Jun 09, 2015 233 234 235 236  Kohn-Sham decomposition of the transition density matrix ========================================================  kuismam committed Jun 25, 2015 237 238 Soon it will be possible to analyse the origin of the transitions the same way as is commonly done in Casida-based codes. The LCAO basis will be transformed to electron-hole basis of the Kohn-Sham system.  tprossi committed Jun 11, 2015 239   jensj committed Jun 26, 2015 240   tprossi committed Jun 11, 2015 241 242 243 244 References ========== .. [#Kuisma2015]  jensj committed Jun 26, 2015 245  M. Kuisma, A. Sakko, T. P. Rossi, A. H. Larsen, J. Enkovaara, L. Lehtovaara, and T. T. Rantala,  tprossi committed Jun 11, 2015 246 247  Localized surface plasmon resonance in silver nanoparticles: Atomistic first-principles time-dependent density functional theory calculations,  tprossi committed Jun 18, 2015 248 249 250 251 252 253 254 255  *Phys. Rev. B* **69**, 245419 (2004). doi:10.1103/PhysRevB.91.115431 _ .. [#Rossi2015] T. P. Rossi, S. Lehtola, A. Sakko, M. J. Puska, and R. M. Nieminen, Nanoplasmonics simulations at the basis set limit through completeness-optimized, local numerical basis sets, *J. Chem. Phys.* **142**, 094114 (2015). doi:10.1063/1.4913739 _