Merge 462b7da into 59c5b25

.flake8

@@ -1,2 +1,2 @@
 [flake8]
-ignore = E501  # lines too long
+max-line-length = 121

.gitignore

@@ -16,3 +16,5 @@ postgres*.log
 _archive/
 .ipynb_checkpoints/
 databases/
+test_workdir/
+_archive/

.style.yapf

@@ -0,0 +1,3 @@
+[style]
+based_on_style = google
+column_limit = 120

.vscode/settings.json

@@ -15,15 +15,19 @@
     },
 
     "python.pythonPath": "/anaconda/envs/aiida-workshop/bin/python",
+    "python.formatting.provider": "yapf",
     "python.linting.pylintEnabled": false,
     "python.linting.flake8Enabled": true,
     "python.linting.enabled": true,
-    "python.testing.pyTestArgs": [
+    "python.testing.pytestArgs": [
         "aiida_crystal17"
     ],
     "python.testing.unittestEnabled": false,
     "python.testing.nosetestsEnabled": false,
-    "python.testing.pyTestEnabled": true,
+    "python.testing.pytestEnabled": false,
     "restructuredtext.preview.sphinx.disabled": true,
-    "restructuredtext.confPath": "${workspaceFolder}/docs/source"
+    "restructuredtext.confPath": "${workspaceFolder}/docs/source",
+    "cSpell.words": [
+        "aiida"
+    ]
 }

MANIFEST.in

@@ -1,4 +1,4 @@
 include LICENSE
 include setup.json
