Commit c8ace994 authored by Tuomas Rossi's avatar Tuomas Rossi

Add tutorial for Kohn-Sham decomposition

parent 51a32207
......@@ -59,8 +59,10 @@ from `t` to `t+\mathrm dt` by solving
This procedure is repeated using 500--2000 time steps of 5--40 as to
calculate the time evolution of the electrons.
Usage
=====
.. _example:
Example usage
=============
First do a standard ground-state calculation with the ``GPAW`` calculator:
......@@ -92,13 +94,8 @@ After the time propagation, the spectrum can be calculated:
.. literalinclude:: lcaotddft.py
:lines: 37-39
Example script
--------------
Here is the previous example as a complete script.
.. literalinclude:: lcaotddft.py
The previous example as a complete script can be downloaded here:
:download:`lcaotddft.py`.
General notes about basis sets
......@@ -238,61 +235,155 @@ For more details, see [#Kuisma2015]_.
Here, we will calculate a small and a large organic molecule with lcao-tddft.
Analysis tools and data recorders
=================================
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
============== ===============================
Kohn--Sham decomposition of the transition density matrix
---------------------------------------------------------
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 an electron-hole basis of the Kohn-Sham system.
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.
.. image:: lcaotddft_Na8/tcm_1.12.png
:scale: 70%
.. image:: lcaotddft_Na8/tcm_2.48.png
:scale: 70%
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
the density to a file on every iteration.
There are certain advantages of writing the density on every iteration
instead of using the predefined frequencies and
on-the-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 most cases.
.. literalinclude:: lcaotddft_induced.py
The density matrix gives access to any other quantities.
For instance, the induced density can be conveniently obtained:
Files with extensions ``.sG`` and ``.asp`` are created, where ``.sG`` files
contain the density on the coarse grid while ``.asp`` files contain
the atomic density matrices. With these, it is possible to reconstruct
the full density.
This can now be fourier transformed at the desired frequency.
Here, we look from the produced spectrum file that plasmonic peak,
and perform Fourier transform at that frequency.
.. literalinclude:: lcaotddft_Na8/ksd_ind.py
.. literalinclude:: lcaotddft_analyse.py
The resulting cube files can be visualized, for example, with
:download:`this script <lcaotddft_Na8/ind_plot.py>` producing
the figures:
Two cube files are created, one for the sin (imag) and cos (real)
transform at the frequency. Usually, one is interested in the absorbing
part, i.e., the imaginary part. Below the plasmon resonance is visualized
in the Na8 wire. In their current form, these cube files contain just
the pseudo part of density.
.. image:: lcaotddft_Na8/ind_1.12.png
:scale: 70%
.. image:: Na8_imag.png
.. image:: lcaotddft_Na8/ind_2.48.png
:scale: 70%
References
==========
.. [#Kuisma2015]
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,
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,
*Phys. Rev. B* **69**, 245419 (2004).
`doi:10.1103/PhysRevB.91.115431 <http://dx.doi.org/10.1103/PhysRevB.91.115431>`_
`doi:10.1103/PhysRevB.91.115431 <https://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,
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>`_
`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>`_
from ase import Atoms
from gpaw import GPAW
from gpaw.poisson import PoissonSolver
from gpaw.poisson_extravacuum import ExtraVacuumPoissonSolver
# Sodium atom chain
atoms = Atoms('Na8',
positions=[[i * 3.0, 0, 0] for i in range(8)],
cell=[33.6, 12.0, 12.0])
atoms.center()
# Use an advanced Poisson solver
eps = 1e-16
ps = ExtraVacuumPoissonSolver(gpts=(512, 256, 256),
poissonsolver_large=PoissonSolver(eps=eps),
coarses=2,
poissonsolver_small=PoissonSolver(eps=eps))
# Ground-state calculation
calc = GPAW(mode='lcao', h=0.3, basis='pvalence.dz', xc='LDA', nbands=6,
setups={'Na': '1'},
poissonsolver=ps,
convergence={'density': 1e-8},
txt='gs.out')
atoms.set_calculator(calc)
energy = atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from ase.io import read
def do(freq):
# Read cube file
cube = read('ind_%.2f.cube' % freq, full_output=True)
d_g = cube['data']
atoms = cube['atoms']
box = np.diag(atoms.get_cell())
ng = d_g.shape
# Take slice of data array
d_yx = d_g[:, :, ng[2] // 2]
x = np.linspace(0, box[0], ng[0])
xlabel = u'x (Å)'
y = np.linspace(0, box[1], ng[1])
ylabel = u'y (Å)'
# Plot
plt.figure(figsize=(8, 3.5))
ax = plt.subplot(1, 1, 1)
X, Y = np.meshgrid(x, y)
dmax = max(d_yx.min(), d_yx.max())
vmax = 0.9 * dmax
vmin = -vmax
plt.pcolormesh(X, Y, d_yx.T, cmap='RdBu_r', vmin=vmin, vmax=vmax)
contours = np.sort(np.outer([-1, 1], [0.02]).ravel() * dmax)
plt.contour(X, Y, d_yx.T, contours, cmap='RdBu_r', vmin=-1e-10, vmax=1e-10)
for atom in atoms:
pos = atom.position
plt.scatter(pos[0], pos[1], s=100, c='k', marker='o')
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.xlim([x[0], x[-1]])
plt.ylim([y[0], y[-1]])
ax.set_aspect('equal')
plt.title('Induced density of Na8 at %.2f eV' % freq)
plt.tight_layout()
plt.savefig('ind_%.2f.png' % freq)
do(1.12)
do(2.48)
import numpy as np
from ase.io import write
from gpaw import GPAW
from gpaw.tddft.units import au_to_eV
from gpaw.lcaotddft.ksdecomposition import KohnShamDecomposition
from gpaw.lcaotddft.densitymatrix import DensityMatrix
from gpaw.lcaotddft.frequencydensitymatrix import FrequencyDensityMatrix
# Load the objects
calc = GPAW('unocc.gpw', txt=None)
calc.initialize_positions() # Initialize in order to calculate density
ksd = KohnShamDecomposition(calc, 'ksd.ulm')
dmat = DensityMatrix(calc)
fdm = FrequencyDensityMatrix(calc, dmat, 'fdm.ulm')
def do(w):
# Select the frequency and the density matrix
rho_uMM = fdm.FReDrho_wuMM[w]
freq = fdm.freq_w[w]
print('Frequency: %.2f eV' % (freq.freq * au_to_eV))
print('Folding: %s' % freq.folding)
# Transform the LCAO density matrix to KS basis
rho_up = ksd.transform(rho_uMM)
# Induced density
rho_g = ksd.get_density([rho_up[0].imag])
# Save as a cube file
write('ind_%.2f.cube' % (freq.freq * au_to_eV), calc.atoms, data=rho_g)
# Calculate dipole moment for reference
dm_v = ksd.density.finegd.calculate_dipole_moment(rho_g, center=True)
absorption = 2 * freq.freq / np.pi * dm_v[0] / au_to_eV * 1e5
print('Total absorption: %.2f eV^-1' % absorption)
do(0)
do(1)
from gpaw import GPAW
from gpaw.lcaotddft.ksdecomposition import KohnShamDecomposition
# Calculate ground state with full unoccupied space
calc = GPAW('gs.gpw', nbands='nao', fixdensity=True,
txt='unocc.out')
atoms = calc.get_atoms()
energy = atoms.get_potential_energy()
calc.write('unocc.gpw', mode='all')
# Construct KS electron-hole basis
ksd = KohnShamDecomposition(calc)
ksd.initialize(calc)
ksd.write('ksd.ulm')
def agts(queue):
gs = queue.add('gs.py', ncpus=8, walltime=10)
td0 = queue.add('td.py', deps=[gs], ncpus=8, walltime=30)
td = queue.add('tdc.py', deps=[td0], ncpus=8, walltime=30)
queue.add('td_replay.py', deps=[td], ncpus=8, walltime=30)
spec = queue.add('spectrum.py', deps=[td], ncpus=1, walltime=2)
fdm = queue.add('td_fdm_replay.py', deps=[td], ncpus=1, walltime=5)
ksd = queue.add('ksd_init.py', deps=[gs], ncpus=1, walltime=5)
queue.add('tcm_plot.py', deps=[ksd, fdm, spec], ncpus=1, walltime=2,
creates=['tcm_1.12.png', 'tcm_2.48.png'])
ind = queue.add('ksd_ind.py', deps=[ksd, fdm], ncpus=1, walltime=2)
queue.add('ind_plot.py', deps=[ind], ncpus=1, walltime=2,
creates=['ind_1.12.png', 'ind_2.48.png'])
from gpaw.tddft.spectrum import photoabsorption_spectrum
photoabsorption_spectrum('dm.dat', 'spec.dat',
folding='Gauss', width=0.1,
e_min=0.0, e_max=10.0, delta_e=0.01)
import numpy as np
from matplotlib import pyplot as plt
from gpaw import GPAW
from gpaw.tddft.units import au_to_eV
from gpaw.lcaotddft.ksdecomposition import KohnShamDecomposition
from gpaw.lcaotddft.densitymatrix import DensityMatrix
from gpaw.lcaotddft.frequencydensitymatrix import FrequencyDensityMatrix
# Load the objects
calc = GPAW('unocc.gpw', txt=None)
ksd = KohnShamDecomposition(calc, 'ksd.ulm')
dmat = DensityMatrix(calc)
fdm = FrequencyDensityMatrix(calc, dmat, 'fdm.ulm')
def do(w):
# Select the frequency and the density matrix
rho_uMM = fdm.FReDrho_wuMM[w]
freq = fdm.freq_w[w]
print('Frequency: %.2f eV' % (freq.freq * au_to_eV))
print('Folding: %s' % freq.folding)
# Transform the LCAO density matrix to KS basis
rho_up = ksd.transform(rho_uMM)
# Photoabsorption decomposition
dmrho_vp = ksd.get_dipole_moment_contributions(rho_up)
weight_p = 2 * freq.freq / np.pi * dmrho_vp[0].imag / au_to_eV * 1e5
print('Total absorption: %.2f eV^-1' % np.sum(weight_p))
# Print contributions as a table
print(ksd.get_contributions_table(weight_p, minweight=0.1))
# Plot the decomposition as a TCM
r = ksd.plot_TCM(weight_p,
occ_energy_min=-3, occ_energy_max=0.1,
unocc_energy_min=-0.1, unocc_energy_max=3,
delta_energy=0.01, sigma=0.1
)
(ax_tcm, ax_occ_dos, ax_unocc_dos, ax_spec) = r
# Plot diagonal line at the analysis frequency
x = np.array([-4, 1])
y = x + freq.freq * au_to_eV
ax_tcm.plot(x, y, c='k')
# Show the plot
plt.savefig('tcm_%.2f.png' % (freq.freq * au_to_eV))
do(0)
do(1)
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
from gpaw.lcaotddft.wfwriter import WaveFunctionWriter
# Read the ground-state file
td_calc = LCAOTDDFT('gs.gpw', txt='td.out')
# Attach any data recording or analysis tools
DipoleMomentWriter(td_calc, 'dm.dat')
WaveFunctionWriter(td_calc, 'wfw.ulm')
# Kick and propagate
td_calc.absorption_kick([1e-5, 0., 0.])
td_calc.propagate(20, 1500)
# Save the state for restarting later
td_calc.write('td.gpw', mode='all')
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.densitymatrix import DensityMatrix
from gpaw.lcaotddft.frequencydensitymatrix import FrequencyDensityMatrix
from gpaw.tddft.folding import frequencies
# Read the ground-state file
td_calc = LCAOTDDFT('gs.gpw', propagator=dict(name='wfw.ulm', update='none'))
# Attach analysis tools
dmat = DensityMatrix(td_calc)
freqs = frequencies([1.12, 2.48], 'Gauss', 0.1)
fdm = FrequencyDensityMatrix(td_calc, dmat, frequencies=freqs)
# Propagate according to the saved trajectory
td_calc.autopropagate()
# Store the density matrix
fdm.write('fdm.ulm')
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
# Read the ground-state file
td_calc = LCAOTDDFT('gs.gpw', propagator=dict(name='wfw.ulm', update='all'))
# Attach analysis tools
DipoleMomentWriter(td_calc, 'dm_replayed.dat')
# Propagate according to the saved trajectory
td_calc.autopropagate()
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
from gpaw.lcaotddft.wfwriter import WaveFunctionWriter
# Read the restart file
td_calc = LCAOTDDFT('td.gpw', txt='tdc.out')
# Attach the data recording for appending the new data
DipoleMomentWriter(td_calc, 'dm.dat')
WaveFunctionWriter(td_calc, 'wfw.ulm')
# Continue propagation
td_calc.propagate(20, 500)
# Save the state for restarting later
td_calc.write('tdc.gpw', mode='all')
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment