mdparser.py 24 KB
Newer Older
1
from __future__ import absolute_import
2
3
from builtins import next
from builtins import range
4
import numpy as np
Lauri Himanen's avatar
Lauri Himanen committed
5
6
from nomadcore.simple_parser import SimpleMatcher as SM
from nomadcore.baseclasses import MainHierarchicalParser
7
from .commonparser import CP2KCommonParser
Lauri Himanen's avatar
Lauri Himanen committed
8
9
10
import cp2kparser.generic.configurationreading
import cp2kparser.generic.csvparsing
from nomadcore.caching_backend import CachingLevel
Lauri Himanen's avatar
Lauri Himanen committed
11
from nomadcore.unit_conversion.unit_conversion import convert_unit
Lauri Himanen's avatar
Lauri Himanen committed
12
13
14
15
16
17
18
19
20
21
22
23
24
25
import logging
logger = logging.getLogger("nomad")


#===============================================================================
class CP2KMDParser(MainHierarchicalParser):
    """Used to parse the CP2K calculation with run types:
        -MD
        -MOLECULAR_DYNAMICS
    """
    def __init__(self, file_path, parser_context):
        """
        """
        super(CP2KMDParser, self).__init__(file_path, parser_context)
26
        self.setup_common_matcher(CP2KCommonParser(parser_context))
Lauri Himanen's avatar
Lauri Himanen committed
27
        self.traj_iterator = None
28
        self.vel_iterator = None
29
        self.energy_iterator = None
30
        self.cell_iterator = None
31
32
33
34
35
        self.n_steps = None
        self.output_freq = None
        self.coord_freq = None
        self.velo_freq = None
        self.energy_freq = None
36
        self.cell_freq = None
37
        self.md_quicksteps = []
38
        self.ensemble = None
Lauri Himanen's avatar
Lauri Himanen committed
39
40
41
42
43
44

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

        #=======================================================================
        # Cache levels
