Adding ENVIRON namelists

aiida_quantumespresso/calculations/__init__.py

@@ -238,6 +238,11 @@ def prepare_for_submission(self, folder):
         if environ_namelist is not None:
             if not isinstance(environ_namelist, dict):
                 raise exceptions.InputValidationError('ENVIRON namelist should be specified as a dictionary')
+
+            ## Check for BOUNDARY and ELECTROSTATIC Keys as well
+            boundary_namelist = settings.pop('BOUNDARY', None)
+            electrostatic_namelist = settings.pop('ELECTROSTATIC', None)
+
             # We first add the environ flag to the command-line options (if not already present)
             try:
                 if '-environ' not in settings['CMDLINE']:
@@ -255,6 +260,18 @@ def prepare_for_submission(self, folder):
                     handle.write(convert_input_to_namelist_entry(key, value, mapping=mapping_species))
                 handle.write('/\n')
 
+                if boundary_namelist is not None:
+                    handle.write('&BOUNDARY\n')
+                    for key, value in sorted(boundary_namelist.items()):
+                        handle.write(convert_input_to_namelist_entry(key, value, mapping=mapping_species))
+                    handle.write('/\n')
+
+                if electrostatic_namelist is not None:
+                    handle.write('&ELECTROSTATIC\n')
+                    for key, value in sorted(electrostatic_namelist.items()):
+                        handle.write(convert_input_to_namelist_entry(key, value, mapping=mapping_species))
+                    handle.write('/\n')
+
         # Check for the deprecated 'ALSO_BANDS' setting and if present fire a deprecation log message
         also_bands = settings.pop('ALSO_BANDS', None)
         if also_bands:

aiida_quantumespresso/parsers/parse_raw/pw.py

@@ -743,6 +743,15 @@ def parse_stdout(stdout, input_parameters, parser_options=None, parsed_xml=None)
                     parsed_data['fermi_energy' + units_suffix] = default_energy_units
                 except Exception:
                     logs.warning.append('Error while parsing Fermi energy from the output file.')
+            elif 'Gaussian-smeared nuclei' in line:
+                try:
+                    value_potential_shift = float(line.split()[-2])
+                    unit_potential_shift = line.split()[-1]
+                    parsed_data['environ_potential_shift'] = value_potential_shift
+                    parsed_data['environ_unit_potential_shift'] = unit_potential_shift
+                except Exception:
+                    logs.warning.append('Not an Environ calculation with an open boundary condition.')
+                    pass
 
             elif 'Forces acting on atoms' in line:
                 try:

docs/source/user_guide/get_started/examples/environ_calculation.py

