mdparser.py 24.7 KB
Newer Older
1
# Copyright 2015-2018 Lauri Himanen, Fawzi Mohamed, Ankit Kariryaa
2
#
3
4
5
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
6
#
7
#     http://www.apache.org/licenses/LICENSE-2.0
8
#
9
10
11
12
13
14
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

15
from __future__ import absolute_import
16
17
from builtins import next
from builtins import range
18
import numpy as np
Lauri Himanen's avatar
Lauri Himanen committed
19
20
from nomadcore.simple_parser import SimpleMatcher as SM
from nomadcore.baseclasses import MainHierarchicalParser
21
22
import nomadcore.configurationreading
import nomadcore.csvparsing
23
from .commonparser import CP2KCommonParser
Lauri Himanen's avatar
Lauri Himanen committed
24
from nomadcore.caching_backend import CachingLevel
Lauri Himanen's avatar
Lauri Himanen committed
25
from nomadcore.unit_conversion.unit_conversion import convert_unit
Lauri Himanen's avatar
Lauri Himanen committed
26
27
28
29
30
31
32
33
34
import logging
logger = logging.getLogger("nomad")


class CP2KMDParser(MainHierarchicalParser):
    """Used to parse the CP2K calculation with run types:
        -MD
        -MOLECULAR_DYNAMICS
    """
35
    def __init__(self, parser_context):
Lauri Himanen's avatar
Lauri Himanen committed
36
37
        """
        """
38
        super(CP2KMDParser, self).__init__(parser_context)
39
        self.setup_common_matcher(CP2KCommonParser(parser_context))
Lauri Himanen's avatar
Lauri Himanen committed
40
        self.traj_iterator = None
41
        self.vel_iterator = None
42
        self.energy_iterator = None
43
        self.cell_iterator = None
44
45
46
47
48
        self.n_steps = None
        self.output_freq = None
        self.coord_freq = None
        self.velo_freq = None
        self.energy_freq = None
49
        self.cell_freq = None
50
        self.md_quicksteps = []
51
        self.ensemble = None
52
        self.cp2k_ensemble = None
Lauri Himanen's avatar
Lauri Himanen committed
53
54
55
56
57
58

        #=======================================================================
        # Globally cached values

        #=======================================================================
        # Cache levels
