Feature/general pdos

aiida_nanotech_empa/plugins/overlap.py

@@ -7,10 +7,14 @@ def define(cls, spec):
         super().define(spec)
         spec.input("parameters", valid_type=orm.Dict, help="Overlap input parameters")
         spec.input(
-            "parent_slab_folder", valid_type=orm.RemoteData, help="slab scf folder"
+            "parent_all_folder",
+            valid_type=orm.RemoteData,
+            help="whole system scf folder",
         )
         spec.input(
-            "parent_mol_folder", valid_type=orm.RemoteData, help="molecule scf folder"
+            "parent_fragment_folder",
+            valid_type=orm.RemoteData,
+            help="fragment scf folder",
         )
         spec.input("settings", valid_type=orm.Dict, help="special settings")
 
@@ -56,10 +60,10 @@ def prepare_for_submission(self, folder):
         calcinfo.retrieve_list = settings.pop("additional_retrieve_list", [])
 
         # Symlinks.
-        if "parent_slab_folder" in self.inputs:
-            comp_uuid = self.inputs.parent_slab_folder.computer.uuid
-            remote_path = self.inputs.parent_slab_folder.get_remote_path()
-            copy_info = (comp_uuid, remote_path, "parent_slab_folder/")
+        if "parent_all_folder" in self.inputs:
+            comp_uuid = self.inputs.parent_all_folder.computer.uuid
+            remote_path = self.inputs.parent_all_folder.get_remote_path()
+            copy_info = (comp_uuid, remote_path, "parent_all_folder/")
             if (
                 self.inputs.code.computer.uuid == comp_uuid
             ):  # if running on the same computer - make a symlink
@@ -68,10 +72,10 @@ def prepare_for_submission(self, folder):
             else:
                 calcinfo.remote_copy_list.append(copy_info)
 
-        if "parent_mol_folder" in self.inputs:
-            comp_uuid = self.inputs.parent_mol_folder.computer.uuid
-            remote_path = self.inputs.parent_mol_folder.get_remote_path()
-            copy_info = (comp_uuid, remote_path, "parent_mol_folder/")
+        if "parent_fragment_folder" in self.inputs:
+            comp_uuid = self.inputs.parent_fragment_folder.computer.uuid
+            remote_path = self.inputs.parent_fragment_folder.get_remote_path()
+            copy_info = (comp_uuid, remote_path, "parent_fragment_folder/")
             if (
                 self.inputs.code.computer.uuid == comp_uuid
             ):  # if running on the same computer - make a symlink

aiida_nanotech_empa/workflows/cp2k/diag_workchain.py

@@ -88,7 +88,10 @@ def setup(self):
         # Get initial magnetization.
         magnetization_per_site = [0 for i in range(len(self.inputs.structure.sites))]
         if "uks" in self.ctx.dft_params and self.ctx.dft_params["uks"]:
