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, 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

Now 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, 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

data = np.load('dm.npz')

time = data['time_t'] * au_to_fs  # Convert from atomic units to femtoseconds
dipmom = data['dm_tv'][:, 2] * 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 obatain 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

The calculation of density matrices takes some time, even if it is orders of magnitudes faster than the time propagation. We set up a ConvolutionDensityMatricesFromWaveFunctions that reads the wave functions, transforms them into density matrices and performs the convolution with the Gausssian pulse. Then we set up a DipoleCalculator and 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.

import numpy as np
from gpaw.lcaotddft.laser import GaussianPulse
from rhodent.density_matrices import ConvolutionDensityMatricesFromWaveFunctions
from rhodent.calculators.dipole import DipoleCalculator

pulse = GaussianPulse(1e-5, 10e3, 3.8, 0.3, 'cos')

# Set up the density matrices object
density_matrices = ConvolutionDensityMatricesFromWaveFunctions(
    ksd='ksd.ulm',
    wfs_fname='wfs.ulm',
    pulses=[pulse],
    perturbation={'name': 'SincPulse', 'strength': 1e-6, 'cutoff_freq': 8, 'time0': 5, 'relative_t0': True},
    times=np.arange(0, 30000, 100))

# Set up the dipole calculator
calc = DipoleCalculator(density_matrices, None, energies_occ=[], energies_unocc=[])

dipmom_from_calc = []
# Iterate over the results. Since we have only specified one time, the loop will only iterate once
for metadata, res in calc.icalculate_gather_on_root(v=2):
    dipmom_from_calc.append(res['dm'][0])

if density_matrices.loop_comm.rank == 0 and density_matrices.calc_comm.rank == 0:
    dipmom_from_calc = np.array(dipmom_from_calc)
    np.savetxt('dm_convolution.dat',
               np.vstack((density_matrices.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
data_direct = np.load('dm_gauss.npz')

time = data_direct['time_t'] * au_to_fs  # Convert from atomic units to femtoseconds
dipmom = data_direct['dm_tv'][:, 2] * 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 use our ConvolutionDensityMatricesFromWaveFunctions with the appropriate parameters instead of TimeDensityMatricesFromWaveFunctions to perform the convolution. Again, the agreement is good!

import numpy as np
from gpaw.lcaotddft.laser import GaussianPulse
from rhodent.density_matrices import ConvolutionDensityMatricesFromWaveFunctions
from rhodent.voronoi import VoronoiReader
from rhodent.writers.hcdist import calculate_hcdist_and_save_by_filename

pulse = GaussianPulse(1e-6, 10e3, 3.8, 0.3, 'cos')

# Set up the density matrices object
density_matrices = ConvolutionDensityMatricesFromWaveFunctions(
    ksd='ksd.ulm',
    wfs_fname='wfs.ulm',
    pulses=[pulse],
    perturbation={'name': 'SincPulse', 'strength': 1e-6, 'cutoff_freq': 8, 'time0': 5, 'relative_t0': True},
    times=[30000])

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

# Calculate the hot carrier distributions and save them to a file
calculate_hcdist_and_save_by_filename('hcdist_convolution.dat', density_matrices=density_matrices, voronoi=voronoi, energies_occ=np.arange(-5, 1, 0.01), energies_unocc=np.arange(-1, 5, 0.01), sigma=0.07, 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.density_matrices import ConvolutionDensityMatricesFromWaveFunctions
from rhodent.voronoi import VoronoiReader
from rhodent.writers.hcdist import calculate_hcsweeptime_and_save_by_filename

pulse = GaussianPulse(1e-5, 10e3, 3.8, 0.3, 'cos')

# Set up the density matrices object
density_matrices = ConvolutionDensityMatricesFromWaveFunctions(
    ksd='ksd.ulm',
    wfs_fname='wfs.ulm',
    pulses=pulse,
    perturbation={'name': 'SincPulse', 'strength': 1e-6, 'cutoff_freq': 8, 'time0': 5, 'relative_t0': True},
    times=np.arange(0, 30e3, 100))

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

# Calculate the hot carrier distributions and save them to a file
calculate_hcsweeptime_and_save_by_filename('hcdist_convolution_time.dat', density_matrices=density_matrices, voronoi=voronoi)

!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.010000  1.74890287e-20  1.74890287e-20  1.90934760e-22  2.82953480e-22
         0.110000  3.65596345e-18  3.65596345e-18  3.21160307e-20  6.43713037e-20
         0.210000  1.36184834e-17  1.36184834e-17  1.10310503e-19  2.58821552e-19
         0.310000  1.96907207e-17  1.96907207e-17  1.40509699e-19  4.03688628e-19

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 using the calculate_hcsweeppulse_and_save_by_filename

import numpy as np
from gpaw.lcaotddft.laser import GaussianPulse
from rhodent.density_matrices import ConvolutionDensityMatricesFromWaveFunctions
from rhodent.voronoi import VoronoiReader
from rhodent.writers.hcdist import calculate_hcsweeppulse_and_save_by_filename

pulses = [GaussianPulse(1e-5, 10e3, frequency, 0.3, 'cos') for frequency in np.arange(3.5, 4.2001, 0.1)]

# Set up the density matrices object
density_matrices = ConvolutionDensityMatricesFromWaveFunctions(
    ksd='ksd.ulm',
    wfs_fname='wfs.ulm',
    pulses=pulses,
    perturbation={'name': 'SincPulse', 'strength': 1e-6, 'cutoff_freq': 8, 'time0': t, 'relative_t0': True},
    times=[30000])

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

# Calculate the hot carrier distributions and save them to a file
calculate_hcsweeppulse_and_save_by_filename('hcdist_convolution_pulses.dat', density_matrices=density_matrices, voronoi=voronoi)

!head hcdist_gauss_pulses.dat

# Hot carrier (H=hole, E=electron) numbers
# Averaged for 1 times between 30.0fs-30.0fs
# Atomic projections:
#      0: [201, 202]
# Pulse freq (eV)   Pulse FWHM (fs)         Total H         Total E   Proj H  0       Proj E  0
         3.500000          5.166569  9.23968202e-05  9.23968202e-05  2.81704628e-07  9.64467041e-07
         3.600000          5.166569  1.20374444e-04  1.20374444e-04  3.61656036e-07  1.27872746e-06
         3.700000          5.166569  1.35025282e-04  1.35025282e-04  4.01493548e-07  1.50933146e-06
         3.800000          5.166569  1.29463977e-04  1.29463977e-04  3.79063480e-07  1.54539065e-06
         3.900000          5.166569  1.06959168e-04  1.06959168e-04  3.02644221e-07  1.37692576e-06

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