First working version of Bands work chain

aiida_cp2k/utils/__init__.py

@@ -23,3 +23,4 @@
 from .workchains import check_resize_unit_cell
 from .workchains import resize_unit_cell
 from .workchains import HARTREE2EV
+from .workchains import seekpath_structure_analysis, update_input_dict_for_bands_calc

aiida_cp2k/utils/workchains.py

@@ -13,6 +13,40 @@
 
 HARTREE2EV = 27.211399
 
+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
@@ -154,3 +188,46 @@ 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):
+    total = 0
+    for s in structure.get_ase().get_chemical_symbols():
+        total += val_elec[s]
+    added_mos = total // 10  # 20% of conduction band
+    if added_mos == 0:
+        added_mos = 1
+    return added_mos
+
+
+def update_input_dict_for_bands_calc(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 p in path:
+        p1 = p[0] + ' ' + " ".join(str(x) for x in coords[p[0]])
+        p2 = p[1] + ' ' + " ".join(str(x) for x in coords[p[1]])
+        kpath.append({'_': "", 'UNITS': 'B_VECTOR', 'NPOINTS': 10, 'SPECIAL_POINT': [p1, p2]})
+
+    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

@@ -9,3 +9,4 @@
 
 from .base import Cp2kBaseWorkChain
 from .multistage import Cp2kMultistageWorkChain
+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_calc
+
+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_calc(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)

examples/workchains/fixme_bands.py

@@ -0,0 +1,184 @@
+# -*- coding: utf-8 -*-
+# pylint: disable=invalid-name
+###############################################################################
+# 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.           #
+###############################################################################
+"""Run simple Band Structure calculation"""
+
+from __future__ import print_function
+from __future__ import absolute_import
+
+import os
+import sys
+import click
+import numpy as np
+from ase.atoms import Atoms
+
+from aiida.engine import run
+from aiida.orm import (Code, Dict, SinglefileData, StructureData)
+from aiida.common import NotExistent
+from aiida.plugins import WorkflowFactory
+
+Cp2kBandsWorkChain = WorkflowFactory('cp2k.bands')
+
+
+def example_bands(cp2k_code):
+    """Run simple Band Structure calculation"""
+
+    print("Computing Band Structure of Si...")
+
+    pwd = os.path.dirname(os.path.realpath(__file__))
+
+    # structure
+    positions = [
+        [0.0000000000, 0.0000000000, 2.6954627656],
+        [4.0431941484, 4.0431941484, 4.0431941484],
+    ]
+    cell = [
+        [0.0, 2.69546276561, 2.69546276561],
+        [2.69546276561, 0.0, 2.69546276561],
+        [2.69546276561, 2.69546276561, 0.0],
+    ]
+    atoms = Atoms('Si2', positions=positions, cell=cell, pbc=[True, True, True])
+    structure = StructureData(ase=atoms)
+
+    # basis set
+    basis_file = SinglefileData(file=os.path.join(pwd, "..", "files", "BASIS_MOLOPT"))
+
+    # pseudopotentials
+    pseudo_file = SinglefileData(file=os.path.join(pwd, "..", "files", "GTH_POTENTIALS"))
+
+    # parameters
+    parameters = Dict(
+        dict={
+            'FORCE_EVAL': {
+                'METHOD': 'Quickstep',
+                'DFT': {
+                    'CHARGE': 0,
+                    'KPOINTS': {
+                        'SCHEME MONKHORST-PACK': '1 1 1',
+                        'SYMMETRY': 'OFF',
+                        'WAVEFUNCTIONS': 'REAL',
+                        'FULL_GRID': '.TRUE.',
+                        'PARALLEL_GROUP_SIZE': 0,
+                    },
+                    'MGRID': {
+                        'CUTOFF': 600,
+                        'NGRIDS': 4,
+                        'REL_CUTOFF': 50,
+                    },
+                    'UKS': False,
+                    'BASIS_SET_FILE_NAME': 'BASIS_MOLOPT',
+                    'POTENTIAL_FILE_NAME': 'GTH_POTENTIALS',
+                    'QS': {
+                        'METHOD': 'GPW',
+                        'EXTRAPOLATION': 'USE_GUESS'
+                    },
+                    'POISSON': {
+                        'PERIODIC': 'XYZ',
+                    },
+                    'SCF': {
+                        'EPS_SCF': 1.0E-4,
+                        'SMEAR': {
+                            'METHOD': 'FERMI_DIRAC',
+                            'ELECTRONIC_TEMPERATURE': 300
+                        },
+                        'DIAGONALIZATION': {
+                            'ALGORITHM': 'STANDARD',
+                            'EPS_ADAPT': 0.01
+                        },
+                        'MIXING': {
+                            'METHOD': 'BROYDEN_MIXING',
+                            'ALPHA': 0.2,
+                            'BETA': 1.5,
+                            'NBROYDEN': 8
+                        },
+                    },
+                    'XC': {
+                        'XC_FUNCTIONAL': {
+                            '_': 'PBE',
+                        },
+                    },
+                },
+                'SUBSYS': {
+                    'KIND': [{
+                        '_': 'Si',
+                        'BASIS_SET': 'DZVP-MOLOPT-SR-GTH',
+                        'POTENTIAL': 'GTH-LDA-q4'
+                    },]
+                },
+                'PRINT': {  # this is to print forces (may be necessary for problems
+                    # detection)
+                    'FORCES': {
+                        '_': 'ON',
+                    }
+                },
+            },
+            'GLOBAL': {
+                "EXTENDED_FFT_LENGTHS": True,  # Needed for large systems
+            }
+        })
+
+    # Construct process builder
+    builder = Cp2kBandsWorkChain.get_builder()
+    builder.structure = structure
+    builder.parameters = parameters
+    builder.cp2k_base.cp2k.code = cp2k_code
+    builder.cp2k_base.cp2k.file = {
+        'basis': basis_file,
+        'pseudo': pseudo_file,
+    }
+
+    builder.cp2k_base.cp2k.metadata.options.resources = {
+        "num_machines": 1,
+        "num_mpiprocs_per_machine": 1,
+    }
+
+    builder.cp2k_base.cp2k.metadata.options.max_wallclock_seconds = 1 * 3 * 60
+
+    print("submitted calculation...")
+    calc = run(builder)
+
+    bands = calc['output_bands']
+
+    # check bands
+    expected_gamma_kpoint = np.array([-5.71237757, 6.5718575, 6.5718575, 6.5718575, 8.88653953])
+
+    if bands.get_kpoints().shape == (66, 3):
+        print("OK, got expected kpoints set size.")
+    else:
+        print("Got unexpected kpoints set.")
+        sys.exit(3)
+
+    if bands.get_bands().shape == (66, 5):
+        print("OK, got expected bands set size.")
+    else:
+        print("Got unexpected bands set.")
+        sys.exit(3)
+
+    if abs(max(bands.get_bands()[0] - expected_gamma_kpoint)) < 1e-7:
+        print("Ok, got expected energy levels at GAMMA point.")
+    else:
+        print("Got unexpected energy levels at GAMMA point.")
+        sys.exit(3)
+
+    sys.exit(0)
+
+
+@click.command('cli')
+@click.argument('codelabel')
+def cli(codelabel):
+    """Click interface"""
+    try:
+        code = Code.get_from_string(codelabel)
+    except NotExistent:
+        print("The code '{}' does not exist".format(codelabel))
+        sys.exit(1)
+    example_bands(code)
+
+
+if __name__ == '__main__':
+    cli()  # pylint: disable=no-value-for-parameter

setup.json

@@ -30,7 +30,8 @@
         ],
         "aiida.workflows": [
             "cp2k.base = aiida_cp2k.workchains:Cp2kBaseWorkChain",
-            "cp2k.multistage= aiida_cp2k.workchains:Cp2kMultistageWorkChain"
+            "cp2k.multistage = aiida_cp2k.workchains:Cp2kMultistageWorkChain",
+            "cp2k.bands = aiida_cp2k.workchains:Cp2kBandsWorkChain"
         ]
     },
     "setup_requires": ["reentry"],