59
        self.caching_levels.update({
60
61
62
63
64
            'x_cp2k_section_md_settings': CachingLevel.ForwardAndCache,
            'x_cp2k_section_md_step': CachingLevel.ForwardAndCache,
            'x_cp2k_section_quickstep_calculation': CachingLevel.ForwardAndCache,
            "x_cp2k_section_scf_iteration": CachingLevel.ForwardAndCache,
            "x_cp2k_section_stress_tensor": CachingLevel.ForwardAndCache,
Lauri Himanen's avatar
Lauri Himanen committed
65
66
67
68
69
70
        })

        #=======================================================================
        # SimpleMatchers
        self.md = SM(
            " MD\| Molecular Dynamics Protocol",
71
72
            forwardMatch=True,
            sections=["section_frame_sequence", "x_cp2k_section_md"],
Lauri Himanen's avatar
Lauri Himanen committed
73
            subMatchers=[
74
75
                SM( " MD\| Molecular Dynamics Protocol",
                    forwardMatch=True,
Lauri Himanen's avatar
Lauri Himanen committed
76
                    sections=["section_sampling_method", "x_cp2k_section_md_settings"],
77
                    subMatchers=[
78
79
80
81
82
83
84
85
                        SM( " MD\| Ensemble Type\s+(?P<x_cp2k_md_ensemble_type>{})".format(self.regexs.word)),
                        SM( " MD\| Number of Time Steps\s+(?P<x_cp2k_md_number_of_time_steps>{})".format(self.regexs.int)),
                        SM( " MD\| Time Step \[fs\]\s+(?P<x_cp2k_md_time_step__fs>{})".format(self.regexs.float)),
                        SM( " MD\| Temperature \[K\]\s+(?P<x_cp2k_md_target_temperature>{})".format(self.regexs.float)),
                        SM( " MD\| Temperature tolerance \[K\]\s+(?P<x_cp2k_md_target_temperature_tolerance>{})".format(self.regexs.float)),
                        SM( " MD\| Pressure \[Bar\]\s+(?P<x_cp2k_md_target_pressure>{})".format(self.regexs.float)),
                        SM( " MD\| Barostat time constant \[  fs\]\s+(?P<x_cp2k_md_barostat_time_constant>{})".format(self.regexs.float)),
                        SM( " MD\| Print MD information every\s+(?P<x_cp2k_md_print_frequency>{}) step\(s\)".format(self.regexs.int)),
86
                        SM( " MD\| File type     Print frequency\[steps\]                             File names"),
87
88
89
90
91
                        SM( " MD\| Coordinates\s+(?P<x_cp2k_md_coordinates_print_frequency>{})\s+(?P<x_cp2k_md_coordinates_filename>{})".format(self.regexs.int, self.regexs.word)),
                        SM( " MD\| Simulation Cel\s+(?P<x_cp2k_md_simulation_cell_print_frequency>{})\s+(?P<x_cp2k_md_simulation_cell_filename>{})".format(self.regexs.int, self.regexs.word)),
                        SM( " MD\| Velocities\s+(?P<x_cp2k_md_velocities_print_frequency>{})\s+(?P<x_cp2k_md_velocities_filename>{})".format(self.regexs.int, self.regexs.word)),
                        SM( " MD\| Energies\s+(?P<x_cp2k_md_energies_print_frequency>{})\s+(?P<x_cp2k_md_energies_filename>{})".format(self.regexs.int, self.regexs.word)),
                        SM( " MD\| Dump\s+(?P<x_cp2k_md_dump_print_frequency>{})\s+(?P<x_cp2k_md_dump_filename>{})".format(self.regexs.int, self.regexs.word)),
92
93
94
95
96
97
98
99
100
                    ]
                ),
                SM( " ************************** Velocities initialization **************************".replace("*", "\*"),
                    name="md_initialization_step",
                    endReStr=" INITIAL CELL ANGLS",
                    sections=["x_cp2k_section_md_step"],
                    subMatchers=[
                        self.cm.quickstep_calculation(),
                        SM( " ******************************** GO CP2K GO! **********************************".replace("*", "\*")),
101
102
103
104
105
106
107
108
                        SM( " INITIAL POTENTIAL ENERGY\[hartree\]     =\s+(?P<x_cp2k_md_potential_energy_instantaneous__hartree>{})".format(self.regexs.float)),
                        SM( " INITIAL KINETIC ENERGY\[hartree\]       =\s+(?P<x_cp2k_md_kinetic_energy_instantaneous__hartree>{})".format(self.regexs.float)),
                        SM( " INITIAL TEMPERATURE\[K\]                =\s+(?P<x_cp2k_md_temperature_instantaneous>{})".format(self.regexs.float)),
                        SM( " INITIAL BAROSTAT TEMP\[K\]              =\s+(?P<x_cp2k_md_barostat_temperature_instantaneous>{})".format(self.regexs.float)),
                        SM( " INITIAL PRESSURE\[bar\]                 =\s+(?P<x_cp2k_md_pressure_instantaneous__bar>{})".format(self.regexs.float)),
                        SM( " INITIAL VOLUME\[bohr\^3\]                =\s+(?P<x_cp2k_md_volume_instantaneous__bohr3>{})".format(self.regexs.float)),
                        SM( " INITIAL CELL LNTHS\[bohr\]              =\s+(?P<x_cp2k_md_cell_length_a_instantaneous__bohr>{0})\s+(?P<x_cp2k_md_cell_length_b_instantaneous__bohr>{0})\s+(?P<x_cp2k_md_cell_length_c_instantaneous__bohr>{0})".format(self.regexs.float)),
                        SM( " INITIAL CELL ANGLS\[deg\]               =\s+(?P<x_cp2k_md_cell_angle_a_instantaneous__deg>{0})\s+(?P<x_cp2k_md_cell_angle_b_instantaneous__deg>{0})\s+(?P<x_cp2k_md_cell_angle_c_instantaneous__deg>{0})".format(self.regexs.float)),
109
                    ],
110
                    adHoc=self.adHoc_save_md_quickstep()
111
112
113
114
                ),
                SM( " SCF WAVEFUNCTION OPTIMIZATION",
                    endReStr=" TEMPERATURE \[K\]              =",
                    name="md_step",
115
                    forwardMatch=True,
116
117
118
119
120
                    repeats=True,
                    sections=["x_cp2k_section_md_step"],
                    subMatchers=[
                        self.cm.quickstep_calculation(),
                        SM( " ENSEMBLE TYPE                ="),
121
                        SM( " STEP NUMBER                  =\s+(?P<x_cp2k_md_step_number>{})".format(self.regexs.int)),
122
123
124
125
126
127
128
129
130
131
132
133
134
135
                        SM( " TIME \[fs\]                    =\s+(?P<x_cp2k_md_time__fs>{})".format(self.regexs.float)),
                        SM( " CONSERVED QUANTITY \[hartree\] =\s+(?P<x_cp2k_md_conserved_quantity__hartree>{})".format(self.regexs.float)),
                        SM( " CPU TIME \[s\]                 =\s+(?P<x_cp2k_md_cpu_time_instantaneous>{})\s+(?P<x_cp2k_md_cpu_time_average>{})".format(self.regexs.float, self.regexs.float)),
                        SM( " ENERGY DRIFT PER ATOM \[K\]    =\s+(?P<x_cp2k_md_energy_drift_instantaneous>{})\s+(?P<x_cp2k_md_energy_drift_average>{})".format(self.regexs.float, self.regexs.float)),
                        SM( " POTENTIAL ENERGY\[hartree\]    =\s+(?P<x_cp2k_md_potential_energy_instantaneous__hartree>{})\s+(?P<x_cp2k_md_potential_energy_average__hartree>{})".format(self.regexs.float, self.regexs.float)),
                        SM( " KINETIC ENERGY \[hartree\]     =\s+(?P<x_cp2k_md_kinetic_energy_instantaneous__hartree>{})\s+(?P<x_cp2k_md_kinetic_energy_average__hartree>{})".format(self.regexs.float, self.regexs.float)),
                        SM( " TEMPERATURE \[K\]              =\s+(?P<x_cp2k_md_temperature_instantaneous>{})\s+(?P<x_cp2k_md_temperature_average>{})".format(self.regexs.float, self.regexs.float)),
                        SM( " PRESSURE \[bar\]               =\s+(?P<x_cp2k_md_pressure_instantaneous__bar>{0})\s+(?P<x_cp2k_md_pressure_average__bar>{0})".format(self.regexs.float)),
                        SM( " BAROSTAT TEMP\[K\]             =\s+(?P<x_cp2k_md_barostat_temperature_instantaneous>{0})\s+(?P<x_cp2k_md_barostat_temperature_average>{0})".format(self.regexs.float)),
                        SM( " VOLUME\[bohr\^3\]               =\s+(?P<x_cp2k_md_volume_instantaneous__bohr3>{0})\s+(?P<x_cp2k_md_volume_average__bohr3>{0})".format(self.regexs.float)),
                        SM( " CELL LNTHS\[bohr\]             =\s+(?P<x_cp2k_md_cell_length_a_instantaneous__bohr>{0})\s+(?P<x_cp2k_md_cell_length_b_instantaneous__bohr>{0})\s+(?P<x_cp2k_md_cell_length_c_instantaneous__bohr>{0})".format(self.regexs.float)),
                        SM( " AVE. CELL LNTHS\[bohr\]        =\s+(?P<x_cp2k_md_cell_length_a_average__bohr>{0})\s+(?P<x_cp2k_md_cell_length_b_average__bohr>{0})\s+(?P<x_cp2k_md_cell_length_c_average__bohr>{0})".format(self.regexs.float)),
                        SM( " CELL ANGLS\[deg\]              =\s+(?P<x_cp2k_md_cell_angle_a_instantaneous__deg>{0})\s+(?P<x_cp2k_md_cell_angle_b_instantaneous__deg>{0})\s+(?P<x_cp2k_md_cell_angle_c_instantaneous__deg>{0})".format(self.regexs.float)),
                        SM( " AVE. CELL ANGLS\[deg\]         =\s+(?P<x_cp2k_md_cell_angle_a_average__deg>{0})\s+(?P<x_cp2k_md_cell_angle_b_average__deg>{0})\s+(?P<x_cp2k_md_cell_angle_c_average__deg>{0})".format(self.regexs.float)),
136
                    ],
137
                    adHoc=self.adHoc_save_md_quickstep()
138
                ),
Lauri Himanen's avatar
Lauri Himanen committed
139
140
141
142
143
144
145
146
            ]
        )

        # Compose root matcher according to the run type. This way the
        # unnecessary regex parsers will not be compiled and searched. Saves
        # computational time.
        self.root_matcher = SM("",
            forwardMatch=True,
Lauri Himanen's avatar
Lauri Himanen committed
147
            sections=["section_run"],
Lauri Himanen's avatar
Lauri Himanen committed
148
149
150
151
152
            subMatchers=[
                SM( "",
                    forwardMatch=True,
                    sections=["section_method"],
                    subMatchers=[
153
                        self.cm.restart(),
Lauri Himanen's avatar
Lauri Himanen committed
154
155
156
157
158
                        self.cm.header(),
                        self.cm.quickstep_header(),
                    ],
                ),
                self.md,
159
                self.cm.footer(),
Lauri Himanen's avatar
Lauri Himanen committed
160
161
162
163
164
            ]
        )

    #===========================================================================
    # onClose triggers
