Coverage for sparc/sparc_parsers/aimd.py: 90%
106 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-18 16:19 +0000
1"""
2Created on Thu Oct 18 14:16:21 2018
4Ben Comer (Georgia Tech)
6This file has been heavily modified since SPARC 0.1
8TODO: more descriptions about this file io parser
9"""
10import re
11from warnings import warn
13import numpy as np
14from ase.units import AUT, Angstrom, Bohr, GPa, Hartree, fs
16# Safe wrappers for both string and fd
17from ase.utils import reader, writer
19from ..api import SparcAPI
20from .utils import strip_comments
23@reader
24def _read_aimd(fileobj):
25 """Parse the aimd information Each geopt is similar to the static
26 block, except that the field name is started by ':' The
27 relaxations are separated by ':MDSTEP:' seperators
29 """
30 contents = fileobj.read()
31 # label = get_label(fileobj, ".ion")
32 # The geopt comments are simply discarded
33 stripped, comments = strip_comments(contents)
34 # Do not include the description lines
35 data = [line for line in stripped if ":Desc" not in line]
37 # find the index for all atom type lines. They should be at the
38 # top of their block
39 step_bounds = [i for i, x in enumerate(data) if ":MDSTEP:" in x] + [len(data)]
40 raw_aimd_blocks = [
41 data[start:end] for start, end in zip(step_bounds[:-1], step_bounds[1:])
42 ]
43 aimd_steps = [_read_aimd_step(step) for step in raw_aimd_blocks]
45 return {"aimd": aimd_steps}
48def _read_aimd_step(raw_aimd_text):
49 """Parse a geopt step and compose the data dict
51 Arguments
52 raw_aimd_text: list of lines within the step
54 Most values are just presented in their output format,
55 higher level function calling _read_aimd_step and _read_aimd
56 should implement how to use the values,
57 e.g. E_tot = E_tot_per_atom * N_atoms
59 """
60 header_pattern = re.compile(r"^:\w+:")
61 header, body = raw_aimd_text[0], raw_aimd_text[1:]
62 if ":MDSTEP:" not in header:
63 raise ValueError("Wrong aimd format! The :MDSTEP: label is missing.")
64 # Geopt file uses 1-indexed step names, convert to 0-indexed
65 step = int(header.split(":MDSTEP:")[-1]) - 1
66 print("Step ", step)
67 bounds = [i for i, x in enumerate(body) if header_pattern.match(x)] + [len(body)]
68 blocks = [body[start:end] for start, end in zip(bounds[:-1], bounds[1:])]
69 data = {}
70 for block in blocks:
71 header_block, body_block = block[0], block[1:]
72 header_name = header_block.split(":")[1]
73 header_data = header_block.split(":")[-1].strip()
74 if len(header_data) > 0:
75 block_raw_data = [header_data] + body_block
76 else:
77 block_raw_data = body_block
78 # import pdb; pdb.set_trace()
79 raw_value = np.genfromtxt(block_raw_data, dtype=float)
80 # The type definitions from MD may be treated from API again?
81 if header_name == "R":
82 name = "positions"
83 value = raw_value.reshape((-1, 3)) * Bohr
84 elif header_name == "V":
85 name = "velocities"
86 value = raw_value.reshape((-1, 3)) * Bohr / AUT
87 elif header_name == "F":
88 name = "forces"
89 value = raw_value.reshape((-1, 3)) * Hartree / Bohr
90 elif header_name == "MDTM":
91 # This is not the md integration time!
92 name = "md_walltime"
93 value = float(raw_value)
94 elif header_name == "TEL":
95 name = "electron temp"
96 value = float(raw_value)
97 elif header_name == "TIO":
98 name = "ion temp"
99 value = float(raw_value)
100 elif header_name == "TEN":
101 # Note it's the total energy per atom!
102 # TODO: shall we convert to ase fashion?
103 name = "total energy per atom"
104 value = float(raw_value) * Hartree
105 elif header_name == "KEN":
106 # Note it's the total energy per atom!
107 # TODO: shall we convert to ase fashion?
108 name = "kinetic energy per atom"
109 value = float(raw_value) * Hartree
110 elif header_name == "KENIG":
111 # Note it's the total energy per atom!
112 # TODO: shall we convert to ase fashion?
113 name = "kinetic energy (ideal gas) per atom"
114 value = float(raw_value) * Hartree
115 elif header_name == "FEN":
116 # Note it's the total energy per atom!
117 # TODO: shall we convert to ase fashion?
118 name = "free energy per atom"
119 value = float(raw_value) * Hartree
120 elif header_name == "UEN":
121 # Note it's the total energy per atom!
122 # TODO: shall we convert to ase fashion?
123 name = "internal energy per atom"
124 value = float(raw_value) * Hartree
125 elif header_name == "TSEN":
126 # Note it's the total energy per atom!
127 # TODO: shall we convert to ase fashion?
128 name = "entropy*T per atom"
129 value = float(raw_value) * Hartree
130 elif header_name == "STRESS":
131 # Same rule as STRESS in geopt
132 # no conversion to Voigt form yet
133 dim = np.atleast_1d(raw_value).shape[0]
134 if dim == 3:
135 name = "stress"
136 value = raw_value * GPa
137 elif dim == 2:
138 name = "stress_2d"
139 value = raw_value * Hartree / Bohr**2
140 elif dim == 1:
141 name = "stress_1d"
142 value = float(raw_value * Hartree / Bohr)
143 else:
144 raise ValueError("Incorrect stress matrix dimension!")
145 elif header_name == "STRIO":
146 # Don't do the volume conversion now
147 name = "stress (ion-kinetic)"
148 value = raw_value * GPa
149 elif header_name == "PRES":
150 # Don't do the volume conversion now
151 name = "pressure"
152 value = raw_value * GPa
153 elif header_name == "PRESIO":
154 # Don't do the volume conversion now
155 name = "pressure (ion-kinetic)"
156 value = raw_value * GPa
157 elif header_name == "PRESIG":
158 # Don't do the volume conversion now
159 name = "pressure (ideal gas)"
160 value = raw_value * GPa
161 elif header_name in ("AVGV", "MAXV", "MIND", "TWIST"):
162 warn(f"MD output keyword {header_name} will not be parsed.")
163 value = None
164 else:
165 warn(f"MD output keyword {header_name} not known to SPARC. " "Ignore.")
166 value = None
167 if value is not None:
168 data[name] = value
169 data["step"] = step
170 return data
173@writer
174def _write_aimd(
175 fileobj,
176 data_dict,
177):
178 raise NotImplementedError("Writing aimd file from SPARC-X-API " "not supported!")