workflow.py 18 KB
Newer Older
1
#
Markus Scheidgen's avatar
Markus Scheidgen committed
2
3
4
5
6
# Copyright The NOMAD Authors.
#
# This file is part of NOMAD. See https://nomad-lab.eu for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
7
8
9
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
Markus Scheidgen's avatar
Markus Scheidgen committed
10
#     http://www.apache.org/licenses/LICENSE-2.0
11
12
#
# Unless required by applicable law or agreed to in writing, software
Markus Scheidgen's avatar
Markus Scheidgen committed
13
# distributed under the License is distributed on an "AS IS" BASIS,
14
15
16
# 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.
Markus Scheidgen's avatar
Markus Scheidgen committed
17
#
18
19
20
21

import numpy as np

from nomad.normalizing.normalizer import Normalizer
22
from nomad.datamodel import EntryArchive
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
23
from nomad.datamodel.metainfo.workflow import (
24
25
    Workflow, SinglePoint, GeometryOptimization, MolecularDynamics, Phonon, Elastic,
    Thermodynamics)
26
27


28
29
def resolve_difference(values):
    delta_values = None
30

31
32
33
34
    values = [v for v in values if v is not None]
    for n in range(-1, -len(values), -1):
        delta_values = abs(values[n] - values[n - 1])
        if delta_values != 0.0:
35
36
            break

37
    return delta_values
38
39


Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
40
41
class TaskNormalizer(Normalizer):
    def __init__(self, entry_archive, workflow_index):
42
        super().__init__(entry_archive)
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
43
44
45
        workflow_index = workflow_index if len(entry_archive.workflow) < workflow_index else -1
        self.workflow = entry_archive.workflow[workflow_index]
        run = self.workflow.run_ref
46
        self.run = run if run else entry_archive.run[-1]
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
47

48

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
49
class SinglePointNormalizer(TaskNormalizer):
50
    def normalize(self):
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
51
52
        super().normalize()
        self.section = self.workflow.single_point
53
        if not self.section:
54
            self.section = self.workflow.m_create(SinglePoint)
55

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
56
        if not self.section.method:
57
            try:
58
                method = self.run.method[-1]
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
59
                self.section.method = method.electronic.method
60
61
62
            except Exception:
                pass

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
63
        scc = self.run.calculation
64
65
66
        if not scc:
            return

67
68
        if not self.section.n_scf_steps:
            self.section.n_scf_steps = len(scc[-1].scf_iteration)
69

70
        energies = [scf.energy.total.value for scf in scc[-1].scf_iteration if scf.energy is not None and scf.energy.total is not None]
71
        delta_energy = resolve_difference(energies)
72
73
74
75
76
        if not self.section.final_scf_energy_difference and delta_energy is not None:
            self.section.final_scf_energy_difference = delta_energy

        if not self.section.is_converged and delta_energy is not None:
            try:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
77
                threshold = self.run.method[-1].scf.threshold_energy_change
78
                self.section.is_converged = bool(delta_energy <= threshold)
79
80
81
82
            except Exception:
                pass

        if not self.section.with_density_of_states:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
83
            self.section.with_density_of_states = len(scc[-1].dos_electronic) > 0
84
85

        if not self.section.with_bandstructure:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
86
            self.section.with_bandstructure = len(scc[-1].band_structure_electronic) > 0
87
88

        if not self.section.with_eigenvalues:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
89
            self.section.with_eigenvalues = len(scc[-1].eigenvalues) > 0
90
91

        if not self.section.with_volumetric_data:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
92
93
            self.section.with_volumetric_data = (
                len(scc[-1].potential) > 0 or len(scc[-1].density_charge) > 0)
94
95

        if not self.section.with_excited_states:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
96
            self.section.with_excited_states = len(scc[-1].excited_states) > 0
97
98


Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
99
class GeometryOptimizationNormalizer(TaskNormalizer):
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
100

101
    def _get_geometry_optimization_type(self):
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
102
        sec_system = self.run.system
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
        if not sec_system:
            return

        def compare_cell(cell1, cell2):
            if (cell1 == cell2).all():
                return None
            else:
                cell1_normed = cell1 / np.linalg.norm(cell1)
                cell2_normed = cell2 / np.linalg.norm(cell2)
                if (cell1_normed == cell2_normed).all():
                    return 'cell_volume'
                else:
                    return 'cell_shape'

        if len(sec_system) < 2:
            return 'static'

        else:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
121
122
123
124
125
126
127
            if sec_system[0].atoms is None or sec_system[-1].atoms is None:
                return 'static'

            cell_init = sec_system[0].atoms.lattice_vectors
            cell_final = sec_system[-1].atoms.lattice_vectors
            if cell_init is None or cell_final is None:
                return 'static'
128

