Dielectric functions

This tutorial demonstrates how to compute the dielectric function tensor \(\epsilon^{\alpha\beta}(\omega) = \epsilon_1^{\alpha\beta}(\omega) + i \epsilon_2^{\alpha\beta}(\omega)\) from a GPUMD molecular-dynamics simulation. GPUMD writes the time derivative of the polarization \(\dot{\mathbf{P}}(t)\) to a file named dpdt.out when using a qNEP model. The imaginary part of the dielectric tensor \(\epsilon_2^{\alpha\beta}(\omega)\) is related to the power spectral density of \(\dot{P}_\alpha(t)\) via

\[\epsilon_2^{\alpha\beta}(\omega) = \frac{1}{\epsilon_0 V k_\mathrm{B} T}\, \int_0^\infty \langle \dot{P}_\alpha(0)\,\dot{P}_\beta(\tau) \rangle\, e^{-i\omega\tau}\, d\tau,\]

where \(V\) is the simulation cell volume, \(T\) is the temperature, and \(\epsilon_0\) is the permittivity of free space. The isotropic average \(\bar{\epsilon}_2(\omega) = \frac{1}{3}\sum_\alpha \epsilon_2^{\alpha\alpha}(\omega)\) reduces to the scalar dielectric function for cubic or isotropic systems. The real part \(\epsilon_1^{\alpha\beta}(\omega)\) is obtained component-wise via the Kramers-Kronig relation

\[\epsilon_1^{\alpha\beta}(\omega) = \frac{2}{\pi}\, \mathcal{P}\int_0^\infty \frac{\omega'\,\epsilon_2^{\alpha\beta}(\omega')}{\omega'^2 - \omega^2}\, d\omega'.\]

The example data are qNEP molecular-dynamics simulations of BaTiO\(_3\) at 100, 221, 331, and 450 K. These are reruns using input files similar to those of Z. Fan et al., J. Chem. Theory Comput. 22, 4787 (2026) using the same structures, but with adapted run parameters.

The data required for running this tutorial notebook can be obtained from Zenodo.

import os
import zipfile
import urllib.request

url = 'https://zenodo.org/records/21198312/files/dielectric_functions.zip'
if not os.path.exists('dielectric_functions'):
    urllib.request.urlretrieve(url, 'dielectric_functions.zip')
    with zipfile.ZipFile('dielectric_functions.zip') as zf:
        zf.extractall('.')
    os.remove('dielectric_functions.zip')
    print('Downloaded and extracted dielectric_functions/')
else:
    print('dielectric_functions/ already present')
os.chdir('dielectric_functions')

dielectric_functions/ already present

Reading the data

Each temperature has its own subdirectory (T100, T221, T331, T450) containing a 5 ns production run, with model.xyz (the simulation cell) and dpdt.out (the polarization time derivative, dumped every 20 fs). The read_dpdt function reads a dpdt.out file into a pandas.DataFrame. The columns dPx, dPy, dPz are the Cartesian components of \(\dot{\mathbf{P}}(t)\) in units of \(e\)Å/fs, and Px, Py, Pz are the corresponding polarization components in units of \(e\)Å. The time step is read directly from the time column as the difference between consecutive rows.

from ase.io import read
from calorine.gpumd import read_dpdt

temperatures = [100, 221, 331, 450]  # K

dfs_raw = {T: read_dpdt(f'T{T}/dpdt.out') for T in temperatures}
volumes = {T: read(f'T{T}/model.xyz').get_volume() for T in temperatures}

for T in temperatures:
    print(f'{T:4} K: volume = {volumes[T]:.1f} Angstrom^3, {len(dfs_raw[T])} frames')

dfs_raw[450].head()

 100 K: volume = 527788.9 Angstrom^3, 250000 frames
 221 K: volume = 529222.8 Angstrom^3, 250000 frames
 331 K: volume = 530789.4 Angstrom^3, 250000 frames
 450 K: volume = 531196.7 Angstrom^3, 250000 frames

	time	dPx	dPy	dPz	Px	Py	Pz
0	0.02	-1.926350	-4.03456	-0.377649	-38.5271	-80.69120	-7.55298
1	0.04	0.434974	-1.83067	-0.611361	-29.8276	-117.30500	-19.78020
2	0.06	0.117376	1.55974	-0.388794	-27.4801	-86.10990	-27.55610
3	0.08	-2.777010	4.71552	0.811376	-83.0202	8.20046	-11.32860
4	0.10	0.286943	5.17358	-1.085420	-77.2814	111.67200	-33.03690