165
    def onClose_x_cp2k_section_md_settings(self, backend, gIndex, section):
Lauri Himanen's avatar
Lauri Himanen committed
166

167
        # Ensemble
168
169
170
        cp2k_ensemble = section.get_latest_value("x_cp2k_md_ensemble_type")
        self.cp2k_ensemble = cp2k_ensemble
        if cp2k_ensemble is not None:
171
172
173
            sampling_map = {
                "NVE": "NVE",
                "NVT": "NVT",
Lauri Himanen's avatar
Lauri Himanen committed
174
175
                "NPT_F": "NPT",
                "NPT_I": "NPT",
176
            }
177
            sampling = sampling_map.get(cp2k_ensemble)
178
            if sampling is not None:
179
                self.ensemble = sampling
180
                backend.addValue("ensemble_type", sampling)
181
182
            else:
                logger.error("Unknown ensemble type {}.".format(cp2k_ensemble))
Lauri Himanen's avatar
Lauri Himanen committed
183

184
185
        # Sampling type
        backend.addValue("sampling_method", "molecular_dynamics")
Lauri Himanen's avatar
Lauri Himanen committed
186

187
188
189
190
191
        # Print frequencies
        self.output_freq = section.get_latest_value("x_cp2k_md_print_frequency")
        self.energy_freq = section.get_latest_value("x_cp2k_md_energies_print_frequency")
        self.coord_freq = section.get_latest_value("x_cp2k_md_coordinates_print_frequency")
        self.velo_freq = section.get_latest_value("x_cp2k_md_velocities_print_frequency")