45
        self.caching_levels.update({
46
47
48
49
50
            '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
51
52
53
54
55
56
        })

        #=======================================================================
        # SimpleMatchers
        self.md = SM(
            " MD\| Molecular Dynamics Protocol",
57
58
            forwardMatch=True,
            sections=["section_frame_sequence", "x_cp2k_section_md"],
Lauri Himanen's avatar
Lauri Himanen committed
59
            subMatchers=[
60
61
                SM( " MD\| Molecular Dynamics Protocol",
                    forwardMatch=True,
Lauri Himanen's avatar
Lauri Himanen committed
62
                    sections=["section_sampling_method", "x_cp2k_section_md_settings"],
63
                    subMatchers=[
64
65
66
67
68
69
70
71
                        SM( " MD\| Ensemble Type\s+(?P<x_cp2k_md_ensemble_type>{})".format(self.regexs.regex_word)),
                        SM( " MD\| Number of Time Steps\s+(?P<x_cp2k_md_number_of_time_steps>{})".format(self.regexs.regex_i)),
                        SM( " MD\| Time Step \[fs\]\s+(?P<x_cp2k_md_time_step__fs>{})".format(self.regexs.regex_f)),
                        SM( " MD\| Temperature \[K\]\s+(?P<x_cp2k_md_target_temperature>{})".format(self.regexs.regex_f)),
                        SM( " MD\| Temperature tolerance \[K\]\s+(?P<x_cp2k_md_target_temperature_tolerance>{})".format(self.regexs.regex_f)),
                        SM( " MD\| Pressure \[Bar\]\s+(?P<x_cp2k_md_target_pressure>{})".format(self.regexs.regex_f)),
                        SM( " MD\| Barostat time constant \[  fs\]\s+(?P<x_cp2k_md_barostat_time_constant>{})".format(self.regexs.regex_f)),
                        SM( " MD\| Print MD information every\s+(?P<x_cp2k_md_print_frequency>{}) step\(s\)".format(self.regexs.regex_i)),
72
                        SM( " MD\| File type     Print frequency\[steps\]                             File names"),
73
74
75
76
77
                        SM( " MD\| Coordinates\s+(?P<x_cp2k_md_coordinates_print_frequency>{})\s+(?P<x_cp2k_md_coordinates_filename>{})".format(self.regexs.regex_i, self.regexs.regex_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.regex_i, self.regexs.regex_word)),
                        SM( " MD\| Velocities\s+(?P<x_cp2k_md_velocities_print_frequency>{})\s+(?P<x_cp2k_md_velocities_filename>{})".format(self.regexs.regex_i, self.regexs.regex_word)),
                        SM( " MD\| Energies\s+(?P<x_cp2k_md_energies_print_frequency>{})\s+(?P<x_cp2k_md_energies_filename>{})".format(self.regexs.regex_i, self.regexs.regex_word)),
                        SM( " MD\| Dump\s+(?P<x_cp2k_md_dump_print_frequency>{})\s+(?P<x_cp2k_md_dump_filename>{})".format(self.regexs.regex_i, self.regexs.regex_word)),
78
79
80
81
82
83
84
85
86
                    ]
                ),
                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("*", "\*")),
87
88
89
90
91
92
93
94
                        SM( " INITIAL POTENTIAL ENERGY\[hartree\]     =\s+(?P<x_cp2k_md_potential_energy_instantaneous__hartree>{})".format(self.regexs.regex_f)),
                        SM( " INITIAL KINETIC ENERGY\[hartree\]       =\s+(?P<x_cp2k_md_kinetic_energy_instantaneous__hartree>{})".format(self.regexs.regex_f)),
                        SM( " INITIAL TEMPERATURE\[K\]                =\s+(?P<x_cp2k_md_temperature_instantaneous>{})".format(self.regexs.regex_f)),
                        SM( " INITIAL BAROSTAT TEMP\[K\]              =\s+(?P<x_cp2k_md_barostat_temperature_instantaneous>{})".format(self.regexs.regex_f)),
                        SM( " INITIAL PRESSURE\[bar\]                 =\s+(?P<x_cp2k_md_pressure_instantaneous__bar>{})".format(self.regexs.regex_f)),
                        SM( " INITIAL VOLUME\[bohr\^3\]                =\s+(?P<x_cp2k_md_volume_instantaneous__bohr3>{})".format(self.regexs.regex_f)),
                        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.regex_f)),
                        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.regex_f)),
95
                    ],
96
                    adHoc=self.adHoc_save_md_quickstep()
97
98
99
100
                ),
                SM( " SCF WAVEFUNCTION OPTIMIZATION",
                    endReStr=" TEMPERATURE \[K\]              =",
                    name="md_step",
101
                    forwardMatch=True,
102
103
104
105
106
107
                    repeats=True,
                    sections=["x_cp2k_section_md_step"],
                    subMatchers=[
                        self.cm.quickstep_calculation(),
                        SM( " ENSEMBLE TYPE                ="),
                        SM( " STEP NUMBER                  ="),
108
109
110
111
112
113
114
115
116
117
118
119
120
121
                        SM( " TIME \[fs\]                    =\s+(?P<x_cp2k_md_time__fs>{})".format(self.regexs.regex_f)),
                        SM( " CONSERVED QUANTITY \[hartree\] =\s+(?P<x_cp2k_md_conserved_quantity__hartree>{})".format(self.regexs.regex_f)),
                        SM( " CPU TIME \[s\]                 =\s+(?P<x_cp2k_md_cpu_time_instantaneous>{})\s+(?P<x_cp2k_md_cpu_time_average>{})".format(self.regexs.regex_f, self.regexs.regex_f)),
                        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.regex_f, self.regexs.regex_f)),
                        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.regex_f, self.regexs.regex_f)),
                        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.regex_f, self.regexs.regex_f)),
                        SM( " TEMPERATURE \[K\]              =\s+(?P<x_cp2k_md_temperature_instantaneous>{})\s+(?P<x_cp2k_md_temperature_average>{})".format(self.regexs.regex_f, self.regexs.regex_f)),
                        SM( " PRESSURE \[bar\]               =\s+(?P<x_cp2k_md_pressure_instantaneous__bar>{0})\s+(?P<x_cp2k_md_pressure_average__bar>{0})".format(self.regexs.regex_f)),
                        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.regex_f)),
                        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.regex_f)),
                        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.regex_f)),
                        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.regex_f)),
                        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.regex_f)),
                        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.regex_f)),