129
            cell_relaxation = compare_cell(cell_init.magnitude, cell_final.magnitude)
130
131
132
133

            if cell_relaxation is not None:
                return cell_relaxation

134
135
            atom_pos_init = sec_system[0].atoms.positions
            atom_pos_final = sec_system[-1].atoms.positions
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
136
137
138

            if atom_pos_init is None or atom_pos_final is None:
                return 'static'
139

140
            if (atom_pos_init.magnitude == atom_pos_final.magnitude).all():
141
142
143
144
                return 'static'

            return 'ionic'

145
    def normalize(self):
146
        super().normalize()
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
147
148
149
        self.section = self.workflow.geometry_optimization
        if self.section is None:
            self.section = self.workflow.m_create(GeometryOptimization)
150

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
151
        if not self.section.type:
152
153
            try:
                geometry_optimization_type = self._get_geometry_optimization_type()
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
154
                self.section.type = geometry_optimization_type
155
156
            except Exception:
                pass
157

158
        if not self.section.optimization_steps:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
159
            scc = self.run.calculation
160
            self.section.optimization_steps = len(scc)
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
161

162
        if not self.section.final_energy_difference:
163
            energies = []
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
164
165
166
            for scc in self.run.calculation:
                if scc.energy is not None and scc.energy.total is not None:
                    energies.append(scc.energy.total.value)
167

168
            delta_energy = resolve_difference(energies)
169
            if delta_energy is not None:
170
                self.section.final_energy_difference = delta_energy
171

172
        if not self.section.final_force_maximum:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
173
174
175
            scc = self.run.calculation
            if len(scc) > 0:
                if scc[-1].forces is not None and scc[-1].forces.total is not None:
176
177
178
179
                    forces = scc[-1].forces.total.value
                    if forces is not None:
                        max_force = np.max(np.linalg.norm(forces.magnitude, axis=1))
                        self.section.final_force_maximum = max_force * forces.units
180

181
182
183
        # Store the energies as an explicit list. If a step within the
        # trajectory does not contain an energy the rest of the energies in the
        # trajectory are not included.
184
        trajectory = self.workflow.calculations_ref
185
186
187
188
189
        if trajectory:
            n_steps = len(trajectory)
            energies = []
            invalid = False
            for step in range(n_steps):
190
191
192
                try:
                    energy = trajectory[step].energy.total.value
                except (IndexError, AttributeError):
193
194
195
196
197
198
199
200
                    invalid = True
                    break
                energies.append(energy.magnitude)
            if invalid:
                self.logger.warning("energy not reported for an scc that is part of a geometry optimization")
            if energies:
                self.section.energies = energies

201
202
        if not self.section.final_displacement_maximum:
            try:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
203
                system = self.run.system
204
                displacements = [np.max(np.abs(
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
205
                    system[n].atoms.positions - system[n - 1].atoms.positions)) for n in range(1, len(system))]
206
207
208
209
210
211
212
213
                self.section.final_displacement_maximum = resolve_difference(displacements)
            except Exception:
                pass

        if not self.section.is_converged_geometry:
            # we can have several criteria for convergence: energy, force, displacement
            criteria = []
            try:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
214
                criteria.append(self.section.final_energy_difference <= self.section.convergence_tolerance_energy_difference)
215
216
217
218
            except Exception:
                pass

            try:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
219
                criteria.append(self.section.final_force_maximum <= self.section.convergence_tolerance_force_maximum)
220
221
222
223
            except Exception:
                pass

            try:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
224
                criteria.append(self.section.final_displacement_maximum <= self.section.convergence_tolerance_displacement_maximum)
225
226
227
228
229
230
231
            except Exception:
                pass

            # converged when either criterion is met
            if criteria:
                self.section.is_converged_geometry = True in criteria

232

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
233
class PhononNormalizer(TaskNormalizer):
234
    def _get_n_imaginary_frequencies(self):
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
235
        scc = self.run.calculation
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
236
237
        if not scc:
            return
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
238
        sec_band = scc[0].band_structure_phonon
239
240
        if not sec_band:
            return
241
        result = 0
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
242
        for band_segment in sec_band[0].segment:
243
            freq = band_segment.energies.magnitude
244
245
246
247
            result += np.count_nonzero(np.array(freq) < 0)
        return result

    def normalize(self):
248
        super().normalize()
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
249
        self.section = self.workflow.phonon
250

251
        if not self.section:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
252
            self.section = self.workflow.m_create(Phonon)
253
254
255
256
257
258

        if not self.section.n_imaginary_frequencies:
            # get number from bands (not complete as this is not the whole mesh)
            self.section.n_imaginary_frequencies = self._get_n_imaginary_frequencies()


Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
259
class ElasticNormalizer(TaskNormalizer):
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
260
    def _resolve_mechanical_stability(self):
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
261
262
        spacegroup, c = None, None
        try:
263
264
            spacegroup = self.run.system[-1].symmetry[-1].space_group_number
            c = self.section.elastic_constants_matrix_second_order
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
265
266
        except Exception:
            return False
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
267

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
268
269
        if c is None or spacegroup is None:
            return False
270

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
        # see Phys. Rev B 90, 224104 (2014)
        res = False
        if spacegroup <= 2:  # Triclinic
            res = np.count_nonzero(c < 0)
        elif spacegroup <= 15:  # Monoclinic
            res = np.count_nonzero(c < 0)
        elif spacegroup <= 74:  # Orthorhombic
            res =\
                c[0][0] > 0 and c[0][0] * c[1][1] > c[0][1] ** 2 and\
                c[0][0] * c[1][1] * c[2][2] + 2 * c[0][1] * c[0][2] * c[1][2] -\
                c[0][0] * c[1][2] ** 2 - c[1][1] * c[0][2] ** 2 - c[2][2] * c[0][1] ** 2 > 0 and\
                c[3][3] > 0 and c[4][4] > 0 and c[5][5] > 0
        elif spacegroup <= 88:  # Tetragonal II
            res =\
                c[0][0] > abs(c[0][1]) and\
                2 * c[0][2] ** 2 < c[2][2] * (c[0][0] + c[0][1])
        elif spacegroup <= 142:  # Tetragonal I
            res =\
                c[0][0] > abs(c[0][1]) and\
                2 * c[0][2] ** 2 < c[2][2] * (c[0][0] + c[0][1]) and\
                c[3][3] > 0 and c[5][5] > 0
        elif spacegroup <= 148:  # rhombohedral II
            res =\
                c[0][0] > abs(c[0][1]) and c[3][3] > 0 and\
                c[0][2] ** 2 < (0.5 * c[2][2] * (c[0][0] + c[0][1])) and\
                c[0][3] ** 2 + c[0][4] ** 2 < 0.5 * c[3][3] * (c[0][0] - c[0][1])
        elif spacegroup <= 167:  # rhombohedral I
            res =\
                c[0][0] > abs(c[0][1]) and c[3][3] > 0 and\
                c[0][2] ** 2 < 0.5 * c[2][2] * (c[0][0] + c[0][1]) and\
                c[0][3] ** 2 < 0.5 * c[3][3] * (c[0][0] - c[0][1])
        elif spacegroup <= 194:  # hexagonal I
            res =\
                c[0][0] > abs(c[0][1]) and\
                2 * c[0][2] ** 2 < c[2][2] * (c[0][0] + c[0][1]) and\
                c[3][3] > 0 and c[5][5] > 0
        else:  # cubic
            res = c[0][0] - c[0][1] > 0 and c[0][0] + 2 * c[0][1] > 0 and c[3][3] > 0

        return res

    def _get_maximum_fit_error(self):
        max_error = 0.0
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
314
315
316
        if len(self.run.calculation) == 0:
            return max_error

317
318
319
        for diagram in self.section.strain_diagrams:
            if diagram.type == 'cross-validation':
                error = np.amax(diagram.value)
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
320
321
322
323
324
                max_error = error if error > max_error else max_error

        return max_error

    def normalize(self):
325
        super().normalize()
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
326
        self.section = self.workflow.elastic
327

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
328
        if not self.section:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
329
            self.section = self.workflow.m_create(Elastic)
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
330

331
332
333
334
335
336
        if self.section.is_mechanically_stable is None:
            self.section.is_mechanically_stable = bool(self._resolve_mechanical_stability())

        if self.section.fitting_error_maximum is None:
            self.section.fitting_error_maximum = self._get_maximum_fit_error()

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
337

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
338
class MolecularDynamicsNormalizer(TaskNormalizer):
339
    def _is_with_thermodynamics(self):
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
340
341
342
343
        try:
            return len(self.run.calculation[-1].thermodynamics[-1].values()) > 0
        except Exception:
            return False
344
345

    def _is_with_trajectory(self):
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
346
347
348
349
        try:
            return self.run.system[-1].atoms.positions is not None
        except Exception:
            return False
350
351

    def normalize(self):
352
        super().normalize()
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
353
        self.section = self.workflow.molecular_dynamics
354
355

        if not self.section:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
356
            self.section = self.workflow.m_create(MolecularDynamics)
357
358
359
360
361
362

        if self.section.with_thermodynamics is None:
            self.section.with_thermodynamics = self._is_with_thermodynamics()

        if self.section.with_trajectory is None:
            self.section.with_trajectory = self._is_with_trajectory()
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
363
364