192
        self.cell_freq = section.get_latest_value("x_cp2k_md_simulation_cell_print_frequency")
Lauri Himanen's avatar
Lauri Himanen committed
193

194
        # Number of steps
195
        self.n_steps = section.get_latest_value("x_cp2k_md_number_of_time_steps")
Lauri Himanen's avatar
Lauri Himanen committed
196

197
198
        # Files
        coord_filename = section.get_latest_value("x_cp2k_md_coordinates_filename")
199
        vel_filename = section.get_latest_value("x_cp2k_md_velocities_filename")
200
        energies_filename = section.get_latest_value("x_cp2k_md_energies_filename")
201
        cell_filename = section.get_latest_value("x_cp2k_md_simulation_cell_filename")
202
203
        coord_filepath = self.file_service.set_file_id(coord_filename, "coordinates")
        vel_filepath = self.file_service.set_file_id(vel_filename, "velocities")
204
        cell_filepath = self.file_service.set_file_id(cell_filename, "cell")
205
        energies_filepath = self.file_service.set_file_id(energies_filename, "energies")
Lauri Himanen's avatar
Lauri Himanen committed
206

207
        # Setup trajectory iterator
Lauri Himanen's avatar
Lauri Himanen committed
208
        traj_format = self.cache_service["trajectory_format"]