Computing the dielectric function

get_dielectric_function computes both parts of the dielectric tensor from all three Cartesian components of \(\dot{\mathbf{P}}(t)\), returning the six unique Voigt components of the imaginary part \(\epsilon_2^{\alpha\beta}(\omega)\) together with the real part \(\epsilon_1^{\alpha\beta}(\omega)\) obtained via the Kramers-Kronig relation. The required physical parameters are

• volume: the simulation cell volume in Å\(^3\),

• temperature: the temperature in K.

The keyword column='dP' (the default) selects the dPx/dPy/dPz columns from read_dpdt. The time step is extracted automatically from the time column. The optional parameter t_sigma sets the width of a Gaussian window applied to the autocorrelation function before the Fourier transform (in ps); smaller values produce smoother spectra while larger values yield sharper features. window_size sets how much of the autocorrelation function (in ps, the same unit as t_sigma) is kept before the Fourier transform. If left unset, it defaults to 5 * t_sigma when t_sigma is given, or to the full trajectory otherwise.

A performance note: window_size directly sets the size of the frequency grid handed to the Kramers-Kronig step, and the default kk_method='vectorized' scales as \(\mathcal{O}(n^2)\) in both time and memory. With our 5 ns/250,000-frame runs, leaving window_size unset together with t_sigma=2 ps would already give a modest ~10,000-point grid (5 * t_sigma), but a longer t_sigma or an explicit large window_size can easily produce a grid whose vectorized Kramers-Kronig step would need tens of gigabytes. get_dielectric_function estimates this upfront and raises a clear error (governed by kk_max_memory_gb, default 4 GB) rather than letting the process hang or get killed due to an out-of-memory (OOM) error. Since the t_sigma=2 ps Gaussian window already damps the autocorrelation function well within a much shorter lag, we explicitly pass a window_size=100 ps window below (rather than relying on the 5 * t_sigma = 10 ps default) to be conservative. This gives a visually identical spectrum to a much larger window while keeping the computation fast and memory-friendly. We are, however, still benefiting from the oversampling of a long trajectory; see section on convergence with run length below.

%%time
from calorine.tools import get_dielectric_function

t_sigma = 2  # Gaussian window width in ps; smaller values give smoother spectra
window_size = 100  # ps; see note above on why this is larger than the 5*t_sigma default

dfs = {}
for T in temperatures:
    dfs[T] = get_dielectric_function(
        dfs_raw[T], volumes[T], T, column='dP', t_sigma=t_sigma, window_size=window_size)
dfs[450].head()

CPU times: user 11min 59s, sys: 969 ms, total: 12min
Wall time: 1min 33s

	angular_frequency	wavenumber_invcm	epsilon_imag_xx	epsilon_imag_yy	epsilon_imag_zz	epsilon_imag_yz	epsilon_imag_xz	epsilon_imag_xy	conductivity_xx	conductivity_yy	conductivity_zz	conductivity_yz	conductivity_xz	conductivity_xy	epsilon_real_xx	epsilon_real_yy	epsilon_real_zz	epsilon_real_yz	epsilon_real_xz	epsilon_real_xy
0	0.031419	0.166799	22854.007479	10772.389987	8401.738100	-1167.499899	-458.920797	1462.165912	6357.763952	2996.774758	2337.282320	-324.787186	-127.667329	406.760422	14262.693743	7828.720257	6471.713500	-626.836486	-262.160865	757.196779
1	0.062838	0.333597	11521.382994	5440.903286	4247.011204	-587.167153	-230.380961	733.964324	6410.274745	3027.213394	2362.954922	-326.688452	-128.179513	408.363559	8522.998682	5125.966239	4364.511263	-333.296968	-146.694889	389.278746
2	0.094257	0.500396	7785.542951	3688.004484	2882.565802	-395.227682	-154.602177	492.492169	6497.588360	3077.901587	2405.705821	-329.845561	-129.026493	411.019682	5490.832410	3697.607153	3250.531107	-178.350842	-85.802425	195.233906
3	0.125676	0.667195	5948.625492	2829.694873	2215.641090	-300.370982	-117.004847	372.681416	6619.400428	3148.775036	2465.479731	-334.241214	-130.198470	414.705469	4530.191194	3245.540609	2897.849206	-129.275370	-66.551626	133.849501
4	0.157095	0.833994	4870.973835	2329.156222	1827.668828	-244.330262	-94.670503	301.512234	6775.289167	3239.743726	2542.199000	-339.851585	-131.682094	419.388943	4094.135852	3040.804083	2737.994672	-107.016662	-57.860336	106.091472

