parser_quantum_espresso.py 139 KB
Newer Older
1
2
3
import setup_paths
from nomadcore.simple_parser import mainFunction, SimpleMatcher as SM
from nomadcore.local_meta_info import loadJsonFile, InfoKindEl
Henning Glawe's avatar
Henning Glawe committed
4
5
6
import os
import sys
import json
7
import re
8
import logging
Henning Glawe's avatar
Henning Glawe committed
9
10
import nomadcore.unit_conversion.unit_conversion as unit_conversion
import math
11
import numpy as np
12
import QuantumEspressoCommon as QeC
Henning Glawe's avatar
Henning Glawe committed
13
from nomadcore.parser_backend import valueForStrValue
14
from QuantumEspressoCommon import RE_f, RE_i, cRE_f, cRE_i
Henning Glawe's avatar
Henning Glawe committed
15
from QuantumEspressoXC import translate_qe_xc_num
16
from nomadcore.parser_backend import valueForStrValue
17
18
19
20


LOGGER = logging.getLogger(__name__)

21

22
23
24
25
26
27
# Lookup table mapping string to bool flag
QE_SPIN_NONCOLLINEAR = {
    'Noncollinear': True,
    'Non magnetic': False,
}

28
29
30
31
32
33
# Lookup table mapping diagonalization scheme match to string
QE_DIAGONALIZATION = {
    'Davidson diagonalization with overlap': 'davidson',
    'CG style diagonalization': 'conjugate_gradient',
}

34
35
36
37
38
39
# Lookup table for VCSMD notation of atom position units
QE_VCSMD_ATPOS_UNITS = {
    'cryst coord': 'crystal',
    'cart coord (alat unit)': 'alat',
}

40
41
42
43
44
45
46
47
QE_MD_RELAX_SAMPLING_METHOD = {
    'langevin_overdamped_dynamics': 'langevin_dynamics',
    'BFGS': 'geometry_optimization',
    'damped_dynamics': 'geometry_optimization',
    'molecular_dynamics': 'molecular_dynamics',
    'vcsmd': 'molecular_dynamics',
    'vcsmd_wentzcovitch_damped_minization': 'geometry_optimization',
}
48

49
class QuantumEspressoParserPWSCF(QeC.ParserQuantumEspresso):
50
51
    """main place to keep the parser status, open ancillary files,..."""
    def __init__(self):
52
53
        QeC.ParserQuantumEspresso.__init__(
            self, re_program_name=re.compile(r"^PWSCF$"))
54

Henning Glawe's avatar
Henning Glawe committed
55
56
57
    def initialize_values(self):
        """allows to reset values if the same superContext is used to parse
        different files"""
58
        self.sectionIdx = {}
59
        self.openSectionIdx = {}
Henning Glawe's avatar
Henning Glawe committed
60
        self.tmp = {}
Henning Glawe's avatar
Henning Glawe committed
61
        self.alat = None
Henning Glawe's avatar
Henning Glawe committed
62
        self.section = {}
Henning Glawe's avatar
Henning Glawe committed
63
64
65
66

    def startedParsing(self, path, parser):
        """called when parsing starts"""
        self.parser = parser
Henning Glawe's avatar
Henning Glawe committed
67
        # reset values if same superContext is used to parse different files
Henning Glawe's avatar
Henning Glawe committed
68
69
        self.initialize_values()

70
71
    def onClose_section_basis_set_cell_dependent(
            self, backend, gIndex, section):
72
73
74
        if section['basis_set_planewave_cutoff'] is None:
            LOGGER.error("basis_set_planewave_cutoff is not set")
            return
75
76
77
78
79
        pwc_ry = unit_conversion.convert_unit(
            section['basis_set_planewave_cutoff'][-1],
            'J', 'rydberg')
        backend.addValue('basis_set_cell_dependent_name', 'PW_%.1f' % (pwc_ry))
        backend.addValue('basis_set_cell_dependent_kind', 'plane_waves')
80
81
        method_basis_set_gIndex = backend.openSection('section_method_basis_set')
        backend.addValue('mapping_section_method_basis_set_cell_associated', gIndex)
82
        backend.addValue('method_basis_set_kind', self.tmp['method_basis_set_kind'])
83
        backend.closeSection('section_method_basis_set', method_basis_set_gIndex)
84

85
86
    def onOpen_section_method(
            self, backend, gIndex, section):
87
        self.sectionIdx['section_method'] = gIndex
88
        self.cache_t_pseudopotential = {}
89
        self.cache_t_pp_report = {}
90
        self.cache_t_pp_renorm_wfc = {}
91
92
        self.cache_t_method = section
        self.atom_kind_idx = -1
93

94
95
    def onClose_x_qe_t_section_pseudopotential(
            self, backend, gIndex, section):
96
97
98
        if section['x_qe_t_pp_label'] is None:
            LOGGER.error("x_qe_t_pp_label is not set")
            return
99
100
        self.cache_t_pseudopotential[section['x_qe_t_pp_label'][0]] = section

101
102
    def onClose_section_method_atom_kind(
            self, backend, gIndex, section):
103
        self.atom_kind_idx += 1
104
105
        pp_label = section['x_qe_pp_label'][0]
        pp = self.cache_t_pseudopotential[pp_label]
Henning Glawe's avatar
Henning Glawe committed
106
        for key, value in sorted(pp.simpleValues.items()):
107
            if key == 'x_qe_t_pp_label' or key == 'x_qe_pp_valence' or value is None:
108
109
110
                continue
            target = re.sub(r'^x_qe_t_',r'x_qe_',key)
            if target == key:
111
                raise Exception('found non-temporary key in pseudopotential cache: "%s"' % (key))
112
            backend.addValue(target, value[-1])
113
114
115
116
117
118
119
120
        if pp['x_qe_t_pp_idx'] is not None:
            pp_num = pp['x_qe_t_pp_idx'][-1]
            pp_report = self.cache_t_pp_report.get(pp_num, None)
            if pp_report is not None:
                backend.addValue('x_qe_pp_report_version', pp_report['x_qe_t_pp_report_version'][-1])
                backend.addValue('x_qe_pp_report_contents', "\n".join(pp_report['x_qe_t_pp_report_line']))
        else:
            LOGGER.error("x_qe_t_pp_idx is not set")
121
122
123
        if pp['x_qe_t_pp_filename'] is not None:
            fpath = pp['x_qe_t_pp_filename'][-1]
            basename = re.sub(r".*\/", r"", fpath)
124
125
126
            if len(basename) > 0:
                backend.addValue('method_atom_kind_pseudopotential_name',
                                 basename)
127
128
129
            renorm_info = self.cache_t_pp_renorm_wfc.get(basename, None)
            if renorm_info is not None:
                backend.addValue('x_qe_pp_renormalized_wfc', renorm_info)
130
131
132
133
134
135
136
        # set charge/magnetization integration radius if applicable
        integration_radius=self.cache_t_method['x_qe_t_species_integration_radius']
        if integration_radius:
            backend.addValue('x_qe_species_integration_radius',
                             unit_conversion.convert_unit(
                                 integration_radius[self.atom_kind_idx],
                                 'usrAlat'))
137
138
139
140
141
142
        # DFT-D data
        cache_dft_d = self.tmp.get('dispersion_correction', {})
        dft_d = cache_dft_d.get(pp_label, None)
        if dft_d is not None:
            for k, v in dft_d.items():
                backend.addValue(k, v)