209
        if traj_format is not None and coord_filepath is not None:
Lauri Himanen's avatar
Lauri Himanen committed
210
211
212

            # Use special parsing for CP2K pdb files because they don't follow the proper syntax
            if traj_format == "PDB":
213
                self.traj_iterator = nomadcore.csvparsing.iread(coord_filepath, columns=[3, 4, 5], start="CRYST", end="END")
Lauri Himanen's avatar
Lauri Himanen committed
214
215
            else:
                try:
216
                    self.traj_iterator = nomadcore.configurationreading.iread(coord_filepath)
Lauri Himanen's avatar
Lauri Himanen committed
217
218
219
                except ValueError:
                    pass

220
221
222
223
        # Setup velocity iterator
        vel_format = self.cache_service["velocity_format"]
        if vel_format is not None and vel_filepath is not None:
            try:
224
                self.vel_iterator = nomadcore.configurationreading.iread(vel_filepath)
225
226
227
            except ValueError:
                pass

228
229
        # Setup energy file iterator
        if energies_filepath is not None:
230
            self.energy_iterator = nomadcore.csvparsing.iread(energies_filepath, columns=[0, 1, 2, 3, 4, 5, 6], comments="#")
Lauri Himanen's avatar
Lauri Himanen committed
231

232
233
        # Setup cell file iterator
        if cell_filepath is not None:
234
            self.cell_iterator = nomadcore.csvparsing.iread(cell_filepath, columns=[2, 3, 4, 5, 6, 7, 8, 9, 10], comments="#")
235

236
    def onClose_x_cp2k_section_md(self, backend, gIndex, section):
Lauri Himanen's avatar
Lauri Himanen committed
237

238
239
240
241
        # Determine the highest print frequency and use that as the number of
        # single configuration calculations
        freqs = {
            "output": [self.output_freq, True],
242
            "trajectory": [self.coord_freq, True],
243
244
            "velocities": [self.velo_freq, True],
            "energies": [self.energy_freq, True],
245
            "cell": [self.cell_freq, True],
246
        }
Lauri Himanen's avatar
Lauri Himanen committed
247

248
249
250
        # See if the files actually exist
        traj_file = self.file_service.get_file_by_id("coordinates")
        if traj_file is None:
251
            freqs["trajectory"][1] = False
252
253
254
255
256
257
        velocities_file = self.file_service.get_file_by_id("velocities")
        if velocities_file is None:
            freqs["velocities"][1] = False
        energies_file = self.file_service.get_file_by_id("energies")
        if energies_file is None:
            freqs["energies"][1] = False
258
259
260
        cell_file = self.file_service.get_file_by_id("cell")
        if cell_file is None:
            freqs["cell"][1] = False
261

262
263
264
265
266
267
268
269
        # See if we can determine the units
        traj_unit = self.cache_service["trajectory_unit"]
        if traj_unit is None:
            freqs["coordinates"][1] = False
        vel_unit = self.cache_service["velocity_unit"]
        if vel_unit is None:
            freqs["velocities"][1] = False

270
271
272
273
274
275
        # Trajectory print settings
        add_last_traj = False
        add_last_traj_setting = self.cache_service["traj_add_last"]
        if add_last_traj_setting == "NUMERIC" or add_last_traj_setting == "SYMBOLIC":
            add_last_traj = True

276
277
278
279
280
281
        # Velocities print settings
        add_last_vel = False
        add_last_vel_setting = self.cache_service["vel_add_last"]
        if add_last_vel_setting == "NUMERIC" or add_last_vel_setting == "SYMBOLIC":
            add_last_vel = True

282
        last_step = self.n_steps
283
284
285
        md_steps = section["x_cp2k_section_md_step"]

        frame_sequence_time = []
286
        ener_frames = []
287

Lauri Himanen's avatar
Lauri Himanen committed
288
        single_conf_gids = []
289
        i_md_step = 0
290

