First working version

aiida_cp2k/utils/__init__.py

@@ -19,3 +19,4 @@
 from .workchains import check_resize_unit_cell
 from .workchains import resize_unit_cell
 from .workchains import HARTREE2EV, HARTREE2KJMOL
+from .workchains import seekpath_structure_analysis, update_input_dict_for_bands

aiida_cp2k/utils/workchains.py

@@ -13,10 +13,47 @@
 HARTREE2EV = 27.211399
 HARTREE2KJMOL = 2625.500
 
+VAL_ELEC = {
+    "H": 1,
+    "He": 2,
+    "Li": 3,
+    "Be": 4,
+    "B": 3,
+    "C": 4,
+    "N": 5,
+    "O": 6,
+    "F": 7,
+    "Ne": 8,
+    "Na": 9,
+    "Mg": 2,
+    "Al": 3,
+    "Si": 4,
+    "P": 5,
+    "S": 6,
+    "Cl": 7,
+    "Ar": 8,
+    "K": 9,
+    "Ca": 10,
+    "Sc": 11,
+    "Ti": 12,
+    "V": 13,
+    "Cr": 14,
+    "Mn": 15,
+    "Fe": 16,
+    "Co": 17,
+    "Ni": 18,
+    "Cu": 19,
+    "Zn": 12,
+    "Zr": 12,
+}
+
 
 def merge_dict(dct, merge_dct):
-    """ Taken from https://gist.github.com/angstwad/bf22d1822c38a92ec0a9
-    Recursive dict merge. Inspired by :meth:``dict.update()``, instead of
+    """Recursive dict merge.
+
+    Taken from https://gist.github.com/angstwad/bf22d1822c38a92ec0a9
+
+    Inspired by :meth:``dict.update()``, instead of
     updating only top-level keys, merge_dict recurses down into dicts nested
     to an arbitrary depth, updating keys. The ``merge_dct`` is merged into
     ``dct``.
@@ -33,16 +70,15 @@ def merge_dict(dct, merge_dct):
 
 
 @calcfunction
-def merge_Dict(d1, d2):  #pylint: disable=invalid-name
-    """ Make all the data in the second Dict overwrite the corrisponding data in the first Dict"""
-    d1_dict = d1.get_dict()
-    d2_dict = d2.get_dict()
-    merge_dict(d1_dict, d2_dict)
-    return Dict(dict=d1_dict)
+def merge_Dict(dict1, dict2):  # pylint: disable=invalid-name
+    """Make all the data in the second Dict overwrite the corrisponding data in the first Dict."""
+    result = dict1.get_dict()
+    merge_dict(result, dict2.get_dict())
+    return Dict(dict=result)
 
 
 def get_kinds_section(structure, protocol_settings):
-    """ Write the &KIND sections given the structure and the settings_dict"""
+    """Write the &KIND sections given the structure and the settings_dict"""
     kinds = []
     all_atoms = set(structure.get_ase().get_chemical_symbols())
     for atom in all_atoms:
@@ -154,3 +190,50 @@ def resize_unit_cell(struct, resize):
     """Resize the StructureData according to the resize Dict"""
     resize_tuple = tuple([resize[x] for x in ['nx', 'ny', 'nz']])
     return StructureData(ase=struct.get_ase().repeat(resize_tuple))
+
+
+def add_condband(structure):
+    """Add 20% of conduction bands to the CP2K input. If 20 % is 0, then add only one."""
+    total = 0
+    for symbol in structure.get_ase().get_chemical_symbols():
+        total += VAL_ELEC[symbol]
+    added_mos = total // 10  # 20% of conduction band
+    if added_mos == 0:
+        added_mos = 1
+    return added_mos
+
+
+def update_input_dict_for_bands(input_dict, seekpath, structure):
+    """Insert kpoint path into the input dictonary of CP2K."""
+
+    i_dict = input_dict.get_dict()
+
+    path = seekpath.dict['path']
+    coords = seekpath.dict['point_coords']
+
+    kpath = []
+    for pnt in path:
+        pnt1 = pnt[0] + ' ' + " ".join(str(x) for x in coords[pnt[0]])
+        pnt1 = pnt[1] + ' ' + " ".join(str(x) for x in coords[pnt[1]])
+        kpath.append({'_': "", 'UNITS': 'B_VECTOR', 'NPOINTS': 10, 'SPECIAL_POINT': [pnt1, pnt1]})
+
+    kpath_dict = {'FORCE_EVAL': {'DFT': {'PRINT': {'BAND_STRUCTURE': {'KPOINT_SET': kpath}}}}}
+    merge_dict(i_dict, kpath_dict)
+
+    added_mos = {'FORCE_EVAL': {'DFT': {'SCF': {'ADDED_MOS': add_condband(structure)}}}}
+    merge_dict(i_dict, added_mos)
+
+    return Dict(dict=i_dict)
+
+
+@calcfunction
+def seekpath_structure_analysis(structure, parameters):
+    """This calcfunction will take a structure and pass it through SeeKpath to get the
+    primitive cell and the path of high symmetry k-points through its Brillouin zone.
+
+    Note that the returned primitive cell may differ from the original structure in
+    which case the k-points are only congruent with the primitive cell.
+    """
+
+    from aiida.tools import get_kpoints_path
+    return get_kpoints_path(structure, **parameters.get_dict())

