1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
|
.. _lrtddft:
=====================
Linear response TDDFT
=====================
Ground state
============
The linear response TDDFT calculation needs a converged ground state
calculation with a set of unoccupied states.
Note, that the eigensolver 'rmm-diis' should not be used
for the calculation of unoccupied states, better use 'dav' or 'cg':
.. literalinclude:: Be_gs_8bands.py
Calculating the Omega Matrix
============================
The next step is to calculate the Omega Matrix from the ground state orbitals:
.. literalinclude:: Be_8bands_lrtddft.py
alternatively one can also restrict the number of transitions by their energy:
.. literalinclude:: Be_8bands_lrtddft_dE.py
Note, that parallelization over spin does not work here. As a workaround,
domain decomposition only (``parallel={'domain': world.size}``,
see :ref:`manual_parsize_domain`)
has to be used for spin polarised
calculations in parallel.
Extracting the spectrum
=======================
The dipole spectrum can be evaluated from the Omega matrix and written to
a file:
.. literalinclude:: Be_spectrum.py
The spectrum may be also extracted and plotted in energy or
wavelength directly::
import pylab as plt
from gpaw.lrtddft import LrTDDFT
from gpaw.lrtddft.spectrum import get_folded_spectrum
lr = LrTDDFT('lr.dat.gz')
plt.subplot(121)
# spectrum in energy
x, y = get_folded_spectrum(lr, width=0.05)
plt.plot(x, y[:, 0])
plt.xlabel('energy [eV]')
plt.ylabel('folded oscillator strength')
plt.subplot(122)
# spectrum in wavelengths
x, y = get_folded_spectrum(lr, energyunit='nm', width=10)
plt.plot(x, y[:, 0])
plt.xlabel('wavelength [nm]')
plt.show()
Testing convergence
===================
You can test the convergence of the Kohn-Sham transition basis size by restricting
the basis in the diagonalisation step, e.g.::
from gpaw.lrtddft import LrTDDFT
lr = LrTDDFT.read('lr.dat.gz')
lr.diagonalize(restrict={'energy_range':2.})
This can be automated by using the check_convergence function::
from gpaw.lrtddft import LrTDDFT
from gpaw.lrtddft.convergence import check_convergence
lr = LrTDDFT.read('lr.dat.gz')
check_convergence(lr,
'linear_response',
'my plot title',
dE=.2,
emax=6.)
which will create a directory 'linear_response'. In this directory there will be a
file 'conv.gpl' for gnuplot that compares the spectra varying the basis size.
Analysing the transitions
=========================
The single transitions (or a list of transitions) can be analysed as follows
(output printed)::
from gpaw.lrtddft import LrTDDFT
lr = LrTDDFT.read('lr.dat.gz')
lr.diagonalize()
# analyse transition 1
lr.analyse(1)
# analyse transition 0-10
lr.analyse(range(11))
Relaxation in the excited state
===============================
Despite that we do not have analytical gradients in linear response TDDFT,
we may use finite differences for relaxation in the excited state.
This example shows how to relax the sodium dimer
in the B excited state:
.. literalinclude:: Na2_relax_excited.py
The example runs on a single core. If started on 8 cores, it will split
``world`` into 4 independent parts (2 cores each) that can calculate
4 displacements in parallel at the same time.
Quick reference
===============
Parameters for LrTDDFT:
================ ============== =================== ========================================
keyword type default value description
================ ============== =================== ========================================
``calculator`` ``GPAW`` Calculator object of ground state
calculation
``nspins`` ``int`` 1 number of excited state spins, i.e.
singlet-triplet transitions are
calculated with ``nspins=2``. Effective
only if ground state is spin-compensated
``xc`` ``string`` xc of calculator Exchange-correlation for LrTDDFT, can
differ from ground state value
``restrict`` ``dict`` {} Restrictions ``eps``, ``istart``, ``jend``
and ``energy_range`` collected as dict.
``eps`` ``float`` 0.001 Minimal occupation difference for a transition
``istart`` ``int`` 0 first occupied state to consider
``jend`` ``int`` number of bands last unoccupied state to consider
``energy_range`` ``float`` None Energy range to consider in the involved
Kohn-Sham orbitals (replaces [istart,jend])
================ ============== =================== ========================================
|