@@ -0,0 +1,146 @@
+# -*- coding: utf-8 -*-
+
+from aiida import orm
+from aiida.engine import run, submit
+import numpy as np
+from ase import Atoms
+from ase import units
+"""
+Calculation that recreates the environ calculation of
+CO on Pt(111) with a unit charge
+
+More details about the solvation models and potential corrections
+can be found here:
+
+O. Andreussi, I. Dabo and N. Marzari, J. Chem. Phys. 136,064102 (2012).
+I. Dabo et al. Phys.Rev. B 77, 115139 (2008)
+
+This system is taken directly from Example04 in the
+examples folder from Environ
+"""
+
+
+def calculator():
+    """
+    This is the SCF calculator, which is identical
+    and unchanged from what one would do in the case
+    of a standard SCF calculation - nothing environ
+    related. Remember to add in the tot_charge here
+    """
+    param_dict = {
+        'CONTROL': {
+            'calculation': 'scf',
+            'tprnfor': True,
+        },
+        'SYSTEM': {
+            'ecutwfc': 35,
+            'ecutrho': 280,
+            'occupations': 'smearing',
+            'smearing': 'cold',
+            'degauss': 0.03,
+            'tot_charge': 1,
+        },
+        'ELECTRONS': {
+            'conv_thr': 5e-9,
+            'electron_maxstep': 200,
+            'mixing_beta': 0.2,
+            'mixing_ndim': 15,
+            'diagonalization': 'david',
+        },
+    }
+
+    return param_dict
+
+
+def environ_calculator():
+    """
+    Environ parameters
+    """
+    param_dict = {
+        'ENVIRON': {
+            'verbose': 0,
+            'environ_thr': 1.0,
+            'environ_type': 'input',
+            'env_static_permittivity': 1,
+            'env_surface_tension': 0.0,
+            'env_pressure': 0.0,
+            'env_electrostatic': True,
+        },
+        'BOUNDARY': {
+            'solvent_mode': 'full',
+        },
+        'ELECTROSTATIC': {
+            'pbc_correction': 'parabolic',
+            'pbc_dim': 2,
+            'pbc_axis': 3,
+            'tol': 5e-13,
+        }
+    }
+    return param_dict
+
+
+def runner(code, structure):
+
+    ## Base calculation for testing
+    PwBaseWorkChain = WorkflowFactory('quantumespresso.pw.base')
+    builder = PwBaseWorkChain.get_builder()
+    builder.pw.structure = structure
+
+    builder.metadata.label = 'Single point environ'
+    builder.metadata.description = 'Testing calculation with environ for Pt(111) + CO'
+
+    KpointsData = DataFactory('array.kpoints')
+    kpoints = KpointsData()
+    kpoints.set_kpoints_mesh([1, 1, 1])  # Not physical just a test
+    builder.kpoints = kpoints
+
+    family = load_group('SSSP/1.1/PBE/efficiency')
+    builder.pw.pseudos = family.get_pseudos(structure=structure)
+
+    calculation = calculator()
+    environ_calculation = environ_calculator()
+    ## these are the main cube files that could potentially be parsed
+    ## if the verbosity is set to 2 or higher
+    # environ_calculation['additional_retrieve_list'] = ['epsilon.cube', \
+    #     'vreference.cube', 'velectrostatic.cube', 'vsoftcavity.cube', \
+    #     'electrons.cube', 'charges.cube', 'smeared_ions.cube']
+
+    builder.pw.parameters = orm.Dict(dict=calculation)
+    builder.pw.settings = orm.Dict(dict=environ_calculation)
+    builder.pw.metadata.options.resources = {'num_machines': 1}
+    builder.pw.metadata.options.max_wallclock_seconds = 5 * 60
+    builder.pw.code = code
+
+    calculation = submit(builder)
+
+
+if __name__ == '__main__':
+    code = load_code('pw_6-7@juwels')
+
+    StructureData = DataFactory('structure')
+    ## these are the original coordinates for the Pt-CO system
+    positions = [
+        [5.335084148, 4.646723426, 12.901029877],
+        [5.335009643, 4.619623254, 15.079854269],
+        [8.061327071, 0.098057998, 8.992142901],
+        [2.608989366, 0.098058283, 8.992140585],
+        [0.000036609, 4.720846294, 8.968756935],
+        [5.335159557, 4.721612729, 9.380196435],
+        [0.000041121, 7.802951963, 4.604626508],
+        [5.335161233, 7.697749113, 4.753489408],
+        [2.697860636, 3.152173889, 4.688412329],
+        [7.972463687, 3.152174491, 4.688415209],
+    ]
+
+    ## setting up the system with ASE
+    ## notice the units that are being used
+    atoms = Atoms('COPt8')
+    atoms.set_positions(np.array(positions) * units.Bohr)
+    atoms.set_pbc([True, True, True])
+    a = 10.6881 * units.Bohr
+    b = 0.866025 * a * units.Bohr
+    c = 3.95422 * a * units.Bohr
+    atoms.set_cell([a, b, c])
+    structure = StructureData(ase=atoms)
+
+    runner(code=code, structure=structure)

docs/source/user_guide/get_started/examples/pw_tutorial.rst

@@ -623,6 +623,32 @@ More options are available, and can be explored by expanding
 ``builder.metadata.options.`` + ``TAB``.
 
 