Kramers-Kronig integration method

The Kramers-Kronig step is performed internally by get_dielectric_function (controlled by return_real_part, which defaults to True). Two integration methods are available via the kk_method keyword:

• 'vectorized' (default): exact trapezoid rule, \(\mathcal{O}(n^2)\) time and memory.

• 'fft': \(\mathcal{O}(n \log n)\) Hilbert-transform approximation; faster for large arrays but can deviate near sharp features.

The section below compares both methods on the 100 K data, which has the sharpest phonon peaks and provides the most stringent test.

Visualization

from matplotlib import pyplot as plt

def isotropic_average(df, prefix):
    return (df[f'{prefix}_xx'] + df[f'{prefix}_yy'] + df[f'{prefix}_zz']) / 3


fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 3.6), dpi=140, sharex=True)

cmap = plt.colormaps['coolwarm']
for T in temperatures:
    df = dfs[T]
    kwargs = dict(label=f'{T} K', color=cmap((T - 100) / (450 - 100)), alpha=0.7)
    ax1.plot(df.wavenumber_invcm, isotropic_average(df, 'epsilon_real'), **kwargs)
    ax2.plot(df.wavenumber_invcm, isotropic_average(df, 'epsilon_imag'), **kwargs)

ax1.axhline(0, color='0.5', linewidth=1, linestyle='--')
ax1.set_ylabel(r'$\widebar{\epsilon}_1(\omega)$')
ax1.legend(frameon=False, fontsize='small')
ax1.set_xlim(0, 830)
ax1.set_ylim(-400, 1400)

ax2.set_xlabel(r'Wavenumber (cm$^{-1}$)')
ax2.set_ylabel(r'$\widebar{\epsilon}_2(\omega)$')
ax2.set_yscale('log')
ax2.set_ylim(0.1, 3e3)

fig.tight_layout()
fig.subplots_adjust(hspace=0)

Convergence with run length

Each of the runs above is 5 ns long, which allows substantial oversampling and thus reduces numerical noise by better statistics. To see how much this matters, we recompute the spectrum for one temperature (331 K) using only the first 0.1, 0.2, 0.5, 1, or 5 ns of the same trajectory; get_dielectric_function takes the time series DataFrame directly, so a shorter “run” is a row-count truncation of it. To isolate the effect of run length from the choice of autocorrelation window, we keep the window itself fixed at 50 ps for every truncation; since window_size is now an absolute time in ps rather than a fraction of the trajectory length, this requires no rescaling. What changes between curves is thus purely how much data went into estimating the autocorrelation function at each lag, not the frequency resolution.

T_conv = 331
lengths_ns = [0.1, 0.2, 0.5, 1, 5]
window_ps = 50  # fixed autocorrelation window, independent of run length

fig, ax = plt.subplots(figsize=(4, 2.8), dpi=140)

cmap = plt.colormaps['viridis']
for k, length_ns in enumerate(lengths_ns):
    df = dfs_raw[T_conv]
    df = df[df.time <= 1e3 * length_ns]
    actual_length_ns = df.time.max() / 1e3
    result = get_dielectric_function(
        df, volumes[T_conv], T_conv, column='dP',
        t_sigma=t_sigma, window_size=window_ps, return_real_part=False)
    ax.plot(result.wavenumber_invcm, isotropic_average(result, 'epsilon_imag'),
            color=cmap(k / (len(lengths_ns) - 1)), label=f'{actual_length_ns} ns')

ax.set_xlabel(r'Wavenumber (cm$^{-1}$)')
ax.set_ylabel(r'$\widebar{\epsilon}_2(\omega)$')
ax.set_yscale('log')
ax.set_xlim(0, 830)
ax.legend(title='Run length', frameon=False, fontsize='small', title_fontsize='small')

fig.tight_layout()

The overall shape is already recognizable from just 0.1 ns, but the shorter runs are visibly noisier. Longer trajectories average over more statistically independent time origins, giving a smoother, more reliable estimate of the autocorrelation function at each lag. 5 ns is comfortably enough to suppress this noise for this system.

Accuracy: FFT vs. vectorized

dfs[100] was computed with the default kk_method='vectorized'. The next cell below recomputes the 100 K spectrum with kk_method='fft' for comparison. The two methods agree well across most of the spectrum, but the FFT approximation can deviate near sharp features and at larger frequencies.

%%time

