Coverage for calorine/gpumd/io.py: 99%
283 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-10-08 12:52 +0000
1from warnings import warn
2from collections.abc import Iterable
3from pathlib import Path
5import numpy as np
6from ase import Atoms
7from ase.io import read, write
8from ase.units import fs
9from pandas import DataFrame
11from ._header import (
12 _insert_time_column, _offset_time_column, _read_gpumd_header,
13 _read_gpumd_header_blocks, _read_gpumd_table, _require_header_key,
14 _require_single_block, _resolve_normalization, _translate_header_columns,
15)
18def read_kappa(filename: str) -> DataFrame:
19 """Parses a file in ``kappa.out`` format from GPUMD and returns the
20 content as a data frame. More information concerning file format,
21 content and units can be found `here
22 <https://gpumd.org/gpumd/output_files/kappa_out.html>`__.
24 Parameters
25 ----------
26 filename
27 Input file name.
28 """
29 data = np.loadtxt(filename, ndmin=2)
30 tags = _require_column_count(data.shape[1], 'kx_in kx_out ky_in ky_out kz_tot'.split())
31 df = DataFrame(data=data, columns=tags)
32 df['kx_tot'] = df.kx_in + df.kx_out
33 df['ky_tot'] = df.ky_in + df.ky_out
34 return df
37def _require_column_count(ncols: int, tags: list[str]) -> list[str]:
38 """Returns ``tags``, raising when the file does not hold one column per tag."""
39 if ncols != len(tags):
40 raise ValueError(
41 f'Input file contains {ncols} data columns.'
42 f' Expected {len(tags)} columns.'
43 )
44 return tags
47def _reject_reserved_column_tags(tags: list[str], filename: str) -> None:
48 """Raises when the header of ``filename`` names a column twice, or names
49 one that ``read_thermo`` derives itself."""
50 duplicates = sorted({tag for tag in tags if tags.count(tag) > 1})
51 if duplicates:
52 raise ValueError(
53 f'The header of `{filename}` names the columns'
54 f' {", ".join(duplicates)} more than once.'
55 )
56 clashes = sorted(set(tags).intersection(_THERMO_DERIVED_COLUMNS))
57 if clashes:
58 raise ValueError(
59 f'The header of `{filename}` names the columns'
60 f' {", ".join(clashes)}, which `read_thermo` derives itself.'
61 )
64# Maps the raw base column tokens GPUMD writes in the `msd.out` header (see
65# https://gpumd.org/gpumd/output_files/msd_out.html) onto the column names
66# calorine has historically used. When MSD is computed over multiple groups,
67# each base token is suffixed with `_<group index>` (e.g. `msdx_0`), which
68# `_translate_header_columns` splits off before the lookup and reattaches
69# after it.
70# The keys are GPUMD names, and the values are calorine names.
71_MSD_HEADER_BASE_COLUMN_TAGS = {
72 'time_ps': 'time',
73 'msdx': 'msd_x',
74 'msdy': 'msd_y',
75 'msdz': 'msd_z',
76 'sdcx': 'sdc_x',
77 'sdcy': 'sdc_y',
78 'sdcz': 'sdc_z',
79}
82def read_msd(filename: str) -> DataFrame:
83 """Parses a file in ``msd.out`` format from GPUMD and returns the
84 content as a data frame. More information concerning file format,
85 content and units can be found `here
86 <https://gpumd.org/gpumd/output_files/msd_out.html>`__.
88 GPUMD writes one header per ``run`` and appends to the same file.
89 Each block is a separate correlation function, so a file holding more
90 than one is rejected rather than combined.
92 Parameters
93 ----------
94 filename
95 Input file name.
96 """
97 def legacy_tags(ncols):
98 ngroups = (ncols - 1) // 6
99 if ngroups * 6 + 1 != ncols:
100 raise ValueError(
101 f'Input file contains {ncols} data columns.'
102 f' Expected {1+ngroups*6} columns (1+6*ngroups).'
103 )
104 fields = 'msd_x msd_y msd_z sdc_x sdc_y sdc_z'.split()
105 if ngroups == 1:
106 return ['time'] + fields
107 return ['time'] + [f'{field}_{g}' for g in range(ngroups) for field in fields]
109 header = _read_gpumd_header(filename, _MSD_HEADER_BASE_COLUMN_TAGS,
110 group_suffixes=True)
111 _require_single_block(header.blocks, filename, 'compute_msd')
112 return _read_gpumd_table(filename, header, legacy_tags)
115# The keys are GPUMD names, and the values are calorine names.
116_SDC_HEADER_COLUMN_TAGS = {
117 'time_ps': 'time',
118 'vacx': 'vac_x',
119 'vacy': 'vac_y',
120 'vacz': 'vac_z',
121 'sdcx': 'sdc_x',
122 'sdcy': 'sdc_y',
123 'sdcz': 'sdc_z',
124}
127def read_sdc(filename: str) -> DataFrame:
128 """Parses a file in ``sdc.out`` format from GPUMD and returns the
129 content as a data frame. More information concerning file format,
130 content and units can be found `here
131 <https://gpumd.org/gpumd/output_files/sdc_out.html>`__.
133 GPUMD writes one header per ``run`` and appends to the same file.
134 Each block is a separate correlation function, so a file holding more
135 than one is rejected rather than combined.
137 Parameters
138 ----------
139 filename
140 Input file name.
141 """
142 def legacy_tags(ncols):
143 return _require_column_count(
144 ncols, 'time vac_x vac_y vac_z sdc_x sdc_y sdc_z'.split())
146 header = _read_gpumd_header(filename, _SDC_HEADER_COLUMN_TAGS)
147 _require_single_block(header.blocks, filename, 'compute_sdc')
148 return _read_gpumd_table(filename, header, legacy_tags)
151# The keys are GPUMD names, and the values are calorine names. The nine
152# `sacf_*` and nine `visc_*` columns keep the names GPUMD gives them, and are
153# carried through by `_translate_header_columns`.
154_VISCOSITY_HEADER_COLUMN_TAGS = {'time_ps': 'time'}
157def read_viscosity(filename: str) -> DataFrame:
158 """Parses a file in ``viscosity.out`` format from GPUMD and returns the
159 content as a data frame. More information concerning file format,
160 content and units can be found `here
161 <https://gpumd.org/gpumd/output_files/viscosity_out.html>`__.
163 GPUMD writes one header per ``run`` and appends to the same file.
164 Each block is a separate correlation function, so a file holding more
165 than one is rejected rather than combined.
167 Parameters
168 ----------
169 filename
170 Input file name.
171 """
172 def legacy_tags(ncols):
173 return _require_column_count(ncols, (
174 'time sacf_xx sacf_yy sacf_zz sacf_xy sacf_xz sacf_yz sacf_yx sacf_zx sacf_zy'
175 ' visc_xx visc_yy visc_zz visc_xy visc_xz visc_yz visc_yx visc_zx visc_zy'
176 ).split())
178 header = _read_gpumd_header(filename, _VISCOSITY_HEADER_COLUMN_TAGS)
179 _require_single_block(header.blocks, filename, 'compute_viscosity')
180 return _read_gpumd_table(filename, header, legacy_tags)
183def read_hac(filename: str,
184 exclude_currents: bool = True,
185 exclude_in_out: bool = True) -> DataFrame:
186 """Parses a file in ``hac.out`` format from GPUMD and returns the
187 content as a data frame. More information concerning file format,
188 content and units can be found `here
189 <https://gpumd.org/gpumd/output_files/hac_out.html>`__.
191 Parameters
192 ----------
193 filename
194 Input file name.
195 exclude_currents
196 Do not include currents in output to save memory.
197 exclude_in_out
198 Do not include `in` and `out` parts of conductivity in output to save memory.
199 """
200 data = np.loadtxt(filename, ndmin=2)
201 tags = 'time'
202 tags += ' jin_jtot_x jout_jtot_x jin_jtot_y jout_jtot_y jtot_jtot_z'
203 tags += ' kx_in kx_out ky_in ky_out kz_tot'
204 tags = _require_column_count(data.shape[1], tags.split())
205 df = DataFrame(data=data, columns=tags)
206 df['kx_tot'] = df.kx_in + df.kx_out
207 df['ky_tot'] = df.ky_in + df.ky_out
208 df['jjx_tot'] = df.jin_jtot_x + df.jout_jtot_x
209 df['jjy_tot'] = df.jin_jtot_y + df.jout_jtot_y
210 df['jjz_tot'] = df.jtot_jtot_z
211 del df['jtot_jtot_z']
212 if exclude_in_out:
213 # remove columns with in/out data to save space
214 for col in df:
215 if 'in' in col or 'out' in col:
216 del df[col]
217 if exclude_currents:
218 # remove columns with currents to save space
219 for col in df:
220 if col.startswith('j'):
221 del df[col]
222 return df
225# The keys are GPUMD names, and the values are calorine names.
226_SHC_CORRELATION_COLUMN_TAGS = {'time_ps': 'time', 'ki': 'ki', 'ko': 'ko'}
227_SHC_SPECTRAL_COLUMN_TAGS = {'omega_THz': 'omega', 'shc_i': 'shc_i', 'shc_o': 'shc_o'}
230def _read_shc_header(filename: str) -> tuple[int, int, int, int, list[str], list[str]]:
231 """Parses the leading ``#``-prefixed header block of ``shc.out`` and
232 returns the number of correlation rows, of frequency rows, of output
233 groups and of data rows, followed by the correlation and spectral column
234 tags. Every data row has exactly three numeric columns regardless of
235 content, so only the header's row-count fields can disambiguate the
236 correlation and spectral blocks (and split repeated per-group blocks).
237 """
238 blocks, n_rows = _read_gpumd_header_blocks(filename)
239 if not blocks:
240 raise ValueError(
241 f'`{filename}` has no header. `read_shc` requires the header to'
242 ' determine how the correlation and spectral blocks are laid out;'
243 ' header-less `shc.out` files cannot be parsed unambiguously.'
244 )
245 _require_single_block(blocks, filename, 'compute_shc')
246 raw = blocks[0].lines
247 num_correlation_rows = _require_header_key(raw, 'num_correlation_rows', filename)
248 num_frequency_rows = _require_header_key(raw, 'num_frequency_rows', filename)
249 columns_correlation = _require_header_key(raw, 'columns_correlation', filename)
250 columns_shc = _require_header_key(raw, 'columns_shc', filename)
252 return (
253 int(num_correlation_rows[0]),
254 int(num_frequency_rows[0]),
255 int(raw['num_output_groups'][0]) if 'num_output_groups' in raw else 1,
256 n_rows,
257 _translate_header_columns(columns_correlation, _SHC_CORRELATION_COLUMN_TAGS),
258 _translate_header_columns(columns_shc, _SHC_SPECTRAL_COLUMN_TAGS),
259 )
262def read_shc(filename: str) -> tuple[DataFrame, DataFrame]:
263 """Parses a file in ``shc.out`` format from GPUMD (written by the
264 ``compute_shc`` keyword) and returns the correlation and spectral blocks
265 as two data frames.
267 When the calculation covers more than one output group, the
268 ``ki``/``ko``/``shc_i``/``shc_o`` columns are suffixed with the group
269 index GPUMD assigns, which runs from 1 (``ki_1``, ``ko_1``, ``ki_2``,
270 ...); the shared ``time``/``omega`` grid is not suffixed.
272 GPUMD writes one header per ``run`` and appends to the same file.
273 Each block is a separate correlation function, so a file holding more
274 than one is rejected rather than combined.
276 This reader requires the header, which alone says how the correlation and
277 spectral blocks are laid out. Its siblings fall back to reading a
278 header-less file by column count.
280 Parameters
281 ----------
282 filename
283 Input file name.
285 Returns
286 -------
287 correlation: DataFrame
288 DataFrame with columns ``time``, ``ki``, ``ko`` (suffixed with the
289 group index, e.g. ``ki_1``, when the calculation covers multiple
290 output groups).
292 shc: DataFrame
293 DataFrame with columns ``omega``, ``shc_i``, ``shc_o`` (suffixed
294 with the group index when the calculation covers multiple output
295 groups).
296 """
297 (n_corr, n_freq, num_groups, n_rows,
298 correlation_tags, shc_tags) = _read_shc_header(filename)
299 rows_per_group = n_corr + n_freq
301 time_tag, ki_tag, ko_tag = correlation_tags
302 omega_tag, shc_i_tag, shc_o_tag = shc_tags
304 suffixes = [''] if num_groups == 1 else [f'_{g + 1}' for g in range(num_groups)]
306 if n_rows == 0:
307 # `np.loadtxt` warns on a file that holds no data rows, which is what
308 # a `run` interrupted before its first dump leaves behind.
309 correlation = [time_tag] + [f'{tag}{s}' for s in suffixes for tag in (ki_tag, ko_tag)]
310 shc = [omega_tag] + [f'{tag}{s}' for s in suffixes
311 for tag in (shc_i_tag, shc_o_tag)]
312 return (DataFrame(columns=correlation, dtype=float),
313 DataFrame(columns=shc, dtype=float))
314 data = np.loadtxt(filename, ndmin=2)
315 if data.shape[1] != 3:
316 raise ValueError(
317 f'Input file contains {data.shape[1]} data columns. Expected 3 columns.'
318 )
319 expected_rows = rows_per_group * num_groups
320 if data.shape[0] != expected_rows:
321 raise ValueError(
322 f'Header declares {n_corr} correlation rows + {n_freq} frequency rows'
323 f' per group, times {num_groups} group(s) ({expected_rows} rows total),'
324 f' but `{filename}` contains {data.shape[0]} data rows.'
325 )
327 correlation = {time_tag: data[0:n_corr, 0]}
328 shc = {omega_tag: data[n_corr:rows_per_group, 0]}
329 for g, suffix in enumerate(suffixes):
330 offset = g * rows_per_group
331 corr_block = data[offset:offset + n_corr]
332 shc_block = data[offset + n_corr:offset + rows_per_group]
333 correlation[f'{ki_tag}{suffix}'] = corr_block[:, 1]
334 correlation[f'{ko_tag}{suffix}'] = corr_block[:, 2]
335 shc[f'{shc_i_tag}{suffix}'] = shc_block[:, 1]
336 shc[f'{shc_o_tag}{suffix}'] = shc_block[:, 2]
338 return DataFrame(correlation), DataFrame(shc)
341# Maps the raw column tokens GPUMD writes in the `thermo.out` header (see
342# https://gpumd.org/gpumd/output_files/thermo_out.html) onto the column
343# names calorine has historically used. `T`/`KE` and their quantum
344# thermostat counterparts `T_target`/`KE_quantum` (written when
345# `integrate.type >= 31`) both map onto the same calorine tags, since they
346# occupy the same column position and mean the same thing to a consumer.
347# The keys are GPUMD names, and the values are calorine names.
348_THERMO_HEADER_COLUMN_TAGS = {
349 'T': 'temperature',
350 'T_target': 'temperature',
351 'KE': 'kinetic_energy',
352 'KE_quantum': 'kinetic_energy',
353 'PE': 'potential_energy',
354 'sxx': 'stress_xx',
355 'syy': 'stress_yy',
356 'szz': 'stress_zz',
357 'syz': 'stress_yz',
358 'sxz': 'stress_xz',
359 'sxy': 'stress_xy',
360 'ax': 'cell_xx',
361 'ay': 'cell_xy',
362 'az': 'cell_xz',
363 'bx': 'cell_yx',
364 'by': 'cell_yy',
365 'bz': 'cell_yz',
366 'cx': 'cell_zx',
367 'cy': 'cell_zy',
368 'cz': 'cell_zz',
369}
371# The GPUMD labels of the columns that `normalize` divides by the atom count.
372# Each label is translated through _THERMO_HEADER_COLUMN_TAGS, and a column the file does not carry
373# is skipped.
374_THERMO_NORMALIZABLE_COLUMNS = ('KE', 'KE_quantum', 'PE')
377# The columns :func:`read_thermo` adds to those the file carries.
378_THERMO_DERIVED_COLUMNS = ('time', 'volume', 'pressure',
379 'cell_length_1', 'cell_length_2', 'cell_length_3',
380 'cell_angle_12', 'cell_angle_13', 'cell_angle_23')
383def read_thermo(filename: str, normalize: bool | None = None,
384 natoms: int = None) -> DataFrame:
385 """Parses a file in ``thermo.out`` format from GPUMD and returns the
386 content as a data frame. More information concerning file format,
387 content and units can be found `here
388 <https://gpumd.org/gpumd/output_files/thermo_out.html>`__. Additionally,
389 the time (in ps), the pressure (in GPa), the volume (in Å:sup:`3`), the
390 lengths of the three cell vectors (``cell_length_1`` to
391 ``cell_length_3``, in Å) and the angles between them (``cell_angle_12``,
392 ``cell_angle_13`` and ``cell_angle_23``, in degrees) are included.
393 A header that names one of those derived columns itself, or names any
394 column twice, raises.
396 GPUMD writes one header per ``run`` and appends to the same file.
397 The blocks are laid end to end, each contributing its own ``dt_output``
398 and ``num_atoms``, so ``time`` runs across the whole file and each block
399 is normalized by its own atom count.
400 The blocks have to name the same columns.
402 Parameters
403 ----------
404 filename
405 Input file name.
406 normalize
407 Divide the energies by the number of atoms, taken from the
408 ``num_atoms`` line of each header block.
409 natoms
410 Number of atoms to divide the energies by, in place of the count in
411 the header.
412 Passing it implies :attr:`normalize`, and naming it together with
413 ``normalize=False`` raises.
415 Returns
416 -------
417 DataFrame
418 DataFrame with columns determined by the ``thermo.out`` format.
419 """
420 def legacy_tags(ncols):
421 layouts = {
422 9: 'temperature kinetic_energy potential_energy'
423 ' stress_xx stress_yy stress_zz'
424 ' cell_xx cell_yy cell_zz',
425 12: 'temperature kinetic_energy potential_energy'
426 ' stress_xx stress_yy stress_zz stress_yz stress_xz stress_xy'
427 ' cell_xx cell_yy cell_zz',
428 15: 'temperature kinetic_energy potential_energy'
429 ' stress_xx stress_yy stress_zz'
430 ' cell_xx cell_xy cell_xz cell_yx cell_yy cell_yz cell_zx cell_zy cell_zz',
431 18: 'temperature kinetic_energy potential_energy'
432 ' stress_xx stress_yy stress_zz stress_yz stress_xz stress_xy'
433 ' cell_xx cell_xy cell_xz cell_yx cell_yy cell_yz cell_zx cell_zy cell_zz',
434 }
435 if ncols not in layouts:
436 raise ValueError(
437 f'Input file contains {ncols} data columns.'
438 ' Expected 9, 12, 15 or 18 columns.'
439 )
440 return layouts[ncols].split()
442 header = _read_gpumd_header(filename, _THERMO_HEADER_COLUMN_TAGS,
443 required=('num_atoms', 'dt_output'))
444 if header.tags is not None:
445 _reject_reserved_column_tags(header.tags, filename)
446 df = _read_gpumd_table(filename, header, legacy_tags)
447 if header.tags is None and len(df) == 0: 447 ↛ 448line 447 didn't jump to line 448 because the condition on line 447 was never true
448 return df
450 required = ['kinetic_energy', 'potential_energy',
451 'stress_xx', 'stress_yy', 'stress_zz', 'cell_xx', 'cell_yy', 'cell_zz']
452 missing = [tag for tag in required if tag not in df]
453 if missing:
454 raise ValueError(
455 f'`{filename}` is missing the columns {", ".join(missing)}.'
456 )
458 divisor = _resolve_normalization(
459 normalize, natoms, header.blocks, len(df), filename, 'natoms')
460 if divisor is not None:
461 normalizable_columns = dict.fromkeys(
462 _THERMO_HEADER_COLUMN_TAGS.get(label, label) for label in _THERMO_NORMALIZABLE_COLUMNS)
463 for column in normalizable_columns:
464 if column in df: 464 ↛ 463line 464 didn't jump to line 463 because the condition on line 464 was always true
465 df[column] /= divisor
467 df = _insert_time_column(df, header.blocks)
469 if 'cell_xy' in df:
470 components = ['cell_xx', 'cell_xy', 'cell_xz',
471 'cell_yx', 'cell_yy', 'cell_yz',
472 'cell_zx', 'cell_zy', 'cell_zz']
473 cell = df[components].to_numpy().reshape(-1, 3, 3)
474 else:
475 cell = np.zeros((len(df), 3, 3))
476 cell[:, [0, 1, 2], [0, 1, 2]] = df[['cell_xx', 'cell_yy', 'cell_zz']].to_numpy()
478 lengths = np.linalg.norm(cell, axis=2)
479 cosines = np.einsum('nij,nkj->nik', cell, cell) / (
480 lengths[:, :, None] * lengths[:, None, :])
481 angles = np.degrees(np.arccos(np.clip(cosines, -1, 1)))
483 df['volume'] = np.abs(np.linalg.det(cell))
484 df['pressure'] = (df.stress_xx + df.stress_yy + df.stress_zz) / 3
485 for i in range(3):
486 df[f'cell_length_{i + 1}'] = lengths[:, i]
487 for i, j in ((0, 1), (0, 2), (1, 2)):
488 df[f'cell_angle_{i + 1}{j + 1}'] = angles[:, i, j]
490 return df
493def read_xyz(filename: str) -> Atoms:
494 """
495 Reads the structure input file (``model.xyz``) for GPUMD and returns the
496 structure.
498 This is a wrapper function around :func:`ase.io.read_xyz` since the ASE implementation does
499 not read velocities properly. Specifically, the velocity unit is converted from GPUMD units
500 (Å/fs) to ASE units (1/sqrt(u/eV)).
502 Parameters
503 ----------
504 filename
505 Name of file from which to read the structure.
507 Returns
508 -------
509 Structure as ASE Atoms object with additional per-atom arrays
510 representing atomic masses, velocities etc.
511 """
512 structure = read(filename, format='extxyz')
513 if structure.has('vel'):
514 gpumd_to_ase_velocity = 1 / fs
515 structure.set_velocities(structure.get_array('vel') * gpumd_to_ase_velocity)
516 return structure
519def read_runfile(filename: str) -> list[tuple[str, list]]:
520 """
521 Parses a GPUMD input file in ``run.in`` format and returns the
522 content in the form a list of keyword-value pairs.
524 The values of a few keywords are cast to a number, among them the dump interval that
525 opens a ``dump_xyz`` line.
526 The remaining fields of such a line, meaning the file name and any flags, are kept as
527 strings, and all other values fall back to strings as well.
529 Parameters
530 ----------
531 filename
532 Input file name.
534 Returns
535 -------
536 List of keyword-value pairs.
537 """
538 data = []
539 with open(filename, 'r') as f:
540 for k, line in enumerate(f.readlines()):
541 flds = line.split()
542 if len(flds) == 0:
543 continue
544 elif len(flds) == 1:
545 raise ValueError(f'Line {k} contains only one field:\n{line}')
546 keyword = flds[0]
547 values = tuple(flds[1:])
548 if keyword in ['time_step', 'velocity']:
549 values = float(values[0])
550 elif keyword in ['dump_thermo', 'dump_restart', 'run']:
551 values = int(values[0])
552 elif keyword == 'dump_xyz':
553 # The first field is the dump interval, the remaining ones are the file name
554 # followed by optional flags such as ``precision double force``.
555 values = (int(values[0]), *values[1:])
556 elif len(values) == 1:
557 values = values[0]
558 data.append((keyword, values))
559 return data
562def write_runfile(
563 file: Path, parameters: list[tuple[str, int | float | tuple[str, float]]]
564):
565 """Write a file in run.in format to define input parameters for MD simulation.
567 Parameters
568 ----------
569 file
570 Path to file to be written.
572 parameters
573 Defines all command-parameter(s) pairs used in run.in file
574 (see GPUMD documentation for a complete list).
575 Values can be either floats, integers, or lists/tuples.
576 """
578 with open(file, 'w') as f:
579 # Write all keywords with parameter(s)
580 for key, val in parameters:
581 f.write(f'{key} ')
582 if isinstance(val, Iterable) and not isinstance(val, str):
583 for v in val:
584 f.write(f'{v} ')
585 else:
586 f.write(f'{val}')
587 f.write('\n')
590def write_xyz(filename: str, structure: Atoms, groupings: list[list[list[int]]] = None):
591 """
592 Writes a structure into GPUMD input format (`model.xyz`).
594 This is a wrapper function around :func:`ase.io.write_xyz` since the ASE implementation does
595 not write velocities properly. Specifically, the velocity unit is converted from ASE units
596 (1/sqrt(u/eV)) to GPUMD units (Å/fs).
598 Parameters
599 ----------
600 filename
601 Name of file to which the structure should be written.
602 structure
603 Input structure.
604 groupings
605 Groups into which the individual atoms should be divided in the form of
606 a list of list of lists. Specifically, the outer list corresponds to
607 the grouping methods, of which there can be three at the most, which
608 contains a list of groups in the form of lists of site indices. The
609 sum of the lengths of the latter must be the same as the total number
610 of atoms.
612 Raises
613 ------
614 ValueError
615 Raised if parameters are incompatible.
616 """
617 # Make a local copy of the atoms object
618 _structure = structure.copy()
620 # Check velocties parameter
621 velocities = _structure.get_velocities()
622 if velocities is None or np.max(np.abs(velocities)) < 1e-6:
623 has_velocity = 0
624 else:
625 has_velocity = 1
627 # Check groupings parameter
628 if groupings is None:
629 number_of_grouping_methods = 0
630 else:
631 number_of_grouping_methods = len(groupings)
632 if number_of_grouping_methods > 3:
633 raise ValueError('There can be no more than 3 grouping methods!')
634 for g, grouping in enumerate(groupings):
635 all_indices = [i for group in grouping for i in group]
636 if len(all_indices) != len(_structure) or set(all_indices) != set(
637 range(len(_structure))
638 ):
639 raise ValueError(
640 f'The indices listed in grouping method {g} are'
641 ' not compatible with the input structure!'
642 )
644 # Allowed keyword=value pairs. Use ASEs extyz write functionality.
645 # pbc="pbc_a pbc_b pbc_c"
646 # lattice="ax ay az bx by bz cx cy cz"
647 # properties=property_name:data_type:number_of_columns
648 # species:S:1
649 # pos:R:3
650 # mass:R:1
651 # vel:R:3
652 # group:I:number_of_grouping_methods
653 if _structure.has('mass'):
654 # If structure already has masses set, use those
655 warn('Structure already has array "mass"; will use existing values.')
656 else:
657 _structure.new_array('mass', _structure.get_masses())
659 if has_velocity:
660 ase_to_gpumd_velocity = fs
661 _structure.new_array('vel', _structure.get_velocities() * ase_to_gpumd_velocity)
662 if groupings is not None:
663 group_indices = np.array(
664 [
665 [
666 [
667 group_index
668 for group_index, group in enumerate(grouping)
669 if structure_idx in group
670 ]
671 for grouping in groupings
672 ]
673 for structure_idx in range(len(_structure))
674 ]
675 ).squeeze() # pythoniccc
676 _structure.new_array('group', group_indices)
678 write(filename=filename, images=_structure, write_info=True, format='extxyz')
681def read_mcmd(filename: str, accumulate: bool = True) -> DataFrame:
682 """Parses a Monte Carlo output file in ``mcmd.out`` format
683 and returns the content in the form of a DataFrame.
685 Parameters
686 ----------
687 filename
688 Path to file to be parsed.
689 accumulate
690 If ``True`` the MD steps between subsequent Monte Carlo
691 runs in the same output file will be accumulated.
693 Returns
694 -------
695 DataFrame containing acceptance ratios and concentrations (if available),
696 as well as key Monte Carlo parameters.
697 """
698 with open(filename, 'r') as f:
699 lines = f.readlines()
700 data = []
701 offset = 0
702 step = 0
703 accummulated_step = 0
704 for line in lines:
705 if line.startswith('# mc'):
706 flds = line.split()
707 mc_type = flds[2]
708 md_steps = int(flds[3])
709 mc_trials = int(flds[4])
710 temperature_initial = float(flds[5])
711 temperature_final = float(flds[6])
712 if mc_type.endswith('sgc'):
713 ntypes = int(flds[7])
714 species = [flds[8+2*k] for k in range(ntypes)]
715 phis = {f'phi_{flds[8+2*k]}': float(flds[9+2*k]) for k in range(ntypes)}
716 kappa = float(flds[8+2*ntypes]) if mc_type == 'vcsgc' else np.nan
717 elif line.startswith('# num_MD_steps'):
718 continue
719 else:
720 flds = line.split()
721 previous_step = step
722 step = int(flds[0])
723 if step <= previous_step and accumulate:
724 offset += previous_step
725 accummulated_step = step + offset
726 record = dict(
727 step=accummulated_step,
728 mc_type=mc_type,
729 md_steps=md_steps,
730 mc_trials=mc_trials,
731 temperature_initial=temperature_initial,
732 temperature_final=temperature_final,
733 acceptance_ratio=float(flds[1]),
734 )
735 if mc_type.endswith('sgc'):
736 record.update(phis)
737 if mc_type == 'vcsgc':
738 record['kappa'] = kappa
739 concentrations = {f'conc_{s}': float(flds[k])
740 for k, s in enumerate(species, start=2)}
741 record.update(concentrations)
742 data.append(record)
743 df = DataFrame.from_dict(data)
744 return df
747# The keys are GPUMD names, and the values are calorine names.
748_DPDT_HEADER_COLUMN_TAGS = {
749 'time_fs': 'time',
750 'dpdt_x': 'dPx', 'dpdt_y': 'dPy', 'dpdt_z': 'dPz',
751 'P_x': 'Px', 'P_y': 'Py', 'P_z': 'Pz',
752}
755def read_dpdt(fname: str) -> DataFrame:
756 """Read a GPUMD ``dpdt.out`` file.
758 The time column is converted from fs (as written by GPUMD) to ps.
760 GPUMD writes one header per ``run`` and appends to the same file.
761 The blocks are laid end to end, each contributing its own
762 ``dt_output``, so ``time`` runs across the whole file.
764 Parameters
765 ----------
766 fname
767 Path to the ``dpdt.out`` file.
769 Returns
770 -------
771 DataFrame
772 DataFrame with columns ``time`` (ps), ``dPx``, ``dPy``, ``dPz``
773 (time derivatives of the polarization in e·Å/fs), and ``Px``, ``Py``,
774 ``Pz`` (polarization components in e·Å).
775 """
776 header = _read_gpumd_header(fname, _DPDT_HEADER_COLUMN_TAGS)
777 if len(header.blocks) > 1:
778 # `time` comes from the file, so the spacing is needed only to lay a
779 # second block after the first.
780 for block in header.blocks:
781 _require_header_key(block.lines, 'dt_output', fname)
782 df = _read_gpumd_table(fname, header,
783 lambda ncols: 'time dPx dPy dPz Px Py Pz'.split())
784 df['time'] *= 1e-3
785 return _offset_time_column(df, header.blocks)
788# The keys are GPUMD names, and the values are calorine names.
789_DIPOLE_HEADER_COLUMN_TAGS = {
790 'step': 'step',
791 'dipole_x': 'mu_x', 'dipole_y': 'mu_y', 'dipole_z': 'mu_z',
792}
795def read_dipole(fname: str) -> DataFrame:
796 r"""Read a GPUMD ``dipole.out`` file written by the ``dump_dipole`` keyword.
798 GPUMD writes one header per ``run`` and appends to the same file.
799 The blocks are laid end to end, each contributing its own
800 ``dt_output``, so ``time`` runs across the whole file.
801 The ``step`` column is the counter GPUMD wrote, which restarts at every
802 ``run``.
804 Parameters
805 ----------
806 fname
807 Path to the ``dipole.out`` file.
809 Returns
810 -------
811 DataFrame
812 DataFrame with columns ``step`` (int), ``mu_x``, ``mu_y``, ``mu_z``
813 (dipole moment :math:`\mu` for molecules, or polarization **P** for extended systems,
814 in e·Å). If a header is present in the ``dipole.out`` file, the ``time`` in ps is
815 also included as a column.
816 """
817 header = _read_gpumd_header(fname, _DIPOLE_HEADER_COLUMN_TAGS,
818 required=('num_atoms', 'dt_output'))
819 df = _read_gpumd_table(fname, header,
820 lambda ncols: 'step mu_x mu_y mu_z'.split())
821 if len(df):
822 df['step'] = df['step'].astype(int)
823 return _insert_time_column(df, header.blocks)
826# The keys are GPUMD names, and the values are calorine names.
827_POLARIZABILITY_HEADER_COLUMN_TAGS = {
828 'step': 'step',
829 'pol_xx': 'xx', 'pol_yy': 'yy', 'pol_zz': 'zz',
830 'pol_xy': 'xy', 'pol_yz': 'yz',
831 'pol_zx': 'xz',
832 # GPUMD's header calls the `xz` component `zx`, while calorine calls it
833 # `xz`.
834}
837def read_polarizability(fname: str, normalize: bool | None = None,
838 scale: float = None) -> DataFrame:
839 r"""Read a GPUMD ``polarizability.out`` file written by ``dump_polarizability``.
841 GPUMD writes one header per ``run`` and appends to the same file.
842 The blocks are laid end to end, each contributing its own ``dt_output``
843 and ``num_atoms``, so ``time`` runs across the whole file. The ``step``
844 column is the counter GPUMD wrote, which restarts at every ``run``.
846 Parameters
847 ----------
848 fname
849 Path to the ``polarizability.out`` file.
850 normalize
851 Divide the six susceptibility columns by the number of atoms, taken
852 from the ``num_atoms`` line of each header block.
853 scale
854 Divisor applied to the six susceptibility columns in place of the
855 count in the header.
856 Passing it implies :attr:`normalize`, and naming it together with
857 ``normalize=False`` raises.
858 GPUMD writes the total supercell susceptibility
859 :math:`\chi_\mathrm{cell}`, so the divisor should match the
860 normalization constant used as the training target scale of the TNEP
861 model, typically the number of atoms.
862 Dividing by it recovers the intensive, per-atom quantity that
863 :func:`~calorine.tools.get_raman_spectrum` expects.
864 The ``step`` column is not affected.
866 Returns
867 -------
868 DataFrame
869 DataFrame with columns ``step`` (int) and the six independent components
870 ``xx``, ``yy``, ``zz``, ``xy``, ``yz``, ``xz`` of the polarizability
871 :math:`\alpha` (molecules) or susceptibility :math:`\chi` (extended
872 systems), in the same units as the TNEP training data (typically Å^3 or
873 bohr^3 per atom when :attr:`scale` equals the number of atoms).
874 The off-diagonal order follows the GPUMD ``polarizability.out`` file
875 (xy, yz, xz). If a header is present in the ``polarizability.out`` file, the
876 ``time`` in ps is also included as a column.
877 """
878 header = _read_gpumd_header(fname, _POLARIZABILITY_HEADER_COLUMN_TAGS,
879 required=('num_atoms', 'dt_output'))
880 df = _read_gpumd_table(fname, header,
881 lambda ncols: 'step xx yy zz xy yz xz'.split())
882 if len(df):
883 df['step'] = df['step'].astype(int)
884 scale = _resolve_normalization(
885 normalize, scale, header.blocks, len(df), fname, 'scale')
886 if scale is not None:
887 cols = [col for col in _POLARIZABILITY_HEADER_COLUMN_TAGS.values()
888 if col != 'step' and col in df.columns]
889 df[cols] = df[cols].to_numpy() / scale[:, None]
890 return _insert_time_column(df, header.blocks)