+Environ calculations
+--------------------
+
+Environ calculations can be run by supplying the required namelists to ``builder.settings``. For example, to run a calculation with a continuum dielectric of water above a 2D slab with the surface normal in the z-direction create the following dictionary.
+
+::
+
+    environ_param_dict = {
+        'ENVIRON':{
+            'environ_type': 'input',
+            'env_static_permittivity': 78.36,
+            'env_electrostatic': True,
+        },
+        'BOUNDARY':{
+            'solvent_mode':'full',
+        },
+        'ELECTROSTATIC':{
+            'pbc_correction':'parabolic',
+            'pbc_dim': 2,
+            'pbc_axis': 3,
+        }
+    }
+
+This dictionary can be passed on to the builder as ::
+
+    builder.pw.settings = orm.Dict(dict=environ_param_dict)
 
 
 Launching the calculation

tests/calculations/test_pw.py

@@ -250,3 +250,33 @@ def test_pw_parallelization_duplicate_cmdline_flag(fixture_sandbox, generate_cal
     with pytest.raises(InputValidationError) as exc:
         generate_calc_job(fixture_sandbox, entry_point_name, inputs)
     assert 'Conflicting' in str(exc.value)
+
+
+def test_environ_namelists(fixture_sandbox, generate_calc_job, generate_inputs_pw, file_regression):
+    """Test that Environ namelists are  created."""
+    entry_point_name = 'quantumespresso.pw'
+
+    inputs = generate_inputs_pw()
+    inputs['settings'] = orm.Dict(
+        dict={
+            'ENVIRON': {
+                'electrolyte_linearized': True,
+                'environ_type': 'input',
+            },
+            'BOUNDARY': {
+                'solvent_mode': 'electronic',
+                'electrolyte_mode': 'electronic',
+            },
+            'ELECTROSTATIC': {
+                'pbc_correction': 'parabolic',
+            }
+        },
+    )
+    generate_calc_job(fixture_sandbox, entry_point_name, inputs)
+
+    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', 'environ.in'])
+    file_regression.check(input_written, encoding='utf-8', extension='.in')

tests/calculations/test_pw/test_environ_namelists.in

@@ -0,0 +1,27 @@
+&CONTROL
+  calculation = 'scf'
+  outdir = './out/'
+  prefix = 'aiida'
+  pseudo_dir = './pseudo/'
+  verbosity = 'high'
+/
+&SYSTEM
+  ecutrho =   2.4000000000d+02
+  ecutwfc =   3.0000000000d+01
+  ibrav = 0
+  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
+CELL_PARAMETERS angstrom
+      2.7150000000       2.7150000000       0.0000000000
+      2.7150000000       0.0000000000       2.7150000000
+      0.0000000000       2.7150000000       2.7150000000

tests/parsers/test_pw.py

@@ -907,3 +907,35 @@ def test_pw_vcrelax_failed_not_converged_nstep(
     assert 'output_parameters' in results
     assert 'output_structure' in results
     assert 'output_trajectory' in results
+
+
+def test_environ(
+    fixture_localhost,
+    generate_calc_job_node,
+    generate_parser,
+    generate_inputs,
+    data_regression,
+):
+    """Test a simple Environ calculation."""
+    name = 'default'
+    entry_point_calc_job = 'quantumespresso.pw'
+    entry_point_parser = 'quantumespresso.pw'
+    environ_settings = {'ENVIRON': {'environ_type': 'water'}}
+    node = generate_calc_job_node(
+        entry_point_calc_job, fixture_localhost, name, generate_inputs(settings=environ_settings)
+    )
+    parser = generate_parser(entry_point_parser)
+    results, calcfunction = parser.parse_from_node(node, store_provenance=False)
+
+    assert calcfunction.is_finished, calcfunction.exception
+    assert calcfunction.is_finished_ok, calcfunction.exit_message
+    assert not orm.Log.objects.get_logs_for(node), [log.message for log in orm.Log.objects.get_logs_for(node)]
+    assert 'output_kpoints' in results
+    assert 'output_parameters' in results
+    assert 'output_trajectory' in results
+
+    data_regression.check({
+        'output_kpoints': results['output_kpoints'].attributes,
+        'output_parameters': results['output_parameters'].get_dict(),
+        'output_trajectory': results['output_trajectory'].attributes,
+    })