from warnings import warn
from collections.abc import Iterable
from pathlib import Path
import numpy as np
from ase import Atoms
from ase.io import read, write
from ase.units import fs
from pandas import DataFrame
from ._header import (
_insert_time_column, _offset_time_column, _read_gpumd_header,
_read_gpumd_header_blocks, _read_gpumd_table, _require_header_key,
_require_single_block, _resolve_normalization, _translate_header_columns,
)
[docs]
def read_kappa(filename: str) -> DataFrame:
"""Parses a file in ``kappa.out`` format from GPUMD and returns the
content as a data frame. More information concerning file format,
content and units can be found `here
<https://gpumd.org/gpumd/output_files/kappa_out.html>`__.
Parameters
----------
filename
Input file name.
"""
data = np.loadtxt(filename, ndmin=2)
tags = _require_column_count(data.shape[1], 'kx_in kx_out ky_in ky_out kz_tot'.split())
df = DataFrame(data=data, columns=tags)
df['kx_tot'] = df.kx_in + df.kx_out
df['ky_tot'] = df.ky_in + df.ky_out
return df
def _require_column_count(ncols: int, tags: list[str]) -> list[str]:
"""Returns ``tags``, raising when the file does not hold one column per tag."""
if ncols != len(tags):
raise ValueError(
f'Input file contains {ncols} data columns.'
f' Expected {len(tags)} columns.'
)
return tags
def _reject_reserved_column_tags(tags: list[str], filename: str) -> None:
"""Raises when the header of ``filename`` names a column twice, or names
one that ``read_thermo`` derives itself."""
duplicates = sorted({tag for tag in tags if tags.count(tag) > 1})
if duplicates:
raise ValueError(
f'The header of `{filename}` names the columns'
f' {", ".join(duplicates)} more than once.'
)
clashes = sorted(set(tags).intersection(_THERMO_DERIVED_COLUMNS))
if clashes:
raise ValueError(
f'The header of `{filename}` names the columns'
f' {", ".join(clashes)}, which `read_thermo` derives itself.'
)
# Maps the raw base column tokens GPUMD writes in the `msd.out` header (see
# https://gpumd.org/gpumd/output_files/msd_out.html) onto the column names
# calorine has historically used. When MSD is computed over multiple groups,
# each base token is suffixed with `_<group index>` (e.g. `msdx_0`), which
# `_translate_header_columns` splits off before the lookup and reattaches
# after it.
# The keys are GPUMD names, and the values are calorine names.
_MSD_HEADER_BASE_COLUMN_TAGS = {
'time_ps': 'time',
'msdx': 'msd_x',
'msdy': 'msd_y',
'msdz': 'msd_z',
'sdcx': 'sdc_x',
'sdcy': 'sdc_y',
'sdcz': 'sdc_z',
}
[docs]
def read_msd(filename: str) -> DataFrame:
"""Parses a file in ``msd.out`` format from GPUMD and returns the
content as a data frame. More information concerning file format,
content and units can be found `here
<https://gpumd.org/gpumd/output_files/msd_out.html>`__.
GPUMD writes one header per ``run`` and appends to the same file.
Each block is a separate correlation function, so a file holding more
than one is rejected rather than combined.
Parameters
----------
filename
Input file name.
"""
def legacy_tags(ncols):
ngroups = (ncols - 1) // 6
if ngroups * 6 + 1 != ncols:
raise ValueError(
f'Input file contains {ncols} data columns.'
f' Expected {1+ngroups*6} columns (1+6*ngroups).'
)
fields = 'msd_x msd_y msd_z sdc_x sdc_y sdc_z'.split()
if ngroups == 1:
return ['time'] + fields
return ['time'] + [f'{field}_{g}' for g in range(ngroups) for field in fields]
header = _read_gpumd_header(filename, _MSD_HEADER_BASE_COLUMN_TAGS)
_require_single_block(header.blocks, filename, 'compute_msd')
return _read_gpumd_table(filename, header, legacy_tags)
# The keys are GPUMD names, and the values are calorine names.
_SDC_HEADER_COLUMN_TAGS = {
'time_ps': 'time',
'vacx': 'vac_x',
'vacy': 'vac_y',
'vacz': 'vac_z',
'sdcx': 'sdc_x',
'sdcy': 'sdc_y',
'sdcz': 'sdc_z',
}
[docs]
def read_sdc(filename: str) -> DataFrame:
"""Parses a file in ``sdc.out`` format from GPUMD and returns the
content as a data frame. More information concerning file format,
content and units can be found `here
<https://gpumd.org/gpumd/output_files/sdc_out.html>`__.
GPUMD writes one header per ``run`` and appends to the same file.
Each block is a separate correlation function, so a file holding more
than one is rejected rather than combined.
Parameters
----------
filename
Input file name.
"""
def legacy_tags(ncols):
return _require_column_count(
ncols, 'time vac_x vac_y vac_z sdc_x sdc_y sdc_z'.split())
header = _read_gpumd_header(filename, _SDC_HEADER_COLUMN_TAGS)
_require_single_block(header.blocks, filename, 'compute_sdc')
return _read_gpumd_table(filename, header, legacy_tags)
# The keys are GPUMD names, and the values are calorine names. The nine
# `sacf_*` and nine `visc_*` columns keep the names GPUMD gives them, and are
# carried through by `_translate_header_columns`.
_VISCOSITY_HEADER_COLUMN_TAGS = {'time_ps': 'time'}
[docs]
def read_viscosity(filename: str) -> DataFrame:
"""Parses a file in ``viscosity.out`` format from GPUMD and returns the
content as a data frame. More information concerning file format,
content and units can be found `here
<https://gpumd.org/gpumd/output_files/viscosity_out.html>`__.
GPUMD writes one header per ``run`` and appends to the same file.
Each block is a separate correlation function, so a file holding more
than one is rejected rather than combined.
Parameters
----------
filename
Input file name.
"""
def legacy_tags(ncols):
return _require_column_count(ncols, (
'time sacf_xx sacf_yy sacf_zz sacf_xy sacf_xz sacf_yz sacf_yx sacf_zx sacf_zy'
' visc_xx visc_yy visc_zz visc_xy visc_xz visc_yz visc_yx visc_zx visc_zy'
).split())
header = _read_gpumd_header(filename, _VISCOSITY_HEADER_COLUMN_TAGS)
_require_single_block(header.blocks, filename, 'compute_viscosity')
return _read_gpumd_table(filename, header, legacy_tags)
[docs]
def read_hac(filename: str,
exclude_currents: bool = True,
exclude_in_out: bool = True) -> DataFrame:
"""Parses a file in ``hac.out`` format from GPUMD and returns the
content as a data frame. More information concerning file format,
content and units can be found `here
<https://gpumd.org/gpumd/output_files/hac_out.html>`__.
Parameters
----------
filename
Input file name.
exclude_currents
Do not include currents in output to save memory.
exclude_in_out
Do not include `in` and `out` parts of conductivity in output to save memory.
"""
data = np.loadtxt(filename, ndmin=2)
tags = 'time'
tags += ' jin_jtot_x jout_jtot_x jin_jtot_y jout_jtot_y jtot_jtot_z'
tags += ' kx_in kx_out ky_in ky_out kz_tot'
tags = _require_column_count(data.shape[1], tags.split())
df = DataFrame(data=data, columns=tags)
df['kx_tot'] = df.kx_in + df.kx_out
df['ky_tot'] = df.ky_in + df.ky_out
df['jjx_tot'] = df.jin_jtot_x + df.jout_jtot_x
df['jjy_tot'] = df.jin_jtot_y + df.jout_jtot_y
df['jjz_tot'] = df.jtot_jtot_z
del df['jtot_jtot_z']
if exclude_in_out:
# remove columns with in/out data to save space
for col in df:
if 'in' in col or 'out' in col:
del df[col]
if exclude_currents:
# remove columns with currents to save space
for col in df:
if col.startswith('j'):
del df[col]
return df
# The keys are GPUMD names, and the values are calorine names.
_SHC_CORRELATION_COLUMN_TAGS = {'time_ps': 'time', 'ki': 'ki', 'ko': 'ko'}
_SHC_SPECTRAL_COLUMN_TAGS = {'omega_THz': 'omega', 'shc_i': 'shc_i', 'shc_o': 'shc_o'}
def _read_shc_header(filename: str) -> tuple[int, int, int, int, list[str], list[str]]:
"""Parses the leading ``#``-prefixed header block of ``shc.out`` and
returns the number of correlation rows, of frequency rows, of output
groups and of data rows, followed by the correlation and spectral column
tags. Every data row has exactly three numeric columns regardless of
content, so only the header's row-count fields can disambiguate the
correlation and spectral blocks (and split repeated per-group blocks).
"""
blocks, n_rows = _read_gpumd_header_blocks(filename)
if not blocks:
raise ValueError(
f'`{filename}` has no header. `read_shc` requires the header to'
' determine how the correlation and spectral blocks are laid out;'
' header-less `shc.out` files cannot be parsed unambiguously.'
)
_require_single_block(blocks, filename, 'compute_shc')
raw = blocks[0].lines
num_correlation_rows = _require_header_key(raw, 'num_correlation_rows', filename)
num_frequency_rows = _require_header_key(raw, 'num_frequency_rows', filename)
columns_correlation = _require_header_key(raw, 'columns_correlation', filename)
columns_shc = _require_header_key(raw, 'columns_shc', filename)
return (
int(num_correlation_rows[0]),
int(num_frequency_rows[0]),
int(raw['num_output_groups'][0]) if 'num_output_groups' in raw else 1,
n_rows,
_translate_header_columns(columns_correlation, _SHC_CORRELATION_COLUMN_TAGS),
_translate_header_columns(columns_shc, _SHC_SPECTRAL_COLUMN_TAGS),
)
[docs]
def read_shc(filename: str) -> tuple[DataFrame, DataFrame]:
"""Parses a file in ``shc.out`` format from GPUMD (written by the
``compute_shc`` keyword) and returns the correlation and spectral blocks
as two data frames.
When the calculation covers more than one output group, the
``ki``/``ko``/``shc_i``/``shc_o`` columns are suffixed with the group
index GPUMD assigns, which runs from 1 (``ki_1``, ``ko_1``, ``ki_2``,
...); the shared ``time``/``omega`` grid is not suffixed.
GPUMD writes one header per ``run`` and appends to the same file.
Each block is a separate correlation function, so a file holding more
than one is rejected rather than combined.
This reader requires the header, which alone says how the correlation and
spectral blocks are laid out. Its siblings fall back to reading a
header-less file by column count.
Parameters
----------
filename
Input file name.
Returns
-------
correlation: DataFrame
DataFrame with columns ``time``, ``ki``, ``ko`` (suffixed with the
group index, e.g. ``ki_1``, when the calculation covers multiple
output groups).
shc: DataFrame
DataFrame with columns ``omega``, ``shc_i``, ``shc_o`` (suffixed
with the group index when the calculation covers multiple output
groups).
"""
(n_corr, n_freq, num_groups, n_rows,
correlation_tags, shc_tags) = _read_shc_header(filename)
rows_per_group = n_corr + n_freq
time_tag, ki_tag, ko_tag = correlation_tags
omega_tag, shc_i_tag, shc_o_tag = shc_tags
suffixes = [''] if num_groups == 1 else [f'_{g + 1}' for g in range(num_groups)]
if n_rows == 0:
# `np.loadtxt` warns on a file that holds no data rows, which is what
# a `run` interrupted before its first dump leaves behind.
correlation = [time_tag] + [f'{tag}{s}' for s in suffixes for tag in (ki_tag, ko_tag)]
shc = [omega_tag] + [f'{tag}{s}' for s in suffixes
for tag in (shc_i_tag, shc_o_tag)]
return (DataFrame(columns=correlation, dtype=float),
DataFrame(columns=shc, dtype=float))
data = np.loadtxt(filename, ndmin=2)
if data.shape[1] != 3:
raise ValueError(
f'Input file contains {data.shape[1]} data columns. Expected 3 columns.'
)
expected_rows = rows_per_group * num_groups
if data.shape[0] != expected_rows:
raise ValueError(
f'Header declares {n_corr} correlation rows + {n_freq} frequency rows'
f' per group, times {num_groups} group(s) ({expected_rows} rows total),'
f' but `{filename}` contains {data.shape[0]} data rows.'
)
correlation = {time_tag: data[0:n_corr, 0]}
shc = {omega_tag: data[n_corr:rows_per_group, 0]}
for g, suffix in enumerate(suffixes):
offset = g * rows_per_group
corr_block = data[offset:offset + n_corr]
shc_block = data[offset + n_corr:offset + rows_per_group]
correlation[f'{ki_tag}{suffix}'] = corr_block[:, 1]
correlation[f'{ko_tag}{suffix}'] = corr_block[:, 2]
shc[f'{shc_i_tag}{suffix}'] = shc_block[:, 1]
shc[f'{shc_o_tag}{suffix}'] = shc_block[:, 2]
return DataFrame(correlation), DataFrame(shc)
# Maps the raw column tokens GPUMD writes in the `thermo.out` header (see
# https://gpumd.org/gpumd/output_files/thermo_out.html) onto the column
# names calorine has historically used. `T`/`KE` and their quantum
# thermostat counterparts `T_target`/`KE_quantum` (written when
# `integrate.type >= 31`) both map onto the same calorine tags, since they
# occupy the same column position and mean the same thing to a consumer.
# The keys are GPUMD names, and the values are calorine names.
_THERMO_HEADER_COLUMN_TAGS = {
'T': 'temperature',
'T_target': 'temperature',
'KE': 'kinetic_energy',
'KE_quantum': 'kinetic_energy',
'PE': 'potential_energy',
'sxx': 'stress_xx',
'syy': 'stress_yy',
'szz': 'stress_zz',
'syz': 'stress_yz',
'sxz': 'stress_xz',
'sxy': 'stress_xy',
'ax': 'cell_xx',
'ay': 'cell_xy',
'az': 'cell_xz',
'bx': 'cell_yx',
'by': 'cell_yy',
'bz': 'cell_yz',
'cx': 'cell_zx',
'cy': 'cell_zy',
'cz': 'cell_zz',
}
# The columns :func:`read_thermo` adds to those the file carries.
_THERMO_DERIVED_COLUMNS = ('time', 'volume', 'pressure',
'cell_length_1', 'cell_length_2', 'cell_length_3',
'cell_angle_12', 'cell_angle_13', 'cell_angle_23')
[docs]
def read_thermo(filename: str, normalize: bool | None = None,
natoms: int = None) -> DataFrame:
"""Parses a file in ``thermo.out`` format from GPUMD and returns the
content as a data frame. More information concerning file format,
content and units can be found `here
<https://gpumd.org/gpumd/output_files/thermo_out.html>`__. Additionally,
the time (in ps), the pressure (in GPa), the volume (in Å:sup:`3`), the
lengths of the three cell vectors (``cell_length_1`` to
``cell_length_3``, in Å) and the angles between them (``cell_angle_12``,
``cell_angle_13`` and ``cell_angle_23``, in degrees) are included.
A header that names one of those derived columns itself, or names any
column twice, raises.
GPUMD writes one header per ``run`` and appends to the same file.
The blocks are laid end to end, each contributing its own ``dt_output``
and ``num_atoms``, so ``time`` runs across the whole file and each block
is normalized by its own atom count.
Parameters
----------
filename
Input file name.
normalize
Divide the energies by the number of atoms, taken from the
``num_atoms`` line of each header block.
natoms
Number of atoms to divide the energies by, in place of the count in
the header.
Passing it implies :attr:`normalize`, and naming it together with
``normalize=False`` raises.
Returns
-------
DataFrame
DataFrame with columns determined by the ``thermo.out`` format.
"""
def legacy_tags(ncols):
layouts = {
9: 'temperature kinetic_energy potential_energy'
' stress_xx stress_yy stress_zz'
' cell_xx cell_yy cell_zz',
12: 'temperature kinetic_energy potential_energy'
' stress_xx stress_yy stress_zz stress_yz stress_xz stress_xy'
' cell_xx cell_yy cell_zz',
15: 'temperature kinetic_energy potential_energy'
' stress_xx stress_yy stress_zz'
' cell_xx cell_xy cell_xz cell_yx cell_yy cell_yz cell_zx cell_zy cell_zz',
18: 'temperature kinetic_energy potential_energy'
' stress_xx stress_yy stress_zz stress_yz stress_xz stress_xy'
' cell_xx cell_xy cell_xz cell_yx cell_yy cell_yz cell_zx cell_zy cell_zz',
}
if ncols not in layouts:
raise ValueError(
f'Input file contains {ncols} data columns.'
' Expected 9, 12, 15 or 18 columns.'
)
return layouts[ncols].split()
header = _read_gpumd_header(filename, _THERMO_HEADER_COLUMN_TAGS,
required=('num_atoms', 'dt_output'))
if header.tags is not None:
_reject_reserved_column_tags(header.tags, filename)
df = _read_gpumd_table(filename, header, legacy_tags)
if header.tags is None and len(df) == 0:
return df
required = ['kinetic_energy', 'potential_energy',
'stress_xx', 'stress_yy', 'stress_zz', 'cell_xx', 'cell_yy', 'cell_zz']
missing = [tag for tag in required if tag not in df]
if missing:
raise ValueError(
f'`{filename}` is missing the columns {", ".join(missing)}.'
)
divisor = _resolve_normalization(
normalize, natoms, header.blocks, len(df), filename, 'natoms')
if divisor is not None:
df.kinetic_energy /= divisor
df.potential_energy /= divisor
df = _insert_time_column(df, header.blocks)
if 'cell_xy' in df:
components = ['cell_xx', 'cell_xy', 'cell_xz',
'cell_yx', 'cell_yy', 'cell_yz',
'cell_zx', 'cell_zy', 'cell_zz']
cell = df[components].to_numpy().reshape(-1, 3, 3)
else:
cell = np.zeros((len(df), 3, 3))
cell[:, [0, 1, 2], [0, 1, 2]] = df[['cell_xx', 'cell_yy', 'cell_zz']].to_numpy()
lengths = np.linalg.norm(cell, axis=2)
cosines = np.einsum('nij,nkj->nik', cell, cell) / (
lengths[:, :, None] * lengths[:, None, :])
angles = np.degrees(np.arccos(np.clip(cosines, -1, 1)))
df['volume'] = np.abs(np.linalg.det(cell))
df['pressure'] = (df.stress_xx + df.stress_yy + df.stress_zz) / 3
for i in range(3):
df[f'cell_length_{i + 1}'] = lengths[:, i]
for i, j in ((0, 1), (0, 2), (1, 2)):
df[f'cell_angle_{i + 1}{j + 1}'] = angles[:, i, j]
return df
[docs]
def read_xyz(filename: str) -> Atoms:
"""
Reads the structure input file (``model.xyz``) for GPUMD and returns the
structure.
This is a wrapper function around :func:`ase.io.read_xyz` since the ASE implementation does
not read velocities properly. Specifically, the velocity unit is converted from GPUMD units
(Å/fs) to ASE units (1/sqrt(u/eV)).
Parameters
----------
filename
Name of file from which to read the structure.
Returns
-------
Structure as ASE Atoms object with additional per-atom arrays
representing atomic masses, velocities etc.
"""
structure = read(filename, format='extxyz')
if structure.has('vel'):
gpumd_to_ase_velocity = 1 / fs
structure.set_velocities(structure.get_array('vel') * gpumd_to_ase_velocity)
return structure
[docs]
def read_runfile(filename: str) -> list[tuple[str, list]]:
"""
Parses a GPUMD input file in ``run.in`` format and returns the
content in the form a list of keyword-value pairs.
The values of a few keywords are cast to a number, among them the dump interval that
opens a ``dump_xyz`` line.
The remaining fields of such a line, meaning the file name and any flags, are kept as
strings, and all other values fall back to strings as well.
Parameters
----------
filename
Input file name.
Returns
-------
List of keyword-value pairs.
"""
data = []
with open(filename, 'r') as f:
for k, line in enumerate(f.readlines()):
flds = line.split()
if len(flds) == 0:
continue
elif len(flds) == 1:
raise ValueError(f'Line {k} contains only one field:\n{line}')
keyword = flds[0]
values = tuple(flds[1:])
if keyword in ['time_step', 'velocity']:
values = float(values[0])
elif keyword in ['dump_thermo', 'dump_restart', 'run']:
values = int(values[0])
elif keyword == 'dump_xyz':
# The first field is the dump interval, the remaining ones are the file name
# followed by optional flags such as ``precision double force``.
values = (int(values[0]), *values[1:])
elif len(values) == 1:
values = values[0]
data.append((keyword, values))
return data
[docs]
def write_runfile(
file: Path, parameters: list[tuple[str, int | float | tuple[str, float]]]
):
"""Write a file in run.in format to define input parameters for MD simulation.
Parameters
----------
file
Path to file to be written.
parameters
Defines all command-parameter(s) pairs used in run.in file
(see GPUMD documentation for a complete list).
Values can be either floats, integers, or lists/tuples.
"""
with open(file, 'w') as f:
# Write all keywords with parameter(s)
for key, val in parameters:
f.write(f'{key} ')
if isinstance(val, Iterable) and not isinstance(val, str):
for v in val:
f.write(f'{v} ')
else:
f.write(f'{val}')
f.write('\n')
[docs]
def write_xyz(filename: str, structure: Atoms, groupings: list[list[list[int]]] = None):
"""
Writes a structure into GPUMD input format (`model.xyz`).
This is a wrapper function around :func:`ase.io.write_xyz` since the ASE implementation does
not write velocities properly. Specifically, the velocity unit is converted from ASE units
(1/sqrt(u/eV)) to GPUMD units (Å/fs).
Parameters
----------
filename
Name of file to which the structure should be written.
structure
Input structure.
groupings
Groups into which the individual atoms should be divided in the form of
a list of list of lists. Specifically, the outer list corresponds to
the grouping methods, of which there can be three at the most, which
contains a list of groups in the form of lists of site indices. The
sum of the lengths of the latter must be the same as the total number
of atoms.
Raises
------
ValueError
Raised if parameters are incompatible.
"""
# Make a local copy of the atoms object
_structure = structure.copy()
# Check velocties parameter
velocities = _structure.get_velocities()
if velocities is None or np.max(np.abs(velocities)) < 1e-6:
has_velocity = 0
else:
has_velocity = 1
# Check groupings parameter
if groupings is None:
number_of_grouping_methods = 0
else:
number_of_grouping_methods = len(groupings)
if number_of_grouping_methods > 3:
raise ValueError('There can be no more than 3 grouping methods!')
for g, grouping in enumerate(groupings):
all_indices = [i for group in grouping for i in group]
if len(all_indices) != len(_structure) or set(all_indices) != set(
range(len(_structure))
):
raise ValueError(
f'The indices listed in grouping method {g} are'
' not compatible with the input structure!'
)
# Allowed keyword=value pairs. Use ASEs extyz write functionality.
# pbc="pbc_a pbc_b pbc_c"
# lattice="ax ay az bx by bz cx cy cz"
# properties=property_name:data_type:number_of_columns
# species:S:1
# pos:R:3
# mass:R:1
# vel:R:3
# group:I:number_of_grouping_methods
if _structure.has('mass'):
# If structure already has masses set, use those
warn('Structure already has array "mass"; will use existing values.')
else:
_structure.new_array('mass', _structure.get_masses())
if has_velocity:
ase_to_gpumd_velocity = fs
_structure.new_array('vel', _structure.get_velocities() * ase_to_gpumd_velocity)
if groupings is not None:
group_indices = np.array(
[
[
[
group_index
for group_index, group in enumerate(grouping)
if structure_idx in group
]
for grouping in groupings
]
for structure_idx in range(len(_structure))
]
).squeeze() # pythoniccc
_structure.new_array('group', group_indices)
write(filename=filename, images=_structure, write_info=True, format='extxyz')
[docs]
def read_mcmd(filename: str, accumulate: bool = True) -> DataFrame:
"""Parses a Monte Carlo output file in ``mcmd.out`` format
and returns the content in the form of a DataFrame.
Parameters
----------
filename
Path to file to be parsed.
accumulate
If ``True`` the MD steps between subsequent Monte Carlo
runs in the same output file will be accumulated.
Returns
-------
DataFrame containing acceptance ratios and concentrations (if available),
as well as key Monte Carlo parameters.
"""
with open(filename, 'r') as f:
lines = f.readlines()
data = []
offset = 0
step = 0
accummulated_step = 0
for line in lines:
if line.startswith('# mc'):
flds = line.split()
mc_type = flds[2]
md_steps = int(flds[3])
mc_trials = int(flds[4])
temperature_initial = float(flds[5])
temperature_final = float(flds[6])
if mc_type.endswith('sgc'):
ntypes = int(flds[7])
species = [flds[8+2*k] for k in range(ntypes)]
phis = {f'phi_{flds[8+2*k]}': float(flds[9+2*k]) for k in range(ntypes)}
kappa = float(flds[8+2*ntypes]) if mc_type == 'vcsgc' else np.nan
elif line.startswith('# num_MD_steps'):
continue
else:
flds = line.split()
previous_step = step
step = int(flds[0])
if step <= previous_step and accumulate:
offset += previous_step
accummulated_step = step + offset
record = dict(
step=accummulated_step,
mc_type=mc_type,
md_steps=md_steps,
mc_trials=mc_trials,
temperature_initial=temperature_initial,
temperature_final=temperature_final,
acceptance_ratio=float(flds[1]),
)
if mc_type.endswith('sgc'):
record.update(phis)
if mc_type == 'vcsgc':
record['kappa'] = kappa
concentrations = {f'conc_{s}': float(flds[k])
for k, s in enumerate(species, start=2)}
record.update(concentrations)
data.append(record)
df = DataFrame.from_dict(data)
return df
# The keys are GPUMD names, and the values are calorine names.
_DPDT_HEADER_COLUMN_TAGS = {
'time_fs': 'time',
'dpdt_x': 'dPx', 'dpdt_y': 'dPy', 'dpdt_z': 'dPz',
'P_x': 'Px', 'P_y': 'Py', 'P_z': 'Pz',
}
[docs]
def read_dpdt(fname: str) -> DataFrame:
"""Read a GPUMD ``dpdt.out`` file.
The time column is converted from fs (as written by GPUMD) to ps.
GPUMD writes one header per ``run`` and appends to the same file.
The blocks are laid end to end, each contributing its own
``dt_output``, so ``time`` runs across the whole file.
Parameters
----------
fname
Path to the ``dpdt.out`` file.
Returns
-------
DataFrame
DataFrame with columns ``time`` (ps), ``dPx``, ``dPy``, ``dPz``
(time derivatives of the polarization in e·Å/fs), and ``Px``, ``Py``,
``Pz`` (polarization components in e·Å).
"""
header = _read_gpumd_header(fname, _DPDT_HEADER_COLUMN_TAGS)
if len(header.blocks) > 1:
# `time` comes from the file, so the spacing is needed only to lay a
# second block after the first.
for block in header.blocks:
_require_header_key(block.lines, 'dt_output', fname)
df = _read_gpumd_table(fname, header,
lambda ncols: 'time dPx dPy dPz Px Py Pz'.split())
df['time'] *= 1e-3
return _offset_time_column(df, header.blocks)
# The keys are GPUMD names, and the values are calorine names.
_DIPOLE_HEADER_COLUMN_TAGS = {
'step': 'step',
'dipole_x': 'mu_x', 'dipole_y': 'mu_y', 'dipole_z': 'mu_z',
}
[docs]
def read_dipole(fname: str) -> DataFrame:
r"""Read a GPUMD ``dipole.out`` file written by the ``dump_dipole`` keyword.
GPUMD writes one header per ``run`` and appends to the same file.
The blocks are laid end to end, each contributing its own
``dt_output``, so ``time`` runs across the whole file.
The ``step`` column is the counter GPUMD wrote, which restarts at every
``run``.
Parameters
----------
fname
Path to the ``dipole.out`` file.
Returns
-------
DataFrame
DataFrame with columns ``step`` (int), ``mu_x``, ``mu_y``, ``mu_z``
(dipole moment :math:`\mu` for molecules, or polarization **P** for extended systems,
in e·Å). If a header is present in the ``dipole.out`` file, the ``time`` in ps is
also included as a column.
"""
header = _read_gpumd_header(fname, _DIPOLE_HEADER_COLUMN_TAGS,
required=('num_atoms', 'dt_output'))
df = _read_gpumd_table(fname, header,
lambda ncols: 'step mu_x mu_y mu_z'.split())
if len(df):
df['step'] = df['step'].astype(int)
return _insert_time_column(df, header.blocks)
# The keys are GPUMD names, and the values are calorine names.
_POLARIZABILITY_HEADER_COLUMN_TAGS = {
'step': 'step',
'pol_xx': 'xx', 'pol_yy': 'yy', 'pol_zz': 'zz',
'pol_xy': 'xy', 'pol_yz': 'yz',
'pol_zx': 'xz',
# GPUMD's header calls the `xz` component `zx`, while calorine calls it
# `xz`.
}
[docs]
def read_polarizability(fname: str, normalize: bool | None = None,
scale: float = None) -> DataFrame:
r"""Read a GPUMD ``polarizability.out`` file written by ``dump_polarizability``.
GPUMD writes one header per ``run`` and appends to the same file.
The blocks are laid end to end, each contributing its own ``dt_output``
and ``num_atoms``, so ``time`` runs across the whole file. The ``step``
column is the counter GPUMD wrote, which restarts at every ``run``.
Parameters
----------
fname
Path to the ``polarizability.out`` file.
normalize
Divide the six susceptibility columns by the number of atoms, taken
from the ``num_atoms`` line of each header block.
scale
Divisor applied to the six susceptibility columns in place of the
count in the header.
Passing it implies :attr:`normalize`, and naming it together with
``normalize=False`` raises.
GPUMD writes the total supercell susceptibility
:math:`\chi_\mathrm{cell}`, so the divisor should match the
normalization constant used as the training target scale of the TNEP
model, typically the number of atoms.
Dividing by it recovers the intensive, per-atom quantity that
:func:`~calorine.tools.get_raman_spectrum` expects.
The ``step`` column is not affected.
Returns
-------
DataFrame
DataFrame with columns ``step`` (int) and the six independent components
``xx``, ``yy``, ``zz``, ``xy``, ``yz``, ``xz`` of the polarizability
:math:`\alpha` (molecules) or susceptibility :math:`\chi` (extended
systems), in the same units as the TNEP training data (typically Å^3 or
bohr^3 per atom when :attr:`scale` equals the number of atoms).
The off-diagonal order follows the GPUMD ``polarizability.out`` file
(xy, yz, xz). If a header is present in the ``polarizability.out`` file, the
``time`` in ps is also included as a column.
"""
header = _read_gpumd_header(fname, _POLARIZABILITY_HEADER_COLUMN_TAGS,
required=('num_atoms', 'dt_output'))
df = _read_gpumd_table(fname, header,
lambda ncols: 'step xx yy zz xy yz xz'.split())
if len(df):
df['step'] = df['step'].astype(int)
scale = _resolve_normalization(
normalize, scale, header.blocks, len(df), fname, 'scale')
if scale is not None:
cols = [col for col in _POLARIZABILITY_HEADER_COLUMN_TAGS.values()
if col != 'step' and col in df.columns]
df[cols] = df[cols].to_numpy() / scale[:, None]
return _insert_time_column(df, header.blocks)
