Convolution trick for TDDFT response

In this tutorial, we start from the system in Observables, and show that the responses to any arbitrary pulses are related to each other, as long as the response is in the linear regime. We can take advantage of this relation by performing one expensive TDDFT calculation with a pulse that has a broad spectral response (like a delta-kick or sinc-shaped pulse), and in post-processing obtain the response to any other pulse.

Preparation

Perform the ground state calculation, the calculation of all unoccupied states, and the Voronoi decomposition calculation, from Observables. We will need the following files:

• gs.gpw

• gs_unocc.gpw

• ksd.ulm

• voronoi.ulm

To show that we can obtain the correct results for the Gaussian pulse response with the convolution trick, we will compare to data files:

• dm_gauss.dat

• hcdist_gauss_proj.dat

• hcdist_gauss_time.dat

No we load the ground state file (the first file, with only the occupied states) and perform time propagation under some external field. We choose a spatially uniform (we are in the long-wavelength limit) field with a temporal shape that is a sinc shape. The sinc shape is chosen such that it starts relatively smoothly from 0, and that its Fourier transform goes to zero above 8eV. In frequency space it thus has similar behavior to the delta-kick while keeping down high-frequency response.

It is instructive to first plot the pulse

import numpy as np
from gpaw.tddft.units import eV_to_au, fs_to_au, au_to_fs, au_to_eV
from matplotlib import pyplot as plt
from gpaw.lcaotddft.laser import SincPulse

def create_pulse():
    strength = 1e-6  # Amplitude
    time0 = 5  # Time of maximum in number of oscillations
    cutoff_freq = 8  # Cut-off frequency in eV

    pulse = SincPulse(strength, time0, cutoff_freq, relative_t0=True)
    return pulse

fig, axes = plt.subplots(1, 2, figsize=(6, 2), dpi=200)

time = np.linspace(0, 30, 500)  # Time grid in femtoseconds for plotting
frequency = 2 * np.pi * np.fft.rfftfreq(2*len(time), (time[1] - time[0]) * fs_to_au) * au_to_eV  # Frequency grid in eV for plotting
pulse = create_pulse()

axes[0].plot(time, pulse.strength(time * fs_to_au))
axes[0].set_xlabel('Time (fs)')
axes[0].set_ylabel('Pulse amplitude')

fourier = np.fft.rfft(pulse.strength(time * fs_to_au), n=2*len(time))
axes[1].plot(frequency, fourier.real, label='Re')
axes[1].plot(frequency, fourier.imag, label='Im')
axes[1].plot(frequency, np.abs(fourier), label='Abs')
axes[1].set_xlabel('Frequency (eV)')
axes[1].legend(loc='upper left')
axes[1].set_xlim([0, 10])

fig.tight_layout()

Now we perform the time propagation. This is out most time consuming task. We refer to the GPAW documentation on how to run using MPI and restart files.

• We attach a DipoleMomentWriter that writes the dipole moment at every timestep to dm.dat

• We attach a WaveFunctionWriter that writes the LCAO wave function coefficients to wfs.ulm. To save space, we write it every 10 steps

import numpy as np
from ase.units import Hartree, Bohr
from gpaw.external import ConstantElectricField
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.laser import GaussianPulse
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
from gpaw.lcaotddft.wfwriter import WaveFunctionWriter

# Temporal shape of the time-dependent potential
pulse = create_pulse()
# Spatial shape of the time-dependent potential
ext = ConstantElectricField(Hartree / Bohr, [0., 0., 1.])
# Time-dependent potential
td_potential = {'ext': ext, 'laser': pulse}

# Set up the time-propagation calculation
td_calc = LCAOTDDFT('gs.gpw',
                    td_potential=td_potential,
                    parallel={'sl_auto': True,
                              'augment_grids': True,
                              'domain': 1},  # Adjust to higher value if out of memory
                    txt='td.out')

# Attach the data recording and analysis tools
DipoleMomentWriter(td_calc, 'dm.dat')
WaveFunctionWriter(td_calc, 'wfs.ulm', interval=10)

# Propagate
td_calc.propagate(20, 1500)

Analyzing the induced dipole moment

Let us plot the induced dipole moment by reading the dm.dat file. The system responds to the applied pulse by oscillating strongly. After the pulse has reached its maximum, however, the dipole moment of the system keeps oscillating. The amplitude of the oscillations slowly decreases over the remainder of the calculation.

from ase.units import Bohr

# Format of the file is (time, norm, dmx, dmy, dmz)
data = np.loadtxt('dm.dat')

time = data[:, 0] * au_to_fs  # Convert from atomic units to femtoseconds
dipmom = data[:, 4] * Bohr  # Convert from atomic units to eÅ

pulse = create_pulse()

fig, ax = plt.subplots(1, 1, figsize=(6, 2), dpi=200)
twinax = ax.twinx()

twinax.plot(time, pulse.strength(time * fs_to_au), color='forestgreen', alpha=0.6)
ax.plot(time, dipmom - dipmom[0], color='cornflowerblue')

ax.set_xlabel('Time (fs)')
ax.set_ylabel('Dipole moment (eÅ)')
twinax.set_ylabel('Pulse amplitude', color='forestgreen')

fig.tight_layout()

Obtaining the respose to other pulses

Suppose that we are interested by the response to another pulse. If we keep the strength of the pulse weak, then there is a trick that we can use to obtain the response without performing another costly time propagation. As an instructive example, consider the dipole moment.

In the linear response regime, we can express the dipole moment \(\boldsymbol{d}\) as the polarizability times the external perturbation in the electric field. In the frequency domain

\[\boldsymbol{d}(\omega) = \boldsymbol{\alpha}(\omega) \boldsymbol{E}_\text{ext}(\omega),\]

or, equivalently, in the time domain

\[\boldsymbol{d}(t) = \int_0^t \boldsymbol{\alpha}(t - t') \boldsymbol{E}_\text{ext}(t') \mathrm{d}t'.\]

Therefore, if we know the dipole response to one perturbation (e.g. \(\boldsymbol{d}_\mathrm{sinc}\) to a sinc perturbation \(\boldsymbol{E}_\text{sinc}\) as earlier), then we know the polarizability (at least at frequencies where the perturbation is non-zero). Hence, we can reconstruct the response to any other perturbation (e.g. \(\boldsymbol{d}_\mathrm{gauss}\) to a Gaussian laser field \(\boldsymbol{E}_\text{gauss}\)). This is easiest obtained in frequency space

\[\begin{split}\begin{align} \boldsymbol{d}_\text{gauss}(\omega) &= \boldsymbol{\alpha}(\omega) \boldsymbol{E}_\text{gauss}(\omega) \\ &= \boldsymbol{d}_\text{sinc}(\omega) \frac{\boldsymbol{E}_\text{gauss}(\omega)}{\boldsymbol{E}_\text{sinc}(\omega)}. \end{align}\end{split}\]

The same principle applies to any other quantity which is linear in the perturbation, including the electron-hole part of the KS density matrix

\[\rho_{ia}^\text{(gauss)}(\omega) = \rho_{ia}^\text{(sinc)}(\omega)\frac{\boldsymbol{E}_\text{gauss}(\omega)}{\boldsymbol{E}_\text{sinc}(\omega)}.\]

The Gaussian pulse has the following shape

import numpy as np
from gpaw.tddft.units import eV_to_au, fs_to_au, au_to_fs
from matplotlib import pyplot as plt
from gpaw.lcaotddft.laser import GaussianPulse

def create_pulse():
    strength = 1e-5  # Amplitude
    time0 = 10e3  # Time of maximum in attoseconds
    frequency = 3.8  # Frequency in eV
    sigma = 0.3  # Width parameter in eV

    # Transform the sigma parameter of the Gaussian to a full-width at half-maximum
    fwhm_eV = np.sqrt(8 * np.log(2)) * sigma
    fwhm_fs = np.sqrt(8 * np.log(2)) / (sigma * eV_to_au) * au_to_fs

    print(f'Full width at half-maximum in time {fwhm_fs:.2f}fs')
    print(f'Full width at half-maximum in frequency {fwhm_eV:.2f}eV')

    # Temporal shape of the time-dependent potential
    pulse = GaussianPulse(strength, time0, frequency, sigma, 'cos')

    return pulse

fig, axes = plt.subplots(1, 2, figsize=(6, 2), dpi=200)

time = np.linspace(0, 30, 500)  # Time grid in femtoseconds for plotting
frequency = np.linspace(0, 6, 500)  # Frequency grid in eV for plotting
pulse = create_pulse()

axes[0].plot(time, pulse.strength(time * fs_to_au))
axes[0].set_xlabel('Time (fs)')
axes[0].set_ylabel('Pulse amplitude')

axes[1].plot(frequency, pulse.fourier(frequency * eV_to_au).real, label='Re')
axes[1].plot(frequency, pulse.fourier(frequency * eV_to_au).imag, label='Im')
axes[1].plot(frequency, np.abs(pulse.fourier(frequency * eV_to_au)), label='Abs')
axes[1].set_xlabel('Frequency (eV)')
axes[1].legend(loc='upper left')

fig.tight_layout()

Full width at half-maximum in time 5.17fs
Full width at half-maximum in frequency 0.71eV

Dipole moment

We set up a ResponseFromWaveFunctions object to tell rhodent what information about the response we have from the TDDFT calculation. We initialize this object with the filename of the KohnShamDecomposition file, the file name of the wave function dump file, and importantly, the type of perturbation that the system was subjected to. It is important that we provide the same perturbation as in the TDDFT calculation, which is in this case a Sinc-pulse.

To compute the dipole moment, we initialize a DipoleCalculator with a list of times, and the pulse that we want to know the response for, which is a Gaussian pulse. Because the Gaussian does not match the Sinc-pulse, the calculator knows that it needs to apply the convolution trick.

Finally we set perform the calculation and (gather the results on the root rank when in MPI) using icalculate_gather_on_root. The parameter v=2 means that we extract only the z component of the dipole moment. We write the results to the file dm_convolution.dat.

The calculation of density matrices, and forwards and backwards Fourier transformation takes a significant amount of time, even if it is orders of magnitudes faster than the time propagation.

import numpy as np
from gpaw.lcaotddft.laser import GaussianPulse
from rhodent.response import ResponseFromWaveFunctions
from rhodent.calculators import DipoleCalculator

# Set up the response object
response = ResponseFromWaveFunctions(
    ksd='ksd.ulm',
    wfs_fname='wfs.ulm',
    perturbation={'name': 'SincPulse', 'strength': 1e-6, 'cutoff_freq': 8, 'time0': 5, 'relative_t0': True})

# Set up the dipole calculator
pulse = GaussianPulse(1e-5, 10e3, 3.8, 0.3, 'cos')
calc = DipoleCalculator(response,
                        times=np.arange(0, 30000, 100),
                        pulses=[pulse])

dipmom_from_calc = []
# Iterate over the results.
for metadata, res in calc.icalculate_gather_on_root(v=2):
    dipmom_from_calc.append(res['dm'][0])

if calc.loop_comm.rank == 0 and calc.calc_comm.rank == 0:
    dipmom_from_calc = np.array(dipmom_from_calc)
    np.savetxt('dm_convolution.dat',
               np.vstack((calc.times * 1e-3, dipmom_from_calc)).T,
               fmt=['%.3f', '%.8e'])

The obtained dipole moment matches perfectly with the one obtained in Observables.

from ase.units import Bohr

# Format of the file is (time, norm, dmx, dmy, dmz)
data_direct = np.loadtxt('dm_gauss.dat')

time = data_direct[:, 0] * au_to_fs  # Convert from atomic units to femtoseconds
dipmom = data_direct[:, 4] * Bohr  # Convert from atomic units to eÅ

data_convolution = np.loadtxt('dm_convolution.dat')

fig, ax = plt.subplots(1, 1, figsize=(6, 2), dpi=200)

ax.plot(data_convolution[:, 0], data_convolution[:, 1], ls='-', color='goldenrod', label='Convolution')
ax.plot(time, dipmom - dipmom[0], color='cornflowerblue', ls='--', label='Direct')

ax.set_xlabel('Time (fs)')
ax.set_ylabel('Dipole moment (eÅ)')
ax.legend()

fig.tight_layout()

Obtaining projected hot-carrier distributions

We obtain projected hot-carrier distributions just like in Observables, with the only difference that we provide a the perturbation parameter to our ResponseFromWaveFunctions and a list of pulses to the calculator. Again, the agreement is good!

import numpy as np
from gpaw.lcaotddft.laser import GaussianPulse
from rhodent.response  import ResponseFromWaveFunctions
from rhodent.calculators import HotCarriersCalculator
from rhodent.voronoi import VoronoiReader

# Set up the density matrices object
response = ResponseFromWaveFunctions(
    ksd='ksd.ulm',
    wfs_fname='wfs.ulm',
    perturbation={'name': 'SincPulse', 'strength': 1e-6, 'cutoff_freq': 8, 'time0': 5, 'relative_t0': True})

# Set up the Voronoi file reader
voronoi = VoronoiReader('voronoi.ulm')

# Set up the hot-carriers calculator that calculates and broadens dipole contributions onto an energy grid
pulse = GaussianPulse(1e-6, 10e3, 3.8, 0.3, 'cos')
calc = HotCarriersCalculator(response,
                             times=[30000],
                             pulses=[pulse],
                             voronoi=voronoi,
                             energies_occ=np.arange(-5, 1, 0.01),
                             energies_unocc=np.arange(-1, 5, 0.01),
                             sigma=0.07)

# Calculate the hot carrier distributions and save them to a file
calc.calculate_hcdist_and_write('hcdist_convolution.dat', average_times=False)

import numpy as np
from matplotlib import pyplot as plt

data = np.loadtxt('hcdist_gauss_proj.dat')
data_convolution = np.loadtxt('hcdist_convolution.dat')

fig, axes = plt.subplots(2, 2, sharey='row', sharex='col', figsize=(6, 3), dpi=200)

axes[0, 0].plot(data_convolution[:, 0], data_convolution[:, 2])
axes[0, 1].plot(data_convolution[:, 1], data_convolution[:, 3], label='Convolution')
axes[1, 0].plot(data_convolution[:, 0], data_convolution[:, 4])
axes[1, 1].plot(data_convolution[:, 1], data_convolution[:, 5])
axes[0, 0].plot(data[:, 0], data[:, 2], 'k', ls='dotted')
axes[0, 1].plot(data[:, 1], data[:, 3], 'k', ls='dotted', label='Direct')
axes[1, 0].plot(data[:, 0], data[:, 4], 'k', ls='dotted')
axes[1, 1].plot(data[:, 1], data[:, 5], 'k', ls='dotted')

axes[0, 0].set_ylabel('Total carriers \n (1/eV)')
axes[1, 0].set_ylabel('Carriers on \n molecule (1/eV)')
axes[-1, 0].set_xlabel('Energy holes (eV)')
axes[-1, 1].set_xlabel('Energy electrons (eV)')

fig.tight_layout()

Obtaining the total number of carriers

Similarly, we obtain the total number of hot carriers in time and compare then to the same quantity obtained in Observables.

import numpy as np
from gpaw.lcaotddft.laser import GaussianPulse
from rhodent.response  import ResponseFromWaveFunctions
from rhodent.calculators import HotCarriersCalculator
from rhodent.voronoi import VoronoiReader

# Set up the density matrices object
response = ResponseFromWaveFunctions(
    ksd='ksd.ulm',
    wfs_fname='wfs.ulm',
    perturbation={'name': 'SincPulse', 'strength': 1e-6, 'cutoff_freq': 8, 'time0': 5, 'relative_t0': True})

# Set up the Voronoi file reader
voronoi = VoronoiReader('voronoi.ulm')

# Set up the hot-carriers calculator that calculates and broadens dipole contributions onto an energy grid
pulse = GaussianPulse(1e-6, 10e3, 3.8, 0.3, 'cos')
calc = HotCarriersCalculator(response,
                             times=np.arange(0, 30e3, 100),
                             pulses=[pulse],
                             voronoi=voronoi)

# Calculate the hot carrier distributions and save them to a file
calc.calculate_hcsweeptime_and_write('hcdist_convolution_time.dat')

!head hcdist_gauss_time.dat

# Hot carrier (H=hole, E=electron) numbers
#
# Atomic projections:
#      0: [201, 202]
#
#       Time (fs)         Total H         Total E       Proj H  0       Proj E  0
         0.020000  2.72573372e-19  2.72573372e-19  2.27957876e-21  3.84953606e-21
         0.220000  5.19251063e-17  5.19251063e-17  3.22714986e-19  7.63459641e-19
         0.420000  4.05332005e-17  4.05332005e-17  1.75037592e-19  7.18875379e-19
         0.620000  3.41373472e-16  3.41373472e-16  1.78031882e-18  3.53014396e-18

import numpy as np
from matplotlib import pyplot as plt

data = np.loadtxt('hcdist_convolution_time.dat')
data_gauss = np.loadtxt('hcdist_gauss_time.dat')

time = data[:, 0]

fig, axes = plt.subplots(3, 1, sharey='row', sharex='col', figsize=(6, 4), dpi=200)

axes[0].plot(time, pulse.strength(time * fs_to_au))
axes[1].plot(time, data[:, 2], label='Convolution')
axes[1].plot(time, data_gauss[:, 2], 'k', ls='dotted', label='Direct')
axes[2].plot([], [], 'k', label='Convolution')
axes[2].plot([], [], 'k', ls='dotted', label='Direct')
axes[2].plot(time, data[:, 4], label='Electrons')
axes[2].plot(time, data[:, 3], label='Holes')
axes[2].plot(time, data_gauss[:, 4], 'k', ls='dotted')
axes[2].plot(time, data_gauss[:, 3], 'k', ls='dotted')

axes[0].set_ylabel('Pulse')
axes[1].set_ylabel('Total carriers')
axes[2].set_ylabel('Carriers on \n molecule')
axes[1].legend(fontsize=8)
axes[2].legend(fontsize=8)
axes[-1].set_xlabel('Time (fs)')

fig.tight_layout()

Response to different pulses

With little added cost, we can also calculate the response to many different pulses at once by setting up a list of pulses and using the calculate_hcsweeppulse_and_write function.

import numpy as np
from gpaw.lcaotddft.laser import GaussianPulse
from rhodent.response  import ResponseFromWaveFunctions
from rhodent.calculators import HotCarriersCalculator
from rhodent.voronoi import VoronoiReader

# Set up the density matrices object
response = ResponseFromWaveFunctions(
    ksd='ksd.ulm',
    wfs_fname='wfs.ulm',
    perturbation={'name': 'SincPulse', 'strength': 1e-6, 'cutoff_freq': 8, 'time0': 5, 'relative_t0': True})

# Set up the Voronoi file reader
voronoi = VoronoiReader('voronoi.ulm')

# Set up the hot-carriers calculator that calculates and broadens dipole contributions onto an energy grid
pulses = [GaussianPulse(1e-5, 10e3, frequency, 0.3, 'cos') for frequency in np.arange(3.5, 4.2001, 0.1)]
calc = HotCarriersCalculator(response,
                             times=[30000],
                             pulses=pulses,
                             voronoi=voronoi)

# Calculate the hot carrier distributions and save them to a file
calc.calculate_hcsweeppulse_and_write('hcdist_convolution_pulses.dat')

!head hcdist_gauss_pulses.dat

head: cannot open 'hcdist_gauss_pulses.dat' for reading: No such file or directory

import numpy as np
from matplotlib import pyplot as plt

data = np.loadtxt('hcdist_convolution_pulses.dat')

pulse_freq = data[:, 0]

fig, axes = plt.subplots(2, 1, sharey='row', sharex='col', figsize=(6, 3), dpi=200)

axes[0].plot(pulse_freq, data[:, 3])
axes[1].plot(pulse_freq, data[:, 4], label='Holes')
axes[1].plot(pulse_freq, data[:, 5], label='Electrons')

axes[0].set_ylabel('Total carriers')
axes[1].set_ylabel('Carriers on \n molecule')
axes[-1].set_xlabel('Pulse frequency (eV)')
axes[1].legend()

fig.tight_layout()