143
144
145
146

    def onClose_x_qe_t_section_pp_report(
            self, backend, gIndex, section):
        self.cache_t_pp_report[section['x_qe_t_pp_report_species'][0]] = section
147

148
149
150
151
152
    def onClose_x_qe_t_section_pp_warning(
            self, backend, gIndex, section):
        self.cache_t_pp_renorm_wfc[
            section['x_qe_t_pp_warning_filename'][-1]] = section['x_qe_t_pp_warning_wfclabel'][-1]

153
154
155
156
157
158
159
160
161
162
163
164
    def onOpen_x_qe_t_section_input_occupations(
            self, backend, gIndex, section):
        self.tmp['occ_spin'] = 'none'
        self.tmp['occ_vals_spin'] = {}
        self.tmp['occ_vals_spin']['none'] = []
        self.tmp['occ_vals'] = self.tmp['occ_vals_spin']['none']

    def adHoc_input_occupations_spin(self, parser):
        self.tmp['occ_spin'] = parser.lastMatch['x_qe_t_input_occupations_spin']
        self.tmp['occ_vals_spin'][self.tmp['occ_spin']] = []
        self.tmp['occ_vals'] = self.tmp['occ_vals_spin'][self.tmp['occ_spin']]

165
166
    def onClose_x_qe_t_section_input_occupations(
            self, backend, gIndex, section):
167
168
169
        if len(self.tmp['occ_vals_spin']['none']) > 0:
            # spin-unpolarized case
            input_occ = np.array([ self.tmp['occ_vals_spin']['none'] ])
170
        else:
171
172
173
174
175
            input_occ = np.array([
                self.tmp['occ_vals_spin']['up'],
                self.tmp['occ_vals_spin']['down']
            ])
        if input_occ.size < 1:
176
            LOGGER.error("closing x_qe_t_section_input_occupations, but no input occupations parsed")
177
178
        else:
            LOGGER.error("FIXME: implement proper output of x_qe_t_section_input_occupations (shape %s)", str(input_occ.shape))
179

180
181
    def onClose_section_method(
            self, backend, gIndex, section):
182
        # set flag if we deal with user-enforced XC functional
183
184
        if section['x_qe_t_xc_functional_shortname_enforced'] is not None:
            backend.addValue('x_qe_xc_functional_user_enforced', True)
185
        # translate XC functional to section_xc_functionals
186
        method_xc_functionals = None
187
188
189
        xc_functionals = None
        if section['x_qe_xc_functional_num'] is not None:
            xc_functional_num = section['x_qe_xc_functional_num'][-1]
190
191
192
193
194
195
196
197
            exx_fraction = None
            if section['x_qe_t_exact_exchange_fraction']:
                # first SCF in EXX is done without HF-X, cache for later
                self.tmp['exx_fraction'] = section['x_qe_t_exact_exchange_fraction'][-1]
                self.tmp['xc_functional_num'] = xc_functional_num
                exx_fraction = 0.0
            elif section['x_qe_exact_exchange_fraction']:
                exx_fraction = section['x_qe_exact_exchange_fraction'][-1]
198
            (method_xc_functionals, xc_functionals) = translate_qe_xc_num(xc_functional_num, exx_fraction)
199
200
        else:
            LOGGER.error("x_qe_xc_functional_num is not set")
201
        if method_xc_functionals is not None:
202
203
204
205
206
            # NOTE: value of XC_functional generated by translate_qe_xc_num
            #       does not fully respect the metaInfo definition
            #       when XC_functional_parameters are involved.
            # Therefore, remove it here
            method_xc_functionals.pop('XC_functional', None)
207
            self.addDict(backend, method_xc_functionals)
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
        if xc_functionals is not None:
            for xc_functional in xc_functionals:
                self.addSectionDict(backend, 'section_XC_functionals', xc_functional)
        else:
            LOGGER.error("error getting xc_functionals")
        if section['x_qe_t_allocated_array_name'] is not None:
            backend.addArrayValues('x_qe_allocated_array_name', np.asarray(section['x_qe_t_allocated_array_name']))
        if section['x_qe_t_allocated_array_size'] is not None:
            backend.addArrayValues('x_qe_allocated_array_size', np.asarray(section['x_qe_t_allocated_array_size']))
        if section['x_qe_t_allocated_array_dimensions'] is not None:
            backend.addArrayValues('x_qe_allocated_array_dimensions', np.asarray(section['x_qe_t_allocated_array_dimensions']))
        if section['x_qe_t_temporary_array_name'] is not None:
            backend.addArrayValues('x_qe_temporary_array_name', np.asarray(section['x_qe_t_temporary_array_name']))
        if section['x_qe_t_temporary_array_size'] is not None:
            backend.addArrayValues('x_qe_temporary_array_size', np.asarray(section['x_qe_t_temporary_array_size']))
        if section['x_qe_t_temporary_array_dimensions'] is not None:
            backend.addArrayValues('x_qe_temporary_array_dimensions', np.asarray(section['x_qe_t_temporary_array_dimensions']))
225
226
227
228
229
230
231
232
        if section['x_qe_t_spin_orbit_magn'] is not None:
            backend.addValue('x_qe_spin_orbit', (section['x_qe_t_spin_orbit_mode'][-1] == 'with'))
            noncollinear = QE_SPIN_NONCOLLINEAR.get(section['x_qe_t_spin_orbit_magn'][-1], None)
            if noncollinear is None:
                LOGGER.error("unimplemented value for 'x_qe_t_spin_orbit_magn': '%s'",
                             str(section['x_qe_t_spin_orbit_magn'][-1]))
            else:
                backend.addValue('x_qe_spin_noncollinear', noncollinear)
233
234
        # TODO check for LDA+U and switch to 'DFT+U' in that clase
        backend.addValue('electronic_structure_method', 'DFT')
235

236
237
238
    def onOpen_section_single_configuration_calculation(
            self, backend, gIndex, section):
        exx_refine = self.tmp.pop('exx_refine', None)
239
        md_relax = self.tmp.get('md_relax', None)
240
        have_new_system = self.tmp.pop('have_new_system', None)
241
242
243
244
        old_scc = self.section.get('single_configuration_calculation', None)
        if old_scc is not None:
            if md_relax is not None:
                LOGGER.info('new scc due to md_relax=="%s"', str(md_relax))
245
246
247
                if self.tmp.get('frames', None) is None:
                    # add previous scc
                    self.tmp['frames'] = [self.sectionIdx['single_configuration_calculation']]
248
249
250
251
            elif exx_refine is not None:
                LOGGER.info('new scc due to exx_refine=="%s"', str(exx_refine))
            else:
                raise Exception('encountered new section_single_configuration_calculation without knowing why')
252
253
        if md_relax:
            if have_new_system:
254
                self.tmp['frames'].append(gIndex)
255
256
257
258
            elif exx_refine:
                LOGGER.info("EXX refinement in MD/relax run")
            else:
                raise Exception('running MD/relax calculation, but no new cell or atom coordinates found')
259
260
        if exx_refine:
            backend.addValue('x_qe_exx_refine', True)
