Allow explicitly specifying 'ibrav'.

If an 'ibrav' value other than zero is specified in the
SYSTEM inputs, the cell of the input structure is converted
into the appropriate A, B, C, cosAB, cosAC, cosBC values. As a
check, the cell is reconstructed using qe_tools. The cells are
compared element-wise, with an absolute tolerance controlled by
the 'IBRAV_CELL_TOLERANCE' setting.

aiida_quantumespresso/calculations/__init__.py

@@ -4,16 +4,20 @@
 
 import abc
 import os
+
 import six
 from six.moves import zip
 
 from aiida import orm
 from aiida.common import datastructures, exceptions
 from aiida.common.lang import classproperty
 
+from qe_tools.parsers.qeinputparser import get_cell_from_parameters
+
 from aiida_quantumespresso.utils.convert import convert_input_to_namelist_entry
 
 from .base import CalcJob
+from .helpers import QEInputValidationError
 
 
 class BasePwCpInputGenerator(CalcJob):
@@ -288,10 +292,16 @@ def _generate_PWCPinputdata(cls, parameters, settings, pseudos, structure, kpoin
 
  # ============ I prepare the input site data =============
  # ------------ CELL_PARAMETERS -----------
- cell_parameters_card = 'CELL_PARAMETERS angstrom\n'
- for vector in structure.cell:
- cell_parameters_card += ('{0:18.10f} {1:18.10f} {2:18.10f}'
- '\n'.format(*vector))
+
+ # Specify cell parameters only if 'ibrav' is either zero or
+ # unspecified.
+ if input_params.get('SYSTEM', {}).get('ibrav', 0) == 0:
+ cell_parameters_card = 'CELL_PARAMETERS angstrom\n'
+ for vector in structure.cell:
+ cell_parameters_card += ('{0:18.10f} {1:18.10f} {2:18.10f}'
+ '\n'.format(*vector))
+ else:
+ cell_parameters_card = ''
 
  # ------------- ATOMIC_SPECIES ------------
  atomic_species_card_list = []
@@ -451,7 +461,14 @@ def _generate_PWCPinputdata(cls, parameters, settings, pseudos, structure, kpoin
  # Set some variables (look out at the case! NAMELISTS should be
  # uppercase, internal flag names must be lowercase)
  input_params.setdefault('SYSTEM', {})
- input_params['SYSTEM']['ibrav'] = 0
+ input_params['SYSTEM'].setdefault('ibrav', 0)
+ ibrav = input_params['SYSTEM']['ibrav']
+ if ibrav != 0:
+ input_params['SYSTEM'].update(_get_parameters_from_cell(
+ ibrav=ibrav,
+ structure=structure,
+ tolerance=settings.pop('IBRAV_CELL_TOLERANCE', 1e-6)
+ ))
  input_params['SYSTEM']['nat'] = len(structure.sites)
  input_params['SYSTEM']['ntyp'] = len(structure.kinds)
 
@@ -641,3 +658,86 @@ def _pop_parser_options(calc_job_instance, settings_dict, ignore_errors=True):
  pass
  else:
  raise exc
+
+
+def _get_parameters_from_cell(ibrav, structure, tolerance): # pylint: disable=too-many-statements,too-many-branches
+ """
+ Get the A, B, C, cosAB, cosAB, cosBC from a given `ibrav` and
+ StructureData. Only the parameters necessary for the given ibrav
+ are returned.
+
+ :raise QEInputValidationError:
+ If an invalid `ibrav` is passed, or the cell reconstructed
+ from the parameters does not match the input cell.
+ """
+ import numpy as np
+ import scipy.linalg as la
+
+ A, B, C = ('a', 'b', 'c') # pylint: disable=invalid-name
+ cosAB, cosAC, cosBC = ('cosab', 'cosac', 'cosbc') # pylint: disable=invalid-name
+
+ cell_input = np.array(structure.get_attribute('cell'))
+ v1, v2, v3 = cell_input # pylint: disable=invalid-name
+
+ if ibrav == 1:
+ parameters = {A: la.norm(v1)}
+ elif ibrav == 2:
+ parameters = {A: np.sqrt(2) * la.norm(v1)}
+ elif ibrav in [-3, 3]:
+ parameters = {A: 2 * la.norm(v1) / 3}
+ elif ibrav == 4:
+ parameters = {A: la.norm(v1), C: la.norm(v3)}
+ elif ibrav in [5, -5]:
+ parameters = {A: la.norm(v1)}
+ parameters[cosAB] = np.dot(v1, v2) / parameters[A]** 2
+ elif ibrav == 6:
+ parameters = {A: la.norm(v1), C: la.norm(v3)}
+ elif ibrav == 7:
+ parameters = {A: np.sqrt(2) * la.norm(v1[:2]), C: 2 * v1[-1]}
+ elif ibrav == 8:
+ parameters = {A: la.norm(v1), B: la.norm(v2), C: la.norm(v3)}
+ elif ibrav in [9, -9]:
+ parameters = {A: 2 * abs(v1[0]), B: 2 * abs(v1[1]), C: la.norm(v3)}
+ elif ibrav == 91:
+ parameters = {A: la.norm(v1), B: la.norm(v2 + v3), C: la.norm(v3 - v2)}
+ elif ibrav == 10:
+ parameters = {A: 2 * v1[0], B: 2 * v2[1], C: 2 * v1[2]}
+ elif ibrav == 11:
+ parameters = {A: 2 * v1[0], B: 2 * v1[1], C: 2 * v1[2]}
+ elif ibrav == 12:
+ parameters = {A: la.norm(v1), B: la.norm(v2), C: la.norm(v3)}
+ parameters[cosAB] = np.dot(v1, v2) / (parameters[A] * parameters[B])
+ elif ibrav == -12:
+ parameters = {A: la.norm(v1), B: la.norm(v2), C: la.norm(v3)}
+ parameters[cosAC] = np.dot(v1, v3) / (parameters[A] * parameters[C])
+ elif ibrav == 13:
+ parameters = {A: 2 * v1[0], B: la.norm(v2), C: 2 * v3[2]}
+ parameters[cosAB] = 2 * np.dot(v1, v2) / (parameters[A] * parameters[B])
+ elif ibrav == -13:
+ parameters = {A: 2 * abs(v1[0]), B: 2 * abs(v1[1]), C: la.norm(v3)}
+ parameters[cosAC] = 2 * np.dot(v1, v2) / (parameters[A] * parameters[C])
+ elif ibrav == 14:
+ parameters = {A: la.norm(v1), B: la.norm(v2), C: la.norm(v3)}
+ parameters[cosAB] = np.dot(v1, v2) / (parameters[A] * parameters[B])
+ parameters[cosAC] = np.dot(v1, v3) / (parameters[A] * parameters[C])
+ parameters[cosBC] = np.dot(v2, v3) / (parameters[B] * parameters[C])
+ else:
+ raise QEInputValidationError("The given 'ibrav' value '{}' is not understood.".format(ibrav))
+
+ # Check that the reconstructed cell matches the input
+ system_dict = {'ibrav': ibrav}
+ system_dict.update(parameters)
+ cell_reconstructed = get_cell_from_parameters(
+ cell_parameters=None, # this is only used for ibrav=0
+ system_dict=system_dict,
+ alat=parameters[A],
+ using_celldm=False
+ )
+ if not np.allclose(cell_reconstructed, cell_input, rtol=0, atol=tolerance):
+ raise QEInputValidationError(
+ "The cell {} constructed with 'ibrav={}' does not match the input cell{}.".format(
+ cell_reconstructed, ibrav, cell_input
+ )
+ )
+
+ return parameters

aiida_quantumespresso/calculations/cp.py

@@ -55,7 +55,6 @@ class CpCalculation(BasePwCpInputGenerator):
  ('CONTROL', 'pseudo_dir'), # set later
  ('CONTROL', 'outdir'), # set later
  ('CONTROL', 'prefix'), # set later
- ('SYSTEM', 'ibrav'), # set later
  ('SYSTEM', 'celldm'),
  ('SYSTEM', 'nat'), # set later
  ('SYSTEM', 'ntyp'), # set later

aiida_quantumespresso/calculations/helpers/__init__.py

@@ -163,7 +163,6 @@ def pw_input_helper(input_params, structure, stop_at_first_error=False, flat_mod
  i.lower() for i in [
  'pseudo_dir',
  'outdir',
- 'ibrav',
  'celldm',
  'nat',
  'ntyp',

aiida_quantumespresso/calculations/pw.py

@@ -30,7 +30,6 @@ class PwCalculation(BasePwCpInputGenerator):
  ('CONTROL', 'pseudo_dir'),
  ('CONTROL', 'outdir'),
  ('CONTROL', 'prefix'),
- ('SYSTEM', 'ibrav'),
  ('SYSTEM', 'celldm'),
  ('SYSTEM', 'nat'),
  ('SYSTEM', 'ntyp'),

aiida_quantumespresso/tools/pwinputparser.py

@@ -100,9 +100,6 @@ def create_builder_from_file(input_folder, input_file_name, code, metadata, pseu
  builder.structure = parsed_file.get_structuredata()
  builder.kpoints = parsed_file.get_kpointsdata()
 
- if parsed_file.namelists['SYSTEM']['ibrav'] != 0:
- raise NotImplementedError('Found ibrav != 0: `aiida-quantumespresso` currently only supports ibrav = 0.')
-
  # Then, strip the namelist items that the plugin doesn't allow or sets later.
  # NOTE: If any of the position or cell units are in alat or crystal
  # units, that will be taken care of by the input parsing tools, and

setup.json

@@ -87,7 +87,9 @@
  "matplotlib<3.0.0; python_version<'3'",
  "packaging",
  "qe-tools==1.1.3",
- "xmlschema==1.0.13"
+ "xmlschema==1.0.13",
+ "numpy",
+ "scipy"
  ],
  "license": "MIT License",
  "name": "aiida_quantumespresso",

tests/calculations/test_pw.py

@@ -52,3 +52,56 @@ def test_pw_default(
  # Checks on the files written to the sandbox folder as raw input
  assert sorted(fixture_sandbox.get_content_list()) == sorted(['aiida.in', 'pseudo', 'out'])
  file_regression.check(input_written, encoding='utf-8', extension='.in')
+
+
+def test_pw_ibrav(
+ aiida_profile, fixture_sandbox, generate_calc_job, fixture_code, generate_kpoints_mesh, generate_upf_data,
+ file_regression
+):
+ """Test a default `PwCalculation`."""
+ entry_point_name = 'quantumespresso.pw'
+
+ parameters = {'CONTROL': {'calculation': 'scf'}, 'SYSTEM': {'ecutrho': 240.0, 'ecutwfc': 30.0, 'ibrav': 2}}
+
+ # The structure needs to be rotated in the same way QE does it for ibrav=2.
+ param = 5.43
+ cell = [[-param / 2., 0, param / 2.], [0, param / 2., param / 2.], [-param / 2., param / 2., 0]]
+ structure = orm.StructureData(cell=cell)
+ structure.append_atom(position=(0., 0., 0.), symbols='Si', name='Si')
+ structure.append_atom(position=(param / 4., param / 4., param / 4.), symbols='Si', name='Si')
+
+ upf = generate_upf_data('Si')
+ inputs = {
+ 'code': fixture_code(entry_point_name),
+ 'structure': structure,
+ 'kpoints': generate_kpoints_mesh(2),
+ 'parameters': orm.Dict(dict=parameters),
+ 'pseudos': {
+ 'Si': upf
+ },
+ 'metadata': {
+ 'options': get_default_options()
+ }
+ }
+
+ calc_info = generate_calc_job(fixture_sandbox, entry_point_name, inputs)
+
+ cmdline_params = ['-in', 'aiida.in']
+ local_copy_list = [(upf.uuid, upf.filename, u'./pseudo/Si.upf')]
+ retrieve_list = ['aiida.out', './out/aiida.save/data-file-schema.xml', './out/aiida.save/data-file.xml']
+ retrieve_temporary_list = [['./out/aiida.save/K*[0-9]/eigenval*.xml', '.', 2]]
+
+ # Check the attributes of the returned `CalcInfo`
+ assert isinstance(calc_info, datastructures.CalcInfo)
+ assert sorted(calc_info.cmdline_params) == sorted(cmdline_params)
+ assert sorted(calc_info.local_copy_list) == sorted(local_copy_list)
+ assert sorted(calc_info.retrieve_list) == sorted(retrieve_list)
+ assert sorted(calc_info.retrieve_temporary_list) == sorted(retrieve_temporary_list)
+ assert sorted(calc_info.remote_symlink_list) == sorted([])
+
+ with fixture_sandbox.open('aiida.in') as handle:
+ input_written = handle.read()
+
+ # Checks on the files written to the sandbox folder as raw input
+ assert sorted(fixture_sandbox.get_content_list()) == sorted(['aiida.in', 'pseudo', 'out'])
+ file_regression.check(input_written, encoding='utf-8', extension='.in')

tests/calculations/test_pw/test_pw_ibrav.in

@@ -0,0 +1,24 @@
+&CONTROL
+ calculation = 'scf'
+ outdir = './out/'
+ prefix = 'aiida'
+ pseudo_dir = './pseudo/'
+ verbosity = 'high'
+/
+&SYSTEM
+ a = 5.4300000000d+00
+ ecutrho = 2.4000000000d+02
+ ecutwfc = 3.0000000000d+01
+ ibrav = 2
+ nat = 2
+ ntyp = 1
+/
+&ELECTRONS
+/
+ATOMIC_SPECIES
+Si 28.0855 Si.upf
+ATOMIC_POSITIONS angstrom
+Si 0.0000000000 0.0000000000 0.0000000000 
+Si 1.3575000000 1.3575000000 1.3575000000 
+K_POINTS automatic
+2 2 2 0 0 0

tests/tools/test_immigrate.py

@@ -3,6 +3,8 @@
 from __future__ import absolute_import
 import os
 
+import numpy as np
+
 from aiida_quantumespresso.tools.pwinputparser import create_builder_from_file
 
 
@@ -48,10 +50,68 @@ def test_create_builder(aiida_profile, fixture_sandbox, fixture_code, generate_u
  },
  'SYSTEM': {
  'ecutrho': 240.0,
- 'ecutwfc': 30.0
+ 'ecutwfc': 30.0,
+ 'ibrav': 0,
+ }
+ }
+ assert 'kpoints' in builder
+ assert 'structure' in builder
+
+ generate_calc_job(fixture_sandbox, entry_point_name, builder)
+
+
+def test_create_builder_nonzero_ibrav(
+ aiida_profile, fixture_sandbox, fixture_code, generate_upf_data, generate_calc_job
+):
+ """Test the `create_builder_from_file` method that parses an existing `pw.x` folder into a process builder.
+
+ The input file used is the one generated for `tests.calculations.test_pw.test_pw_ibrav`.
+ """
+ entry_point_name = 'quantumespresso.pw'
+ code = fixture_code(entry_point_name)
+
+ metadata = {
+ 'options': {
+ 'resources': {
+ 'num_machines': 1,
+ 'num_mpiprocs_per_machine': 32,
+ },
+ 'max_memory_kb': 1000,
+ 'max_wallclock_seconds': 60 * 60 * 12,
+ 'withmpi': True,
+ }
+ }
+
+ in_foldername = os.path.join('tests', 'calculations', 'test_pw')
+ in_folderpath = os.path.abspath(in_foldername)
+
+ upf_foldername = os.path.join('tests', 'fixtures', 'pseudos')
+ upf_folderpath = os.path.abspath(upf_foldername)
+ si_upf = generate_upf_data('Si')
+ si_upf.store()
+
+ builder = create_builder_from_file(
+ in_folderpath, 'test_pw_ibrav.in', code, metadata, upf_folderpath, use_first=True
+ )
+
+ assert builder['code'] == code
+ assert builder['metadata'] == metadata
+ assert builder['pseudos']['Si'].id == si_upf.id
+ assert builder['parameters'].get_dict() == {
+ 'CONTROL': {
+ 'calculation': 'scf',
+ 'verbosity': 'high'
+ },
+ 'SYSTEM': {
+ 'ecutrho': 240.0,
+ 'ecutwfc': 30.0,
+ 'ibrav': 2,
  }
  }
  assert 'kpoints' in builder
  assert 'structure' in builder
+ param = 5.43
+ np.testing.assert_allclose([[-param / 2., 0, param / 2.], [0, param / 2., param / 2.], [-param / 2., param / 2., 0]
+ ], builder.structure.get_attribute('cell'))
 
  generate_calc_job(fixture_sandbox, entry_point_name, builder)