text: # -*- coding: utf-8 -*-
  """A collection of function that are used to parse the output of Quantum Espresso PW.
  
  The function that needs to be called from outside is parse_raw_output(). The functions mostly work without aiida
  specific functionalities. The parsing will try to convert whatever it can in some dictionary, which by operative
  decision doesn't have much structure encoded, [the values are simple ]
  """
  import re
  import numpy
  from qe_tools import CONSTANTS
  
  from aiida_quantumespresso.parsers import QEOutputParsingError
  from aiida_quantumespresso.parsers.parse_raw import convert_qe_time_to_sec
  from aiida_quantumespresso.utils.mapping import get_logging_container
  
  lattice_tolerance = 1.e-5
  units_suffix = '_units'
  default_charge_units = 'e'
  default_dipole_units = 'Debye'
  default_energy_units = 'eV'
  default_force_units = 'ev / angstrom'
  default_k_points_units = '1 / angstrom'
  default_length_units = 'Angstrom'
  default_magnetization_units = 'Bohrmag / cell'
  default_polarization_units = 'C / m^2'
  default_stress_units = 'GPascal'
  
  
  def reduce_symmetries(parsed_parameters, parsed_structure, logger):
      """Reduce the symmetry information parsed from the output to save space.
  
      In the standard output, each symmetry operation print two rotation matrices:
  
          * S_cryst^T: matrix in crystal coordinates, transposed
          * S_cart: matrix in cartesian coordinates,
  
      The XML files only print one matrix:
  
          * S_cryst: matrix in crystal coordinates
  
      The raw parsed symmetry information from the XML is large and will load the database heavily if stored as
      is for each calculation. Instead, we will map these dictionaries onto a static dictionary of rotation
      matrices generated by the _get_qe_symmetry_list static method. This dictionary will return the rotation
      matrices in cartesian coordinates, i.e. S_cart. In order to compare the raw matrices from the XML to these
      static matrices we have to convert S_cryst into S_cart. We derive here how that is done:
  
          S_cryst * v_cryst = v_cryst'
  
      where v_cryst' is the rotated vector v_cryst under S_cryst
      We define `cell` where cell vectors are rows. Converting a vector from crystal to cartesian
      coordinates is defined as:
  
          cell^T * v_cryst = v_cart
  
      The inverse of this operation is defined as
  
          v_cryst = cell^Tinv * v_cart
  
      Replacing the last equation into the first we find:
  
          S_cryst * cell^Tinv * v_cart = cell^Tinv * v_cart'
  
      Multiply on the left with cell^T gives:
  
          cell^T * S_cryst * cell^Tinv * v_cart = v_cart'
  
      which can be rewritten as:
  
          S_cart * v_cart = v_cart'
  
      where:
  
          S_cart = cell^T * S_cryst * cell^Tinv
  
      We compute here the transpose and its inverse of the structure cell basis, which is needed to transform
      the parsed rotation matrices, which are in crystal coordinates, to cartesian coordinates, which are the
      matrices that are returned by the _get_qe_symmetry_list staticmethod
      """
      from aiida_quantumespresso.parsers.parse_raw.pw import get_symmetry_mapping
      from aiida_quantumespresso.utils.linalg import are_matrices_equal
  
      cell = parsed_structure['cell']['lattice_vectors']
      cell_T = numpy.transpose(cell)
      cell_Tinv = numpy.linalg.inv(cell_T)
      possible_symmetries = get_symmetry_mapping()
  
      for symmetry_type in ['symmetries', 'lattice_symmetries']:  # crystal vs. lattice symmetries
          if symmetry_type in list(parsed_parameters.keys()):
              try:
                  old_symmetries = parsed_parameters[symmetry_type]
                  new_symmetries = []
                  for this_sym in old_symmetries:
                      name = this_sym['name'].strip()
                      for i, this in enumerate(possible_symmetries):
                          # Since we do an exact comparison we strip the string name from whitespace
                          # and as soon as it is matched, we break to prevent it from matching another
                          if name == this['name'].strip():
                              index = i
                              break
                      else:
                          index = None
                          logger.error(f'Symmetry {name} not found')
  
                      new_dict = {}
                      if index is not None:
                          # The raw parsed rotation matrix is in crystal coordinates, whereas the mapped rotation
                          # in possible_symmetries is in cartesian coordinates. To allow them to be compared
                          # to make sure we matched the correct rotation symmetry, we first convert the parsed matrix
                          # to cartesian coordinates. For explanation of the method, see comment above.
                          rotation_cryst = this_sym['rotation']
                          rotation_cart_new = possible_symmetries[index]['matrix']
                          rotation_cart_old = numpy.dot(cell_T, numpy.dot(rotation_cryst, cell_Tinv))
  
                          inversion = possible_symmetries[index]['inversion']
                          if not are_matrices_equal(rotation_cart_old, rotation_cart_new, swap_sign_matrix_b=inversion):
                              logger.error(
                                  'Mapped rotation matrix {} does not match the original rotation {}'.format(
                                      rotation_cart_new, rotation_cart_old
                                  )
                              )
                              new_dict['all_symmetries'] = this_sym
                          else:
                              # Note: here I lose the information about equivalent ions and fractional_translation
                              # since I don't copy them to new_dict (but they can be reconstructed).
                              new_dict['t_rev'] = this_sym['t_rev']
                              new_dict['symmetry_number'] = index
                      else:
                          new_dict['all_symmetries'] = this_sym
  
                      new_symmetries.append(new_dict)
  
                  parsed_parameters[symmetry_type] = new_symmetries  # and overwrite the old one
              except KeyError:
                  logger.warning(f"key '{symmetry_type}' is not present in raw output dictionary")
          else:
              # backwards-compatiblity: 'lattice_symmetries' is not created in older versions of the parser
              if symmetry_type != 'lattice_symmetries':
                  logger.warning(f"key '{symmetry_type}' is not present in raw output dictionary")
  
  
  def get_symmetry_mapping():
      """Hard coded names and rotation matrices + inversion from QE v 5.0.2 function for Parser class usage only.
  
      :return: a list of dictionaries, each containing name (string),
          inversion (boolean) and matrix (list of lists)
      """
      sin3 = 0.866025403784438597
      cos3 = 0.5
      msin3 = -0.866025403784438597
      mcos3 = -0.5
  
      # 32 rotations that are checked + inversion taken from symm_base.f90 from the QE source code
      # They are in Fortran format and therefore transposed with respect to the default python format
      transposed_matrices_cartesian = [
          [[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]],
          [[-1., 0., 0.], [0., -1., 0.], [0., 0., 1.]],
          [[-1., 0., 0.], [0., 1., 0.], [0., 0., -1.]],
          [[1., 0., 0.], [0., -1., 0.], [0., 0., -1.]],
          [[0., 1., 0.], [1., 0., 0.], [0., 0., -1.]],
          [[0., -1., 0.], [-1., 0., 0.], [0., 0., -1.]],
          [[0., -1., 0.], [1., 0., 0.], [0., 0., 1.]],
          [[0., 1., 0.], [-1., 0., 0.], [0., 0., 1.]],
          [[0., 0., 1.], [0., -1., 0.], [1., 0., 0.]],
          [[0., 0., -1.], [0., -1., 0.], [-1., 0., 0.]],
          [[0., 0., -1.], [0., 1., 0.], [1., 0., 0.]],
          [[0., 0., 1.], [0., 1., 0.], [-1., 0., 0.]],
          [[-1., 0., 0.], [0., 0., 1.], [0., 1., 0.]],
          [[-1., 0., 0.], [0., 0., -1.], [0., -1., 0.]],
          [[1., 0., 0.], [0., 0., -1.], [0., 1., 0.]],
          [[1., 0., 0.], [0., 0., 1.], [0., -1., 0.]],
          [[0., 0., 1.], [1., 0., 0.], [0., 1., 0.]],
          [[0., 0., -1.], [-1., 0., 0.], [0., 1., 0.]],
          [[0., 0., -1.], [1., 0., 0.], [0., -1., 0.]],
          [[0., 0., 1.], [-1., 0., 0.], [0., -1., 0.]],
          [[0., 1., 0.], [0., 0., 1.], [1., 0., 0.]],
          [[0., -1., 0.], [0., 0., -1.], [1., 0., 0.]],
          [[0., -1., 0.], [0., 0., 1.], [-1., 0., 0.]],
          [[0., 1., 0.], [0., 0., -1.], [-1., 0., 0.]],
          [[cos3, sin3, 0.], [msin3, cos3, 0.], [0., 0., 1.]],
          [[cos3, msin3, 0.], [sin3, cos3, 0.], [0., 0., 1.]],
          [[mcos3, sin3, 0.], [msin3, mcos3, 0.], [0., 0., 1.]],
          [[mcos3, msin3, 0.], [sin3, mcos3, 0.], [0., 0., 1.]],
          [[cos3, msin3, 0.], [msin3, mcos3, 0.], [0., 0., -1.]],
          [[cos3, sin3, 0.], [sin3, mcos3, 0.], [0., 0., -1.]],
          [[mcos3, msin3, 0.], [msin3, cos3, 0.], [0., 0., -1.]],
          [[mcos3, sin3, 0.], [sin3, cos3, 0.], [0., 0., -1.]],
      ]
  
      # Names for the 32 matrices, with and without inversion
      matrices_name = [
          'identity                                     ', '180 deg rotation - cart. axis [0,0,1]        ',
          '180 deg rotation - cart. axis [0,1,0]        ', '180 deg rotation - cart. axis [1,0,0]        ',
          '180 deg rotation - cart. axis [1,1,0]        ', '180 deg rotation - cart. axis [1,-1,0]       ',
          ' 90 deg rotation - cart. axis [0,0,-1]       ', ' 90 deg rotation - cart. axis [0,0,1]        ',
          '180 deg rotation - cart. axis [1,0,1]        ', '180 deg rotation - cart. axis [-1,0,1]       ',
          ' 90 deg rotation - cart. axis [0,1,0]        ', ' 90 deg rotation - cart. axis [0,-1,0]       ',
          '180 deg rotation - cart. axis [0,1,1]        ', '180 deg rotation - cart. axis [0,1,-1]       ',
          ' 90 deg rotation - cart. axis [-1,0,0]       ', ' 90 deg rotation - cart. axis [1,0,0]        ',
          '120 deg rotation - cart. axis [-1,-1,-1]     ', '120 deg rotation - cart. axis [-1,1,1]       ',
          '120 deg rotation - cart. axis [1,1,-1]       ', '120 deg rotation - cart. axis [1,-1,1]       ',
          '120 deg rotation - cart. axis [1,1,1]        ', '120 deg rotation - cart. axis [-1,1,-1]      ',
          '120 deg rotation - cart. axis [1,-1,-1]      ', '120 deg rotation - cart. axis [-1,-1,1]      ',
          ' 60 deg rotation - cryst. axis [0,0,1]       ', ' 60 deg rotation - cryst. axis [0,0,-1]      ',
          '120 deg rotation - cryst. axis [0,0,1]       ', '120 deg rotation - cryst. axis [0,0,-1]      ',
          '180 deg rotation - cryst. axis [1,-1,0]      ', '180 deg rotation - cryst. axis [2,1,0]       ',
          '180 deg rotation - cryst. axis [0,1,0]       ', '180 deg rotation - cryst. axis [1,1,0]       ',
          'inversion                                    ', 'inv. 180 deg rotation - cart. axis [0,0,1]   ',
          'inv. 180 deg rotation - cart. axis [0,1,0]   ', 'inv. 180 deg rotation - cart. axis [1,0,0]   ',
          'inv. 180 deg rotation - cart. axis [1,1,0]   ', 'inv. 180 deg rotation - cart. axis [1,-1,0]  ',
          'inv.  90 deg rotation - cart. axis [0,0,-1]  ', 'inv.  90 deg rotation - cart. axis [0,0,1]   ',
          'inv. 180 deg rotation - cart. axis [1,0,1]   ', 'inv. 180 deg rotation - cart. axis [-1,0,1]  ',
          'inv.  90 deg rotation - cart. axis [0,1,0]   ', 'inv.  90 deg rotation - cart. axis [0,-1,0]  ',
          'inv. 180 deg rotation - cart. axis [0,1,1]   ', 'inv. 180 deg rotation - cart. axis [0,1,-1]  ',
          'inv.  90 deg rotation - cart. axis [-1,0,0]  ', 'inv.  90 deg rotation - cart. axis [1,0,0]   ',
          'inv. 120 deg rotation - cart. axis [-1,-1,-1]', 'inv. 120 deg rotation - cart. axis [-1,1,1]  ',
          'inv. 120 deg rotation - cart. axis [1,1,-1]  ', 'inv. 120 deg rotation - cart. axis [1,-1,1]  ',
          'inv. 120 deg rotation - cart. axis [1,1,1]   ', 'inv. 120 deg rotation - cart. axis [-1,1,-1] ',
          'inv. 120 deg rotation - cart. axis [1,-1,-1] ', 'inv. 120 deg rotation - cart. axis [-1,-1,1] ',
          'inv.  60 deg rotation - cryst. axis [0,0,1]  ', 'inv.  60 deg rotation - cryst. axis [0,0,-1] ',
          'inv. 120 deg rotation - cryst. axis [0,0,1]  ', 'inv. 120 deg rotation - cryst. axis [0,0,-1] ',
          'inv. 180 deg rotation - cryst. axis [1,-1,0] ', 'inv. 180 deg rotation - cryst. axis [2,1,0]  ',
          'inv. 180 deg rotation - cryst. axis [0,1,0]  ', 'inv. 180 deg rotation - cryst. axis [1,1,0]  '
      ]
  
      rotations = []
  
      for key, value in zip(matrices_name[:len(transposed_matrices_cartesian)], transposed_matrices_cartesian):
          rotations.append({'name': key, 'matrix': numpy.transpose(value), 'inversion': False})
  
      for key, value in zip(matrices_name[len(transposed_matrices_cartesian):], transposed_matrices_cartesian):
          rotations.append({'name': key, 'matrix': numpy.transpose(value), 'inversion': True})
  
      return rotations
  
  
  def detect_important_message(logs, line):
  
      message_map = {
          'error': {
              'Maximum CPU time exceeded': 'ERROR_OUT_OF_WALLTIME',
              'convergence NOT achieved after': 'ERROR_ELECTRONIC_CONVERGENCE_NOT_REACHED',
              'history already reset at previous step: stopping': 'ERROR_IONIC_CYCLE_BFGS_HISTORY_FAILURE',
              'problems computing cholesky': 'ERROR_DIAGONALIZATION_CHOLESKY_DECOMPOSITION',
              'charge is wrong': 'ERROR_CHARGE_IS_WRONG',
              'not orthogonal operation': 'ERROR_SYMMETRY_NON_ORTHOGONAL_OPERATION',
              'problems computing cholesky': 'ERROR_COMPUTING_CHOLESKY',
              'dexx is negative': 'ERROR_DEXX_IS_NEGATIVE',
              'some nodes have no k-points': 'ERROR_NPOOLS_TOO_HIGH',
          },
          'warning': {
              'Warning:': None,
              'DEPRECATED:': None,
          }
      }
  
      # Match any known error and warning messages
      for marker, message in message_map['error'].items():
          if marker in line:
              if message is None:
                  message = line
              logs.error.append(message)
  
      for marker, message in message_map['warning'].items():
          if marker in line:
              if message is None:
                  message = line
              logs.warning.append(message)
  
  
  def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None):
      """Parses the stdout content of a Quantum ESPRESSO `pw.x` calculation.
  
      :param stdout: the stdout content as a string
      :param input_parameters: dictionary with the input parameters
      :param parser_options: the parser options from the settings input parameter node
      :param parsed_xml: dictionary with data parsed from the XML output file
      :returns: tuple of two dictionaries, with the parsed data and log messages, respectively
      """
      if parser_options is None:
          parser_options = {}
  
      if parsed_xml is None:
          parsed_xml = {}
  
      # Separate the input string into separate lines
      data_lines = stdout.split('\n')
  
      logs = get_logging_container()
  
      parsed_data = {}
      vdw_correction = False
      bands_data = parsed_xml.pop('bands', {})
      structure_data = parsed_xml.pop('structure', {})
      trajectory_data = {}
  
      maximum_ionic_steps = None
      marker_bfgs_converged = False
  
      # First check whether the `JOB DONE` message was written, otherwise the job was interrupted
      for line in data_lines:
          if 'JOB DONE' in line:
              break
      else:
          logs.error.append('ERROR_OUTPUT_STDOUT_INCOMPLETE')
  
      # Determine whether the input switched on an electric field
      lelfield = input_parameters.get('CONTROL', {}).get('lelfield', False)
  
      # Find some useful quantities.
      if not parsed_xml.get('number_of_bands', None):
          try:
              for line in stdout.split('\n'):
                  if 'lattice parameter (alat)' in line:
                      alat = float(line.split('=')[1].split('a.u')[0])
                  elif 'number of atoms/cell' in line:
                      nat = int(line.split('=')[1])
                  elif 'number of atomic types' in line:
                      ntyp = int(line.split('=')[1])
                  elif 'unit-cell volume' in line:
                      if '(a.u.)^3' in line:
                          volume = float(line.split('=')[1].split('(a.u.)^3')[0])
                      else:
                          # occurs in v5.3.0
                          volume = float(line.split('=')[1].split('a.u.^3')[0])
                  elif 'number of Kohn-Sham states' in line:
                      nbnd = int(line.split('=')[1])
                  elif 'number of k points' in line:
                      nk = int(line.split('=')[1].split()[0])
                      if input_parameters.get('SYSTEM', {}).get('nspin', 1) > 1:
                          # QE counts twice each k-point in spin-polarized calculations
                          nk /= 2
                  elif 'Dense  grid' in line:
                      FFT_grid = [int(g) for g in line.split('(')[1].split(')')[0].split(',')]
                  elif 'Smooth grid' in line:
                      smooth_FFT_grid = [int(g) for g in line.split('(')[1].split(')')[0].split(',')]
                      break
              alat *= CONSTANTS.bohr_to_ang
              volume *= CONSTANTS.bohr_to_ang**3
              parsed_data['lattice_parameter_initial'] = alat
              parsed_data['number_of_bands'] = nbnd
              try:
                  parsed_data['number_of_k_points'] = nk
                  parsed_data['fft_grid'] = FFT_grid
                  parsed_data['smooth_fft_grid'] = smooth_FFT_grid
              except NameError:  # these are not crucial, so parsing does not fail if they are not found
                  pass
          except NameError:  # nat or other variables where not found, and thus not initialized
  
              # Try to get some error messages
              lines = stdout.split('\n')
  
              for line_number, line in enumerate(lines):
                  # Compare the line to the known set of error and warning messages and add them to the log container
                  detect_important_message(logs, line)
  
              if len(logs.error) or len(logs.warning) > 0:
                  parsed_data['trajectory'] = trajectory_data
                  return parsed_data, logs
  
              # did not find any error message -> raise an Error and do not return anything
              raise QEOutputParsingError('Parser cannot load basic info.')
      else:
          nat = structure_data['number_of_atoms']
          ntyp = structure_data['number_of_species']
          alat = structure_data['lattice_parameter_xml']
          volume = structure_data['cell']['volume']
      # NOTE: lattice_parameter_xml is the lattice parameter of the xml file
      # in the units used by the code. lattice_parameter instead in angstroms.
  
      # Save these two quantities in the parsed_data, because they will be
      # useful for queries (maybe), and structure_data will not be stored as a Dict
      parsed_data['number_of_atoms'] = nat
      parsed_data['number_of_species'] = ntyp
      parsed_data['volume'] = volume
  
      c_bands_error = False
  
      # now grep quantities that can be considered isolated informations.
      for count, line in enumerate(data_lines):
  
          # Compare the line to the known set of error and warning messages and add them to the log container
          detect_important_message(logs, line)
  
          # to be used for later
          if 'Non-local correlation energy' in line:
              vdw_correction = True
  
          elif 'Cartesian axes' in line:
              # this is the part when initial positions and chemical
              # symbols are printed (they do not change during a run)
              i = count + 1
              while i < count + 10 and not ('site n.' in data_lines[i] and 'atom' in data_lines[i]):
                  i += 1
              if 'site n.' in data_lines[i] and 'atom' in data_lines[i]:
                  trajectory_data['atomic_species_name'] = [data_lines[i + 1 + j].split()[1] for j in range(nat)]
  
          # parse the initialization time (take only first occurence)
          elif ('init_wall_time_seconds' not in parsed_data and 'total cpu time spent up to now is' in line):
              init_time = float(line.split('total cpu time spent up to now is')[1].split('secs')[0])
              parsed_data['init_wall_time_seconds'] = init_time
  
          # parse dynamical RAM estimates
          elif 'Estimated max dynamical RAM per process' in line:
              value = line.split('>')[-1]
              match = re.match(r'\s+([+-]?\d+(\.\d*)?|\.\d+([eE][+-]?\d+)?)\s*(Mb|MB|GB)', value)
              if match:
                  try:
                      parsed_data['estimated_ram_per_process'] = float(match.group(1))
                      parsed_data[f'estimated_ram_per_process{units_suffix}'] = match.group(4)
                  except (IndexError, ValueError):
                      pass
  
          # parse dynamical RAM estimates
          elif 'Estimated total dynamical RAM' in line:
              value = line.split('>')[-1]
              match = re.match(r'\s+([+-]?\d+(\.\d*)?|\.\d+([eE][+-]?\d+)?)\s*(Mb|MB|GB)', value)
              if match:
                  try:
                      parsed_data['estimated_ram_total'] = float(match.group(1))
                      parsed_data[f'estimated_ram_total{units_suffix}'] = match.group(4)
                  except (IndexError, ValueError):
                      pass
  
          # parse the global file, for informations that are written only once
          elif 'PWSCF' in line and 'WALL' in line:
              try:
                  time = line.split('CPU')[1].split('WALL')[0]
                  parsed_data['wall_time'] = time
              except Exception:
                  logs.warning.append('Error while parsing wall time.')
              try:
                  parsed_data['wall_time_seconds'] = convert_qe_time_to_sec(time)
              except ValueError:
                  raise QEOutputParsingError('Unable to convert wall_time in seconds.')
  
          # for later control on relaxation-dynamics convergence
          elif 'nstep' in line and '=' in line:
              maximum_ionic_steps = int(line.split()[2])
  
          elif 'bfgs converged in' in line:
              marker_bfgs_converged = True
  
          elif 'number of bfgs steps' in line:
  
              try:
                  parsed_data['number_ionic_steps'] += 1
              except KeyError:
                  parsed_data['number_ionic_steps'] = 1
  
          elif 'A final scf calculation at the relaxed structure' in line:
              parsed_data['final_scf'] = True
  
          elif 'point group' in line:
              if 'k-point group' not in line:
                  try:
                      # Split line in components delimited by either space(s) or
                      # parenthesis and filter out empty strings
                      line_elems = [_f for _f in re.split(r' +|\(|\)', line) if _f]
  
                      pg_international = line_elems[-1]
                      pg_schoenflies = line_elems[-2]
  
                      parsed_data['pointgroup_international'] = pg_international
                      parsed_data['pointgroup_schoenflies'] = pg_schoenflies
  
                  except Exception:
                      warning = f'Problem parsing point group, I found: {line.strip()}'
                      logs.warning.append(warning)
  
          # special parsing of c_bands error
          elif 'c_bands' in line and 'eigenvalues not converged' in line:
              c_bands_error = True
  
          elif 'iteration #' in line:
              if 'Calculation restarted' not in line and 'Calculation stopped' not in line:
                  try:
                      parsed_data['total_number_of_scf_iterations'] += 1
                  except KeyError:
                      parsed_data['total_number_of_scf_iterations'] = 1
  
              if c_bands_error:
                  # if there is another iteration, c_bands is not necessarily a problem
                  # I put a warning only if c_bands error appears in the last iteration
                  c_bands_error = False
  
      if c_bands_error:
          logs.warning.append('c_bands: at least 1 eigenvalues not converged')
  
      # I split the output text in the atomic SCF calculations.
      # the initial part should be things already contained in the xml.
      # (cell, initial positions, kpoints, ...) and I skip them.
      # In case, parse for them before this point.
      # Put everything in a trajectory_data dictionary
      relax_steps = stdout.split('Self-consistent Calculation')[1:]
      relax_steps = [i.split('\n') for i in relax_steps]
  
      # now I create a bunch of arrays for every step.
  
      for data_step in relax_steps:
          trajectory_frame = {}
  
          for count, line in enumerate(data_step):
  
              if 'CELL_PARAMETERS' in line:
                  try:
                      a1 = [float(s) for s in data_step[count + 1].split()]
                      a2 = [float(s) for s in data_step[count + 2].split()]
                      a3 = [float(s) for s in data_step[count + 3].split()]
                      # try except indexerror for not enough lines
                      lattice = line.split('(')[1].split(')')[0].split('=')
                      if lattice[0].lower() not in ['alat', 'bohr', 'angstrom']:
                          raise QEOutputParsingError(
                              'Error while parsing cell_parameters: ' + f'unsupported units {lattice[0]}'
                          )
  
                      if 'alat' in lattice[0].lower():
                          a1 = [alat * CONSTANTS.bohr_to_ang * float(s) for s in a1]
                          a2 = [alat * CONSTANTS.bohr_to_ang * float(s) for s in a2]
                          a3 = [alat * CONSTANTS.bohr_to_ang * float(s) for s in a3]
                          lattice_parameter_b = float(lattice[1])
                          if abs(lattice_parameter_b - alat) > lattice_tolerance:
                              raise QEOutputParsingError(
                                  'Lattice parameters mismatch! ' + f'{lattice_parameter_b} vs {alat}'
                              )
                      elif 'bohr' in lattice[0].lower():
                          lattice_parameter_b *= CONSTANTS.bohr_to_ang
                          a1 = [CONSTANTS.bohr_to_ang * float(s) for s in a1]
                          a2 = [CONSTANTS.bohr_to_ang * float(s) for s in a2]
                          a3 = [CONSTANTS.bohr_to_ang * float(s) for s in a3]
                      trajectory_data.setdefault('lattice_vectors_relax', []).append([a1, a2, a3])
  
                  except Exception:
                      logs.warning.append('Error while parsing relaxation cell parameters.')
  
              elif 'ATOMIC_POSITIONS' in line:
                  try:
                      this_key = 'atomic_positions_relax'
                      # the inizialization of tau prevent parsed_data to be associated
                      # to the pointer of the previous iteration
                      metric = line.split('(')[1].split(')')[0]
                      if metric == 'crystal':
                          this_key = 'atomic_fractionals_relax'
                      elif metric not in ['alat', 'bohr', 'angstrom']:
                          raise QEOutputParsingError('Error while parsing atomic_positions: units not supported.')
                      # TODO: check how to map the atoms in the original scheme
                      positions = []
                      for i in range(nat):
                          line2 = data_step[count + 1 + i].split()
                          tau = [float(s) for s in line2[1:4]]
                          if metric == 'alat':
                              tau = [alat * float(s) for s in tau]
                          elif metric == 'bohr':
                              tau = [CONSTANTS.bohr_to_ang * float(s) for s in tau]
                          positions.append(tau)
                      trajectory_data.setdefault(this_key, []).append(positions)
                  except Exception:
                      logs.warning.append('Error while parsing relaxation atomic positions.')
  
              # NOTE: in the above, the chemical symbols are not those of AiiDA
              # since the AiiDA structure is different. So, I assume now that the
              # order of atoms is the same of the input atomic structure.
  
              # Computed dipole correction in slab geometries.
              # save dipole in debye units, only at last iteration of scf cycle
              elif 'Computed dipole along edir' in line:
                  j = count + 3
                  line2 = data_step[j]
                  try:
                      units = line2.split()[-1]
                      if default_dipole_units.lower() not in units.lower():  # only debye
                          raise QEOutputParsingError(
                              f'Error parsing the dipole correction. Units {units} are not supported.'
                          )
                      value = float(line2.split()[-2])
                  except IndexError:  # on units
                      pass
                  # save only the last dipole correction
                  while 'Computed dipole along edir' not in line2:
                      j += 1
                      try:
                          line2 = data_step[j]
                      except IndexError:  # The dipole is also written at the beginning of a new bfgs iteration
                          break
                      if 'End of self-consistent calculation' in line2:
                          trajectory_data.setdefault('dipole', []).append(value)
                          parsed_data['dipole' + units_suffix] = default_dipole_units
                          break
  
              # saving the SCF convergence accuracy for each SCF cycle
              # If for some step this line is not printed, the later check with the scf_accuracy array length should catch it
              elif 'estimated scf accuracy' in line:
                  try:
                      value = float(line.split()[-2]) * CONSTANTS.ry_to_ev
                      trajectory_data.setdefault('scf_accuracy', []).append(value)
                  except Exception:
                      logs.warning.append('Error while parsing scf accuracy.')
  
              elif 'convergence has been achieved in' in line or 'convergence NOT achieved after' in line:
                  try:
                      value = int(line.split('iterations')[0].split()[-1])
                      trajectory_data.setdefault('scf_iterations', []).append(value)
                  except Exception:
                      logs.warning.append('Error while parsing scf iterations.')
  
              elif 'Calculation stopped in scf loop at iteration' in line:
                  try:
                      value = int(line.split()[-1])
                      trajectory_data.setdefault('scf_iterations', []).append(value)
                  except Exception:
                      logs.warning.append('Error while parsing scf iterations.')
  
              elif 'End of self-consistent calculation' in line:
                  # parse energy threshold for diagonalization algorithm
                  try:
                      j = 0
                      while True:
                          j -= 1
                          line2 = data_step[count + j]
                          if 'ethr' in line2:
                              value = float(line2.split('=')[1].split(',')[0])
                              break
                      trajectory_data.setdefault('energy_threshold', []).append(value)
                  except Exception:
                      logs.warning.append('Error while parsing ethr.')
  
                  # parse final magnetic moments, if present
                  try:
                      j = 0
                      while True:
                          j -= 1
                          line2 = data_step[count + j]
                          if 'Magnetic moment per site' in line2:
                              break
                          if 'iteration' in line2:
                              raise QEOutputParsingError
                      mag_moments = []
                      charges = []
                      while True:
                          j += 1
                          line2 = data_step[count + j]
                          if 'atom:' in line2:
                              mag_moments.append(float(line2.split('magn:')[1].split()[0]))
                              charges.append(float(line2.split('charge:')[1].split()[0]))
                          # Alternate parsing for new format introduced in v6.8
                          elif 'atom' in line2:
                              mag_moments.append(float(line2.split('magn=')[1].split()[0]))
                              charges.append(float(line2.split('charge=')[1].split()[0]))
                          if len(mag_moments) == nat:
                              break
                      trajectory_data.setdefault('atomic_magnetic_moments', []).append(mag_moments)
                      trajectory_data.setdefault('atomic_charges', []).append(charges)
                      parsed_data['atomic_magnetic_moments' + units_suffix] = default_magnetization_units
                      parsed_data['atomic_charges' + units_suffix] = default_charge_units
                  except QEOutputParsingError:
                      pass
  
              # grep energy and possibly, magnetization
              elif '!' in line:
                  try:
  
                      En = float(line.split('=')[1].split('Ry')[0]) * CONSTANTS.ry_to_ev
  
                      # Up till v6.5, the line after total energy would be the Harris-Foulkes estimate, followed by the
                      # estimated SCF accuracy. However, pw.x v6.6 removed the HF estimate line.
                      marker = 'estimated scf accuracy'
                      for i in range(5):
                          subline = data_step[count + i]
                          if marker in subline:
                              try:
                                  E_acc = float(subline.split('<')[1].split('Ry')[0]) * CONSTANTS.ry_to_ev
                              except Exception:
                                  pass
                              else:
                                  break
                      else:
                          raise KeyError(f'could not find and parse the line with `{marker}`')
  
                      for key, value in [['energy', En], ['energy_accuracy', E_acc]]:
                          trajectory_data.setdefault(key, []).append(value)
                          parsed_data[key + units_suffix] = default_energy_units
                      # TODO: decide units for magnetization. now bohr mag/cell
                      j = 0
                      while True:
                          j += 1
                          line2 = data_step[count + j]
  
                          for string, key in [
                              ['one-electron contribution', 'energy_one_electron'],
                              ['hartree contribution', 'energy_hartree'],
                              ['xc contribution', 'energy_xc'],
                              ['ewald contribution', 'energy_ewald'],
                              ['smearing contrib.', 'energy_smearing'],
                              ['one-center paw contrib.', 'energy_one_center_paw'],
                              ['est. exchange err', 'energy_est_exchange'],
                              ['Fock energy', 'energy_fock'],
                              ['Hubbard energy', 'energy_hubbard'],
                              # Add also ENVIRON specific contribution to the total energy
                              ['solvation energy', 'energy_solvation'],
                              ['cavitation energy', 'energy_cavitation'],
                              ['PV energy', 'energy_pv'],
                              ['periodic energy correct.', 'energy_pbc_correction'],
                              ['ionic charge energy', 'energy_ionic_charge'],
                              ['external charges energy', 'energy_external_charges']
                          ]:
                              if string in line2:
                                  value = grep_energy_from_line(line2)
                                  trajectory_data.setdefault(key, []).append(value)
                                  parsed_data[key + units_suffix] = default_energy_units
                          # magnetizations
                          if 'total magnetization' in line2:
                              this_m = line2.split('=')[1].split('Bohr')[0]
                              try:  # magnetization might be a scalar
                                  value = float(this_m)
                              except ValueError:  # but can also be a three vector component in non-collinear calcs
                                  value = [float(i) for i in this_m.split()]
                              trajectory_data.setdefault('total_magnetization', []).append(value)
                              parsed_data['total_magnetization' + units_suffix] = default_magnetization_units
  
                          elif 'absolute magnetization' in line2:
                              value = float(line2.split('=')[1].split('Bohr')[0])
                              trajectory_data.setdefault('absolute_magnetization', []).append(value)
                              parsed_data['absolute_magnetization' + units_suffix] = default_magnetization_units
                          # exit loop
                          elif 'convergence' in line2:
                              break
  
                      if vdw_correction:
                          j = 0
                          while True:
                              j += -1
                              line2 = data_step[count + j]
                              if 'Non-local correlation energy' in line2:
                                  value = grep_energy_from_line(line2)
                                  trajectory_data.setdefault('energy_vdw', []).append(value)
                                  break
                          parsed_data['energy_vdw' + units_suffix] = default_energy_units
                  except Exception as exception:
                      import traceback
                      traceback.print_exc()
                      print(exception)
                      logs.warning.append('Error while parsing for energy terms.')
  
              elif 'the Fermi energy is' in line:
                  try:
                      value = float(line.split('is')[1].split('ev')[0])
                      trajectory_data.setdefault('fermi_energy', []).append(value)
                      parsed_data['fermi_energy' + units_suffix] = default_energy_units
                  except Exception:
                      logs.warning.append('Error while parsing Fermi energy from the output file.')
  
              elif 'Forces acting on atoms' in line:
                  try:
                      forces = []
                      j = 0
                      while True:
                          j += 1
                          line2 = data_step[count + j]
                          if 'atom ' in line2:
                              line2 = line2.split('=')[1].split()
                              # CONVERT FORCES IN eV/Ang
                              vec = [float(s) * CONSTANTS.ry_to_ev / CONSTANTS.bohr_to_ang for s in line2]
                              forces.append(vec)
                          if len(forces) == nat:
                              break
                      trajectory_data.setdefault('forces', []).append(forces)
                      parsed_data['forces' + units_suffix] = default_force_units
                  except Exception:
                      logs.warning.append('Error while parsing forces.')
  
              # TODO: adding the parsing support for the decomposition of the forces
  
              elif 'Total force =' in line:
                  try:  # note that I can't check the units: not written in output!
                      value = float(line.split('=')[1].split('Total')[0]) * CONSTANTS.ry_to_ev / CONSTANTS.bohr_to_ang
                      trajectory_data.setdefault('total_force', []).append(value)
                      parsed_data['total_force' + units_suffix] = default_force_units
                  except Exception:
                      logs.warning.append('Error while parsing total force.')
  
              elif ('entering subroutine stress ...'
                    in line) or ('Computing stress (Cartesian axis) and pressure' in line):
                  try:
                      stress = []
                      count2 = None
                      for k in range(15):  # Up to 15 lines later - more than 10 are needed if vdW is turned on
                          if 'P=' in data_step[count + k + 1]:
                              count2 = count + k + 1
                      if count2 is None:
                          logs.warning.append(
                              'Error while parsing stress tensor: '
                              '"P=" not found within 15 lines from the start of the stress block'
                          )
                      else:
                          if '(Ry/bohr**3)' not in data_step[count2]:
                              raise QEOutputParsingError('Error while parsing stress: unexpected units.')
                          for k in range(3):
                              line2 = data_step[count2 + k + 1].split()
                              vec = [float(s) * 10**(-9) * CONSTANTS.ry_si / (CONSTANTS.bohr_si)**3 for s in line2[0:3]]
                              stress.append(vec)
                          trajectory_data.setdefault('stress', []).append(stress)
                          parsed_data['stress' + units_suffix] = default_stress_units
                  except Exception:
                      import traceback
                      logs.warning.append(f'Error while parsing stress tensor: {traceback.format_exc()}')
  
              # Electronic and ionic dipoles when 'lelfield' was set to True in input parameters
              elif lelfield is True:
  
                  if 'Electronic Dipole per cell' in line:
                      electronic_dipole = float(line.split()[-1])
                      trajectory_frame.setdefault('electronic_dipole_cell_average', []).append(electronic_dipole)
  
                  elif 'Ionic Dipole per cell' in line:
                      ionic_dipole = float(line.split()[-1])
                      trajectory_frame.setdefault('ionic_dipole_cell_average', []).append(ionic_dipole)
  
                  elif 'Electronic Dipole on Cartesian axes' in line:
                      electronic_dipole = [float(data_step[count + i + 1].split()[1]) for i in range(3)]
                      trajectory_frame.setdefault('electronic_dipole_cartesian_axes', []).append(electronic_dipole)
  
                  elif 'Ionic Dipole on Cartesian axes' in line:
                      ionic_dipole = [float(data_step[count + i + 1].split()[1]) for i in range(3)]
                      trajectory_frame.setdefault('ionic_dipole_cartesian_axes', []).append(ionic_dipole)
  
          # End of trajectory frame, only keep last entries for dipole related values
          if lelfield is True:
  
              # For every property only get the last entry if possible
              try:
                  ed_cell = trajectory_frame['electronic_dipole_cell_average'].pop()
              except IndexError:
                  ed_cell = None
  
              try:
                  ed_axes = trajectory_frame['electronic_dipole_cartesian_axes'].pop()
              except IndexError:
                  ed_axes = None
  
              try:
                  id_cell = trajectory_frame['ionic_dipole_cell_average'].pop()
              except IndexError:
                  id_cell = None
  
              try:
                  id_axes = trajectory_frame['ionic_dipole_cartesian_axes'].pop()
              except IndexError:
                  id_axes = None
  
              # Only add them if all four properties were successfully parsed
              if all([value is not None for value in [ed_cell, ed_axes, id_cell, id_axes]]):
                  trajectory_data.setdefault('electronic_dipole_cell_average', []).append(ed_cell)
                  trajectory_data.setdefault('electronic_dipole_cartesian_axes', []).append(ed_axes)
                  trajectory_data.setdefault('ionic_dipole_cell_average', []).append(id_cell)
                  trajectory_data.setdefault('ionic_dipole_cartesian_axes', []).append(id_axes)
  
      # check consistency of scf_accuracy and scf_iterations
      if 'scf_accuracy' in trajectory_data:
          if 'scf_iterations' in trajectory_data:
              if len(trajectory_data['scf_accuracy']) != sum(trajectory_data['scf_iterations']):
                  logs.warning.append(
                      'the length of scf_accuracy does not match the sum of the elements of scf_iterations.'
                  )
          else:
              logs.warning.append('"the scf_accuracy array was parsed but the scf_iterations was not.')
  
      # If specified in the parser options, parse the atomic occupations
      parse_atomic_occupations = parser_options.get('parse_atomic_occupations', False)
  
      if parse_atomic_occupations:
  
          atomic_occupations = {}
          hubbard_blocks = stdout.split('LDA+U parameters')
  
          for line in hubbard_blocks[-1].split('\n'):
  
              if 'Tr[ns(na)]' in line:
  
                  values = line.split('=')
                  atomic_index = values[0].split()[1]
                  occupations = values[1].split()
  
                  if len(occupations) == 1:
                      atomic_occupations[atomic_index] = {'total': occupations[0]}
                  elif len(occupations) == 3:
                      atomic_occupations[atomic_index] = {
                          'up': occupations[0],
                          'down': occupations[1],
                          'total': occupations[2]
                      }
                  else:
                      continue
  
          parsed_data['atomic_occupations'] = atomic_occupations
  
      # Ionic calculations and BFGS algorithm did not print that calculation is converged
      if 'atomic_positions_relax' in trajectory_data and not marker_bfgs_converged:
          logs.error.append('ERROR_IONIC_CONVERGENCE_NOT_REACHED')
  
      # Ionic calculation that hit the maximum number of ionic steps. Note: does not necessarily mean that convergence was
      # not reached as it could have occurred in the last step.
      if maximum_ionic_steps is not None and maximum_ionic_steps == parsed_data.get('number_ionic_steps', None):
          logs.warning.append('ERROR_MAXIMUM_IONIC_STEPS_REACHED')
  
      # Remove duplicate log messages by turning it into a set. Then convert back to list as that is what is expected
      logs.error = list(set(logs.error))
      logs.warning = list(set(logs.warning))
  
      parsed_data['bands'] = bands_data
      parsed_data['structure'] = structure_data
      parsed_data['trajectory'] = trajectory_data
  
      return parsed_data, logs
  
  
  def grep_energy_from_line(line):
      try:
          return float(line.split('=')[1].split('Ry')[0]) * CONSTANTS.ry_to_ev
      except Exception:
          raise QEOutputParsingError('Error while parsing energy')
  