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