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

3 4 5
================================
Time-propagation TDDFT with LCAO
================================
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 66
.. _example:

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

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

70
.. literalinclude:: lcaotddft.py
71
   :lines: 3-20
kuismam's avatar
kuismam committed
72

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

75 76 77 78 79 80 81 82
* 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).
Tuomas Rossi's avatar
Tuomas Rossi committed
83 84
* One should use multipole-corrected Poisson solvers or
  other advanced Poisson solvers in any TDDFT run
85 86
  in order to guarantee the convergence of the potential with respect to
  the vacuum size.
87 88
  See the documentation on :ref:`advancedpoisson`.

89
Next we kick the system in the z direction and propagate 3000 steps of 10 as:
kuismam's avatar
kuismam committed
90

91
.. literalinclude:: lcaotddft.py
92
   :lines: 22-34
kuismam's avatar
kuismam committed
93

94
After the time propagation, the spectrum can be calculated:
kuismam's avatar
kuismam committed
95

96
.. literalinclude:: lcaotddft.py
97 98 99 100
   :lines: 36-38

This example input script can be downloaded :download:`here <lcaotddft.py>`.

kuismam's avatar
kuismam committed
101

102 103 104 105 106 107
Extending the functionality
---------------------------

The real-time propagation LCAOTDDFT provides very modular framework
for calculating many properties during the propagation.
See :ref:`analysis` for tutorial how to extend the analysis capabilities.
kuismam's avatar
kuismam committed
108 109


Tuomas Rossi's avatar
Tuomas Rossi committed
110 111
.. _note basis sets:

kuismam's avatar
kuismam committed
112 113 114
General notes about basis sets
==============================

115
In time-propagation LCAO-TDDFT, it is much more important to think
116
about the basis sets compared to ground-state LCAO calculations.  It
117
is required that the basis set can represent both the occupied
118
(holes) and relevant unoccupied states (electrons) adequately.  Custom
119 120
basis sets for the time propagation should be generated according to
one's need, and then benchmarked.
kuismam's avatar
kuismam committed
121

122
**Irrespective of the basis sets you choose always benchmark LCAO
123 124 125 126 127 128
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.

129 130 131
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`.
132

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

134 135 136 137 138
.. _pvalence basis sets:

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

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
139 140 141 142 143 144 145 146 147 148
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
149

Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
150
See :ref:`installation of paw datasets` for more information on basis set
Jens Jørgen Mortensen's avatar
Jens Jørgen Mortensen committed
151 152 153
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**.
154 155 156 157 158 159 160


.. _coopt basis sets:

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

jensj's avatar
jensj committed
161 162 163
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,
164
silver, and gold clusters.
165

jensj's avatar
jensj committed
166 167
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
168
obtained with::
169

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

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

kuismam's avatar
kuismam committed
176

177 178
.. _parallelization:

kuismam's avatar
kuismam committed
179 180 181
Parallelization
===============

182 183 184 185 186
For maximum performance on large systems, it is advicable to use
ScaLAPACK and large band parallelization with ``augment_grids`` enabled.
This can be achieved with parallelization settings like
``parallel={'sl_auto': True, 'domain': 8, 'augment_grids': True}``
(see :ref:`manual_parallel`),
Tuomas Rossi's avatar
Tuomas Rossi committed
187
which will use groups of 8 tasks for domain parallelization and the rest for
188 189 190 191 192 193 194 195 196 197 198 199
band parallelization (for example, with a total of 144 cores this would mean
domain and band parallelizations of 8 and 18, respectively).

Instead of ``sl_auto``, the ScaLAPACK settings can be set by hand
as ``sl_default=(m, n, block)`` (see :ref:`manual_ScaLAPACK`,
in which case it is important that ``m * n``` equals
the total number of cores used by the calculator
and that ``max(m, n) * block < nbands``.

It is also possible to run the code without ScaLAPACK, but it is
very inefficient for large systems as in that case only a single core
is used for linear algebra.
kuismam's avatar
kuismam committed
200