-            magnetization_per_site = self.ctx.dft_params["magnetization_per_site"]
+            try:
+                magnetization_per_site = self.ctx.dft_params["magnetization_per_site"]
+            except KeyError:
+                pass
 
         structure_with_tags, kinds_dict = cp2k_utils.determine_kinds(
             structure, magnetization_per_site

aiida_nanotech_empa/workflows/cp2k/pdos_workchain.py

@@ -3,7 +3,7 @@
 import numpy as np
 from aiida import engine, orm, plugins
 
-from ...utils import common_utils, split_structure
+from ...utils import common_utils
 from . import cp2k_utils
 
 Cp2kDiagWorkChain = plugins.WorkflowFactory("nanotech_empa.cp2k.diag")
@@ -20,12 +20,18 @@ def define(cls, spec):
         spec.input("overlap_code", valid_type=orm.Code)
 
         # Structures.
+        # Whole system
         spec.input(
             "structure",
             valid_type=orm.StructureData,
             help="Coordinates of the whole system.",
         )
-        spec.input("molecule_indices", valid_type=orm.List)
+        # Portion of the system, could also be not derived from the whole system
+        spec.input(
+            "fragment_structure",
+            valid_type=orm.StructureData,
+            help="Coordinates of the fragment system.",
+        )
 
         spec.input("pdos_lists", valid_type=orm.List)
 
@@ -66,62 +72,60 @@ def define(cls, spec):
 
     def setup(self):
         self.report("Setting up workchain")
-        structure_generator = split_structure.split_structure(
-            structure=self.inputs.structure,
-            fixed_atoms=[],
-            magnetization_per_site=self.inputs.dft_params["magnetization_per_site"]
-            if "magnetization_per_site" in self.inputs.dft_params
-            else None,
-            fragments={"molecule": self.inputs.molecule_indices},
-        )
 
-        self.ctx.n_slab_atoms = len(self.inputs.structure.sites)
+        self.ctx.n_all_atoms = len(self.inputs.structure.sites)
         emax = float(self.inputs.overlap_params.get_dict()["--emax1"])
         nlumo = int(self.inputs.overlap_params.get_dict()["--nlumo2"])
 
         dft_parameters = self.inputs.dft_params.get_dict()
         charges = dft_parameters.pop("charges")
         multiplicities = dft_parameters.pop("multiplicities", {})
-
+        magnetization_per_site = dft_parameters.pop("magnetization_per_site", {})
         # Set up DFT parameters of the whole system.
-        slab_info = next(structure_generator)
-        self.ctx.structure = slab_info["structure"]
+        self.ctx.structure = self.inputs.structure
         self.ctx.dft_parameters = copy.deepcopy(dft_parameters)
         self.ctx.dft_parameters["charge"] = charges["all"]
         if "all" in multiplicities:
             self.ctx.dft_parameters["multiplicity"] = multiplicities["all"]
+        if "all" in magnetization_per_site:
+            self.ctx.dft_parameters["magnetization_per_site"] = magnetization_per_site[
+                "all"
+            ]
         self.ctx.dft_parameters["added_mos"] = np.max(
-            [100, int(1.2 * self.ctx.n_slab_atoms * emax / 5.0)]
+            [100, int(1.2 * self.ctx.n_all_atoms * emax / 5.0)]
         )
 
-        # Use the same cutoff for molecule and slab.
+        # Use the same cutoff for fragment and whole system.
         self.ctx.dft_parameters["cutoff"] = cp2k_utils.get_cutoff(self.ctx.structure)
 
-        # Set up molecular DFT parameters.
-        molecule_info = next(structure_generator)
-        self.ctx.molecule_structure = molecule_info["structure"]
-        self.ctx.mol_dft_parameters = copy.deepcopy(self.ctx.dft_parameters)
-        self.ctx.mol_dft_parameters["charge"] = charges["molecule"]
-        if "molecule" in multiplicities:
-            self.ctx.mol_dft_parameters["multiplicity"] = multiplicities["molecule"]
-        self.ctx.mol_dft_parameters["added_mos"] = nlumo + 2
-        self.ctx.mol_dft_parameters[
+        # Set up fragment DFT parameters.
+        self.ctx.fragment_structure = self.inputs.fragment_structure
+        self.ctx.fragment_dft_parameters = copy.deepcopy(self.ctx.dft_parameters)
+        self.ctx.fragment_dft_parameters["charge"] = charges["fragment"]
+        if "fragment" in multiplicities:
+            self.ctx.fragment_dft_parameters["multiplicity"] = multiplicities[
+                "fragment"
+            ]
+        if "fragment" in magnetization_per_site:
+            self.ctx.fragment_dft_parameters[
+                "magnetization_per_site"
+            ] = magnetization_per_site["fragment"]
+        self.ctx.fragment_dft_parameters["added_mos"] = nlumo + 2
+        self.ctx.fragment_dft_parameters["cutoff"] = self.ctx.dft_parameters["cutoff"]
+        self.ctx.fragment_dft_parameters[
             "elpa_switch"
         ] = False  # Elpa can cause problems with small systems
-        self.ctx.mol_dft_parameters["magnetization_per_site"] = molecule_info[
-            "magnetization_per_site"
-        ]
 
     def run_diags(self):
-        # Slab part.
-        self.report("Running Diag Workchain for slab")
+        # Whole system  part.
+        self.report("Running Diag Workchain for whole system")
         builder = Cp2kDiagWorkChain.get_builder()
         builder.cp2k_code = self.inputs.cp2k_code
         builder.structure = self.ctx.structure
         builder.protocol = self.inputs.protocol
         builder.dft_params = orm.Dict(self.ctx.dft_parameters)
         builder.settings = orm.Dict({"additional_retrieve_list": ["*.pdos"]})
-        builder.options = orm.Dict(self.inputs.options["slab"])
+        builder.options = orm.Dict(self.inputs.options["all"])
 
         # Restart WFN.
         if "parent_calc_folder" in self.inputs:
@@ -131,30 +135,38 @@ def run_diags(self):
         if self.inputs.pdos_lists is not None:
             builder.pdos_lists = orm.List([pdos[0] for pdos in self.inputs.pdos_lists])
 
-        self.to_context(slab_diag_scf=self.submit(builder))
+        self.to_context(all_diag_scf=self.submit(builder))
 
-        # Molecule part.
-        self.report("Running Diag Workchain for molecule")
+        # Fragment part.
+        self.report("Running Diag Workchain for fragment")
         builder = Cp2kDiagWorkChain.get_builder()
         builder.cp2k_code = self.inputs.cp2k_code
-        builder.structure = self.ctx.molecule_structure
+        builder.structure = self.ctx.fragment_structure
         builder.protocol = self.inputs.protocol
-        builder.dft_params = orm.Dict(self.ctx.mol_dft_parameters)
-        builder.options = orm.Dict(self.inputs.options["molecule"])
-        self.to_context(mol_diag_scf=self.submit(builder))
+        builder.dft_params = orm.Dict(self.ctx.fragment_dft_parameters)
+        builder.options = orm.Dict(self.inputs.options["fragment"])
+        self.to_context(fragment_diag_scf=self.submit(builder))
 
     def run_overlap(self):
-        for calculation in [self.ctx.slab_diag_scf, self.ctx.mol_diag_scf]:
+        for calculation in [self.ctx.all_diag_scf, self.ctx.fragment_diag_scf]:
             if not common_utils.check_if_calc_ok(self, calculation):
                 return self.exit_codes.ERROR_TERMINATION
         self.report("Running overlap")
         builder = OverlapCalculation.get_builder()
         builder.code = self.inputs.overlap_code
         builder.parameters = self.inputs.overlap_params
-        builder.parent_slab_folder = self.ctx.slab_diag_scf.outputs.remote_folder
-        builder.parent_mol_folder = self.ctx.mol_diag_scf.outputs.remote_folder
+        builder.parent_all_folder = self.ctx.all_diag_scf.outputs.remote_folder
+        builder.parent_fragment_folder = (
+            self.ctx.fragment_diag_scf.outputs.remote_folder
+        )
 
-        n_machines = 4 if self.ctx.n_slab_atoms < 2000 else 8
+        # set n_machines to 1 4 8 if n_all_atoms < 100 2000 or else
+        if self.ctx.n_all_atoms < 100:
+            n_machines = 1
+        elif self.ctx.n_all_atoms < 2000:
+            n_machines = 4
+        else:
+            n_machines = 8
 
         builder.metadata = {
             "label": "overlap",
@@ -174,7 +186,7 @@ def finalize(self):
         ]:
             self.report("Overlap calculation did not finish correctly")
             return self.exit_codes.ERROR_TERMINATION
-        self.out("slab_retrieved", self.ctx.slab_diag_scf.outputs.retrieved)
+        self.out("slab_retrieved", self.ctx.all_diag_scf.outputs.retrieved)
 
         # Add the workchain uuid to the input structure extras.
         common_utils.add_extras(self.inputs.structure, "surfaces", self.node.uuid)

examples/workflows/c2h2_for_pdos.xyz

@@ -0,0 +1,6 @@
+4
+Lattice="8.84399664 0.0 0.0 0.0 10.21216768 0.0 0.0 0.0 20.0" Properties=species:S:1:pos:R:3:masses:R:1 pbc="T T T"
+H        2.89445002       3.11134486       9.51675909       1.00800000
+C        3.40914678       4.04942122       9.51675909      12.01100000
+C        4.06762878       5.24274647       9.51675909      12.01100000
+H        4.58232554       6.18082283       9.51675909       1.00800000

examples/workflows/example_cp2k_geo_opt.py

@@ -18,8 +18,7 @@ def _example_cp2k_geo_opt(cp2k_code, sys_type, uks, n_nodes, n_cores_per_node):
         filters={"label": {"in": GEOS}},
     )
     structures = {}
-    for node_tuple in qb.iterall():
-        node = node_tuple[0]
+    for node in qb.all(flat=True):
         structures[node.label] = node
     for required in GEOS:
         if required in structures:

examples/workflows/example_cp2k_pdos.py

@@ -7,71 +7,74 @@
 Cp2kPdosWorkChain = plugins.WorkflowFactory("nanotech_empa.cp2k.pdos")
 
 DATA_DIR = pathlib.Path(__file__).parent.absolute()
-GEO_FILE = "c2h2_on_au111.xyz"
+GEOS = ["c2h2_on_au111.xyz", "c2h2_for_pdos.xyz"]
 
 
 def _example_cp2k_pdos(
     cp2k_code, overlap_code, sc_diag, force_multiplicity, uks, n_nodes, n_cores_per_node
 ):
     # Check test geometry is already in database.
     qb = orm.QueryBuilder()
-    qb.append(orm.Node, filters={"label": {"in": [GEO_FILE]}})
-    structure = None
-    for node_tuple in qb.iterall():
-        node = node_tuple[0]
-        structure = node
-    if structure is not None:
-        print(f"Found existing structure: {structure.pk}")
-    else:
-        structure = orm.StructureData(ase=ase.io.read(DATA_DIR / GEO_FILE))
-        structure.label = GEO_FILE
-        structure.store()
-        print(f"Created new structure: {structure.pk}")
+    qb.append(orm.Node, filters={"label": {"in": GEOS}})
+    structures = {}
+    for node in qb.all(flat=True):
+        structures[node.label] = node
+    for required in GEOS:
+        if required in structures:
+            print("found existing structure: ", required, structures[required].pk)
+        else:
+            structure = orm.StructureData(ase=ase.io.read(DATA_DIR / required))
+            structure.label = required
+            structure.store()
+            structures[required] = structure
+            print("created new structure: ", required, structure.pk)
 
     builder = Cp2kPdosWorkChain.get_builder()
 
     builder.metadata.label = "CP2K_PDOS"
-    builder.metadata.description = "test description"
     builder.cp2k_code = cp2k_code
-    ase_geom_slab = ase.io.read(DATA_DIR / GEO_FILE)
-    ase_geom_mol = ase_geom_slab[0:4]
-    builder.slabsys_structure = orm.StructureData(ase=ase_geom_slab)
-    builder.mol_structure = orm.StructureData(ase=ase_geom_mol)
-    builder.pdos_lists = orm.List([("1..4", "molecule"), ("1", "cat")])
+    builder.structure = structures["c2h2_on_au111.xyz"]
+    builder.fragment_structure = structures["c2h2_for_pdos.xyz"]
+    builder.pdos_lists = orm.List([("1..4", "C2H2"), ("1", "C_at")])
     builder.protocol = orm.Str("debug")
     if uks:
+        builder.metadata.description = "automatic test PDOS UKS"
         dft_params = {
             "sc_diag": sc_diag,
             "force_multiplicity": force_multiplicity,
             "elpa_switch": False,
             "periodic": "XYZ",
             "uks": uks,
-            "multiplicity": 1,
+            "multiplicities": {"all": 1, "fragment": 1},
+            "magnetization_per_site": {"fragment": [0, 1, -1, 0]},
+            "charges": {"all": 0, "fragment": 0},
             "smear_t": 150,
             "spin_up_guess": [0],
             "spin_dw_guess": [1],
         }
     else:
+        builder.metadata.description = "automatic test PDOS RKS"
         dft_params = {
             "sc_diag": sc_diag,
             "force_multiplicity": force_multiplicity,
             "elpa_switch": False,
             "periodic": "XYZ",
             "uks": uks,
+            "charges": {"all": 0, "fragment": 0},
             "smear_t": 150,
         }
     builder.dft_params = orm.Dict(dft_params)
 
     builder.options = {
-        "slab": {
+        "all": {
             "max_wallclock_seconds": 600,
             "resources": {
                 "num_machines": n_nodes,
                 "num_mpiprocs_per_machine": n_cores_per_node,
                 "num_cores_per_mpiproc": 1,
             },
         },
-        "molecule": {
+        "fragment": {
             "max_wallclock_seconds": 600,
             "resources": {
                 "num_machines": n_nodes,
@@ -84,16 +87,16 @@ def _example_cp2k_pdos(
     builder.overlap_code = overlap_code
     builder.overlap_params = orm.Dict(
         {
-            "--cp2k_input_file1": "parent_slab_folder/aiida.inp",
-            "--basis_set_file1": "parent_slab_folder/BASIS_MOLOPT",
-            "--xyz_file1": "parent_slab_folder/aiida.coords.xyz",
-            "--wfn_file1": "parent_slab_folder/aiida-RESTART.wfn",
+            "--cp2k_input_file1": "parent_all_folder/aiida.inp",
+            "--basis_set_file1": "parent_all_folder/BASIS_MOLOPT",
+            "--xyz_file1": "parent_all_folder/aiida.coords.xyz",
+            "--wfn_file1": "parent_all_folder/aiida-RESTART.wfn",
             "--emin1": "-2",
             "--emax1": "2",
-            "--cp2k_input_file2": "parent_mol_folder/aiida.inp",
-            "--basis_set_file2": "parent_mol_folder/BASIS_MOLOPT",
-            "--xyz_file2": "parent_mol_folder/aiida.coords.xyz",
-            "--wfn_file2": "parent_mol_folder/aiida-RESTART.wfn",
+            "--cp2k_input_file2": "parent_fragment_folder/aiida.inp",
+            "--basis_set_file2": "parent_fragment_folder/BASIS_MOLOPT",
+            "--xyz_file2": "parent_fragment_folder/aiida.coords.xyz",
+            "--wfn_file2": "parent_fragment_folder/aiida-RESTART.wfn",
             "--nhomo2": "2",
             "--nlumo2": "2",
             "--output_file": "./overlap.npz",
@@ -129,8 +132,18 @@ def run_all(cp2k_code, overlap_code, n_nodes, n_cores_per_node):
         sc_diag=False,
         force_multiplicity=True,
         uks=False,
-        n_nodes=n_nodes,
-        n_cores_per_node=n_cores_per_node,
+        n_nodes=1,
+        n_cores_per_node=1,
+    )
+    print("#### sc_diag UKS magnetization guess only on fragment")
+    _example_cp2k_pdos(
+        orm.load_code(cp2k_code),
+        orm.load_code(overlap_code),
+        sc_diag=True,
+        force_multiplicity=True,
+        uks=True,
+        n_nodes=1,
+        n_cores_per_node=1,
     )