122
                    ],
123
                    adHoc=self.adHoc_save_md_quickstep()
124
                ),
Lauri Himanen's avatar
Lauri Himanen committed
125
126
127
128
129
130
131
132
            ]
        )

        # 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
133
            sections=["section_run"],
Lauri Himanen's avatar
Lauri Himanen committed
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
            subMatchers=[
                SM( "",
                    forwardMatch=True,
                    sections=["section_method"],
                    subMatchers=[
                        self.cm.header(),
                        self.cm.quickstep_header(),
                    ],
                ),
                self.md,
            ]
        )

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

151
152
153
154
155
156
        # Ensemble
        sampling = section.get_latest_value("x_cp2k_md_ensemble_type")
        if sampling is not None:
            sampling_map = {
                "NVE": "NVE",
                "NVT": "NVT",
Lauri Himanen's avatar
Lauri Himanen committed
157
158
                "NPT_F": "NPT",
                "NPT_I": "NPT",
159
160
161
            }
            sampling = sampling_map.get(sampling)
            if sampling is not None:
162
                self.ensemble = sampling
163
                backend.addValue("ensemble_type", sampling)
Lauri Himanen's avatar
Lauri Himanen committed
164

165
166
        # Sampling type
        backend.addValue("sampling_method", "molecular_dynamics")
Lauri Himanen's avatar
Lauri Himanen committed
167

168
169
170
171
172
        # 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")
173
        self.cell_freq = section.get_latest_value("x_cp2k_md_simulation_cell_print_frequency")
Lauri Himanen's avatar
Lauri Himanen committed
174

175
176
        # Step number
        self.n_steps = section.get_latest_value("x_cp2k_md_number_of_time_steps")
Lauri Himanen's avatar
Lauri Himanen committed
177

178
179
        # Files
        coord_filename = section.get_latest_value("x_cp2k_md_coordinates_filename")
180
        vel_filename = section.get_latest_value("x_cp2k_md_velocities_filename")
181
        energies_filename = section.get_latest_value("x_cp2k_md_energies_filename")
182
        cell_filename = section.get_latest_value("x_cp2k_md_simulation_cell_filename")
183
184
        coord_filepath = self.file_service.set_file_id(coord_filename, "coordinates")
        vel_filepath = self.file_service.set_file_id(vel_filename, "velocities")
185
        cell_filepath = self.file_service.set_file_id(cell_filename, "cell")
186
        energies_filepath = self.file_service.set_file_id(energies_filename, "energies")
Lauri Himanen's avatar
Lauri Himanen committed
187

188
        # Setup trajectory iterator
Lauri Himanen's avatar
Lauri Himanen committed
189
        traj_format = self.cache_service["trajectory_format"]
190
        if traj_format is not None and coord_filepath is not None:
Lauri Himanen's avatar
Lauri Himanen committed
191
192
193

            # Use special parsing for CP2K pdb files because they don't follow the proper syntax
            if traj_format == "PDB":
194
                self.traj_iterator = cp2kparser.generic.csvparsing.iread(coord_filepath, columns=[3, 4, 5], start="CRYST", end="END")
Lauri Himanen's avatar
Lauri Himanen committed
195
196
            else:
                try:
197
                    self.traj_iterator = cp2kparser.generic.configurationreading.iread(coord_filepath)
Lauri Himanen's avatar
Lauri Himanen committed
198
199
200
                except ValueError:
                    pass

201
202
203
204
205
206
207
208
        # Setup velocity iterator
        vel_format = self.cache_service["velocity_format"]
        if vel_format is not None and vel_filepath is not None:
            try:
                self.vel_iterator = cp2kparser.generic.configurationreading.iread(vel_filepath)
            except ValueError:
                pass

209
210
211
        # Setup energy file iterator
        if energies_filepath is not None:
            self.energy_iterator = cp2kparser.generic.csvparsing.iread(energies_filepath, columns=[0, 1, 2, 3, 4, 5, 6], comments="#")
