lcaotddft.rst 13.7 KB
Newer Older
1 2
.. _lcaotddft:

3
===================================================
4
Time-propagation TDDFT with LCAO: Theory and usage
5
===================================================
6

7
This page documents the use of time-propagation TDDFT in :ref:`LCAO
8
mode <lcao>`. The implementation is described in Ref. [#Kuisma2015]_.
9

10
Real-time propagation of LCAO functions
kuismam's avatar
kuismam committed
11 12
=======================================

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's avatar
kuismam committed
16 17 18

.. math::

19
  \tilde{\psi}_n(\mathbf{r},t) = \sum_{\mu} \tilde{\phi}_{\mu}(\mathbf{r}-\mathbf{R}^\mu) c_{\mu n}(t) .
kuismam's avatar
kuismam committed
20

21
The time-dependent Kohn--Sham equation in the PAW formalism can be written as
kuismam's avatar
kuismam committed
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.

27 28
From this, the following matrix equation can be derived for the
LCAO wave function coefficients:
kuismam's avatar
kuismam committed
29 30 31 32

.. math::
  {\rm i}\mathbf{S} \frac{{\rm d}\mathbf{C}(t)}{{\rm d}t} = \mathbf{H}(t) \mathbf{C}(t).

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's avatar
kuismam committed
39 40 41

.. math::

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's avatar
kuismam committed
43

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's avatar
kuismam committed
47 48 49

.. math::

50
   \mathbf{H}(t+{\rm d}t/2) = (\mathbf{H}(t) + \mathbf{H}'(t+{\rm d}t)) / 2
kuismam's avatar
kuismam committed
51

52 53
With the improved Hamiltonian, the wave functions are again propagated
from `t` to `t+\mathrm dt` by solving
kuismam's avatar
kuismam committed
54

kuismam's avatar
kuismam committed
55 56
.. math::

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

59 60
This procedure is repeated using 500--2000 time steps of 5--40 as to
calculate the time evolution of the electrons.
61

62 63 64 65
.. _example:

Example usage
=============
kuismam's avatar
kuismam committed
66

67
First do a standard ground-state calculation with the ``GPAW`` calculator:
kuismam's avatar
kuismam committed
68

69 70
.. literalinclude:: lcaotddft.py
   :lines: 3-21
kuismam's avatar
kuismam committed
71

72
Some important points are:
kuismam's avatar
kuismam committed
73

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's avatar
kuismam committed
88

89 90
.. literalinclude:: lcaotddft.py
   :lines: 23-35
kuismam's avatar
kuismam committed
91

92
After the time propagation, the spectrum can be calculated:
kuismam's avatar
kuismam committed
93

94 95
.. literalinclude:: lcaotddft.py
   :lines: 37-39
kuismam's avatar
kuismam committed
96

97 98
The previous example as a complete script can be downloaded here:
:download:`lcaotddft.py`.
kuismam's avatar
kuismam committed
99 100 101 102 103


General notes about basis sets
==============================

104
In time-propagation LCAO-TDDFT, it is much more important to think
105
about the basis sets compared to ground-state LCAO calculations.  It
106
is required that the basis set can represent both the occupied
107
(holes) and relevant unoccupied states (electrons) adequately.  Custom
108 109
basis sets for the time propagation should be generated according to
one's need, and then benchmarked.
kuismam's avatar
kuismam committed
110

jensj's avatar
jensj committed
111
**Irrespective of the basis sets you choose, ALWAYS, ALWAYS, benchmark LCAO
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.

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`.
121

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
122

123 124 125 126 127
.. _pvalence basis sets:

p-valence basis sets
--------------------

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
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 {<directory>} --basis --version=pvalence
138

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
139
See :ref:`installation of paw datasets` for more information on basis set
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
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**.
143 144 145 146 147 148 149


.. _coopt basis sets:

Completeness-optimized basis sets
---------------------------------

jensj's avatar
jensj committed
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,
153
silver, and gold clusters.
154

jensj's avatar
jensj committed
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's avatar
Jens Jørgen Mortensen committed
157
obtained with::
158

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
159
    $ gpaw install-data {<directory>} --basis --version=coopt
160

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
161
See :ref:`installation of paw datasets` for basis set installation. Finally,
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
162 163
it is again emphasized that when using the basis sets, **it is essential to
benchmark their suitability for your application**.
164

kuismam's avatar
kuismam committed
165 166 167 168

Parallelization
===============

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's avatar
kuismam committed
171

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's avatar
kuismam committed
179

180

181 182 183 184 185 186
.. TODO

    Timing
    ======
    
    Add ``ParallelTimer`` example
tprossi's avatar
tprossi committed
187

188

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
==============  ===============================
237 238


239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287
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.

First, let's calculate the spectrum:

.. literalinclude:: lcaotddft_Na8/spectrum.py

We notice the main resonances at TODO

Frequency-space density matrix
""""""""""""""""""""""""""""""

We generate the density matrix for the frequencies of interest:

.. literalinclude:: lcaotddft_Na8/td_fdm_replay.py

.. tip::

   Instead of using replaying, you can do the same analysis on-the-fly by
   attaching the analysis tools to the usual time-propagation calculation.

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.

288 289 290 291 292
We obtain the following transition contributions:

.. literalinclude:: lcaotddft_Na8/table_1.12.txt
   :language: none

293 294 295
.. image:: lcaotddft_Na8/tcm_1.12.png
   :scale: 70%

296 297 298
.. literalinclude:: lcaotddft_Na8/table_2.48.txt
   :language: none

299 300 301
.. image:: lcaotddft_Na8/tcm_2.48.png
   :scale: 70%

302

303
Induced density
304
---------------
305

306 307
The density matrix gives access to any other quantities.
For instance, the induced density can be conveniently obtained:
308

309
.. literalinclude:: lcaotddft_Na8/ksd_ind.py
310

311 312 313
The resulting cube files can be visualized, for example, with
:download:`this script <lcaotddft_Na8/ind_plot.py>` producing
the figures:
314

315 316
.. image:: lcaotddft_Na8/ind_1.12.png
   :scale: 70%
317

318 319
.. image:: lcaotddft_Na8/ind_2.48.png
   :scale: 70%
320 321


322 323 324 325 326 327 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
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.
   

371 372 373 374
References
==========

.. [#Kuisma2015]
375 376 377 378 379
   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,
380
   *Phys. Rev. B* **69**, 245419 (2004).
381
   `doi:10.1103/PhysRevB.91.115431 <https://doi.org/10.1103/PhysRevB.91.115431>`_
382 383 384

.. [#Rossi2015]
   T. P. Rossi, S. Lehtola, A. Sakko, M. J. Puska, and R. M. Nieminen,
385 386
   Nanoplasmonics simulations at the basis set limit
   through completeness-optimized, local numerical basis sets,
387
   *J. Chem. Phys.* **142**, 094114 (2015).
388 389 390 391 392 393 394 395 396
   `doi:10.1063/1.4913739 <https://doi.org/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 <https://doi.org/10.1021/acs.jctc.7b00589>`_