# This module contains components adapted from GPAW:
# gpaw.tddft.spectrum
# https://gitlab.com/gpaw/gpaw/-/blob/aca9ed6f520d6a855013247119b50c630f5eebc9/gpaw/tddft/spectrum.py
import re
from collections import namedtuple
from functools import partial
from pathlib import Path
from typing import List, Optional, Tuple, Union
import numpy as np
from ase import Atoms
import ase.io
from . import TimeDomainResponse
from .utilities import Broadening
from ..types import Path_str
Kick = namedtuple('Kick', ['time', 'strength_v'])
[docs]
def read_gpaw_tddft(dipole_files: Union[Path_str, List[Path_str]],
broadening: Broadening,
atoms_file: Optional[Path_str] = None,
name: str = 'GPAWTDDFT') -> TimeDomainResponse:
"""Builds a :class:`~strongcoca.response.TimeDomainResponse` instance
using time-dependent dipole moment data from GPAW real-time
time-dependent density functional theory (RT-TDDFT) calculations.
Parameters
----------
dipole_files
Name(s) of file(s) with time-dependent dipole moment data; each file
is assumed to provide a different Cartesian direction;
if there is only one file, the response is assumed to be diagonal.
broadening
Artificial broadening applied to the response.
atoms_file
Name of file from which to read atomic coordinates; needs to be readable
by `ASE <https://ase-lib.org/ase/io/io.html>`_.
name
Name of the response.
Returns
-------
response
Time domain response.
Raises
------
ValueError
If dipole moment files are not suitable or mutually compatible.
"""
atol_time = 1e-6
allclose_time = partial(np.allclose, rtol=0, atol=atol_time)
if isinstance(dipole_files, (Path, str)):
dipole_files = [dipole_files]
elif len(dipole_files) != 3:
raise ValueError('Three dipole files expected')
# Read files
data_i = []
Data = namedtuple('Data', ['kickv', 'time_t', 'dm_tv'])
for fpath in dipole_files:
time_t, dm_tv, kick_v = _read_and_clean_dipole_moment_data(fpath)
dm_tv -= dm_tv[0]
# Find kick direction
kickv = np.argmax(np.abs(kick_v))
for v in range(3):
if v != kickv and not np.isclose(kick_v[v], 0, atol=1e-24):
raise ValueError(f'Kick in {fpath} not aligned in '
f'Cartesian direction')
data_i.append(Data(kickv, time_t, dm_tv / kick_v[kickv]))
# Check file consistency
time_t = data_i[0].time_t
if not all([allclose_time(d.time_t, time_t) for d in data_i]):
raise ValueError('Inconsistent times in dipole files')
# Assume diagonal response if single dipole file given
if len(data_i) == 1:
data = data_i[0]
dm_t = data.dm_tv[:, data.kickv]
# Build dipole tensor
dm_tvv = np.eye(3)[np.newaxis, :, :] * dm_t[:, np.newaxis, np.newaxis]
else:
# Check that all kick directions are included
kickv_i = [d.kickv for d in data_i]
for v in range(3):
if v not in kickv_i:
raise ValueError(f'Kick direction {"xyz"[v]} missing')
# Build dipole tensor
dm_tvv = np.zeros((len(time_t), 3, 3))
for data in data_i:
dm_tvv[:, data.kickv, :] = data.dm_tv
# Read atoms
if atoms_file is not None:
atoms = ase.io.read(atoms_file)
assert isinstance(atoms, Atoms)
else:
atoms = None
# Initialize response object
response = TimeDomainResponse(time_t, dm_tvv, broadening=broadening,
name=name, units='au')
response._atoms = atoms
return response
def _read_and_clean_dipole_moment_data(filepath: Path_str, atol_time: float = 1e-6) \
-> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Reads and cleans a time-dependent dipole moment file from a real-time
time-dependent density functional theory calculation.
Parameters
----------
filepath
File path.
atol_time
Absolute tolerance for comparing time values.
Returns
-------
time_t
Time values with time zero set to the kick time; shape {t}.
dm_tv
Dipole moment values, data epsilon before kick discarded; shape {t}x{3}.
kick_v
Kick strength; shape {3}.
Raises
------
ValueError
If dipole moment file contains multiple kicks or
if dipole moment file has varying time step.
"""
kick_i, time_t, _, dm_tv = _read_dipole_moment_file(filepath)
# Check kick
if len(kick_i) > 1:
raise ValueError(f'Multiple kicks in {filepath}')
kick = kick_i[0]
# Set time zero to the kick time and discard data before kick
time_t -= kick.time
flt_t = time_t > 0 - atol_time
time_t = time_t[flt_t]
dm_tv = dm_tv[flt_t]
# Check time step
dt_t = time_t[1:] - time_t[:-1]
dt = dt_t[0]
if not np.allclose(dt_t, dt, rtol=0, atol=atol_time):
raise ValueError('Varying time step')
return time_t, dm_tv, kick.strength_v
def _read_dipole_moment_file(filepath: Path_str, remove_duplicates: bool = True) \
-> Tuple[List[Kick], np.ndarray, np.ndarray, np.ndarray]:
"""Reads a time-dependent dipole moment file generated by a real-time
time-dependent density functional theory calculation.
Parameters
----------
filepath
File path.
remove_duplicates
If true, data for overlapping time instances is removed;
the data read first is kept.
Returns
-------
kick_i
List of kicks.
time_t
Time values; shape {t}.
norm_t
Norm values; shape {t}.
dm_tv
Dipole moment values; shape {t}x{3}.
Raises
------
ValueError
If kick lines cannot be read.
"""
# Search and read kicks
kick_re = re.compile(r'Kick = \['
r'(?P<k0>[-+0-9\.e\ ]+), '
r'(?P<k1>[-+0-9\.e\ ]+), '
r'(?P<k2>[-+0-9\.e\ ]+)\]')
time_re = re.compile(r'Time = (?P<time>[-+0-9\.e\ ]+)')
kick_i = []
with open(filepath, 'r') as f:
for line in f:
if line.startswith('# Kick'):
# Kick strength
m = kick_re.search(line)
if m is None:
raise ValueError('Kick line unreadable')
strength_v = np.array([float(m.group(f'k{v}')) for v in range(3)])
# Kick time
m = time_re.search(line)
if m is None:
time = 0.0
else:
time = float(m.group('time'))
kick_i.append(Kick(time, strength_v))
# Read data
time_t, data_ti = _read_td_file(filepath, remove_duplicates)
norm_t = data_ti[:, 0]
dm_tv = data_ti[:, 1:]
return kick_i, time_t, norm_t, dm_tv
def _read_td_file(filepath: Path_str, remove_duplicates: bool = True) \
-> Tuple[np.ndarray, np.ndarray]:
"""Reads a time data series from file.
The data file is assumed to be a plain-text file with the first column
being the time values and the remaining columns the data values.
Parameters
----------
filepath
File path.
remove_duplicates
If true, data for overlapping time instances is removed;
the data read first is kept.
Returns
-------
time_t
Time values; shape `{t}`.
data_ti
Data values; shape {t}x{the number of data columns}.
"""
# Read data
data_tj = np.loadtxt(filepath)
time_t = data_tj[:, 0]
data_ti = data_tj[:, 1:]
# Remove duplicates due to abruptly stopped and restarted calculation
if remove_duplicates:
flt_t = np.ones_like(time_t, dtype=bool)
maxtime = time_t[0]
for t in range(1, len(time_t)):
if time_t[t] > maxtime:
maxtime = time_t[t]
else:
flt_t[t] = False
time_t = time_t[flt_t]
data_ti = data_ti[flt_t]
return time_t, data_ti