Lauri Himanen's avatar
Lauri Himanen committed
212

213
214
215
216
        # Setup cell file iterator
        if cell_filepath is not None:
            self.cell_iterator = cp2kparser.generic.csvparsing.iread(cell_filepath, columns=[2, 3, 4, 5, 6, 7, 8, 9, 10], comments="#")

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

219
220
221
222
        # Determine the highest print frequency and use that as the number of
        # single configuration calculations
        freqs = {
            "output": [self.output_freq, True],
223
            "trajectory": [self.coord_freq, True],
224
225
            "velocities": [self.velo_freq, True],
            "energies": [self.energy_freq, True],
226
            "cell": [self.cell_freq, True],
227
        }
Lauri Himanen's avatar
Lauri Himanen committed
228

229
230
231
        # See if the files actually exist
        traj_file = self.file_service.get_file_by_id("coordinates")
        if traj_file is None:
232
            freqs["trajectory"][1] = False
233
234
235
236
237
238
        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
239
240
241
        cell_file = self.file_service.get_file_by_id("cell")
        if cell_file is None:
            freqs["cell"][1] = False
242

243
244
245
246
247
248
249
250
        # 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

251
252
253
254
255
256
        # 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

257
258
259
260
261
262
        # 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

263
264
265
266
267
268
269
270
        last_step = self.n_steps - 1
        md_steps = section["x_cp2k_section_md_step"]

        frame_sequence_potential_energy = []
        frame_sequence_temperature = []
        frame_sequence_time = []
        frame_sequence_kinetic_energy = []
        frame_sequence_conserved_quantity = []
Lauri Himanen's avatar
Lauri Himanen committed
271
        frame_sequence_pressure = []
272

Lauri Himanen's avatar
Lauri Himanen committed
273
        single_conf_gids = []
274
275
276
277
278
        i_md_step = 0
        for i_step in range(self.n_steps + 1):

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

281
282
283
284
285
            # 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:
                self.cache_service.push_array_values("simulation_cell", unit="angstrom")

286
            # Trajectory
287
288
            if freqs["trajectory"][1] and self.traj_iterator is not None:
                if (i_step + 1) % freqs["trajectory"][0] == 0 or (i_step == last_step and add_last_traj):
289
290
291
292
293
                    try:
                        pos = next(self.traj_iterator)
                    except StopIteration:
                        logger.error("Could not get the next geometries 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 trajectory file.")
                    else:
294
295
296
297
298
299
300
301
302
303
304
                        backend.addArrayValues("atom_positions", pos, unit=traj_unit)

            # Velocities
            if freqs["velocities"][1] and self.vel_iterator is not None:
                if (i_step + 1) % freqs["velocities"][0] == 0 or (i_step == last_step and add_last_vel):
                    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)
305
306
307
308

            # Energy file
            if self.energy_iterator is not None:
                if (i_step + 1) % freqs["energies"][0] == 0:
309
                    line = next(self.energy_iterator)
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326

                    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)
                    frame_sequence_kinetic_energy.append(kinetic_energy)
                    frame_sequence_temperature.append(temperature)
                    frame_sequence_potential_energy.append(potential_energy)
                    frame_sequence_conserved_quantity.append(conserved_quantity)

                    backend.addValue("energy_total", conserved_quantity)
                    backend.addValue("time_calculation", wall_time)

327
328
329
            # Cell file
            if self.cell_iterator is not None:
                if (i_step + 1) % freqs["cell"][0] == 0:
330
                    line = next(self.cell_iterator)
331
                    cell = np.reshape(line, (3, 3))
332
333
                    self.backend.addArrayValues("simulation_cell", cell, unit="angstrom")
                    # self.cache_service["simulation_cell"] = cell
334

335
336
337
338
            # Output file
            if md_steps:
                if (i_step + 1) % freqs["output"][0] == 0:
                    md_step = md_steps[i_md_step]
339
                    quickstep = self.md_quicksteps[i_md_step]