261
262
263
264
265
266
            exx_fraction = self.tmp.pop('exx_fraction', None)
            if exx_fraction:
                # we need to add a section_method including exx
                method_gIndex = backend.openSection('section_method')
                backend.addValue('x_qe_xc_functional_num', self.tmp.pop('xc_functional_num'))
                backend.addValue('x_qe_exact_exchange_fraction',  exx_fraction)
267
268
269
270
271
272
                self.addSectionDict(
                    backend, 'section_method_to_method_refs', {
                        'method_to_method_ref': (method_gIndex-1),
                        'method_to_method_kind': 'starting_point',
                    }
                )
273
                backend.closeSection('section_method', method_gIndex)
274
275
276
277
278
279
            self.addSectionDict(
                backend, 'section_calculation_to_calculation_refs', {
                    'calculation_to_calculation_ref': (gIndex-1),
                    'calculation_to_calculation_kind': 'starting_point',
                }
            )
280
        self.close_header_sections(backend)
281
282
        # reset temporary storage for band structures
        self.tmp['k_energies'] = []
283
        self.tmp['k_occupations'] = []
284
285
        self.tmp['kspin'] = {}
        self.section['single_configuration_calculation'] = section
286
        self.sectionIdx['single_configuration_calculation'] = gIndex
287

Henning Glawe's avatar
Henning Glawe committed
288
289
290
291
    def onClose_section_single_configuration_calculation(
            self, backend, gIndex, section):
        """trigger called when section_single_configuration_calculation
        is closed"""
292
293
        backend.addValue('single_configuration_to_calculation_method_ref', self.sectionIdx['section_method'])
        backend.addValue('single_configuration_calculation_to_system_ref', self.sectionIdx['section_system'])
294
295
        # extract k band structure data if available
        self.create_section_eigenvalues(backend, section)
296
297
298
299
300
301
302
        if section['x_qe_t_energy_decomposition_name'] is not None:
            backend.addArrayValues('x_qe_energy_decomposition_name', np.asarray(
                section['x_qe_t_energy_decomposition_name']))
        if section['x_qe_t_energy_decomposition_value'] is not None:
            backend.addArrayValues('x_qe_energy_decomposition_value', np.asarray(
                section['x_qe_t_energy_decomposition_value']))
        if section['x_qe_t_force_x'] is not None:
Henning Glawe's avatar
Henning Glawe committed
303
304
            # constraints etc. not part of the reported forces, so correct metaInfo is 'atom_forces_raw'
            backend.addArrayValues('atom_forces_raw', np.array([
305
306
                section['x_qe_t_force_x'], section['x_qe_t_force_y'], section['x_qe_t_force_z']
                ]).T)
307
308
309
310
        if section['x_qe_t_dispersion_force_x'] is not None:
            backend.addArrayValues('x_qe_atom_dispersion_force', np.array([
                section['x_qe_t_dispersion_force_x'], section['x_qe_t_dispersion_force_y'], section['x_qe_t_dispersion_force_z']
                ]).T)
311
312
313
314
        if section['x_qe_t_stress_x'] is not None:
            backend.addArrayValues('stress_tensor', np.array([
                section['x_qe_t_stress_x'], section['x_qe_t_stress_y'], section['x_qe_t_stress_z']
                ]).T)
315
        had_energy_reference = (self.section['section_system']['number_of_electrons'] is not None)
316
        HOMO = section['x_qe_t_energy_reference_highest_occupied']
317
        if HOMO is not None:
318
319
            if len(HOMO)>1:
                LOGGER.error('more than one value for HOMO: %s', str(HOMO))
320
            had_energy_reference = True
321
322
            backend.addArrayValues('energy_reference_highest_occupied', np.asarray([HOMO[0]]))
        LUMO = section['x_qe_t_energy_reference_lowest_unoccupied']
323
        if LUMO is not None:
324
325
326
            if len(LUMO)>1:
                LOGGER.error('more than one value for LUMO: %s', str(LUMO))
            backend.addArrayValues('energy_reference_lowest_unoccupied', np.asarray([LUMO[0]]))
Henning Glawe's avatar
Henning Glawe committed
327
        E_Fermi_up = section['x_qe_t_energy_reference_fermi_up']
328
        E_Fermi = section['x_qe_t_energy_reference_fermi']
Henning Glawe's avatar
Henning Glawe committed
329
330
331
332
333
334
335
336
        if E_Fermi_up is not None:
            if len(E_Fermi_up)>1:
                LOGGER.error('more than one value for E_Fermi_up: %s', str(E_Fermi_up))
            backend.addArrayValues('energy_reference_fermi', np.asarray([
                E_Fermi_up[0], section['x_qe_t_energy_reference_fermi_down'][0]
            ]))
            had_energy_reference = True
        elif E_Fermi is not None:
337
338
339
            if len(E_Fermi)>1:
                LOGGER.error('more than one value for E_Fermi: %s', str(E_Fermi))
            backend.addArrayValues('energy_reference_fermi', np.asarray([E_Fermi[0]]))
340
341
            had_energy_reference = True
        if not (had_energy_reference):
342
            LOGGER.error("Neither HOMO, Fermi energy nor number of electrons are defined")
343
344
345
346
347
348
        # setup section_system for next scf
        if section['x_qe_t_md_atom_labels'] or section['x_qe_t_md_vec_a_x']:
            next_system_gIndex = backend.openSection('section_system')
            # we cannot simply do
            #   backend.addValue('x_qe_t_vec_a_x', section['x_qe_t_md_vec_a_x'])
            # as this adds an outer list with one element (WTF)
349
            new_system = {}
350
            if section['x_qe_t_md_vec_a_units']:
351
                # we got new cell vectors
352
                new_system['x_qe_t_vec_a_units'] = section['x_qe_t_md_vec_a_units']
353
354
355
                new_system['x_qe_t_vec_a_x'] = section['x_qe_t_md_vec_a_x']
                new_system['x_qe_t_vec_a_y'] = section['x_qe_t_md_vec_a_y']
                new_system['x_qe_t_vec_a_z'] = section['x_qe_t_md_vec_a_z']
Henning Glawe's avatar
Henning Glawe committed
356
357
358
359
360
361
362
            if section['x_qe_t_md_atom_labels']:
                # we got new atom positions and labels
                new_system['x_qe_t_atom_labels'] = section['x_qe_t_md_atom_labels']
                new_system['x_qe_t_atpos_units'] = section['x_qe_t_md_atom_positions_units']
                new_system['x_qe_t_atpos_x'] = section['x_qe_t_md_atom_positions_x']
                new_system['x_qe_t_atpos_y'] = section['x_qe_t_md_atom_positions_y']
                new_system['x_qe_t_atpos_z'] = section['x_qe_t_md_atom_positions_z']
363
364
365
366
367
368
369
370
            if section['x_qe_t_md_k_info_vec_x']:
                # we got new atom positions and labels
                new_system['x_qe_t_k_info_wk'] = section['x_qe_t_md_k_info_wk']
                new_system['x_qe_t_k_info_vec_x'] = section['x_qe_t_md_k_info_vec_x']
                new_system['x_qe_t_k_info_vec_y'] = section['x_qe_t_md_k_info_vec_y']
                new_system['x_qe_t_k_info_vec_z'] = section['x_qe_t_md_k_info_vec_z']
            if section['x_qe_t_md_k_info_ik']:
                new_system['x_qe_t_k_info_ik'] = section['x_qe_t_md_k_info_ik']