201

202 203 204 205
.. TODO

    Timing
    ======
206

207
    Add ``ParallelTimer`` example
tprossi's avatar
tprossi committed
208

209

210 211 212 213
.. _analysis:

Modular analysis tools
======================
214 215 216

In :ref:`example` it was demonstrated how to calculate photoabsorption
spectrum from the time-dependent dipole moment data collected with
217
``DipoleMomentWriter`` observer.
Tuomas Rossi's avatar
Tuomas Rossi committed
218
However, any (also user-written) analysis tools can be attached
219 220 221 222
as a separate observers in the general time-propagation framework.

There are two ways to perform analysis:
   1. Perform analysis on-the-fly during the propagation::
223 224

         # Read starting point
225
         td_calc = LCAOTDDFT('gs.gpw')
226

227 228
         # Attach analysis tools
         MyObserver(td_calc, ...)
229

230 231 232
         # Kick and propagate
         td_calc.absorption_kick([1e-5, 0., 0.])
         td_calc.propagate(10, 3000)
233

234 235
      For example, the analysis tools can be ``DipoleMomentWriter`` observer
      for spectrum or Fourier transform of density at specific frequencies etc.
236

237 238
   2. Record the wave functions during the first propagation and
      perform the analysis later by replaying the propagation::
239 240

         # Read starting point
241
         td_calc = LCAOTDDFT('gs.gpw')
242

243 244
         # Attach analysis tools
         MyObserver(td_calc, ...)
245

246
         # Replay propagation from a stored file
247 248
         td_calc.replay(name='wf.ulm', update='all')

249 250 251 252 253 254 255 256 257 258 259 260 261 262
      From the perspective of the attached observers the replaying
      is identical to actual propagation.

   The latter method is recommended, because one might not know beforehand
   what to analyze.
   For example, the interesting resonance frequencies are often not know before
   the time-propagation calculation.

In the following we give an example how to utilize the replaying capability
in practice and describe some analysis tools available in GPAW.


Example
-------
263 264 265 266 267 268 269 270 271 272 273 274 275 276

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
277
with ``WaveFunctionWriter`` observer:
278 279 280 281 282

.. literalinclude:: lcaotddft_Na8/td.py

.. tip::

Tuomas Rossi's avatar
Tuomas Rossi committed
283 284
   The time propagation can be continued in the same manner
   from the restart file:
285 286 287

   .. literalinclude:: lcaotddft_Na8/tdc.py

288
The created ``wf.ulm`` file contains the time-dependent wave functions
289 290
`C_{\mu n}(t)` that define the state of the system at each time.
We can use the file to replay the time propagation:
291 292 293

.. literalinclude:: lcaotddft_Na8/td_replay.py

294
The ``update`` keyword in ``replay()`` has following options:
295 296 297 298 299 300 301 302

==============  ===============================
``update``      variables updated during replay
==============  ===============================
``'all'``       Hamiltonian and density
``'density'``   density
``'none'``      nothing
==============  ===============================
303

304 305 306 307 308 309 310 311 312 313 314 315
.. tip::

   The wave functions can be written in separate files
   by using ``split=True``::

      WaveFunctionWriter(td_calc, 'wf.ulm', split=True)

   This creates additional ``wf*.ulm`` files containing the wave functions.
   The replay functionality works as in the above example
   even with splitted files.


316

Tuomas Rossi's avatar
Tuomas Rossi committed
317 318
Kohn--Sham decomposition of density matrix
------------------------------------------
319 320 321