aiida_cp2k/workchains/__init__.py

@@ -8,3 +8,4 @@
 """AiiDA-CP2K workchains"""
 
 from .base import Cp2kBaseWorkChain
+from .bands import Cp2kBandsWorkChain

aiida_cp2k/workchains/bands.py

@@ -0,0 +1,83 @@
+# -*- coding: utf-8 -*-
+###############################################################################
+# Copyright (c), The AiiDA-CP2K authors.                                      #
+# SPDX-License-Identifier: MIT                                                #
+# AiiDA-CP2K is hosted on GitHub at https://github.com/cp2k/aiida-cp2k        #
+# For further information on the license, see the LICENSE.txt file.           #
+###############################################################################
+"""Work chain to compute a band structure using CP2K."""
+from __future__ import absolute_import
+import six
+
+from aiida.common import AttributeDict
+from aiida.engine import WorkChain, ToContext
+from aiida.orm import BandsData, Dict, StructureData
+from aiida.plugins import WorkflowFactory
+
+from aiida_cp2k.utils import seekpath_structure_analysis, update_input_dict_for_bands
+
+Cp2kBaseWorkChain = WorkflowFactory('cp2k.base')  # pylint: disable=invalid-name
+
+
+class Cp2kBandsWorkChain(WorkChain):
+    """Compute Band Structure of a material."""
+
+    @classmethod
+    def define(cls, spec):
+        super(Cp2kBandsWorkChain, cls).define(spec)
+        spec.expose_inputs(Cp2kBaseWorkChain,
+                           namespace='cp2k_base',
+                           exclude=('cp2k.structure', 'cp2k.parameters', 'cp2k.metadata.options.parser_name'))
+        spec.input('structure', valid_type=StructureData)
+        spec.input("parameters", valid_type=Dict)
+        spec.input('cp2k_base.cp2k.metadata.options.parser_name',
+                   valid_type=six.string_types,
+                   default='cp2k_advanced_parser',
+                   non_db=True,
+                   help='Parser of the calculation: the default is cp2k_advanced_parser to get the bands.')
+        spec.outline(
+            cls.setup,
+            cls.run_seekpath,
+            cls.prepare_bands_calculation,
+            cls.run_bands_calculation,
+            cls.return_results,
+        )
+        spec.output('output_bands', valid_type=BandsData)
+
+    def setup(self):
+        """Perform initial setup."""
+        self.ctx.structure = self.inputs.structure
+
+    def run_seekpath(self):
+        """Run Seekpath to get the primitive structure
+        N.B. If, after cell optimization the symmetry change,
+        the primitive cell will be different!"""
+
+        seekpath_parameters = Dict(dict={})
+
+        seekpath_result = seekpath_structure_analysis(self.ctx.structure, seekpath_parameters)
+        self.ctx.seekpath_analysis = seekpath_result['parameters']
+        self.ctx.structure = seekpath_result['primitive_structure']
+
+    def prepare_bands_calculation(self):
+        """Prepare all the neccessary input links to run the calculation."""
+
+        # Add molecular orbitals and kpoints path that was generated by seekpath
+        self.ctx.parameters = update_input_dict_for_bands(self.inputs.parameters, self.ctx.seekpath_analysis,
+                                                          self.ctx.structure)
+
+    def run_bands_calculation(self):
+        """Run cp2k calculation."""
+
+        inputs = AttributeDict(self.exposed_inputs(Cp2kBaseWorkChain, namespace='cp2k_base'))
+        inputs.cp2k.structure = self.ctx.structure
+        inputs.cp2k.parameters = self.ctx.parameters
+
+        # Create the calculation process and launch it
+        running = self.submit(Cp2kBaseWorkChain, **inputs)
+        self.report("Submitted Cp2kBaseWorkChain for band structure calculation.")
+        return ToContext(calculation=running)
+
+    def return_results(self):
+        """Extract output_parameters."""
+        self.out('output_bands', self.ctx.calculation.outputs.output_bands)