340
341
342
343
344
345
346
347
348
                    if quickstep is not None:
                        quickstep.add_latest_value("x_cp2k_atom_forces", "atom_forces")
                        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")
                        quickstep.add_latest_value(["x_cp2k_section_scf_iteration", "x_cp2k_energy_XC_scf_iteration"], "energy_XC_scf_iteration")
                        backend.closeSection("section_scf_iteration", scfGID)
                    i_md_step += 1
Lauri Himanen's avatar
Lauri Himanen committed
349
350
351
                    pressure = md_step.get_latest_value("x_cp2k_md_pressure_instantaneous")
                    if pressure is not None:
                        frame_sequence_pressure.append(pressure)
352
353
354
355
356

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

        # Add the summaries to frame sequence
Lauri Himanen's avatar
Lauri Himanen committed
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
        frame_sequence_potential_energy = convert_unit(np.array(frame_sequence_potential_energy), "hartree")
        frame_sequence_kinetic_energy = convert_unit(np.array(frame_sequence_kinetic_energy), "hartree")
        frame_sequence_conserved_quantity = convert_unit(np.array(frame_sequence_conserved_quantity), "hartree")
        frame_sequence_time = np.array(frame_sequence_time)
        frame_sequence_temperature = np.array(frame_sequence_temperature)
        frame_sequence_pressure = np.array(frame_sequence_pressure)

        backend.addArrayValues("frame_sequence_potential_energy", frame_sequence_potential_energy)
        backend.addArrayValues("frame_sequence_kinetic_energy", frame_sequence_kinetic_energy)
        backend.addArrayValues("frame_sequence_conserved_quantity", frame_sequence_conserved_quantity)
        backend.addArrayValues("frame_sequence_temperature", frame_sequence_temperature)
        backend.addArrayValues("frame_sequence_time", frame_sequence_time, unit="fs")
        if frame_sequence_pressure.size != 0:
            backend.addArrayValues("frame_sequence_pressure", frame_sequence_pressure)

        # Number of frames
        backend.addValue("number_of_frames_in_sequence", self.n_steps + 1)

        # Reference to sampling method
        backend.addValue("frame_sequence_to_sampling_ref", 0)

        # References to local frames
        backend.addArrayValues("frame_sequence_local_frames_ref", np.array(single_conf_gids))

        # Temperature stats
        mean_temp = frame_sequence_temperature.mean()
383
        std_temp = frame_sequence_temperature.std()
Lauri Himanen's avatar
Lauri Himanen committed
384
385
386
387
        backend.addArrayValues("frame_sequence_temperature_stats", np.array([mean_temp, std_temp]))

        # Potential energy stats
        mean_pot = frame_sequence_potential_energy.mean()
388
        std_pot = frame_sequence_potential_energy.std()
Lauri Himanen's avatar
Lauri Himanen committed
389
390
391
392
        backend.addArrayValues("frame_sequence_potential_energy_stats", np.array([mean_pot, std_pot]))

        # Kinetic energy stats
        mean_kin = frame_sequence_kinetic_energy.mean()
393
        std_kin = frame_sequence_kinetic_energy.std()
Lauri Himanen's avatar
Lauri Himanen committed
394
395
396
397
        backend.addArrayValues("frame_sequence_kinetic_energy_stats", np.array([mean_kin, std_kin]))

        # Conserved quantity stats
        mean_cons = frame_sequence_conserved_quantity.mean()
398
        std_cons = frame_sequence_conserved_quantity.std()
Lauri Himanen's avatar
Lauri Himanen committed
399
400
401
402
403
        backend.addArrayValues("frame_sequence_conserved_quantity_stats", np.array([mean_cons, std_cons]))

        # Pressure stats
        if frame_sequence_pressure.size != 0:
            mean_pressure = frame_sequence_pressure.mean()
404
            std_pressure = frame_sequence_pressure.std()
Lauri Himanen's avatar
Lauri Himanen committed
405
            backend.addArrayValues("frame_sequence_pressure_stats", np.array([mean_pressure, std_pressure]))
406
407
408
409
410
411
412
413
414
415
416

    #===========================================================================
    # 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