Coverage for calorine/gpumd/io.py: 100%
267 statements
« prev ^ index » next coverage.py v7.16.0, created at 2026-09-14 16:36 +0000
1from warnings import warn
2from collections.abc import Iterable
3from pathlib import Path
4from typing import List, Tuple, Union
6import numpy as np
7from ase import Atoms
8from ase.io import read, write
9from ase.units import fs
10from pandas import DataFrame
13def read_kappa(filename: str) -> DataFrame:
14 """Parses a file in ``kappa.out`` format from GPUMD and returns the
15 content as a data frame. More information concerning file format,
16 content and units can be found `here
17 <https://gpumd.org/gpumd/output_files/kappa_out.html>`__.
19 Parameters
20 ----------
21 filename
22 Input file name.
23 """
24 data = np.loadtxt(filename)
25 if isinstance(data[0], np.float64):
26 # If only a single row in kappa.out, append a dimension
27 data = data.reshape(1, -1)
28 tags = 'kx_in kx_out ky_in ky_out kz_tot'.split()
29 if len(data[0]) != len(tags):
30 raise ValueError(
31 f'Input file contains {len(data[0])} data columns.'
32 f' Expected {len(tags)} columns.'
33 )
34 df = DataFrame(data=data, columns=tags)
35 df['kx_tot'] = df.kx_in + df.kx_out
36 df['ky_tot'] = df.ky_in + df.ky_out
37 return df
40def read_msd(filename: str) -> DataFrame:
41 """Parses a file in ``msd.out`` format from GPUMD and returns the
42 content as a data frame. More information concerning file format,
43 content and units can be found `here
44 <https://gpumd.org/gpumd/output_files/msd_out.html>`__.
46 Parameters
47 ----------
48 filename
49 Input file name.
50 """
51 data = np.loadtxt(filename)
52 if isinstance(data[0], np.float64):
53 # If only a single row in msd.out, append a dimension
54 data = data.reshape(1, -1)
55 ncols = len(data[0])
56 ngroups = (ncols - 1) // 6
57 if ngroups * 6 + 1 != ncols:
58 raise ValueError(
59 f'Input file contains {ncols} data columns.'
60 f' Expected {1+ngroups*6} columns (1+6*ngroups).'
61 )
62 tags = ['time']
63 flds = 'msd_x msd_y msd_z sdc_x sdc_y sdc_z'.split()
64 if ngroups == 1:
65 tags.extend(flds)
66 else:
67 tags.extend([f'{f}_{n}' for n in range(ngroups) for f in flds])
68 df = DataFrame(data=data, columns=tags)
69 return df
72def read_hac(filename: str,
73 exclude_currents: bool = True,
74 exclude_in_out: bool = True) -> DataFrame:
75 """Parses a file in ``hac.out`` format from GPUMD and returns the
76 content as a data frame. More information concerning file format,
77 content and units can be found `here
78 <https://gpumd.org/gpumd/output_files/hac_out.html>`__.
80 Parameters
81 ----------
82 filename
83 Input file name.
84 exclude_currents
85 Do not include currents in output to save memory.
86 exclude_in_out
87 Do not include `in` and `out` parts of conductivity in output to save memory.
88 """
89 data = np.loadtxt(filename)
90 if isinstance(data[0], np.float64):
91 # If only a single row in hac.out, append a dimension
92 data = data.reshape(1, -1)
93 tags = 'time'
94 tags += ' jin_jtot_x jout_jtot_x jin_jtot_y jout_jtot_y jtot_jtot_z'
95 tags += ' kx_in kx_out ky_in ky_out kz_tot'
96 tags = tags.split()
97 if len(data[0]) != len(tags):
98 raise ValueError(
99 f'Input file contains {len(data[0])} data columns.'
100 f' Expected {len(tags)} columns.'
101 )
102 df = DataFrame(data=data, columns=tags)
103 df['kx_tot'] = df.kx_in + df.kx_out
104 df['ky_tot'] = df.ky_in + df.ky_out
105 df['jjx_tot'] = df.jin_jtot_x + df.jout_jtot_x
106 df['jjy_tot'] = df.jin_jtot_y + df.jout_jtot_y
107 df['jjz_tot'] = df.jtot_jtot_z
108 del df['jtot_jtot_z']
109 if exclude_in_out:
110 # remove columns with in/out data to save space
111 for col in df:
112 if 'in' in col or 'out' in col:
113 del df[col]
114 if exclude_currents:
115 # remove columns with currents to save space
116 for col in df:
117 if col.startswith('j'):
118 del df[col]
119 return df
122def read_thermo(filename: str, natoms: int = 1) -> DataFrame:
123 """Parses a file in ``thermo.out`` format from GPUMD and returns the
124 content as a data frame. More information concerning file format,
125 content and units can be found `here
126 <https://gpumd.org/gpumd/output_files/thermo_out.html>`__.
128 Parameters
129 ----------
130 filename
131 Input file name.
132 natoms
133 Number of atoms; used to normalize energies.
134 """
135 data = np.loadtxt(filename)
136 if len(data) == 0:
137 return DataFrame(data=data)
138 if isinstance(data[0], np.float64):
139 # If only a single row in loss.out, append a dimension
140 data = data.reshape(1, -1)
141 if len(data[0]) == 9:
142 # orthorhombic box
143 tags = 'temperature kinetic_energy potential_energy'
144 tags += ' stress_xx stress_yy stress_zz'
145 tags += ' cell_xx cell_yy cell_zz'
146 elif len(data[0]) == 12:
147 # orthorhombic box with stresses in Voigt notation (v3.3.1+)
148 tags = 'temperature kinetic_energy potential_energy'
149 tags += ' stress_xx stress_yy stress_zz stress_yz stress_xz stress_xy'
150 tags += ' cell_xx cell_yy cell_zz'
151 elif len(data[0]) == 15:
152 # triclinic box
153 tags = 'temperature kinetic_energy potential_energy'
154 tags += ' stress_xx stress_yy stress_zz'
155 tags += (
156 ' cell_xx cell_xy cell_xz cell_yx cell_yy cell_yz cell_zx cell_zy cell_zz'
157 )
158 elif len(data[0]) == 18:
159 # triclinic box with stresses in Voigt notation (v3.3.1+)
160 tags = 'temperature kinetic_energy potential_energy'
161 tags += ' stress_xx stress_yy stress_zz stress_yz stress_xz stress_xy'
162 tags += (
163 ' cell_xx cell_xy cell_xz cell_yx cell_yy cell_yz cell_zx cell_zy cell_zz'
164 )
165 else:
166 raise ValueError(
167 f'Input file contains {len(data[0])} data columns.'
168 ' Expected 9, 12, 15 or 18 columns.'
169 )
170 df = DataFrame(data=data, columns=tags.split())
171 assert natoms > 0, 'natoms must be positive'
172 df.kinetic_energy /= natoms
173 df.potential_energy /= natoms
174 return df
177def read_xyz(filename: str) -> Atoms:
178 """
179 Reads the structure input file (``model.xyz``) for GPUMD and returns the
180 structure.
182 This is a wrapper function around :func:`ase.io.read_xyz` since the ASE implementation does
183 not read velocities properly. Specifically, the velocity unit is converted from GPUMD units
184 (Å/fs) to ASE units (1/sqrt(u/eV)).
186 Parameters
187 ----------
188 filename
189 Name of file from which to read the structure.
191 Returns
192 -------
193 Structure as ASE Atoms object with additional per-atom arrays
194 representing atomic masses, velocities etc.
195 """
196 structure = read(filename, format='extxyz')
197 if structure.has('vel'):
198 gpumd_to_ase_velocity = 1 / fs
199 structure.set_velocities(structure.get_array('vel') * gpumd_to_ase_velocity)
200 return structure
203def read_runfile(filename: str) -> List[Tuple[str, list]]:
204 """
205 Parses a GPUMD input file in ``run.in`` format and returns the
206 content in the form a list of keyword-value pairs.
208 The values of a few keywords are cast to a number, among them the dump interval that
209 opens a ``dump_xyz`` line.
210 The remaining fields of such a line, meaning the file name and any flags, are kept as
211 strings, and all other values fall back to strings as well.
213 Parameters
214 ----------
215 filename
216 Input file name.
218 Returns
219 -------
220 List of keyword-value pairs.
221 """
222 data = []
223 with open(filename, 'r') as f:
224 for k, line in enumerate(f.readlines()):
225 flds = line.split()
226 if len(flds) == 0:
227 continue
228 elif len(flds) == 1:
229 raise ValueError(f'Line {k} contains only one field:\n{line}')
230 keyword = flds[0]
231 values = tuple(flds[1:])
232 if keyword in ['time_step', 'velocity']:
233 values = float(values[0])
234 elif keyword in ['dump_thermo', 'dump_restart', 'run']:
235 values = int(values[0])
236 elif keyword == 'dump_xyz':
237 # The first field is the dump interval, the remaining ones are the file name
238 # followed by optional flags such as ``precision double force``.
239 values = (int(values[0]), *values[1:])
240 elif len(values) == 1:
241 values = values[0]
242 data.append((keyword, values))
243 return data
246def write_runfile(
247 file: Path, parameters: List[Tuple[str, Union[int, float, Tuple[str, float]]]]
248):
249 """Write a file in run.in format to define input parameters for MD simulation.
251 Parameters
252 ----------
253 file
254 Path to file to be written.
256 parameters
257 Defines all command-parameter(s) pairs used in run.in file
258 (see GPUMD documentation for a complete list).
259 Values can be either floats, integers, or lists/tuples.
260 """
262 with open(file, 'w') as f:
263 # Write all keywords with parameter(s)
264 for key, val in parameters:
265 f.write(f'{key} ')
266 if isinstance(val, Iterable) and not isinstance(val, str):
267 for v in val:
268 f.write(f'{v} ')
269 else:
270 f.write(f'{val}')
271 f.write('\n')
274def write_xyz(filename: str, structure: Atoms, groupings: List[List[List[int]]] = None):
275 """
276 Writes a structure into GPUMD input format (`model.xyz`).
278 This is a wrapper function around :func:`ase.io.write_xyz` since the ASE implementation does
279 not write velocities properly. Specifically, the velocity unit is converted from ASE units
280 (1/sqrt(u/eV)) to GPUMD units (Å/fs).
282 Parameters
283 ----------
284 filename
285 Name of file to which the structure should be written.
286 structure
287 Input structure.
288 groupings
289 Groups into which the individual atoms should be divided in the form of
290 a list of list of lists. Specifically, the outer list corresponds to
291 the grouping methods, of which there can be three at the most, which
292 contains a list of groups in the form of lists of site indices. The
293 sum of the lengths of the latter must be the same as the total number
294 of atoms.
296 Raises
297 ------
298 ValueError
299 Raised if parameters are incompatible.
300 """
301 # Make a local copy of the atoms object
302 _structure = structure.copy()
304 # Check velocties parameter
305 velocities = _structure.get_velocities()
306 if velocities is None or np.max(np.abs(velocities)) < 1e-6:
307 has_velocity = 0
308 else:
309 has_velocity = 1
311 # Check groupings parameter
312 if groupings is None:
313 number_of_grouping_methods = 0
314 else:
315 number_of_grouping_methods = len(groupings)
316 if number_of_grouping_methods > 3:
317 raise ValueError('There can be no more than 3 grouping methods!')
318 for g, grouping in enumerate(groupings):
319 all_indices = [i for group in grouping for i in group]
320 if len(all_indices) != len(_structure) or set(all_indices) != set(
321 range(len(_structure))
322 ):
323 raise ValueError(
324 f'The indices listed in grouping method {g} are'
325 ' not compatible with the input structure!'
326 )
328 # Allowed keyword=value pairs. Use ASEs extyz write functionality.
329 # pbc="pbc_a pbc_b pbc_c"
330 # lattice="ax ay az bx by bz cx cy cz"
331 # properties=property_name:data_type:number_of_columns
332 # species:S:1
333 # pos:R:3
334 # mass:R:1
335 # vel:R:3
336 # group:I:number_of_grouping_methods
337 if _structure.has('mass'):
338 # If structure already has masses set, use those
339 warn('Structure already has array "mass"; will use existing values.')
340 else:
341 _structure.new_array('mass', _structure.get_masses())
343 if has_velocity:
344 ase_to_gpumd_velocity = fs
345 _structure.new_array('vel', _structure.get_velocities() * ase_to_gpumd_velocity)
346 if groupings is not None:
347 group_indices = np.array(
348 [
349 [
350 [
351 group_index
352 for group_index, group in enumerate(grouping)
353 if structure_idx in group
354 ]
355 for grouping in groupings
356 ]
357 for structure_idx in range(len(_structure))
358 ]
359 ).squeeze() # pythoniccc
360 _structure.new_array('group', group_indices)
362 write(filename=filename, images=_structure, write_info=True, format='extxyz')
365def read_mcmd(filename: str, accumulate: bool = True) -> DataFrame:
366 """Parses a Monte Carlo output file in ``mcmd.out`` format
367 and returns the content in the form of a DataFrame.
369 Parameters
370 ----------
371 filename
372 Path to file to be parsed.
373 accumulate
374 If ``True`` the MD steps between subsequent Monte Carlo
375 runs in the same output file will be accumulated.
377 Returns
378 -------
379 DataFrame containing acceptance ratios and concentrations (if available),
380 as well as key Monte Carlo parameters.
381 """
382 with open(filename, 'r') as f:
383 lines = f.readlines()
384 data = []
385 offset = 0
386 step = 0
387 accummulated_step = 0
388 for line in lines:
389 if line.startswith('# mc'):
390 flds = line.split()
391 mc_type = flds[2]
392 md_steps = int(flds[3])
393 mc_trials = int(flds[4])
394 temperature_initial = float(flds[5])
395 temperature_final = float(flds[6])
396 if mc_type.endswith('sgc'):
397 ntypes = int(flds[7])
398 species = [flds[8+2*k] for k in range(ntypes)]
399 phis = {f'phi_{flds[8+2*k]}': float(flds[9+2*k]) for k in range(ntypes)}
400 kappa = float(flds[8+2*ntypes]) if mc_type == 'vcsgc' else np.nan
401 elif line.startswith('# num_MD_steps'):
402 continue
403 else:
404 flds = line.split()
405 previous_step = step
406 step = int(flds[0])
407 if step <= previous_step and accumulate:
408 offset += previous_step
409 accummulated_step = step + offset
410 record = dict(
411 step=accummulated_step,
412 mc_type=mc_type,
413 md_steps=md_steps,
414 mc_trials=mc_trials,
415 temperature_initial=temperature_initial,
416 temperature_final=temperature_final,
417 acceptance_ratio=float(flds[1]),
418 )
419 if mc_type.endswith('sgc'):
420 record.update(phis)
421 if mc_type == 'vcsgc':
422 record['kappa'] = kappa
423 concentrations = {f'conc_{s}': float(flds[k])
424 for k, s in enumerate(species, start=2)}
425 record.update(concentrations)
426 data.append(record)
427 df = DataFrame.from_dict(data)
428 return df
431def read_thermodynamic_data(
432 directory_name: str,
433 normalize: bool = False,
434) -> DataFrame:
435 """Parses the data in a GPUMD output directory
436 and returns the content in the form of a :class:`DataFrame`.
437 This function reads the ``thermo.out``, ``run.in``, and ``model.xyz``
438 (optionally) files, and returns the thermodynamic data including the
439 time (in ps), the pressure (in GPa), the side lengths of the simulation
440 cell (in Å), and the volume (in Å:sup:`3` or Å:sup:`3`/atom).
442 Parameters
443 ----------
444 directory_name
445 Path to directory to be parsed.
446 normalize
447 Normalize thermodynamic quantities per atom.
448 This requires the ``model.xyz`` file to be present.
450 Returns
451 -------
452 :class:`DataFrame` containing (augmented) thermodynamic data.
453 """
455 try:
456 params = read_runfile(f'{directory_name}/run.in')
457 except FileNotFoundError:
458 raise FileNotFoundError(f'No `run.in` file found in {directory_name}')
460 if normalize:
461 try:
462 structure = read(f'{directory_name}/model.xyz')
463 except FileNotFoundError:
464 raise FileNotFoundError(f'No `model.xyz` file found in {directory_name}')
465 natoms = len(structure)
466 else:
467 natoms = 1
469 blocks = []
470 time_step = 1.0 # GPUMD default
471 dump_thermo = None
472 for p, v in params:
473 if p == 'time_step':
474 time_step = v
475 elif p == 'dump_thermo':
476 dump_thermo = v
477 elif p == 'run':
478 if dump_thermo is None:
479 continue
480 # We do not require dump_thermo to exist for subsequent
481 # runs if it has been used for atleast one before.
482 # But if there has been no new dump_thermo, we
483 # should not create a block.
484 if (dump_thermo != 'DEFINEDONCE'):
485 blocks.append(dict(
486 nsteps=v,
487 time_step=time_step,
488 dump_thermo=dump_thermo,
489 ))
490 dump_thermo = 'DEFINEDONCE'
492 try:
493 df = read_thermo(f'{directory_name}/thermo.out', natoms=natoms)
494 except FileNotFoundError:
495 raise FileNotFoundError(f'No `thermo.out` file found in {directory_name}')
497 expected_rows = sum([int(round(b['nsteps'] / b['dump_thermo'], 0))
498 for b in blocks if b['dump_thermo'] is not None])
499 if len(df) != expected_rows:
500 warn(f'Number of rows in `thermo.out` file ({len(df)}) does not match'
501 f' expectation based on `run.in` file ({expected_rows}).')
502 if len(df) > expected_rows:
503 raise ValueError(f'Too many rows in `thermo.out` file ({len(df)}) compared to'
504 f' expectation based on `run.in` file ({expected_rows}).')
505 if len(df) == 0:
506 # Could be the case when a run has just started and thermo.out has been created
507 # but not populated yet
508 warn('`thermo.out` is empty')
509 return df
511 # Fewer rows than expected should be ok, since the run may have crashed/not completed yet.
512 times = []
513 offset = 0.0
514 for b in blocks:
515 ns = int(round(b['nsteps'] / b['dump_thermo'], 0))
516 block_times = np.array(range(1, 1 + ns)) \
517 * b['dump_thermo'] * b['time_step'] * 1e-3 # in ps
518 block_times += offset
519 times.extend(block_times)
520 offset = times[-1]
521 df['time'] = times[:len(df)]
523 df['pressure'] = (df.stress_xx + df.stress_yy + df.stress_zz) / 3
524 if 'cell_xy' in df:
525 xx, xy, xz = df.cell_xx.to_numpy(), df.cell_xy.to_numpy(), df.cell_xz.to_numpy()
526 yx, yy, yz = df.cell_yx.to_numpy(), df.cell_yy.to_numpy(), df.cell_yz.to_numpy()
527 zx, zy, zz = df.cell_zx.to_numpy(), df.cell_zy.to_numpy(), df.cell_zz.to_numpy()
528 df['alat'] = np.sqrt(xx ** 2 + xy ** 2 + xz ** 2)
529 df['blat'] = np.sqrt(yx ** 2 + yy ** 2 + yz ** 2)
530 df['clat'] = np.sqrt(zx ** 2 + zy ** 2 + zz ** 2)
531 volume = (xx * yy * zz +
532 xy * yz * zx +
533 xz * yx * zy -
534 xx * yz * zy -
535 xy * yx * zz -
536 xz * yy * zx)
537 else:
538 df['alat'] = df.cell_xx
539 df['blat'] = df.cell_yy
540 df['clat'] = df.cell_zz
541 volume = (df.cell_xx * df.cell_yy * df.cell_zz)
542 df['volume'] = volume
543 if normalize:
544 df.volume /= natoms
546 return df
549def read_dpdt(fname: str) -> DataFrame:
550 """Read a GPUMD ``dpdt.out`` file.
552 The time column is converted from fs (as written by GPUMD) to ps.
554 Parameters
555 ----------
556 fname
557 Path to the ``dpdt.out`` file.
559 Returns
560 -------
561 DataFrame
562 DataFrame with columns ``time`` (ps), ``dPx``, ``dPy``, ``dPz``
563 (time derivatives of the polarization in e·Å/fs), and ``Px``, ``Py``,
564 ``Pz`` (polarization components in e·Å).
565 """
566 df = DataFrame(np.loadtxt(fname), columns='time dPx dPy dPz Px Py Pz'.split())
567 df['time'] *= 1e-3
568 return df
571def read_dipole(fname: str) -> DataFrame:
572 r"""Read a GPUMD ``dipole.out`` file written by the ``dump_dipole`` keyword.
574 Parameters
575 ----------
576 fname
577 Path to the ``dipole.out`` file.
579 Returns
580 -------
581 DataFrame
582 DataFrame with columns ``step`` (int), ``mu_x``, ``mu_y``, ``mu_z``
583 (dipole moment :math:`\mu` for molecules, or polarization **P** for extended systems,
584 in e·Å).
585 """
586 df = DataFrame(np.loadtxt(fname), columns='step mu_x mu_y mu_z'.split())
587 df['step'] = df['step'].astype(int)
588 return df
591def read_polarizability(fname: str, scale: float = None) -> DataFrame:
592 r"""Read a GPUMD ``polarizability.out`` file written by ``dump_polarizability``.
594 Parameters
595 ----------
596 fname
597 Path to the ``polarizability.out`` file.
598 scale
599 Divisor applied to the six susceptibility columns after reading.
600 Should match the normalisation constant used when training the TNEP
601 model. GPUMD writes the *total* supercell susceptibility
602 :math:`\chi_\mathrm{cell}`, so dividing by the same value that was used
603 as the training target scale (typically the number of atoms,
604 ``len(atoms)``) recovers the intensive, per-atom quantity expected by
605 :func:`~calorine.tools.get_raman_spectrum`. The ``step`` column is
606 not affected. ``None`` (default) leaves the data unscaled.
608 Returns
609 -------
610 DataFrame
611 DataFrame with columns ``step`` (int) and the six independent components
612 ``xx``, ``yy``, ``zz``, ``xy``, ``yz``, ``xz`` of the polarizability
613 :math:`\alpha` (molecules) or susceptibility :math:`\chi` (extended
614 systems), in the same units as the TNEP training data (typically Å^3 or
615 bohr^3 per atom when :attr:`scale` equals the number of atoms).
616 The off-diagonal order follows the GPUMD ``polarizability.out`` file
617 (xy, yz, xz).
618 """
619 df = DataFrame(np.loadtxt(fname), columns='step xx yy zz xy yz xz'.split())
620 df['step'] = df['step'].astype(int)
621 if scale is not None:
622 cols = ['xx', 'yy', 'zz', 'xy', 'yz', 'xz']
623 df[cols] = df[cols] / scale
624 return df