👌 Add fixes for SIRIUS FORTRAN XML issue

See

electronic-structure/q-e-sirius#45

Also:

* Add some basic parsing for charge and magnetic moments.
* Add `DIRECT_MINIMIZATION` to the possible namelists.

src/aiida_quantumespresso/calculations/__init__.py

@@ -727,7 +727,7 @@ def _generate_PWCPinputdata(cls, parameters, settings, pseudos, structure, kpoin
  ) from exception
 
  try:
- namelists_toprint = cls._automatic_namelists[calculation_type]
+ namelists_toprint = list(cls._automatic_namelists[calculation_type])
  except KeyError as exception:
  raise exceptions.InputValidationError(
  'Unknown `calculation` value in CONTROL namelist {calculation_type}. Otherwise, specify the list of'

src/aiida_quantumespresso/calculations/pw.py

@@ -13,14 +13,16 @@
 class PwCalculation(BasePwCpInputGenerator):
  """`CalcJob` implementation for the pw.x code of Quantum ESPRESSO."""
 
+ _default_namelists = ('CONTROL', 'SYSTEM', 'ELECTRONS', 'DIRECT_MINIMIZATION')
+
  _automatic_namelists = {
- 'scf': ['CONTROL', 'SYSTEM', 'ELECTRONS'],
- 'nscf': ['CONTROL', 'SYSTEM', 'ELECTRONS'],
- 'bands': ['CONTROL', 'SYSTEM', 'ELECTRONS'],
- 'relax': ['CONTROL', 'SYSTEM', 'ELECTRONS', 'IONS'],
- 'md': ['CONTROL', 'SYSTEM', 'ELECTRONS', 'IONS'],
- 'vc-md': ['CONTROL', 'SYSTEM', 'ELECTRONS', 'IONS', 'CELL'],
- 'vc-relax': ['CONTROL', 'SYSTEM', 'ELECTRONS', 'IONS', 'CELL'],
+ 'scf': _default_namelists,
+ 'nscf': _default_namelists,
+ 'bands': _default_namelists,
+ 'relax': _default_namelists + ('IONS',),
+ 'md': _default_namelists + ('IONS',),
+ 'vc-md': _default_namelists + ('IONS', 'CELL'),
+ 'vc-relax': _default_namelists + ('IONS', 'CELL'),
  }
 
  # Keywords that cannot be set by the user but will be set by the plugin

src/aiida_quantumespresso/parsers/parse_raw/pw.py

@@ -942,6 +942,42 @@ def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None,
  return parsed_data, logs
 
 
+def parse_sirius_stdout(stdout):
+ """Parse the `stdout` of a SIRIUS-enabled Quantum ESPRESSO calculation."""
+
+ def block_search(string: str, block_start: str, block_end: str) -> list:
+
+ start_match = re.search(block_start, string)
+
+ if start_match is None:
+ return []
+
+ start_index = start_match.end()
+ end_index = re.search(block_end, string[start_index:]).start()
+ block_lines = scf[start_index:start_index + end_index].split('\n')
+ return block_lines
+
+ scf_runs = stdout.split('* running SCF ground state *')[1:]
+
+ trajectory_data = {}
+
+ for scf in scf_runs:
+
+ magnetic_moments = []
+
+ for line in block_search(scf, r'\[find\] Charges and magnetic moments', r'\[find\] total charge'):
+ match = re.match(r'\[find\]\s+\d+\s+\[[\s\d\.\,-]+\]\s+(?P<magmom>[\d\.\,-]+)', line)
+ if match is not None:
+ magnetic_moments.append(float(match.group('magmom')))
+
+ if magnetic_moments:
+ trajectory_data.setdefault('atomic_magnetic_moments', []).append(magnetic_moments)
+
+ return {
+ 'trajectory': trajectory_data,
+ }
+
+
 def grep_energy_from_line(line):
  try:
  return float(line.split('=')[1].split('Ry')[0]) * CONSTANTS.ry_to_ev

src/aiida_quantumespresso/parsers/parse_xml/parse.py

@@ -76,11 +76,7 @@ def parse_xml_post_6_2(xml):
  # xml_dictionary['key']['@attr'] returns its attribute 'attr'
  # xml_dictionary['key']['nested_key'] goes one level deeper.
 
- xml_dictionary, errors = xsd.to_dict(xml, validation='lax')
- if errors:
- logs.error.append(f'{len(errors)} XML schema validation error(s) schema: {schema_filepath}:')
- for err in errors:
- logs.error.append(str(err))
+ xml_dictionary = xsd.to_dict(xml, validation='skip')
 
  xml_version = Version(xml_dictionary['general_info']['xml_format']['@VERSION'])
  inputs = xml_dictionary.get('input', {})
@@ -572,6 +568,8 @@ def parse_xml_post_6_2(xml):
 
 def parse_step_to_trajectory(trajectory, data, skip_structure=False):
  """."""
+ from ..pw import fix_sirius_xml_prints
+
  if 'scf_conv' in data and 'n_scf_steps' in data['scf_conv']:
  scf_iterations = data['scf_conv']['n_scf_steps'] # Can be zero in case of initialization-only calculation
  if scf_iterations:
@@ -586,12 +584,12 @@ def parse_step_to_trajectory(trajectory, data, skip_structure=False):
  atomic_structure = data['atomic_structure']
 
  if 'atomic_positions' in atomic_structure:
- positions = np.array([a['$'] for a in atomic_structure['atomic_positions']['atom']])
+ positions = fix_sirius_xml_prints(np.array([a['$'] for a in atomic_structure['atomic_positions']['atom']]))
  trajectory['positions'].append(positions * CONSTANTS.bohr_to_ang)
 
  if 'cell' in atomic_structure:
  cell = atomic_structure['cell']
- cell = np.array([cell['a1'], cell['a2'], cell['a3']])
+ cell = fix_sirius_xml_prints(np.array([cell['a1'], cell['a2'], cell['a3']]))
  trajectory['cells'].append(cell * CONSTANTS.bohr_to_ang)
 
  if 'total_energy' in data:

src/aiida_quantumespresso/parsers/pw.py

@@ -9,11 +9,34 @@
 
 from aiida_quantumespresso.calculations.pw import PwCalculation
 from aiida_quantumespresso.utils.mapping import get_logging_container
+from aiida_quantumespresso.workflows.protocols.utils import recursive_merge
 
 from .base import BaseParser
 from .parse_raw.pw import reduce_symmetries
 
 
+def fix_sirius_xml_prints(array):
+ """Fix an issue where SIRIUS prints very small numbers incorrectly.
+ 
+ In some cases SIRIUS prints very small numbers in scientific notation, but
+ without the capital `E` to indicate the exponent:
+
+ 
+ <forces rank="2" dims=" 3 6">
+ 3.220681500679849E-86 4.347898683069749-103 -3.696810758148386E-49
+
+ This function will fix this by converting any number that cannot be converted into
+ a float to zero.
+ """
+ def try_convert(s):
+ try:
+ return float(s)
+ except (ValueError, TypeError):
+ return 0
+
+ return numpy.vectorize(try_convert)(array)
+
+
 class PwParser(BaseParser):
  """`Parser` implementation for the `PwCalculation` calculation job class."""
 
@@ -40,11 +63,11 @@ def parse(self, **kwargs):
  parser_options = settings.get(self.get_parser_settings_key(), None)
 
  # Verify that the retrieved_temporary_folder is within the arguments if temporary files were specified
- if self.node.base.attributes.get('retrieve_temporary_list', None):
- try:
- dir_with_bands = kwargs['retrieved_temporary_folder']
- except KeyError:
- return self.exit(self.exit_codes.ERROR_NO_RETRIEVED_TEMPORARY_FOLDER)
+ # if self.node.base.attributes.get('retrieve_temporary_list', None):
+ #  try:
+ #  dir_with_bands = kwargs['retrieved_temporary_folder']
+ #  except KeyError:
+ #  return self.exit(self.exit_codes.ERROR_NO_RETRIEVED_TEMPORARY_FOLDER)
 
  # We check if the `CRASH` file was retrieved. If so, we parse its output
  crash_file_filename = self.node.process_class._CRASH_FILE
@@ -404,6 +427,7 @@ def parse_stdout(self, parameters, parser_options=None, crash_file=None):
  :return: tuple of two dictionaries, first with raw parsed data and second with log messages
  """
  from aiida_quantumespresso.parsers.parse_raw.pw import parse_stdout
+ from aiida_quantumespresso.parsers.parse_raw.pw import parse_sirius_stdout
 
  logs = get_logging_container()
  parsed_data = {}
@@ -426,6 +450,8 @@ def parse_stdout(self, parameters, parser_options=None, crash_file=None):
  logs.critical.append(traceback.format_exc())
  self.exit_code_stdout = self.exit_codes.ERROR_UNEXPECTED_PARSER_EXCEPTION.format(exception=exc)
 
+ parsed_data = recursive_merge(parsed_data, parse_sirius_stdout(stdout))
+
  # If the stdout was incomplete, most likely the job was interrupted before it could cleanly finish, so the
  # output files are most likely corrupt and cannot be restarted from
  if 'ERROR_OUTPUT_STDOUT_INCOMPLETE' in logs['error']:
@@ -519,7 +545,13 @@ def build_output_trajectory(parsed_trajectory, structure):
  )
 
  for key, value in parsed_trajectory.items():
- trajectory.set_array(key, numpy.array(value))
+ if key in (
+ 'forces',
+ 'stress'
+ ):
+ trajectory.set_array(key, fix_sirius_xml_prints(numpy.array(value)))
+ else:
+ trajectory.set_array(key, numpy.array(value))
 
  return trajectory
 
@@ -571,14 +603,16 @@ def build_output_bands(self, parsed_bands, parsed_kpoints=None):
 
  # Correct the occupation for nspin=1 calculations where Quantum ESPRESSO populates each band only halfway
  if len(parsed_bands['occupations']) > 1:
- occupations = parsed_bands['occupations']
+ occupations = numpy.array(parsed_bands['occupations'])
  else:
- occupations = 2. * numpy.array(parsed_bands['occupations'][0])
+ occupations = numpy.array(parsed_bands['occupations'][0])
 
- if len(parsed_bands['bands']) > 1:
- bands_energies = parsed_bands['bands']
- else:
- bands_energies = parsed_bands['bands'][0]
+ occupations = fix_sirius_xml_prints(occupations)
+
+ if len(parsed_bands['occupations']) > 1:
+ occupations *= 2.
+
+ bands_energies = parsed_bands['bands'][0] if len(parsed_bands['bands']) == 1 else parsed_bands['bands']
 
  bands = orm.BandsData()
  bands.set_kpointsdata(parsed_kpoints)