df_vec = dfs[100].iloc[5:]
df_fft = get_dielectric_function(dfs_raw[100], volumes[100], 100, column='dP',
                                 t_sigma=t_sigma, window_size=window_size, kk_method='fft').iloc[5:]

CPU times: user 3min 1s, sys: 31 ms, total: 3min 1s
Wall time: 22.7 s

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 3.6), dpi=140, sharex=True)

ax1.plot(df_vec.wavenumber_invcm, isotropic_average(df_vec, 'epsilon_real'), label='Vectorized')
ax1.plot(df_fft.wavenumber_invcm, isotropic_average(df_fft, 'epsilon_real'), '--', label='FFT')
ax1.axhline(0, color='k', linewidth=0.5, linestyle='--')
ax1.set_ylabel(r'$\bar{\epsilon}_1(\omega)$')
ax1.legend(frameon=False, fontsize='small')

diff = isotropic_average(df_fft, 'epsilon_real').to_numpy() - isotropic_average(df_vec, 'epsilon_real').to_numpy()
ax2.plot(df_vec.wavenumber_invcm, diff)
ax1.axhline(0, color='0.5', linewidth=1, linestyle='--')
ax2.set_xlabel(r'Wavenumber (cm$^{-1}$)')
ax2.set_ylabel('FFT − Vectorized')

fig.tight_layout()
fig.subplots_adjust(hspace=0)

Kramers-Kronig performance

The timing comparison below illustrates the performance advantage of the 'fft' method over the 'vectorized' approach, which comes at a small cost in accuracy (see above).

%%time

from calorine.tools import apply_kramers_kronig
import time
import numpy as np
from pandas import DataFrame

sizes = [100, 500, 1000, 5000, 10000]
times_vec = []
times_fft = []

for n in sizes:
    df_test = DataFrame({'angular_frequency': np.linspace(0.1, 100.0, n),
                         'epsilon_imag': np.sin(np.linspace(0.1, 100.0, n))})

    nrep = max(3, 200 // n)
    t0 = time.perf_counter()
    for _ in range(nrep):
        apply_kramers_kronig(df_test, method='vectorized')
    times_vec.append((time.perf_counter() - t0) / nrep * 1000)

    t0 = time.perf_counter()
    for _ in range(500):
        apply_kramers_kronig(df_test, method='fft')
    times_fft.append((time.perf_counter() - t0) / 500 * 1000)

CPU times: user 2.25 s, sys: 543 ms, total: 2.79 s
Wall time: 2.79 s

fig, ax = plt.subplots(figsize=(4, 2.8), dpi=140)
ax.loglog(sizes, times_vec, 'o-', label=r"method='vectorized', $\mathcal{O}(n^2)$")
ax.loglog(sizes, times_fft, 's-', label=r"method='fft', $\mathcal{O}(n \log n)$")
ax.set_xlabel('Array length $n$')
ax.set_ylabel('Wall time (ms)')
ax.legend(frameon=False, fontsize='small')
fig.tight_layout()

Refractive index and absorption coefficient

The complex dielectric function can be converted into the complex index of refraction \(N(\omega) = n(\omega) + i\kappa(\omega)\) via get_refractive_index, and the extinction coefficient \(\kappa(\omega)\) into the absorption coefficient \(\alpha(\omega)\) (in cm\(^{-1}\)) via get_absorption_coefficient. Both functions operate component-wise on the same Voigt-labeled columns produced by get_dielectric_function, so they can be chained directly onto its output. The inverse conversion, from \(n,\kappa\) (or \(n,\alpha\)) back to \(\epsilon_1,\epsilon_2\), is available via get_dielectric_function_from_refractive_index.

from calorine.tools import get_refractive_index, get_absorption_coefficient

df_nk = get_refractive_index(dfs[450])
df_nk = get_absorption_coefficient(df_nk)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 3.6), dpi=140, sharex=True)

ax1.plot(df_nk.wavenumber_invcm, isotropic_average(df_nk, 'refractive_index_real'), label='n')
ax1.plot(df_nk.wavenumber_invcm, isotropic_average(df_nk, 'refractive_index_imag'), label=r'$\kappa$')
ax1.legend(frameon=False, fontsize='small')
ax1.set_ylabel('Refractive index')
ax1.set_xlim(0, 830)

ax2.plot(df_nk.wavenumber_invcm, isotropic_average(df_nk, 'absorption_coefficient'))
ax2.set_xlabel(r'Wavenumber (cm$^{-1}$)')
ax2.set_ylabel(r'$\bar\alpha$ (cm$^{-1}$)')

fig.tight_layout()
fig.subplots_adjust(hspace=0)