Kohn--Sham decomposition is an illustrative way of analyzing electronic
excitations in Kohn--Sham electron-hole basis.
322 323
See Ref. [#Rossi2017]_ for the description of the GPAW implementation
and demonstration.
324 325 326 327

Here we demonstrate how to construct the photoabsorption decomposition
at a specific frequency in Kohn--Sham electon-hole basis.

Tuomas Rossi's avatar
Tuomas Rossi committed
328 329
First, let's calculate and :download:`plot <lcaotddft_Na8/spec_plot.py>`
the spectrum:
330 331 332

.. literalinclude:: lcaotddft_Na8/spectrum.py

Tuomas Rossi's avatar
Tuomas Rossi committed
333 334 335 336
.. image:: lcaotddft_Na8/spec.png
   :scale: 70%

The two main resonances are analyzed in the following.
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

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

We generate the density matrix for the frequencies of interest:

.. literalinclude:: lcaotddft_Na8/td_fdm_replay.py

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.

365 366
We obtain the following transition contributions for the resonances
(presented both as tables and TCMs):
367 368 369 370

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

371 372 373
.. image:: lcaotddft_Na8/tcm_1.12.png
   :scale: 70%

374 375 376
.. literalinclude:: lcaotddft_Na8/table_2.48.txt
   :language: none

377 378 379
.. image:: lcaotddft_Na8/tcm_2.48.png
   :scale: 70%

380

381
Induced density
382
---------------
383

384
The density matrix gives access to any other quantities.
385 386
For instance, the induced density can be conveniently obtained
from the density matrix:
387

388
.. literalinclude:: lcaotddft_Na8/fdm_ind.py
389

390 391 392
The resulting cube files can be visualized, for example, with
:download:`this script <lcaotddft_Na8/ind_plot.py>` producing
the figures:
393

394 395
.. image:: lcaotddft_Na8/ind_1.12.png
   :scale: 70%
396

397 398
.. image:: lcaotddft_Na8/ind_2.48.png
   :scale: 70%
399 400


401 402 403 404 405 406
Advanced tutorials
==================

Plasmon resonance of silver cluster
-----------------------------------

Tuomas Rossi's avatar
Tuomas Rossi committed
407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437
This tutorial demonstrates
:ref:`the efficient parallelization options <parallelization>` and
the importance of :ref:`proper basis sets <note basis sets>`.

We calculate the photoabsorption spectrum of
:download:`an icosahedral Ag55 cluster <lcaotddft_Ag55/Ag55.xyz>`.
We use GLLB-SC potential to significantly improve the description of d states,
:ref:`pvalence basis sets` to improve the description of unoccupied states, and
11-electron Ag setup to reduce computational cost.

The workflow is the same as in the previous examples.
First, we calculate ground state (takes around 20 minutes with 36 cores):

.. literalinclude:: lcaotddft_Ag55/gs.py

Then, we carry the time propagation for 30 femtoseconds in steps of
10 attoseconds (takes around 3.5 hours with 36 cores):

.. literalinclude:: lcaotddft_Ag55/td.py

Finally, we calculate the spectrum (takes a few seconds):

.. literalinclude:: lcaotddft_Ag55/spec.py

The resulting spectrum shows an emerging plasmon resonance at around 4 eV:

.. image:: lcaotddft_Ag55/Ag55_spec.png

For further insight on plasmon resonance in metal nanoparticles,
see [#Kuisma2015]_ and [#Rossi2017]_.

438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456
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.



.. TODO

   Large organic molecule
   ----------------------
457

458 459 460 461 462 463
   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.
464

465
   Here, we will calculate a small and a large organic molecule with lcao-tddft.
466

467

468 469 470 471
References
==========

.. [#Kuisma2015]
472 473 474 475 476
   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,
477
   *Phys. Rev. B* **69**, 245419 (2004).
478
   `doi:10.1103/PhysRevB.91.115431 <https://doi.org/10.1103/PhysRevB.91.115431>`_
479 480 481

.. [#Rossi2015]
   T. P. Rossi, S. Lehtola, A. Sakko, M. J. Puska, and R. M. Nieminen,
482 483
   Nanoplasmonics simulations at the basis set limit
   through completeness-optimized, local numerical basis sets,
484
   *J. Chem. Phys.* **142**, 094114 (2015).
485 486 487 488 489 490 491 492 493
   `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>`_