diff --git a/aiida_quantumespresso/calculations/__init__.py b/aiida_quantumespresso/calculations/__init__.py index d17852038..855d450e7 100644 --- a/aiida_quantumespresso/calculations/__init__.py +++ b/aiida_quantumespresso/calculations/__init__.py @@ -4,6 +4,7 @@ import abc import os + import six from six.moves import zip @@ -11,9 +12,12 @@ 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 diff --git a/aiida_quantumespresso/calculations/cp.py b/aiida_quantumespresso/calculations/cp.py index 833f9b931..3cbd4ffd1 100644 --- a/aiida_quantumespresso/calculations/cp.py +++ b/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 diff --git a/aiida_quantumespresso/calculations/helpers/__init__.py b/aiida_quantumespresso/calculations/helpers/__init__.py index 33ae2a18c..04a8c4028 100644 --- a/aiida_quantumespresso/calculations/helpers/__init__.py +++ b/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', diff --git a/aiida_quantumespresso/calculations/pw.py b/aiida_quantumespresso/calculations/pw.py index a7792c190..dd1c79ea8 100644 --- a/aiida_quantumespresso/calculations/pw.py +++ b/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'), diff --git a/aiida_quantumespresso/tools/pwinputparser.py b/aiida_quantumespresso/tools/pwinputparser.py index 839d86aca..d2351e394 100644 --- a/aiida_quantumespresso/tools/pwinputparser.py +++ b/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 diff --git a/setup.json b/setup.json index f337a3707..c1b23d35e 100644 --- a/setup.json +++ b/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", diff --git a/tests/calculations/test_pw.py b/tests/calculations/test_pw.py index 6a723542c..bef00c3e8 100644 --- a/tests/calculations/test_pw.py +++ b/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') diff --git a/tests/calculations/test_pw/test_pw_ibrav.in b/tests/calculations/test_pw/test_pw_ibrav.in new file mode 100644 index 000000000..dc5efacd7 --- /dev/null +++ b/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 diff --git a/tests/tools/test_immigrate.py b/tests/tools/test_immigrate.py index 87f2416c7..5d5f767a8 100644 --- a/tests/tools/test_immigrate.py +++ b/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)