365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
class ThermodynamicsNormalizer(TaskNormalizer):
    def normalize(self):
        super().normalize()
        self.section = self.workflow.thermodynamics

        if not self.run.calculation or not self.run.calculation[0].thermodynamics:
            return

        if not self.section:
            self.section = self.workflow.m_create(Thermodynamics)

        def set_thermo_property(name):
            values = []
            quantity = None
            for scc in self.run.calculation:
                try:
                    for thermo in scc.thermodynamics:
                        quantity = thermo[name]
                        values.append(quantity.magnitude if hasattr(quantity, 'magnitude') else quantity)
                except Exception:
                    pass
386
387
388
            if len(values) == 0:
                return

389
            unit = quantity.units if hasattr(quantity, 'units') else 1.0
390
391
            setattr(self.section, name, np.array(values) * unit)

392
        if self.section.temperature is None:
393
394
            set_thermo_property('temperature')

395
        if self.section.pressure is None:
396
397
            set_thermo_property('pressure')

398
399
        if self.section.helmholtz_free_energy is None:
            set_thermo_property('helmholtz_free_energy')
400

401
        if self.section.vibrational_free_energy_at_constant_volume is None:
402
403
            set_thermo_property('vibrational_free_energy_at_constant_volume')

404
        if self.section.heat_capacity_c_v is None:
405
406
407
408
409
            set_thermo_property('heat_capacity_c_v')

        # TODO add values for specific energy


410
411
class WorkflowNormalizer(Normalizer):
    '''
412
    This normalizer produces information specific to a workflow.
413
414
415
    '''
    def __init__(self, entry_archive):
        super().__init__(entry_archive)
416
417
        self._elastic_programs = ['elastic']
        self._phonon_programs = ['phonopy']
418
        self._molecular_dynamics_programs = ['lammps']
419

420
    def _resolve_workflow_type(self, run):
421
        # resolve it from parser
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
422
        workflow_type = None
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
423
424
425
426
427
        try:
            program_name = run.program.name
        except Exception:
            return

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
428
429
        if program_name:
            program_name = program_name.lower()
430

431
        if program_name in self._elastic_programs:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
432
            workflow_type = 'elastic'
433

434
        elif program_name in self._molecular_dynamics_programs:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
435
            workflow_type = 'molecular_dynamics'
436

437
        elif program_name in self._phonon_programs:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
438
            workflow_type = 'phonon'
439

440
        # resolve if from scc
441
442
        if workflow_type is None:
            if len(run.calculation) == 1:
443
                workflow_type = 'single_point'
444
445
            else:
                workflow_type = 'geometry_optimization'
446

447
        return workflow_type
448
449

    def normalize(self, logger=None) -> None:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
450
        super().normalize()
451
452

        # Do nothing if section_run is not present
453
        if not self.entry_archive.run:
454
455
            return

456
        if not self.entry_archive.workflow:
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
457
            self.entry_archive.m_create(Workflow)
458

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
459
        for n, sec_workflow in enumerate(self.entry_archive.workflow):
460
461
            # we get reference the referenced run from which information can be extracted
            sec_run = sec_workflow.run_ref
462
            sec_run = sec_run if sec_run else self.entry_archive.run[-1]
463

464
465
466
467
468
469
470
471
472
            scc = sec_run.calculation
            if not sec_workflow.calculation_result_ref:
                if scc:
                    sec_workflow.calculation_result_ref = scc[-1]

            if not sec_workflow.calculations_ref:
                if scc:
                    sec_workflow.calculations_ref = scc

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
473
            if sec_workflow.type is None:
474
                workflow_type = self._resolve_workflow_type(sec_run)
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
475
                sec_workflow.type = workflow_type
476

477
            if sec_workflow.type == 'geometry_optimization':
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
478
                GeometryOptimizationNormalizer(self.entry_archive, n).normalize()
479

480
            elif sec_workflow.type == 'phonon':
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
481
                PhononNormalizer(self.entry_archive, n).normalize()
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
482

483
            elif sec_workflow.type == 'elastic':
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
484
                ElasticNormalizer(self.entry_archive, n).normalize()
485

486
            elif sec_workflow.type == 'molecular_dynamics':
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
487
                MolecularDynamicsNormalizer(self.entry_archive, n).normalize()
488

489
            elif sec_workflow.type == 'single_point':
Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
490
                SinglePointNormalizer(self.entry_archive, n).normalize()
491

492
493
494
            # add thermodynamics data
            ThermodynamicsNormalizer(self.entry_archive, n).normalize()

Alvin Noe Ladines's avatar
Alvin Noe Ladines committed
495
496
497
            # remove the section workflow again, if the parser/normalizer could not produce a result
            if sec_workflow.calculation_result_ref is None:
                self.entry_archive.m_remove_sub_section(EntryArchive.workflow, n)