291
292
293
294
295
        # Check that is the calculation starting from first frame. If not, this
        # is a restarted calculation. In this case it is practically impossible
        # to know from which frame number we should start reading the
        # configurations, because the print settings in the previous runs may
        # have been different from now, and the step numbers are hard to align.
296
297
298
        # In this case we just parse what is found in this output file. The
        # REFTRAJ calculations use a different indexing that starts from 2
        # instead of one.
299
        start_step_number = md_steps[1]["x_cp2k_md_step_number"][0]
300
301
        if (self.cp2k_ensemble.upper() == "REFTRAJ" and start_step_number != 2) or \
           (self.cp2k_ensemble.upper() != "REFTRAJ" and start_step_number != 1):
302
303
304
305
306
            self.traj_iterator = None
            self.cell_iterator = None
            self.vel_iterator = None
            self.energy_iterator = None

307
308
309
310
311
312
313
314
        # The reftraj mode has axactly the requested number of steps, but other
        # modes have one more.
        if self.cp2k_ensemble.upper() == "REFTRAJ":
            n_steps = self.n_steps
        else:
            n_steps = self.n_steps + 1

        for i_step in range(n_steps):
315
316
317

            sectionGID = backend.openSection("section_single_configuration_calculation")
            systemGID = backend.openSection("section_system")
Lauri Himanen's avatar
Lauri Himanen committed
318
            single_conf_gids.append(sectionGID)
319

320
321
322
            # If NPT is run, and the cell file is not available, output the
            # simulation cel only on the first step to section_system
            if i_step == 0 and self.ensemble == "NPT" and self.cell_iterator is None:
323
                self.cache_service.addArrayValues("simulation_cell", unit="angstrom")
324
                self.cache_service.addArrayValues("lattice_vectors", unit="angstrom")
325

326
327
328
            # Unchanging cell
            if self.ensemble != "NPT":
                self.cache_service.addArrayValues("simulation_cell", unit="angstrom")
329
                self.cache_service.addArrayValues("lattice_vectors", unit="angstrom")
330

331
            # Trajectory
332
            if freqs["trajectory"][1] and self.traj_iterator is not None:
333
                if (i_step) % freqs["trajectory"][0] == 0 or (i_step == last_step and add_last_traj):
334
335
336
                    try:
                        pos = next(self.traj_iterator)
                    except StopIteration:
337
                        logger.error("Could not get the next geometries from an external file. It seems that the number of optimization steps in the CP2K outputfile doesn't match the number of steps found in the external trajectory file.")
338
                    else:
339
340
341
342
                        backend.addArrayValues("atom_positions", pos, unit=traj_unit)

            # Velocities
            if freqs["velocities"][1] and self.vel_iterator is not None:
343
                if (i_step) % freqs["velocities"][0] == 0 or (i_step == last_step and add_last_vel):
344
345
346
347
348
349
                    try:
                        vel = next(self.vel_iterator)
                    except StopIteration:
                        logger.error("Could not get the next velociies from an external file. It seems that the number of optimization steps in the CP2K outpufile doesn't match the number of steps found in the external velocities file.")
                    else:
                        backend.addArrayValues("atom_velocities", vel, unit=vel_unit)
350
351
352

            # Energy file
            if self.energy_iterator is not None:
353
                if (i_step) % freqs["energies"][0] == 0:
354
                    line = next(self.energy_iterator)
355
356
357
358
359
360
361
362
363

                    time = line[1]
                    kinetic_energy = line[2]
                    temperature = line[3]
                    potential_energy = line[4]
                    conserved_quantity = line[5]
                    wall_time = line[6]

                    frame_sequence_time.append(time)
Mohamed, Fawzi Roberto (fawzi)'s avatar
Mohamed, Fawzi Roberto (fawzi) committed
364
365
366
367
                    backend.addRealValue('kinetic_energy', kinetic_energy, unit="hartree")
                    backend.addRealValue('instant_temperature', temperature, unit="K")
                    backend.addRealValue('potential_energy', potential_energy, unit="hartree")
                    backend.addRealValue('conserved_quantity', conserved_quantity, unit="hartree")
368

369
370
                    ener_frames.append(i_step)

