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

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

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

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

jensj's avatar
jensj committed
13
In real time LCAO-TDDFT approach, the time-dependent wave functions are
14
represented using localized basis sets as
kuismam's avatar
kuismam committed
15 16 17

.. math::

18
  \tilde{\Psi}_i(\mathbf{r},t) = \sum_{\mu} C_{\mu i}(t) \tilde{\phi}^{\mu}(\mathbf{r}-\mathbf{R}^\mu).
kuismam's avatar
kuismam committed
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's avatar
jensj committed
26
Using these equations, following matrix equation can be derived for LCAO wave
27
function coefficients
kuismam's avatar
kuismam committed
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's avatar
jensj committed
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
35
propagates the system forward using H(t) and solving a linear equation
kuismam's avatar
kuismam committed
36 37 38

.. math::

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

jensj's avatar
jensj committed
41
Using the predicted wave functions at C'(t+dt), the Hamiltonian H'(t+dt) is
42
calculated and the Hamiltonian at middle of the time step is estimated as
kuismam's avatar
kuismam committed
43 44 45

.. math::

46
   \mathbf{H}(t+{\rm d}t/2) = (\mathbf{H}(t) + \mathbf{H}'(t+{\rm d}t)) / 2
kuismam's avatar
kuismam committed
47 48

With the improved Hamiltonian, have functions are again propagated from t to t+dt
kuismam's avatar
kuismam committed
49

kuismam's avatar
kuismam committed
50 51
.. math::

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)
53

jensj's avatar
jensj committed
54
This procedure is repeated using time step of 5-40as and for 500-2000 times to
55
obtain time evolution of electrons.
56

57 58 59
=====
Usage
=====
kuismam's avatar
kuismam committed
60 61 62 63

Create LCAOTDDFT object like a GPAW calculator::

 >>> from gpaw.lcaotddft import LCAOTDDFT
64
 >>> td_calc = LCAOTDDFT(setups={'Na':'1'}, basis='dzp', xc='LDA', h=0.3, nbands=1,
kuismam's avatar
kuismam committed
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).
75
 * One should use multipole corrected poisson solver in any TDDFT run. For more advanced poisson solvers, see :ref:`advancedpoisson`.
kuismam's avatar
kuismam committed
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's avatar
kuismam committed
88
 >>> td_calc.absorption_kick([1e-5, 0.0, 0.0])
kuismam's avatar
kuismam committed
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's avatar
jensj committed
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
106
need, and then benchmarked.
kuismam's avatar
kuismam committed
107

jensj's avatar
jensj committed
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
112
enough to calculate with grid propagation mode. For example, in relevance of benchmarking the basis set, see Fig. 4 and 5 of [#Kuisma2015]_.
113

jensj's avatar
jensj committed
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,
116
:ref:`pvalence basis sets` and :ref:`coopt basis sets`.
117 118 119 120 121 122

.. _pvalence basis sets:

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

jensj's avatar
jensj committed
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
127
states of unoccupied states.
128

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's avatar
jensj committed
133
are not thoroughly tested and **it is essential to benchmark the performance of
134
the basis sets for your application**.
135 136 137 138 139 140

.. _coopt basis sets:

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

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

jensj's avatar
jensj committed
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
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's avatar
jensj committed
152
emphasized that when using the basis sets, **it is essential to benchmark their
153
suitability for your application**.
154

kuismam's avatar
kuismam committed
155 156 157 158

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

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

166 167
Timing
======
168

tprossi's avatar
tprossi committed
169 170
TODO: add ``ParallelTimer`` example

171

172
Advanced tutorial - Plasmon resonance of silver cluster
kuismam's avatar
kuismam committed
173 174
=======================================================

jensj's avatar
jensj committed
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
179 180
set.

jensj's avatar
jensj committed
181
Here is how to generate a double-zeta basis set with 5p orbital in valence for
kuismam's avatar
kuismam committed
182
Silver for GLLB-SC potential. We will use GPAW 0.8 setup definition, since semi-core p states are not relevant here.
183 184

.. literalinclude:: lcaotddft_basis.py
kuismam's avatar
kuismam committed
185

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

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

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's avatar
kuismam committed
212
in Na8 wire. In the current form, these cube files contain just the pseudo part of density.
213 214 215 216

.. image:: Na8_imag.png


kuismam's avatar
kuismam committed
217 218 219 220
Advanced tutorial - large organic molecule
==========================================

General notes
221
-------------
kuismam's avatar
kuismam committed
222

jensj's avatar
jensj committed
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
228
around the fermi level.
kuismam's avatar
kuismam committed
229 230 231

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

232
**TODO**
kuismam's avatar
kuismam committed
233 234 235 236

Kohn-Sham decomposition of the transition density matrix
========================================================

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.
239

jensj's avatar
jensj committed
240

241 242 243 244
References
==========

.. [#Kuisma2015]
jensj's avatar
jensj committed
245
   M. Kuisma, A. Sakko, T. P. Rossi, A. H. Larsen, J. Enkovaara, L. Lehtovaara, and T. T. Rantala,
246 247
   Localized surface plasmon resonance in silver nanoparticles: Atomistic first-principles time-dependent
   density functional theory calculations,
248 249 250 251 252 253 254 255
   *Phys. Rev. B* **69**, 245419 (2004).
   `doi:10.1103/PhysRevB.91.115431 <http://dx.doi.org/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 <http://dx.doi.org/10.1063/1.4913739>`_