-recursive-include aiida_crystal17/tests/ *
+recursive-include aiida_crystal17/tests/raw_files *
 include aiida_crystal17/validation/*.json

aiida_crystal17/__init__.py

@@ -4,4 +4,4 @@
 AiiDA plugin for running the CRYSTAL17 code
 """
 
-__version__ = "0.6.0b3"
+__version__ = "0.9.0b5"

aiida_crystal17/calcfunctions/band_gap.py

@@ -0,0 +1,193 @@
+from collections import namedtuple
+import traceback
+
+import numpy as np
+
+from aiida.orm import ArrayData, Dict, Float, List
+from aiida.engine import calcfunction, ExitCode
+
+BandResult = namedtuple('BandResult',
+                        ['fermi', 'left_edge', 'right_edge', 'non_zero_fermi'])
+
+
+def calculate_band_gap(energies,
+                       densities,
+                       fermi=0,
+                       dtol=1e-8,
+                       try_fshifts=(),
+                       missing_edge=None):
+    """calculate the band gap, given an energy vs density plot
+
+    Parameters
+    ----------
+    energies : list[float]
+    densities : list[float]
+    fermi : float
+    dtol : float
+        tolerance for checking if density is zero
+    try_fshifts : tuple[float]
+        if the density at the fermi energy is non-zero,
+        try shifting the fermi energy by these values, until a non-zero density is found.
+        Useful for dealing with band edges at the fermi energy
+    missing_edge : object
+        the value to return if an edge cannot be determind
+
+    Returns
+    -------
+    BandResult
+
+    """
+    energies = np.array(energies, float)
+    densities = np.abs(np.array(densities, float))
+    if not len(energies) == len(densities):
+        raise AssertionError("the energies and densities arrays are of different lengths")
+    if not fermi < energies.max():
+        raise AssertionError("the energies range does not contain the fermi energy")
+    if not fermi > energies.min():
+        raise AssertionError("the energies range does not contain the fermi energy")
+
+    # sort energies
+    order = np.argsort(energies)
+    energies = energies[order]
+    densities = densities[order]
+
+    # find index closest to fermi
+    fermi_idx = (np.abs(energies - fermi)).argmin()
+
+    # check density isn't non-zero at fermi
+    fermi_non_zero = (densities[fermi_idx] > 0 + dtol)
+
+    # if the density at the fermi is non-zero, try shifting the fermi
+    # (useful to deal with band edges at the fermi)
+    for fshift in try_fshifts:
+        if not fermi_non_zero:
+            break
+        new_index = np.abs(energies - (fermi + fshift)).argmin()
+        if densities[new_index] <= 0 + dtol:
+            fermi_non_zero = False
+            fermi_idx = new_index
+
+    if fermi_non_zero:
+        return BandResult(
+            energies[fermi_idx], missing_edge, missing_edge, True)
+
+    # find left edge
+    found_left = False
+    for left_idx in reversed(range(0, fermi_idx + 1)):
+        if densities[left_idx] > 0 + dtol:
+            found_left = True
+            break
+
+    # find right edge
+    found_right = False
+    for right_idx in range(fermi_idx, len(densities)):
+        if densities[right_idx] > 0 + dtol:
+            found_right = True
+            break
+
+    return BandResult(energies[fermi_idx],
+                      energies[left_idx] if found_left else missing_edge,
+                      energies[right_idx] if found_right else missing_edge,
+                      False)
+
+
+@calcfunction
+def calcfunction_band_gap(doss_results, doss_array, dtol=None, try_fshifts=None):
+    """calculate the band gap, given DoS data computed by CryDossCalculation
+
+    Parameters
+    ----------
+    doss_array : aiida.orm.ArrayData
+    dtol : aiida.orm.Float
+        tolerance for checking if density is zero
+    try_fshifts : aiida.orm.List
+        if the density at the fermi energy is non-zero,
+        try shifting the fermi energy by these values, until a non-zero density is found.
+        Useful for dealing with band edges at the fermi energy
+
+    """
+    if not isinstance(doss_results, Dict):
+        return ExitCode(
+            101, 'doss_results is not of type `aiida.orm.Dict`: {}'.format(doss_results))
+    if "fermi_energy" not in doss_results.get_dict():
+        return ExitCode(
+            102, '`fermi_energy` not in doss_results')
+    if "energy_units" not in doss_results.get_dict():
+        return ExitCode(
+            102, '`energy_units` not in doss_results')
+    if not isinstance(doss_array, ArrayData):
+        return ExitCode(
+            103, 'doss_array is not of type `aiida.orm.ArrayData`: {}'.format(doss_array))
+
+    kwargs = {
+        "fermi": doss_results.get_dict()["fermi_energy"]
+    }
+
+    if dtol is not None:
+        if not isinstance(dtol, Float):
+            return ExitCode(
+                104, 'dtol is not of type `aiida.orm.Float`: {}'.format(dtol))
+        kwargs["dtol"] = dtol.value
+
+    if try_fshifts is not None:
+        if not isinstance(try_fshifts, List):
+            return ExitCode(
+                105, 'try_fshifts is not of type `aiida.orm.List`: {}'.format(try_fshifts))
+        kwargs["try_fshifts"] = try_fshifts.get_list()
+
+    array_names = doss_array.get_arraynames()
+    if "energies" not in array_names:
+        return ExitCode(
+            111, 'doss_array does not contain array `energies`: {}'.format(doss_array))
+    if "total" in array_names:
+        if "total_alpha" in array_names and "total_beta" in array_names:
+            return ExitCode(
+                112, ('doss_array does not contains both array `total` and '
+                      '`total_alpha`, `total_beta`: {}'.format(doss_array)))
+    elif "total_alpha" in array_names and "total_beta" in array_names:
+        if "total" in array_names:
+            return ExitCode(
+                112, ('doss_array does not contains both array `total` and '
+                      '`total_alpha`, `total_beta`: {}'.format(doss_array)))
+    else:
+        return ExitCode(
+            113, 'doss_array does not contain array `total` or `total_alpha` and `total_beta`: {}'.format(doss_array))
+
+    if "total" in array_names:
+        calcs = {'total': doss_array.get_array('total')}
+    else:
+        alpha_density = doss_array.get_array('total_alpha')
+        beta_density = doss_array.get_array('total_beta')
+        total_density = np.abs(alpha_density) + np.abs(beta_density)
+        calcs = {'alpha': alpha_density, 'beta': beta_density, 'total': total_density}
+
+    final_dict = {
+        "energy_units": doss_results.get_dict()["energy_units"]
+    }
+
+    for name, density in calcs.items():
+        try:
+            result = calculate_band_gap(
+                doss_array.get_array('energies'),
+                density,
+                **kwargs)
+        except Exception:
+            traceback.print_exc()
+            return ExitCode(201, 'calculate_band_gap failed')
+        if result.non_zero_fermi:
+            bandgap = 0.
+        elif result.left_edge is None or result.right_edge is None:
+            bandgap = None
+        else:
+            bandgap = result.right_edge - result.left_edge
+        final_dict.update(
+            {
+                name+'_fermi': result.fermi,
+                name+'_left_edge': result.left_edge,
+                name+'_right_edge': result.right_edge,
+                name+'_zero_fermi': not result.non_zero_fermi,
+                name+'_bandgap': bandgap
+            }
+        )
+
+    return {"results": Dict(dict=final_dict)}

aiida_crystal17/calculations/cry_abstract.py

@@ -1,8 +1,10 @@
 """
 Plugin to create a CRYSTAL17 output file from a supplied input file.
 """
+import os
 import six
 
+from aiida.common.datastructures import (CalcInfo, CodeInfo)
 from aiida.engine import CalcJob
 from aiida.plugins import DataFactory
 
@@ -23,9 +25,7 @@ def define(cls, spec):
         super(CryAbstractCalculation, cls).define(spec)
 
         spec.input('metadata.options.input_file_name',
-                   valid_type=six.string_types, default='main.d12')
-        spec.input('metadata.options.external_file_name',
-                   valid_type=six.string_types, default='main.gui')
+                   valid_type=six.string_types, default='INPUT')
         spec.input('metadata.options.output_main_file_name',
                    valid_type=six.string_types, default='main.out')
 
@@ -40,42 +40,58 @@ def define(cls, spec):
             message='The retrieved folder data node could not be accessed.')
         spec.exit_code(
             210, 'ERROR_OUTPUT_FILE_MISSING',
-            message='the main output file was not found')
+            message='the main (stdout) output file was not found')
+        spec.exit_code(
+            211, 'ERROR_TEMP_FOLDER_MISSING',
+            message='the temporary retrieved folder was not found')
 
         # Unrecoverable errors: required retrieved files could not be read, parsed or are otherwise incomplete
         spec.exit_code(
-            300, 'ERROR_OUTPUT_PARSING',
+            300, 'ERROR_PARSING_STDOUT',
             message=('An error was flagged trying to parse the '
-                     'main crystal output file'))
+                     'crystal exec stdout file'))
+        spec.exit_code(  # TODO is this an unrecoverable error?
+            301, 'ERROR_PARSING_OPTIMISATION_GEOMTRIES',
+            message=("An error occurred parsing the 'opta'/'optc' geomerty files"))
 
         spec.exit_code(
             350, 'ERROR_CRYSTAL_INPUT',
-            message='the input file was could not be read by CRYSTAL')
+            message='the input file could not be read by CRYSTAL')
         spec.exit_code(
             351, 'ERROR_WAVEFUNCTION_NOT_FOUND',
             message='CRYSTAL could not find the required wavefunction file')
 
         # Significant errors but calculation can be used to restart
         spec.exit_code(
-            401, 'UNCONVERGED_SCF',
+            400, 'ERROR_OUT_OF_WALLTIME',
+            message='The calculation stopped prematurely because it ran out of walltime.')
+        spec.exit_code(
+            401, 'ERROR_OUT_OF_MEMORY',
+            message='The calculation stopped prematurely because it ran out of memory.')
+        spec.exit_code(
+            402, 'ERROR_OUT_OF_VMEMORY',
+            message='The calculation stopped prematurely because it ran out of virtual memory.')
+
+        spec.exit_code(
+            411, 'UNCONVERGED_SCF',
             message='SCF convergence did not finalise (usually due to reaching step limit)')
         spec.exit_code(
-            402, 'UNCONVERGED_GEOMETRY',
+            412, 'UNCONVERGED_GEOMETRY',
             message='Geometry convergence did not finalise (usually due to reaching step limit)')
         spec.exit_code(
-            410, 'ERROR_SCF_ABNORMAL_END',
-            message='an error was encountered during an SCF computation')
-        spec.exit_code(
-            411, 'BASIS_SET_LINEARLY_DEPENDENT',
+            413, 'BASIS_SET_LINEARLY_DEPENDENT',
             message='an error encountered usually during geometry optimisation')
         spec.exit_code(
-            412, 'ERROR_MPI_ABORT',
+            414, 'ERROR_SCF_ABNORMAL_END',
+            message='an error was encountered during an SCF computation')
+        spec.exit_code(
+            415, 'ERROR_MPI_ABORT',
             message='an unknown error was encountered, causing the MPI to abort')
         spec.exit_code(
             499, 'ERROR_CRYSTAL_RUN',
             message='The main crystal output file flagged an unhandled error')
 
-        # errors in symmetry node consistency check
+        # errors in symmetry node consistency checks
         spec.exit_code(
             510, 'ERROR_SYMMETRY_INCONSISTENCY',
             message=('inconsistency in the input and output symmetry'))
@@ -97,8 +113,43 @@ def define(cls, spec):
                     required=False,
                     help='the symmetry data from the calculation')
 
-        # TODO retrieve .f9 / .f98 from remote folder (for GUESSP or RESTART)
-        # spec.input(
-        #     'parent_folder', valid_type=RemoteData, required=False,
-        #     help=('Use a remote folder as parent folder (for '
-        #           'restarts and similar.'))
+    def create_calc_info(
+            self, tempfolder,
+            local_copy_list=None, remote_copy_list=None, remote_symlink_list=None,
+            retrieve_list=None, retrieve_temporary_list=None):
+        """Prepare CalcInfo object for aiida,
+        to describe how the computation will be executed and recovered
+        """
+        # Prepare CodeInfo object for aiida,
+        # describes how a code has to be executed
+        codeinfo = CodeInfo()
+        codeinfo.code_uuid = self.inputs.code.uuid
+        if self.metadata.options.withmpi:
+            # parallel versions of crystal (Pcrystal, Pproperties & MPPcrystal)
+            # read data specifically from a file called INPUT
+            if self.metadata.options.input_file_name != "INPUT":
+                tempfolder.insert_path(
+                    os.path.join(tempfolder.abspath,
+                                 self.metadata.options.input_file_name),
+                    dest_name="INPUT",
+                )
+        else:
+            codeinfo.stdin_name = self.metadata.options.input_file_name
+        codeinfo.stdout_name = self.metadata.options.output_main_file_name
+        # serial version output to stdout, but parallel version output to stderr!
+        # so we join the files
+        codeinfo.join_files = True
+        codeinfo.cmdline_params = []
+        codeinfo.withmpi = self.metadata.options.withmpi
+
+        # Prepare CalcInfo object for aiida
+        calcinfo = CalcInfo()
+        calcinfo.uuid = self.uuid
+        calcinfo.codes_info = [codeinfo]
+        calcinfo.local_copy_list = local_copy_list or []
+        calcinfo.remote_copy_list = remote_copy_list or []
+        calcinfo.remote_symlink_list = remote_symlink_list or []
+        calcinfo.retrieve_list = retrieve_list or []
+        calcinfo.retrieve_temporary_list = retrieve_temporary_list or []
+
+        return calcinfo