Mohamed, Fawzi Roberto (fawzi)'s avatar
Mohamed, Fawzi Roberto (fawzi) committed
371
                    backend.addRealValue("energy_total", conserved_quantity, unit="hartree")
372
373
                    backend.addValue("time_calculation", wall_time)

374
375
            # Cell file
            if self.cell_iterator is not None:
376
                if (i_step) % freqs["cell"][0] == 0:
377
                    line = next(self.cell_iterator)
378
                    cell = np.reshape(line, (3, 3))
379
                    self.backend.addArrayValues("simulation_cell", cell, unit="angstrom")
380
                    self.backend.addArrayValues("lattice_vectors", cell, unit="angstrom")
381

382
383
            # Output file
            if md_steps:
384
385
386
387
388
389
390
391
392
393
                if (i_step) % freqs["output"][0] == 0:
                    try:
                        md_step = md_steps[i_md_step]

                    except IndexError:
                        # The process has stopped unexpectedly, but still we
                        # should report all the steps so far. So we break the
                        # loop and retain what was parsed.
                        break
                    # print(md_step["x_cp2k_md_step_number"])
394
                    quickstep = self.md_quicksteps[i_md_step]
395
                    if quickstep is not None:
396
                        if quickstep.get_latest_value("x_cp2k_atom_forces") is not None:
397
398
399
400
                            # uglyness
                            fId = quickstep.caching_backend.openSection("section_atom_forces")
                            quickstep.add_latest_value("x_cp2k_atom_forces", "atom_forces")
                            quickstep.caching_backend.closeSection("section_atom_forces", fId)
401
402
403
404
                        quickstep.add_latest_value("x_cp2k_stress_tensor", "stress_tensor")
                        scfGID = backend.openSection("section_scf_iteration")
                        quickstep.add_latest_value(["x_cp2k_section_scf_iteration", "x_cp2k_energy_total_scf_iteration"], "energy_total_scf_iteration")
                        quickstep.add_latest_value(["x_cp2k_section_scf_iteration", "x_cp2k_energy_change_scf_iteration"], "energy_change_scf_iteration")
405
                        quickstep.add_latest_value(["x_cp2k_section_scf_iteration", "x_cp2k_energy_xc_scf_iteration"], "energy_xc_scf_iteration")
406
407
                        backend.closeSection("section_scf_iteration", scfGID)
                    i_md_step += 1
Lauri Himanen's avatar
Lauri Himanen committed
408
409
                    pressure = md_step.get_latest_value("x_cp2k_md_pressure_instantaneous")
                    if pressure is not None:
410
                        backend.addValue("instant_pressure", pressure)
411
412
413
414
415

            backend.closeSection("section_system", systemGID)
            backend.closeSection("section_single_configuration_calculation", sectionGID)

        # Add the summaries to frame sequence
416
417
418
        if len(frame_sequence_time) != 0:
            frame_sequence_time = np.array(frame_sequence_time)
            backend.addArrayValues("frame_sequence_time", frame_sequence_time, unit="fs")
Lauri Himanen's avatar
Lauri Himanen committed
419

420
421
422
        # Number of frames. We open a SCC for each frame, even if there is no
        # information for it.
        backend.addValue("number_of_frames_in_sequence", len(single_conf_gids))
Lauri Himanen's avatar
Lauri Himanen committed
423
424

        # Reference to sampling method
425
        backend.addValue("frame_sequence_to_sampling_method_ref", 0)
Lauri Himanen's avatar
Lauri Himanen committed
426
427

        # References to local frames
428
        backend.addArrayValues("frame_sequence_to_frames_ref", np.array(single_conf_gids))
Lauri Himanen's avatar
Lauri Himanen committed
429

430
431
432
433
434
435
436
437
438
439
    #===========================================================================
    # adHoc functions
    def adHoc_save_md_quickstep(self):
        def wrapper(parser):
            section_managers = parser.backend.sectionManagers
            section_run_manager = section_managers["section_run"]
            section_run = section_run_manager.openSections[0]
            quickstep = section_run.get_latest_value("x_cp2k_section_quickstep_calculation")
            self.md_quicksteps.append(quickstep)
        return wrapper