371
            for target, data in new_system.items():
372
373
374
                for val in data:
                    backend.addValue(target, val)
            backend.closeSection('section_system', next_system_gIndex)
375

376
377
378
    def onClose_section_scf_iteration(
            self, backend, gIndex, section):
        """trigger called when section_scf_iteration is closed"""
379
380
381
382
383
        if section["x_qe_t_iter_mpersite_idx"] is not None:
            backend.addArrayValues("x_qe_iter_mpersite_idx", np.asarray(section["x_qe_t_iter_mpersite_idx"]))
            backend.addArrayValues("x_qe_iter_mpersite_charge", np.asarray(section["x_qe_t_iter_mpersite_charge"]))
            backend.addArrayValues("x_qe_iter_mpersite_magn", np.asarray(section["x_qe_t_iter_mpersite_magn"]))
            backend.addArrayValues("x_qe_iter_mpersite_constr", np.asarray(section["x_qe_t_iter_mpersite_constr"]))
384
        self.tmp['last_iteration'] = section['x_qe_iteration_number'][-1]
385

386
387
388
    def create_section_eigenvalues(self, backend, src_sec):
        if src_sec['x_qe_t_k_x'] is None or len(src_sec['x_qe_t_k_x']) < 1:
            LOGGER.error("no k-points, not creating section_eigenvalues!")
Henning Glawe's avatar
Henning Glawe committed
389
            return
390
        sec_eigenvalues_gIndex = backend.openSection('section_eigenvalues')
Henning Glawe's avatar
Henning Glawe committed
391
        # prepare numpy arrays
392
393
        k_energies = np.array([self.tmp['k_energies']], dtype=np.float64)
        k_energies = unit_conversion.convert_unit(k_energies, 'eV')
394
        k_occupations = np.array([self.tmp['k_occupations']], dtype=np.float64)
395
396
397
        npw = None
        if src_sec['x_qe_t_k_pw'] is not None:
            npw = np.array(src_sec['x_qe_t_k_pw'])
Henning Glawe's avatar
Henning Glawe committed
398
        k_point_cartesian = np.array([
399
            src_sec['x_qe_t_k_x'], src_sec['x_qe_t_k_y'], src_sec['x_qe_t_k_z']
Henning Glawe's avatar
Henning Glawe committed
400
401
402
403
404
405
406
407
        ], dtype=np.float64).T
        # check if we are dealing with spin-polarized data
        #   QE represents this as 2*k-points, with repeating coordinates
        nk = len(k_point_cartesian)
        if (nk & 1):
            # odd number of k points cannot describe spin-polarized case
            nspin = 1
        else:
408
            kd = (k_point_cartesian[0:nk//2,:] - k_point_cartesian[nk//2:,:])
Henning Glawe's avatar
Henning Glawe committed
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
            kd_len = np.sqrt(np.einsum('ij,ij->i', kd, kd))
            # LOGGER.error("kdl: %s", str(kd_len))
            # difference in k coordinates
            # values in 1/m are in the order of 1e+10, so threshold is sufficient
            if np.all(kd_len<0.1):
                # first half of k coodinates same as second half -> spin
                # polarized case
                nspin = 2
            else:
                nspin = 1
        # sanity checks
        len_kspin = len(self.tmp['kspin'])
        if (nspin == 2) and (len_kspin!=2):
            raise Exception("nspin=2 and kspin!=2")
        elif (nspin == 1) and (len_kspin>1):
            raise Exception("nspin=1 and kspin>1")
        elif (len_kspin>2):
            raise Exception("total fuckup: kspin=%d" % (len_kspin))
        # transform spin-polarized case
        if nspin == 2:
429
            k_point_cartesian[0:nk//2,:]
430
            if npw is not None:
431
                npw = npw[0:nk//2]
Henning Glawe's avatar
Henning Glawe committed
432
433
            # put spin channel into first dimension
            k_energies = np.concatenate((
434
435
                k_energies[:,0:nk//2,:],
                k_energies[:,nk//2:,:]), axis=0)
436
            k_occupations = np.concatenate((
437
438
                k_occupations[:,0:nk//2,:],
                k_occupations[:,nk//2:,:]), axis=0)
Henning Glawe's avatar
Henning Glawe committed
439
440
441
        # k-points are in cartesian, but metaInfo specifies crystal
        k_point_crystal = self.bmat_inv.dot(k_point_cartesian.T).T
        # emit data
442
443
        if npw is not None:
            backend.addArrayValues('x_qe_eigenvalues_number_of_planewaves', npw)
444
445
        if k_occupations.shape[2] > 0:
            backend.addArrayValues('eigenvalues_occupation', k_occupations)
Henning Glawe's avatar
Henning Glawe committed
446
447
        backend.addArrayValues('eigenvalues_kpoints', k_point_crystal)
        backend.addArrayValues('eigenvalues_values', k_energies)
448
        backend.closeSection('section_eigenvalues', sec_eigenvalues_gIndex)
Henning Glawe's avatar
Henning Glawe committed
449

Henning Glawe's avatar
Henning Glawe committed
450
    def onClose_section_system(self, backend, gIndex, section):
451
        old_system = self.section.get('section_system', None)
452
453
        if old_system is not None and self.tmp.get('md_relax', None) is None:
            raise Exception('encountered new section_system without knowing why')
454
455
        if section['number_of_atoms']:
            self.number_of_atoms = section['number_of_atoms'][-1]
456

457
        # store direct lattice matrix and inverse for transformation crystal <-> cartesian
458
        if section['x_qe_t_vec_a_x'] is not None:
459
460
            # there may have been more than one instance of the bravais matrix
            # in different units (seen in relax/MD)
461
            self.amat = np.array([
462
                section['x_qe_t_vec_a_x'][-3:], section['x_qe_t_vec_a_y'][-3:], section['x_qe_t_vec_a_z'][-3:],
463
            ], dtype=np.float64).T
464
465
466
467
468
469
470
471
472
473
474
            # explicit unit conversion
            amat_units = section['x_qe_t_vec_a_units'][-1]
            if amat_units == 'a_0' or amat_units == 'alat':
                self.amat = unit_conversion.convert_unit(self.amat, 'usrAlat')
            elif amat_units == 'bohr':
                self.amat = unit_conversion.convert_unit(self.amat, 'bohr')
            elif amat_units == 'angstrom':
                self.amat = unit_conversion.convert_unit(self.amat, 'angstrom')
            else:
                raise RuntimeError("unknown amat_units: %s" % (amat_units))

475
476
477
478
            # store inverse for transformation cartesian -> crystal
            try:
                self.amat_inv = np.linalg.inv(self.amat)
            except np.linalg.linalg.LinAlgError:
479
                raise Exception("error inverting bravais matrix " + str(self.amat))
480
            LOGGER.info('NewCell')
481
482
        elif old_system is not None:
            # we did not get new cell vectors, recycle old ones
483
            LOGGER.info('OldCell')
484
        else:
485
486
            raise Exception("missing bravais vectors")
        backend.addArrayValues('simulation_cell', self.amat)
487
        self.tmp['have_new_system'] = True
488
489

        # store reciprocal lattice matrix and inverse for transformation crystal <-> cartesian
490
491
492
493
        if section['x_qe_t_vec_b_x'] is not None:
            self.bmat = np.array([
                section['x_qe_t_vec_b_x'], section['x_qe_t_vec_b_y'], section['x_qe_t_vec_b_z'],
            ], dtype=np.float64).T
494
495
496
497
498
499
500
501
            # store inverse for transformation cartesian -> crystal
            try:
                self.bmat_inv = np.linalg.inv(self.bmat)
            except np.linalg.linalg.LinAlgError:
                raise Exception("error inverting reciprocal cell matrix")
        elif section['x_qe_t_vec_a_x'] is not None:
            # we got new lattice vectors, but no reciprocal ones, calculate
            # on-the-fly
Henning Glawe's avatar
Henning Glawe committed
502
            LOGGER.info('calculating bmat on the fly from amat')
503
504
505
506
507
508
            abmat = np.zeros((3,3), dtype=np.float64)
            abmat[0] = np.cross(self.amat[1],self.amat[2])
            abmat[1] = np.cross(self.amat[2],self.amat[0])
            abmat[2] = np.cross(self.amat[0],self.amat[1])
            abmat *= 2*math.pi / np.dot(abmat[0],self.amat[0])
            self.bmat = abmat
509
510
511
512
513
514
515
516
            # store inverse for transformation cartesian -> crystal
            try:
                self.bmat_inv = np.linalg.inv(self.bmat)
            except np.linalg.linalg.LinAlgError:
                raise Exception("error inverting reciprocal cell matrix")
        elif old_system is not None:
            # keep what we had
            pass
517
        else:
518
            raise Exception("missing reciprocal cell vectors")
519
        backend.addArrayValues('x_qe_reciprocal_cell', self.bmat)
520
        # atom positions
521
        if section['x_qe_t_atpos_x'] is not None:
522
523
            # there may have been more than one instance of the atom positions
            # in different units (seen in relax/MD)
524
            atpos = np.array([
525
526
527
                section['x_qe_t_atpos_x'][-self.number_of_atoms:],
                section['x_qe_t_atpos_y'][-self.number_of_atoms:],
                section['x_qe_t_atpos_z'][-self.number_of_atoms:]
528
            ], dtype=np.float64).T
529
530
531
            atpos_units = section['x_qe_t_atpos_units'][-1]
            if atpos_units == 'a_0' or atpos_units == 'alat':
                atpos_cart = unit_conversion.convert_unit(atpos, 'usrAlat')
532
533
            elif atpos_units == 'bohr':
                atpos_cart = unit_conversion.convert_unit(atpos, 'bohr')
534
535
            elif atpos_units == 'angstrom':
                atpos_cart = unit_conversion.convert_unit(atpos, 'angstrom')
536
            elif atpos_units == 'crystal' or atpos_units == 'cryst. coord.':
537
                atpos_cart = self.amat.dot(atpos.T).T
538
539
            else:
                raise RuntimeError("unknown atpos_units: %s" % (atpos_units))
540
            LOGGER.info('NewAtpos')
541
        elif old_system is not None:
542
            atpos_cart = old_system['atom_positions'][-1]
543
            LOGGER.info('OldAtpos')
544
545
546
547
548
        else:
            raise Exception("missing atom positions")
        backend.addArrayValues('atom_positions',atpos_cart)

        if section['x_qe_t_atom_labels'] is not None:
549
            backend.addArrayValues('atom_labels',np.asarray(section['x_qe_t_atom_labels'][-self.number_of_atoms:]))
550
        elif old_system is not None:
551
            backend.addArrayValues('atom_labels',old_system['atom_labels'][-1])
552
553
554
        else:
            raise Exception("missing atom labels")

555
        if section['x_qe_t_atom_idx'] is not None:
556
            backend.addArrayValues('x_qe_atom_idx', np.asarray(section['x_qe_t_atom_idx'][-self.number_of_atoms:]))
557
        elif old_system is not None:
558
            backend.addArrayValues('x_qe_atom_idx', old_system['x_qe_atom_idx'][-1])
559
        else:
560
561
562
563
564
565
566
567
568
569
570
            raise Exception("missing x_qe_atom_idx")

        if section['x_qe_t_starting_magnetization_species'] is not None:
            # build dict with per-species magnetization
            sp_magn = {}
            for (label, magn) in zip(
                    section['x_qe_t_starting_magnetization_species'],
                    section['x_qe_t_starting_magnetization_value']):
                sp_magn[label] = magn
            at_magn = []
            # transform to per-atom magnetization
571
            for label in section['atom_labels'][-1]:
572
573
574
                at_magn.append(sp_magn[label])
            backend.addArrayValues('x_qe_atom_starting_magnetization',np.array(at_magn))

575
576
577
578
579
580
581
        if section['x_qe_t_celldm'] is not None:
            celldm_joint = " ".join(section['x_qe_t_celldm'])
            celldm = [None, None, None, None, None, None]
            for match in re.findall(r"celldm\(\s*(\d+)\s*\)\s*=\s*(" + RE_f + r")", celldm_joint):
                celldm[int(match[0])-1] = valueForStrValue(match[1], 'f')
            celldm[0] = self.alat
            backend.addArrayValues('x_qe_celldm', np.array(celldm))
582

583
584
585
586
        if section['x_qe_t_k_info_vec_x'] is not None:
            backend.addArrayValues('x_qe_k_info_vec', np.array([
                section['x_qe_t_k_info_vec_x'], section['x_qe_t_k_info_vec_y'], section['x_qe_t_k_info_vec_z']
            ]).T)
587
588
589
        elif old_system is not None:
            # unless espresso explicitly writes new k-points, sampling is kept fixed
            backend.addArrayValues('x_qe_k_info_vec', old_system['x_qe_k_info_vec'][-1])
590
        else:
591
            LOGGER.error("No K-point info found in output")
592

593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
        if section['x_qe_t_k_info_ik'] is not None:
            backend.addArrayValues('x_qe_k_info_ik', np.array(section['x_qe_t_k_info_ik']))
        elif old_system is not None:
            # unless espresso explicitly writes new k-points, sampling is kept fixed
            backend.addArrayValues('x_qe_k_info_ik', old_system['x_qe_k_info_ik'][-1])
        else:
            LOGGER.error("No K-point index info found in output")

        if section['x_qe_t_k_info_wk'] is not None:
            backend.addArrayValues('x_qe_k_info_wk', np.array(section['x_qe_t_k_info_wk']))
        elif old_system is not None:
            # unless espresso explicitly writes new k-points, sampling is kept fixed
            backend.addArrayValues('x_qe_k_info_wk', old_system['x_qe_k_info_wk'][-1])
        else:
            LOGGER.error("No K-point weight info found in output")


610
611
612
613
        if section['x_qe_t_dense_FFT_grid_x'] is not None:
            backend.addArrayValues('x_qe_dense_FFT_grid', np.array([
                section['x_qe_t_dense_FFT_grid_x'], section['x_qe_t_dense_FFT_grid_y'], section['x_qe_t_dense_FFT_grid_z']
            ]).T)
614
615
616
        elif old_system is not None:
            # unless espresso explicitly writes new FFT grid info, sampling is kept fixed
            backend.addArrayValues('x_qe_dense_FFT_grid', old_system['x_qe_dense_FFT_grid'][-1])
617
        else:
618
619
            LOGGER.warning("No dense FFT grid info found in output")

Henning Glawe's avatar
Henning Glawe committed
620
621
622
623
        if section['x_qe_t_smooth_FFT_grid_x'] is not None:
            backend.addArrayValues('x_qe_smooth_FFT_grid', np.array([
                section['x_qe_t_smooth_FFT_grid_x'], section['x_qe_t_smooth_FFT_grid_y'], section['x_qe_t_smooth_FFT_grid_z']
            ]).T)
624
625
626
        elif old_system is not None and old_system['x_qe_smooth_FFT_grid'] is not None:
            backend.addArrayValues('x_qe_smooth_FFT_grid', old_system['x_qe_smooth_FFT_grid'][-1])

Henning Glawe's avatar
Henning Glawe committed
627
628
629
630
        if section['x_qe_t_vec_supercell_x'] is not None:
            backend.addArrayValues('x_qe_vec_supercell', np.array([
                section['x_qe_t_vec_supercell_x'], section['x_qe_t_vec_supercell_y'], section['x_qe_t_vec_supercell_z']
            ]).T)
631
632

        if section['x_qe_t_number_of_electrons_up'] is not None:
633
            # spin polarized case, with explicit up/down electrons
634
635
636
637
638
639
640
641
642
643
            if len(section['x_qe_t_number_of_electrons_up'])>1:
                LOGGER.error("got multiple nelec_up: %s", str(
                    section['x_qe_t_number_of_electrons_up']))
            backend.addArrayValues('number_of_electrons', np.array([
                section['x_qe_t_number_of_electrons_up'][0],
                section['x_qe_t_number_of_electrons_down'][0]
            ]))
        elif section['x_qe_t_number_of_electrons'] is not None:
            if len(section['x_qe_t_number_of_electrons'])!=1:
                LOGGER.error("got wrong nelec: %s", str(section['x_qe_t_number_of_electrons']))
644
            backend.addArrayValues('number_of_electrons', np.array([
645
                section['x_qe_t_number_of_electrons'][-1]
646
            ]))
647
648
        elif old_system is not None:
            backend.addArrayValues('number_of_electrons', old_system['number_of_electrons'][-1])
649
        else:
650
651
            raise Exception("missing info about number of electrons in system")

652
        backend.addArrayValues('configuration_periodic_dimensions', np.asarray([True, True, True]))
653
654
        self.sectionIdx['section_system'] = gIndex
        self.section['section_system'] = section
Henning Glawe's avatar
Henning Glawe committed
655

656
657
658
659
660
661
    def onOpen_section_run(
            self, backend, gIndex, section):
        """trigger called when section_single_configuration_calculation
        is closed"""
        self.tmp.pop('x_qe_t_profile_caller', None)
        self.tmp.pop('x_qe_t_profile_category', None)
662
663
664
665
        # manually open header sections, closed at the beginning of scf
        for sec in self.header_sections():
            gIndex = backend.openSection(sec)
            self.openSectionIdx[sec] = gIndex
666

667
668
669
670
671
672
673
    def adHoc_final_scf_MD(self, parser):
        """final SCF calculation in VC-relax runs needs open header sections"""
        # manually open header sections, closed at the beginning of scf
        for sec in self.header_sections():
            gIndex = parser.backend.openSection(sec)
            self.openSectionIdx[sec] = gIndex

Henning Glawe's avatar
Henning Glawe committed
674
675
676
677
    def onClose_section_run(
            self, backend, gIndex, section):
        """trigger called when section_single_configuration_calculation
        is closed"""
678
679
680
681
682
683
684
685
686
        if section['x_qe_t_profile_function'] is not None:
            backend.addArrayValues('x_qe_profile_function', np.asarray(
                section['x_qe_t_profile_function']))
            backend.addArrayValues('x_qe_profile_cputime', np.asarray(
                section['x_qe_t_profile_cputime']))
            backend.addArrayValues('x_qe_profile_walltime', np.asarray(
                section['x_qe_t_profile_walltime']))
            backend.addArrayValues('x_qe_profile_ncalls', np.asarray(
                section['x_qe_t_profile_ncalls']))
687
688
689
690
            backend.addArrayValues('x_qe_profile_category', np.asarray(
                section['x_qe_t_profile_category_list']))
            backend.addArrayValues('x_qe_profile_caller', np.asarray(
                section['x_qe_t_profile_caller_list']))
Henning Glawe's avatar
Henning Glawe committed
691

692
693
694
695
696
697
698
699
700
701
        frames = self.tmp.get('frames', None)
        if frames:
            sampling_method_gIndex = backend.openSection('section_sampling_method')
            backend.addValue('sampling_method', QE_MD_RELAX_SAMPLING_METHOD[self.tmp['md_relax']])
            backend.closeSection('section_sampling_method', sampling_method_gIndex)
            frame_sequence_gIndex = backend.openSection('section_frame_sequence')
            backend.addValue('frame_sequence_to_sampling_ref', sampling_method_gIndex)
            backend.addArrayValues('frame_sequence_local_frames_ref', np.array(self.tmp['frames']))
            backend.closeSection('section_frame_sequence', frame_sequence_gIndex)

702
703
704
    def appendToTmp(self, tmpname, value):
        self.tmp[tmpname] += value

705
706
707
    def setTmp(self, tmpname, value):
        self.tmp[tmpname] = value

Henning Glawe's avatar
Henning Glawe committed
708
709
710
711
    def setTmpUnlessExists(self, tmpname, value):
        if self.tmp.get(tmpname,None) is None:
            self.tmp[tmpname] = value

Henning Glawe's avatar
Henning Glawe committed
712
713
714
    def popTmp(self, tmpname, fallback=None):
        return self.tmp.pop(tmpname, fallback)

715
716
717
718
719
720
721
722
    def adHoc_dispersion_correction_values(self, parser):
        if 'dispersion_correction' not in self.tmp:
            self.tmp['dispersion_correction'] = {}
        self.tmp['dispersion_correction'][parser.lastMatch['x_qe_t_species_dispersion_correction_label']] = {
            'x_qe_dispersion_correction_vdw_radius': parser.lastMatch['x_qe_t_species_dispersion_correction_vdw_radius'],
            'x_qe_dispersion_correction_C6': parser.lastMatch['x_qe_t_species_dispersion_correction_C6'],
        }

723
    def adHoc_alat(self, parser):
724
        alat = parser.lastMatch['x_qe_alat']
Henning Glawe's avatar
Henning Glawe committed
725
        self.alat = parser.lastMatch['x_qe_alat']
726
727
        unit_conversion.register_userdefined_quantity('usrAlat', 'm', alat)
        unit_conversion.register_userdefined_quantity('usrTpiba', '1/m', 2*math.pi/alat)
728

729
730
731
732
    def adHoc_pp_renorm(self, parser):
        self.cache_t_pp_renorm_wfc[parser.lastMatch[
            'x_qe_t_pp_renormalized_filename']] = parser.lastMatch['x_qe_t_pp_renormalized_wfc']

733
    def adHoc_bands_spin(self, parser):
734
        k_x = self.section['single_configuration_calculation']['x_qe_t_k_x']
735
736
737
738
739
740
        if k_x is None:
            nk_current = 0
        else:
            nk_current = len(k_x)
        self.tmp['kspin'][parser.lastMatch['x_qe_t_spin_channel'].lower()] = nk_current

741
    def adHoc_profiling_complete(self, parser):
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
        # we have effectively 3 SMs in one, split between the 3 cases
        if parser.lastMatch.get('x_qe_t_profile_caller', None) is not None:
            # caller info
            self.setTmp('x_qe_t_profile_caller', parser.lastMatch['x_qe_t_profile_caller'])
        elif parser.lastMatch.get('x_qe_t_profile_category', None) is not None:
            # category info
            self.setTmp('x_qe_t_profile_category', parser.lastMatch['x_qe_t_profile_category'])
            self.tmp.pop('x_qe_t_profile_caller', None)
        else:
            # actual profiling data
            if parser.lastMatch.get('x_qe_t_profile_cputime', None) is None:
                parser.backend.addValue('x_qe_t_profile_cputime', QeC.NAN)
            if parser.lastMatch.get('x_qe_t_profile_walltime', None) is None:
                parser.backend.addValue('x_qe_t_profile_walltime', QeC.NAN)
            if parser.lastMatch.get('x_qe_t_profile_ncalls', None) is None:
                parser.backend.addValue('x_qe_t_profile_ncalls', QeC.NAN)
            parser.backend.addValue('x_qe_t_profile_caller_list', self.tmp.get('x_qe_t_profile_caller', ''))
            parser.backend.addValue('x_qe_t_profile_category_list', self.tmp.get('x_qe_t_profile_category', ''))
Henning Glawe's avatar
Henning Glawe committed
760

761
    def SMs_summaryf90(self):
762
        return [
763
764
765
            SM(name='ibrav', required=True,
               # first line printed by summary.f90
               startReStr=r"\s*bravais-lattice index\s*=\s*(?P<x_qe_ibrav>" + RE_i + r")\s*$",
766
               subMatchers=[
767
768
769
770
                   # other info coming from/triggered by summary.f90
                   SM(name='alat', required=True,
                      startReStr=r"\s*lattice parameter \((?:a_0|alat)\)\s*=\s*(?P<x_qe_alat__bohr>" + RE_f + r")\s*a\.u\.\s*$",
                      adHoc=self.adHoc_alat,
771
                   ),
772
773
                   SM(name='cell_volume', required=True,
                      startReStr=r"\s*unit-cell volume\s*=\s*(?P<x_qe_cell_volume__bohr3>" + RE_f + r")\s*\(a\.u\.\)\^3\s*$",
774
                   ),
775
776
                   SM(name='nat', required=True,
                      startReStr=r"\s*number of atoms/cell\s*=\s*(?P<number_of_atoms>\d+)\s*$",
777
                   ),
778
779
                   SM(name='nsp', required=True,
                      startReStr=r"\s*number of atomic types\s*=\s*(?P<x_qe_number_of_species>" + RE_i + r")\s*$",
780
                   ),
781
782
783
784
                   SM(name='nelec', required=True,
                      startReStr=(r"\s*number of electrons\s*=\s*(?P<x_qe_t_number_of_electrons>" + RE_f +
                                  r")(?:\s*\(up:\s*(?P<x_qe_t_number_of_electrons_up>" + RE_f +
                                  r")\s*,\s*down:\s*(?P<x_qe_t_number_of_electrons_down>" + RE_f + ")\s*\))?\s*$"),
785
                   ),
786
787
                   SM(name='nbnd', required=True,
                       startReStr=r"\s*number of Kohn-Sham states\s*=\s*(?P<x_qe_number_of_states>" + RE_i + r")\s*$"
788
                   ),
789
790
791
                   SM(name='ecutwfc', required=True, sections=['section_basis_set_cell_dependent'],
                      startReStr=r"\s*kinetic-energy cutoff\s*=\s*(?P<basis_set_planewave_cutoff__rydberg>" + RE_f + r")\s*Ry\s*$",
                      adHoc=lambda p: self.setTmp('method_basis_set_kind', 'wavefunction'),
792
                   ),
793
794
795
                   SM(name='ecut_density', required=True, sections=['section_basis_set_cell_dependent'],
                      startReStr=r"\s*charge density cutoff\s*=\s*(?P<basis_set_planewave_cutoff__rydberg>" + RE_f + r")\s*Ry\s*$",
                      adHoc=lambda p: self.setTmp('method_basis_set_kind', 'density'),
796
                   ),
797
798
                   SM(name='ecutfock',
                      startReStr=r"\s*cutoff for Fock operator\s*=\s*(?P<x_qe_fock_operator_cutoff__rydberg>" + RE_f + r")\s*Ry\s*$"
799
                   ),
800
                   SM(name='convergence_threshold',
801
                      startReStr=r"\s*convergence threshold\s*=\s*(?P<scf_threshold_energy_change__rydberg>" + RE_f + r")\s*$",
802
                   ),
803
804
                   SM(name='mixing_beta',
                      startReStr=r"\s*mixing beta\s*=\s*(?P<x_qe_potential_mixing_beta>" + RE_f + r")\s*$",
805
                   ),
806
807
                   SM(name='mixing_scheme',
                      startReStr=r"\s*number of iterations used\s*=\s*(?P<x_qe_potential_mixing_iterations>\d+)\s*(?P<x_qe_potential_mixing_scheme>.*?)\s+mixing\s*$",
808
                   ),
809
810
                   SM(name='xc_functional', required=True,
                      startReStr=r"\s*Exchange-correlation\s*=\s*(?P<x_qe_xc_functional_shortname>.*?)\s*\((?P<x_qe_xc_functional_num>[^\)]*)\s*\)\s*$"
811
                   ),
812
813
                   SM(name='exx_fraction',
                      startReStr=r"\s*EXX-fraction\s*=\s*(?P<x_qe_t_exact_exchange_fraction>" + RE_f + r")\s*$",
814
                   ),
815
                   SM(name='nstep',
Henning Glawe's avatar
Henning Glawe committed
816
                      startReStr=r"\s*nstep\s*=\s*(?P<x_qe_md_max_steps>" + RE_i + r")\s*$",
817
                   ),
818
819
                   SM(name='spin_orbit_mode',
                      startReStr=r"\s*(?P<x_qe_t_spin_orbit_magn>.*?)\s*calculation\s*(?P<x_qe_t_spin_orbit_mode>with(?:out)?)\s*spin-orbit\s*$",
820
                   ),
821
822
                   SM(name='berry_efield',
                      startReStr=r"\s*Using Berry phase electric field\s*$",
823
                      subMatchers=[
824
825
                          SM(name='berry_efield_direction',
                             startReStr=r"\s*Direction\s*:\s*(?P<x_qe_berry_efield_direction>\d+)\s*$",
826
                          ),
827
828
829
830
831
832
833
834
835
836
837
838
                          SM(name='berry_efield_intensity',
                             # Ry unit is not printed in 4.0
                             startReStr=(r"\s*Intensity \((?:Ry\s*)?a.u.\)\s*:\s*(?P<x_qe_berry_efield_intensity__rydberg>" + RE_f +
                                         r")\s*$"),
                          ),
                          SM(name='berry_efield_strings',
                             startReStr=(r"\s*Strings composed by\s*:\s*(?P<x_qe_berry_efield_strings_nk>" + RE_i +
                                         r")\s*k-points\s*$"),
                          ),
                          SM(name='berry_efield_niter',
                             startReStr=(r"\s*Number of iterative cycles:\s*(?P<x_qe_berry_efield_niter>" + RE_i +
                                         r")\s*$"),
839
840
                          ),
                      ],
841
                      adHoc=lambda p: p.backend.addValue('x_qe_berry_efield', True),
842
                   ),
843
844
                   SM(name='assume_isolated',
                      startReStr=r"\s*Assuming isolated system,\s*(?P<x_qe_isolated_system_method>.*?)\s*method",
845
                   ),
846
847
                   SM(name='celldm', repeats = True,
                      startReStr=r"(?P<x_qe_t_celldm>(?:\s*celldm\(\d+\)\s*=\s*" + RE_f + r")+)\s*$",
848
                   ),
849
                   SM(name='simulation_cell',
850
                      startReStr=r"\s*crystal axes: \(cart. coord.\s*in units of (?P<x_qe_t_vec_a_units>a_0|alat)\s*\)\s*$",
851
852
                      subMatchers=[
                          SM(name='cell_vec_a', repeats=True,
853
                             startReStr=r"\s*a\(\d\)\s*=\s*\(\s*" + QeC.re_vec('x_qe_t_vec_a') + r"\s*\)\s*$",
854
855
                          ),
                      ],
856
                   ),
857
858
859
860
861
862
863
                   SM(name='reciprocal_cell',
                      startReStr=r"\s*reciprocal axes: \(cart. coord. in units 2 pi/(?:alat|a_0)\)\s*$",
                      subMatchers=[
                          SM(name='cell_vec_b', repeats=True,
                             startReStr=r"\s*b\(\d\)\s*=\s*\(\s*" + QeC.re_vec('x_qe_t_vec_b', 'usrTpiba') + r"\s*\)\s*$",
                          ),
                      ],
864
                   ),
865
866
867
868
                   SM(name='pseudopotential', repeats=True,
                      startReStr=(r"\s*PseudoPot\.\s*#\s*(?P<x_qe_t_pp_idx>" + RE_i + r") for\s+(?P<x_qe_t_pp_label>\S+)\s+read from file" +
                                  r"(?::|\s*(?P<x_qe_t_pp_filename>\S+))\s*"),
                      sections=['x_qe_t_section_pseudopotential'],
Henning Glawe's avatar
Henning Glawe committed
869
                      subMatchers=[
870
871
                          SM(name='new_pp_filename',
                             startReStr=r"^\s*(?P<x_qe_t_pp_filename>\S+)\s*$",
872
                          ),
873
874
                          SM(name='pp_md5',
                             startReStr=r"\s*MD5 check sum:\s*(?P<x_qe_t_pp_md5sum>\S+)\s*$",
875
                          ),
876
877
                          SM(name='pp_type_val',
                             startReStr=r"\s*Pseudo is\s*(?P<x_qe_t_pp_type>.*?),\s*Zval\s*=\s*(?P<x_qe_t_pp_valence>" + RE_f + r")\s*$",
Henning Glawe's avatar
Henning Glawe committed
878
                          ),
879
880
                          SM(name='pp_comment',
                             startReStr=r"\s*(?P<x_qe_t_pp_comment>.*?)\s*$",
Henning Glawe's avatar
Henning Glawe committed
881
                          ),
882
883
884
885
                          SM(name='pp_integral_directions',
                             startReStr=(r"\s*Setup to integrate on\s*" +
                                         r"(?P<x_qe_t_pp_integral_ndirections>\d+)\s+directions:\s*" +
                                         r"integral exact up to l =\s*(?P<x_qe_t_pp_integral_lmax_exact>\d+)\s*$"),
Henning Glawe's avatar
Henning Glawe committed
886
                          ),
887
888
                          SM(name='pp_augmentation_shape',
                             startReStr=r"\s*Shape of augmentation charge:\s*(?P<x_qe_t_pp_augmentation_shape>.*?)\s*$",
Henning Glawe's avatar
Henning Glawe committed
889
                          ),
890
891
                          SM(name='pp_dimensions',
                             startReStr=r"\s*Using radial grid of\s*(?P<x_qe_t_pp_ndmx>\d+) points,\s*(?P<x_qe_t_pp_nbeta>\d+) beta functions with\s*:\s*$",
892
                          ),
893
894
                          SM(name='pp_beta', repeats=True,
                             startReStr=r"\s*l\(\s*(?P<x_qe_t_pp_l_idx>\d+)\s*\)\s*=\s*(?P<x_qe_t_pp_l>\d+)\s*$",
895
                          ),
896
897
898
899
900
901
902
                          SM(name='pp_coefficients',
                             startReStr=r"\s*Q\(r\) pseudized with\s*(?P<x_qe_t_pp_ncoefficients>\d+)\s*coefficients,\s*rinner\s*=\s*(?P<x_qe_t_rinner>(?:\s*" + RE_f + r")+)\s*$",
                             subMatchers=[
                                 SM(name='pp_coefficients2',
                                    startReStr=r"\s*(?P<x_qe_t_rinner>(?:\s*" + RE_f + r")+)\s*$",
                                 ),
                             ],
903
                          ),
904
905
                          SM(name='pp_coefficients3',
                             startReStr=r"\s*Q\(r\) pseudized with\s*(?P<x_qe_t_pp_ncoefficients>\d+)\s*coefficients\s*$",
906
                          ),
907
                       ],
908
                   ),
909
910
911
912
913
914
915
916
917
                   SM(name='vdw_kernel_table_file',
                      startReStr=r"\s*vdW kernel table read from file .*?\s*$",
                      adHoc=lambda p: LOGGER.info("do sth with vdW kernel table file"),
                      subMatchers=[
                          SM(name='vdw_kernel_table_md5sum',
                             startReStr=r"\s*MD5 check sum:\s*\S+\s*$",
                          ),
                      ],
                   ),
918
919
                   SM(name='pp_atom_kind_map',
                      startReStr=r"\s*atomic species\s+valence\s+mass\s+pseudopotential\s*$",
920
                      subMatchers=[
921
                          SM(name='atom_kind', repeats=True,
922
                             startReStr=r"\s*(?P<method_atom_kind_label>\S+)\s+(?P<method_atom_kind_explicit_electrons>" + RE_f + r")" +
923
924
925
                                        r"\s+(?P<x_qe_kind_mass>" + RE_f + r")\s+(?P<x_qe_pp_label>[^\s\(]+)\s*" +
                                        r"\(\s*(?P<x_qe_pp_weight>" + RE_f + r")\s*\)\s*$",
                             sections=['section_method_atom_kind'],
926
927
928
                          ),
                      ],
                   ),
929
930
                   SM(name='starting_magnetization',
                      startReStr=r"\s*Starting magnetic structure\s*$",
931
                      subMatchers=[
932
933
                          SM(name='starting_magnetization_header',
                             startReStr=r"\s*atomic species\s+magnetization\s*$",
934
                             subMatchers=[
935
                                 SM(name='starting_magnetization_data', repeats=True,
936
                                    startReStr=(r"\s*(?P<x_qe_t_starting_magnetization_species>\S+)\s*" +
937
                                                r"(?P<x_qe_t_starting_magnetization_value>" + RE_f + r")\s*$"),
938
939
940
941
942
                                 ),
                             ],
                          ),
                      ],
                   ),
943