From ffe1872b180ce840e8e7a9ce7559934dfb5fe908 Mon Sep 17 00:00:00 2001
From: Carlo Antonio Pignedoli <c.pignedoli@gmail.com>
Date: Wed, 10 Jul 2024 07:39:45 +0000
Subject: [PATCH 1/9] work in progress

---
 aiida_nanotech_empa/utils/cycle_tools.py      | 121 +++++
 aiida_nanotech_empa/utils/nmr.py              | 161 +++++++
 .../workflows/gaussian/nics_workchain.py      |  85 ++++
 .../workflows/gaussian/oldnics_workchain.py   | 266 +++++++++++
 examples/workflows/aa                         | 431 ------------------
 examples/workflows/benzene.xyz                |  14 +
 examples/workflows/example_gaussian_opt.py    |  55 +++
 7 files changed, 702 insertions(+), 431 deletions(-)
 create mode 100644 aiida_nanotech_empa/utils/cycle_tools.py
 create mode 100644 aiida_nanotech_empa/utils/nmr.py
 create mode 100644 aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
 create mode 100644 aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py
 delete mode 100755 examples/workflows/aa
 create mode 100644 examples/workflows/benzene.xyz
 create mode 100644 examples/workflows/example_gaussian_opt.py

diff --git a/aiida_nanotech_empa/utils/cycle_tools.py b/aiida_nanotech_empa/utils/cycle_tools.py
new file mode 100644
index 0000000..f4bed8d
--- /dev/null
+++ b/aiida_nanotech_empa/utils/cycle_tools.py
@@ -0,0 +1,121 @@
+import numpy as np
+
+import scipy as sp
+import scipy.linalg
+
+import ase
+import ase.io
+import ase.visualize
+import ase.neighborlist
+
+import os
+import shutil
+
+
+def convert_neighbor_list(nl):
+    new = {}
+    n_vert = np.max(nl)+1
+    
+    for i_v in range(n_vert):
+        new[i_v] = []
+    
+    for i_v, j_v in zip(nl[0], nl[1]):
+        
+        new[i_v].append(j_v)
+        
+    return new
+
+def find_cycles(i_vert, cnl, max_length, cur_path, passed_edges):
+    
+    if len(cur_path)-1 == max_length:
+        return []
+    
+    acc_cycles = []
+    sort_cycles = []
+    
+    neighbs = cnl[i_vert]
+
+    # if we are connected to something that is not the end
+    # then we crossed multiple cycles
+    for n in neighbs:
+        edge = (np.min([i_vert, n]), np.max([i_vert, n]))
+        if edge not in passed_edges:
+
+            if n in cur_path[1:]:
+                # path went too close to itself...
+                return []
+
+    # CHeck if we are at the end
+    for n in neighbs:
+        edge = (np.min([i_vert, n]), np.max([i_vert, n]))
+        if edge not in passed_edges:
+            
+            if n == cur_path[0]:
+                # found cycle
+                return [cur_path]
+    
+    # Continue in all possible directions
+    for n in neighbs:
+        edge = (np.min([i_vert, n]), np.max([i_vert, n]))
+        if edge not in passed_edges:
+
+            cycs = find_cycles(n, cnl, max_length, cur_path + [n], passed_edges + [edge])
+            for cyc in cycs:
+                sorted_cyc = tuple(sorted(cyc))
+                if sorted_cyc not in sort_cycles:
+                    sort_cycles.append(sorted_cyc)
+                    acc_cycles.append(cyc)
+    
+    return acc_cycles
+    
+
+def dumb_cycle_detection(ase_atoms_no_h, max_length):
+    
+    neighbor_list = ase.neighborlist.neighbor_list("ij", ase_atoms_no_h, 2.0)
+    
+    cycles = []
+    sorted_cycles = []
+    n_vert = np.max(neighbor_list)+1
+    
+    cnl = convert_neighbor_list(neighbor_list)
+    
+    for i_vert in range(n_vert):
+        
+        cycs = find_cycles(i_vert, cnl, max_length, [i_vert], [])
+        for cyc in cycs:
+            sorted_cyc = tuple(sorted(cyc))
+            if sorted_cyc not in sorted_cycles:
+                sorted_cycles.append(sorted_cyc)
+                cycles.append(cyc)
+    
+    return cycles
+
+
+def cycle_normal(cycle, h):
+    cycle = np.array(cycle)
+    centroid = np.mean(cycle, axis=0)
+
+    points = cycle - centroid
+    u, s, v = np.linalg.svd(points.T)
+    normal = u[:, -1]
+    normal /= np.linalg.norm(normal)
+    if np.dot(normal, h*np.array([1, 1, 1])) < 0.0:
+        normal *= -1.0
+    return normal
+
+
+def find_cycle_centers_and_normals(ase_atoms_no_h, cycles, h=0.0):
+    """
+    positive h means projection to z axis is positive and vice-versa
+    """
+    if h == 0.0:
+        h = 1.0
+    normals = []
+    centers = []
+    for cyc in cycles:
+        cyc_p = []
+        for i_at in cyc:
+            cyc_p.append(ase_atoms_no_h[i_at].position)
+        normals.append(cycle_normal(cyc_p, h))
+        centers.append(np.mean(cyc_p, axis=0))
+    return np.array(centers), np.array(normals)
diff --git a/aiida_nanotech_empa/utils/nmr.py b/aiida_nanotech_empa/utils/nmr.py
new file mode 100644
index 0000000..4145659
--- /dev/null
+++ b/aiida_nanotech_empa/utils/nmr.py
@@ -0,0 +1,161 @@
+import numpy as np
+
+import scipy as sp
+import scipy.linalg
+
+import ase
+import ase.io
+import ase.visualize
+import ase.neighborlist
+
+import os
+import shutil
+import cclib
+
+from . import cycle_tools
+
+### -----------------------------------------------------------------
+### SETUP
+
+
+def find_ref_points(ase_atoms_no_h, cycles, h = 0.0):
+    """
+    positive h means projection to z axis is positive and vice-versa
+    """
+
+    centers, normals = cycle_tools.find_cycle_centers_and_normals(ase_atoms_no_h, cycles, h)
+    
+    ase_ref_p = ase.Atoms()
+
+    for i_cyc in range(len(cycles)):
+        if h == 0.0:
+            ase_ref_p.append(ase.Atom('X', centers[i_cyc]))
+        else:
+            pos = centers[i_cyc] + np.abs(h)*normals[i_cyc]
+            ase_ref_p.append(ase.Atom('X', pos))
+            
+    return ase_ref_p
+
+def dist(p1, p2):
+    if isinstance(p1, ase.Atom):
+        p1 = p1.position
+    if isinstance(p2, ase.Atom):
+        p2 = p2.position
+    return np.linalg.norm(p2-p1)
+
+def interp_pts(p1, p2, dx):
+    vec = p2 - p1
+    dist = np.sqrt(np.sum(vec**2))
+    num_p = int(np.round(dist/dx))
+    dx_real = dist/num_p
+    dvec = vec/dist*dx_real
+    
+    points = np.outer(np.arange(0, num_p), dvec) + p1
+    return points
+
+def build_path(ref_pts, dx=0.1):
+
+    point_arr = None
+
+    for i_rp in range(len(ref_pts)-1):
+
+        pt1 = ref_pts[i_rp].position
+        pt2 = ref_pts[i_rp+1].position
+
+        points = interp_pts(pt1, pt2, dx)
+
+        if i_rp == len(ref_pts)-2:
+            points = np.concatenate([points, [pt2]], axis=0)
+
+        if point_arr is None:
+            point_arr = points
+        else:
+            point_arr = np.concatenate([point_arr, points], axis=0)
+    
+    ase_arr = ase.Atoms("X%d"%len(point_arr), point_arr)
+    return ase_arr
+
+
+### -----------------------------------------------------------------
+### PROCESS
+
+def load_nics_gaussian(nics_path):
+
+    sigma = []
+
+    def extract_values(line):
+        parts = line.split()
+        return np.array([parts[1], parts[3], parts[5]], dtype=float)
+
+    with open(nics_path) as file:
+        lines = file.readlines()
+        i_line = 0
+        while i_line < len(lines):
+            if "Bq   Isotropic" in lines[i_line]:
+                s = np.zeros((3, 3))
+                s[0] = extract_values(lines[i_line + 1])
+                s[1] = extract_values(lines[i_line + 2])
+                s[2] = extract_values(lines[i_line + 3])
+                sigma.append(s)
+                i_line += 4
+            else:
+                i_line += 1
+
+    return np.array(sigma)
+
+
+def is_number(x):
+    try:
+        float(x)
+        return True
+    except ValueError:
+        return False
+    
+def parse_nmr_cmo_matrix(log_file_str, property_dict):
+    
+    lines = log_file_str.splitlines()
+    
+    # build the object
+    n_atom = property_dict['natom']
+    n_occupied_mo = property_dict['homos'][0] + 1
+
+    nmr_cmo_matrix = np.zeros((n_atom, n_occupied_mo, 3, 3))
+
+    i_line = 0
+    while i_line < len(lines):
+
+        # --------------------------------------------------------------------------
+        # Full Cartesian NMR shielding tensor (ppm) for atom  C(  1):
+        # Canonical MO contributions
+        #
+        #   MO        XX      XY      XZ      YX      YY      YZ      ZX      ZY      ZZ
+        # ===============================================================================
+        #    1.      0.01    0.00    0.00    0.00    0.00    0.00    0.00    0.00   -0.16
+        #    2.      0.01    0.00    0.00    0.00    0.00    0.00    0.00    0.00   -0.13
+        #    3.      0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.24
+        #    4.      0.00    0.00    0.00    0.00   -0.03    0.00    0.00    0.00    0.27
+        #    ...
+
+        if "Full Cartesian NMR shielding tensor (ppm) for atom" in lines[i_line]:
+
+            i_atom = int(lines[i_line].replace('(', ')').split(')')[-2]) - 1
+
+            i_line += 1
+
+            if 'Canonical MO contributions' in lines[i_line]:
+
+                for _i in range(2000):
+                    i_line += 1
+                    if 'Total' in lines[i_line]:
+                        break
+                    split = lines[i_line].split()
+
+                    if len(split) == 10 and is_number(split[-1]):
+                        i_mo = int(split[0][:-1]) - 1
+
+                        arr = np.array([float(x) for x in split[1:]])                    
+                        nmr_cmo_matrix[i_atom, i_mo, :, :] = arr.reshape(3, 3)
+
+        i_line += 1
+        
+    return nmr_cmo_matrix
diff --git a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
new file mode 100644
index 0000000..14020c8
--- /dev/null
+++ b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
@@ -0,0 +1,85 @@
+# pylint: disable=invalid-name
+"""Run simple DFT calculation"""
+
+
+import sys
+
+import ase.io
+import click
+from aiida.common import NotExistent
+from aiida.engine import run_get_node
+from aiida.orm import Code, Dict, StructureData
+from aiida.plugins import CalculationFactory
+
+GaussianCalculation = CalculationFactory("gaussian")
+
+
+def example_dft(gaussian_code):
+    """Run a simple gaussian optimization"""
+
+    # structure
+    structure = StructureData(ase=ase.io.read("./ch4.xyz"))
+
+    num_cores = 1
+    memory_mb = 300
+
+    # Main parameters: geometry optimization
+    parameters = Dict(
+        {
+            "link0_parameters": {
+                "%chk": "aiida.chk",
+                "%mem": "%dMB" % memory_mb,
+                "%nprocshared": num_cores,
+            },
+            "functional": "BLYP",
+            "basis_set": "6-31g",
+            "charge": 0,
+            "multiplicity": 1,
+            "dieze_tag": "#P",
+            "route_parameters": {
+                "scf": {
+                    "cdiis": None,
+                },
+                "nosymm": None,
+                "opt": None,
+            },
+        }
+    )
+
+    # Construct process builder
+
+    builder = GaussianCalculation.get_builder()
+
+    builder.structure = structure
+    builder.parameters = parameters
+    builder.code = gaussian_code
+
+    builder.metadata.options.resources = {
+        "num_machines": 1,
+        "tot_num_mpiprocs": num_cores,
+    }
+
+    # Should ask for extra +25% extra memory
+    builder.metadata.options.max_memory_kb = int(1.25 * memory_mb) * 1024
+    builder.metadata.options.max_wallclock_seconds = 5 * 60
+
+    print("Running calculation...")
+    res, _node = run_get_node(builder)
+
+    print("Final scf energy: %.4f" % res["output_parameters"]["scfenergies"][-1])
+
+
+@click.command("cli")
+@click.argument("codelabel", default="gaussian@localhost")
+def cli(codelabel):
+    """Click interface"""
+    try:
+        code = Code.get_from_string(codelabel)
+    except NotExistent:
+        print(f"The code '{codelabel}' does not exist")
+        sys.exit(1)
+    example_dft(code)
+
+
+if __name__ == "__main__":
+    cli()  # pylint: disable=no-value-for-parameter
diff --git a/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py
new file mode 100644
index 0000000..06fe78e
--- /dev/null
+++ b/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py
@@ -0,0 +1,266 @@
+import numpy as np
+from aiida import engine, orm, plugins
+
+from ...utils import common_utils,nmr,cycle_tools
+#from .delta_scf_workchain import GaussianDeltaScfWorkChain
+#from .natorb_workchain import GaussianNatOrbWorkChain
+#from .relax_workchain import GaussianRelaxWorkChain
+f#rom .scf_workchain import GaussianScfWorkChain
+#from aiida.plugins import CalculationFactory
+
+#GaussianBaseWorkChain = plugins.WorkflowFactory("gaussian.base")
+#GaussianCubesWorkChain = plugins.WorkflowFactory("gaussian.cubes")
+GaussianCalculation = plugins.CalculationFactory("gaussian")
+
+@engine.calcfunction
+def nics_structure(structure=None, h=1.0):
+    
+    ase_geo = structure.get_ase()
+    # orient the geometry
+    pos = ase_geo.positions
+    extents = np.array([
+        np.max(pos[:,0]) - np.min(pos[:,0]),
+        np.max(pos[:,1]) - np.min(pos[:,1]),
+        np.max(pos[:,2]) - np.min(pos[:,2]),
+    ])
+    inds = np.argsort(-extents)
+    ase_atoms = ase.Atoms(numbers=ase_geo.numbers, positions=pos[:, inds])
+    ase_atoms_no_h = ase.Atoms([a for a in ase_atoms if a.symbol != 'H'])        
+    cycles = cycle_tools.dumb_cycle_detection(ase_atoms_no_h, 8)
+    ase_geom = ase_atoms + nmr.find_ref_points(ase_atoms_no_h, cycles, h.value) 
+    return orm.StructureData(ase=ase_geom)   
+
+class GaussianNicsWorkChain(engine.WorkChain):
+    @classmethod
+    def define(cls, spec):
+        super().define(spec)
+
+        spec.input("gaussian_code", valid_type=orm.Code)
+
+        spec.input(
+            "structure",
+            valid_type=orm.StructureData,
+            required=True,
+            help="input geometry",
+        )
+        spec.input(
+            "functional", valid_type=orm.Str, required=False,default=lambda: orm.Str("b3lyp"), help="xc functional"
+        )
+        spec.input(
+            "basis_set", valid_type=orm.Str, required=False, default=lambda: orm.Str("6-311+G(d,p)"), help="basis_set for opt"
+        )
+        spec.input(
+            "multiplicity",
+            valid_type=orm.Int,
+            required=True,
+            help="spin multiplicity",
+        )
+
+        spec.outline(
+            cls.setup
+            cls.submit_nics,
+            cls.finalize,
+        )
+
+        spec.outputs.dynamic = True
+
+        spec.exit_code(
+            390,
+            "ERROR_TERMINATION",
+            message="One or more steps of the work chain failed.",
+        )
+
+    def setup(self):
+        self.ctx.structure_wtih_centers = nics_structure(self.inputs.structure, h=1.0)
+        self.ctx.parameters = {
+                "link0_parameters": {
+                    "%chk": "aiida.chk",
+                    "%mem": "4096MB",
+                    "%nprocshared": 4,
+                },
+                "functional": "BLYP",
+                "basis_set": "6-311+G(d,p)",
+                "multiplicity": self.inputs.multiplicity.value,
+                "dieze_tag": "#P",
+                "route_parameters": {
+                    "scf": {
+                        "maxcycle":2048,
+                        "conver":8,
+                    },
+                    "int": "superfine",
+                    "chpf":{"conver":8}
+                    "nmr": None,
+                },
+            }
+
+    def submit_nics(self):
+        # Multiplicity 0 means RKS calculation.
+        label = f"m{self.inputs.multiplicity.value}_opt"
+
+        # Main parameters: geometry optimization
+        parameters = Dict(
+            {
+                "link0_parameters": {
+                    "%chk": "aiida.chk",
+                    "%mem": "4096MB",
+                    "%nprocshared": 4,
+                },
+                "functional": "BLYP",
+                "basis_set": "6-311+G(d,p)",
+                "multiplicity": self.inputs.multiplicity.value,
+                "dieze_tag": "#P",
+                "route_parameters": {
+                    "scf": {
+                        "maxcycle":2048,
+                        "conver":8,
+                    },
+                    "int": "superfine",
+                    "chpf":{"conver":8}
+                    "nmr": None,
+                },
+            }
+        )
+
+        # Construct process builder
+
+        builder = GaussianCalculation.get_builder()
+
+        builder.structure = structure
+        builder.parameters = parameters
+        builder.code = gaussian_code
+
+        builder.metadata.options.resources = {
+            "num_machines": 1,
+            "tot_num_mpiprocs": num_cores,
+        }
+
+        # Should ask for extra +25% extra memory
+        builder.metadata.options.max_memory_kb = int(1.25 * memory_mb) * 1024
+        builder.metadata.options.max_wallclock_seconds = 5 * 60
+
+        submitted_node = self.submit(builder)
+        submitted_node.description = label
+        self.to_context(**{label: submitted_node})
+
+    def inspect_opt(self):
+
+        label = f"m{self.inputs.multiplicity.value}_opt"
+
+        # check if everything finished nicely
+        if not common_utils.check_if_calc_ok(self, self.ctx[label]):
+            return self.exit_codes.ERROR_TERMINATION
+
+        opt_energy = self.ctx[label].outputs.scf_energy_ev
+        self.out("opt_energy", opt_energy)
+        self.out("opt_structure", self.ctx[label].outputs.output_structure)
+        self.out(
+            "opt_out_params", self.ctx[label].outputs.scf_output_parameters
+        )
+
+        return engine.ExitCode(0)
+
+    def submit_next_steps(self):
+
+        self.report("Submitting NMR cubes")
+
+        parameters = orm.Dict(
+                    {
+                        "link0_parameters": self.ctx.link0.copy(),
+                        "dieze_tag": "#P",
+                        "functional": self.ctx.functional,
+                        "basis_set": self.inputs.basis_set.value,
+                        "charge": 0,
+                        "multiplicity": self.ctx.mult,
+                        "route_parameters": {
+                            "scf": {"maxcycle": 140},
+                            "opt": None,
+                        },
+                    }
+                )
+
+        for mult in self.inputs.multiplicity_list:
+            label = "nmr"
+
+            if mult == self.ctx.gs_mult:
+                continue
+
+            builder = GaussianScfWorkChain.get_builder()
+            builder.gaussian_code = self.inputs.gaussian_code
+            builder.structure = self.ctx.gs_structure
+            builder.functional = self.inputs.functional
+            builder.empirical_dispersion = self.inputs.empirical_dispersion
+            builder.basis_set = 
+            builder.multiplicity = self.inputs.multiplicity
+
+            #if "options" in self.inputs:
+            #    builder.options = self.inputs.options
+
+            submitted_node = self.submit(builder)
+            submitted_node.description = label
+            self.to_context(**{label: submitted_node})
+
+    def inspect_next_steps(self):
+        if not common_utils.check_if_calc_ok(self, self.ctx.gs_cubes):
+            return self.exit_codes.ERROR_TERMINATION
+
+        self.out("gs_cube_images", self.ctx.gs_cubes.outputs.cube_image_folder)
+        self.out("gs_cube_planes", self.ctx.gs_cubes.outputs.cube_planes_array)
+
+        if not common_utils.check_if_calc_ok(self, self.ctx.dscf):
+            return self.exit_codes.ERROR_TERMINATION
+
+        self.out("gs_ionization_potential", self.ctx.dscf.outputs.ionization_potential)
+        self.out("gs_electron_affinity", self.ctx.dscf.outputs.electron_affinity)
+
+        for mult in self.inputs.multiplicity_list:
+            label = f"m{mult}_vert"
+
+            if mult == self.ctx.gs_mult:
+                continue
+
+            # Check if everything finished nicely.
+            if not common_utils.check_if_calc_ok(self, self.ctx[label]):
+                return self.exit_codes.ERROR_TERMINATION
+
+            vert_energy = self.ctx[label].outputs.energy_ev
+            self.out(f"m{mult}_vert_energy", vert_energy)
+            self.out(
+                f"m{mult}_vert_out_params", self.ctx[label].outputs.output_parameters
+            )
+            self.out(
+                f"m{mult}_vert_cube_images", self.ctx[label].outputs.cube_image_folder
+            )
+            self.out(
+                f"m{mult}_vert_cube_planes", self.ctx[label].outputs.cube_planes_array
+            )
+
+        return engine.ExitCode(0)
+
+    def is_gs_oss(self):
+        """Is ground state an open-shell singlet?"""
+        return self.ctx.gs_mult == 1
+
+    def submit_nat_orb(self):
+        self.report("Submitting natural pop analysis")
+
+        builder = GaussianNatOrbWorkChain.get_builder()
+        builder.gaussian_code = self.inputs.gaussian_code
+        builder.parent_calc_folder = self.ctx.gs_scf_remote_folder
+        builder.parent_calc_params = self.ctx.gs_out_params
+        if "options" in self.inputs:
+            builder.options = self.inputs.options
+
+        submitted_node = self.submit(builder)
+        submitted_node.description = "natural orbitals pop"
+        self.to_context(natorb=submitted_node)
+
+    def inspect_nat_orb(self):
+        if not common_utils.check_if_calc_ok(self, self.ctx.natorb):
+            return self.exit_codes.ERROR_TERMINATION
+
+        self.out("gs_natorb_params", self.ctx.natorb.outputs.natorb_proc_parameters)
+
+        return engine.ExitCode(0)
+
+    def finalize(self):
+        self.report("Finalizing...")
diff --git a/examples/workflows/aa b/examples/workflows/aa
deleted file mode 100755
index f9b714b..0000000
--- a/examples/workflows/aa
+++ /dev/null
@@ -1,431 +0,0 @@
-verdi node delete 337
-verdi node delete 344
-verdi node delete 345
-verdi node delete 353
-verdi node delete 354
-verdi node delete 362
-verdi node delete 369
-verdi node delete 370
-verdi node delete 378
-verdi node delete 379
-verdi node delete 387
-verdi node delete 394
-verdi node delete 395
-verdi node delete 403
-verdi node delete 404
-verdi node delete 414
-verdi node delete 422
-verdi node delete 423
-verdi node delete 434
-verdi node delete 442
-verdi node delete 443
-verdi node delete 498
-verdi node delete 501
-verdi node delete 508
-verdi node delete 509
-verdi node delete 517
-verdi node delete 518
-verdi node delete 523
-verdi node delete 537
-verdi node delete 545
-verdi node delete 550
-verdi node delete 555
-verdi node delete 556
-verdi node delete 557
-verdi node delete 558
-verdi node delete 571
-verdi node delete 575
-verdi node delete 576
-verdi node delete 577
-verdi node delete 588
-verdi node delete 593
-verdi node delete 596
-verdi node delete 603
-verdi node delete 610
-verdi node delete 611
-verdi node delete 612
-verdi node delete 621
-verdi node delete 622
-verdi node delete 630
-verdi node delete 631
-verdi node delete 638
-verdi node delete 655
-verdi node delete 663
-verdi node delete 668
-verdi node delete 673
-verdi node delete 674
-verdi node delete 675
-verdi node delete 676
-verdi node delete 689
-verdi node delete 693
-verdi node delete 694
-verdi node delete 695
-verdi node delete 710
-verdi node delete 715
-verdi node delete 718
-verdi node delete 725
-verdi node delete 732
-verdi node delete 733
-verdi node delete 734
-verdi node delete 743
-verdi node delete 744
-verdi node delete 752
-verdi node delete 753
-verdi node delete 760
-verdi node delete 975
-verdi node delete 981
-verdi node delete 988
-verdi node delete 989
-verdi node delete 997
-verdi node delete 998
-verdi node delete 1002
-verdi node delete 1003
-verdi node delete 1012
-verdi node delete 1018
-verdi node delete 1025
-verdi node delete 1026
-verdi node delete 1034
-verdi node delete 1035
-verdi node delete 1039
-verdi node delete 1040
-verdi node delete 1051
-verdi node delete 1056
-verdi node delete 1063
-verdi node delete 1064
-verdi node delete 1072
-verdi node delete 1073
-verdi node delete 1077
-verdi node delete 1080
-verdi node delete 1086
-verdi node delete 1090
-verdi node delete 1097
-verdi node delete 1098
-verdi node delete 1106
-verdi node delete 1107
-verdi node delete 1112
-verdi node delete 1118
-verdi node delete 1121
-verdi node delete 1128
-verdi node delete 1129
-verdi node delete 1137
-verdi node delete 1138
-verdi node delete 1143
-verdi node delete 1153
-verdi node delete 1161
-verdi node delete 1162
-verdi node delete 1174
-verdi node delete 1182
-verdi node delete 1183
-verdi node delete 1198
-verdi node delete 1206
-verdi node delete 1207
-verdi node delete 1222
-verdi node delete 1230
-verdi node delete 1231
-verdi node delete 1243
-verdi node delete 1251
-verdi node delete 1252
-verdi node delete 1263
-verdi node delete 1269
-verdi node delete 1275
-verdi node delete 1281
-verdi node delete 1287
-verdi node delete 1293
-verdi node delete 1299
-verdi node delete 1305
-verdi node delete 1311
-verdi node delete 1319
-verdi node delete 1320
-verdi node delete 1324
-verdi node delete 1332
-verdi node delete 1333
-verdi node delete 1341
-verdi node delete 1349
-verdi node delete 1350
-verdi node delete 1355
-verdi node delete 1363
-verdi node delete 1364
-verdi node delete 1369
-verdi node delete 1377
-verdi node delete 1378
-verdi node delete 1386
-verdi node delete 1394
-verdi node delete 1395
-verdi node delete 1433
-verdi node delete 1441
-verdi node delete 1442
-verdi node delete 1450
-verdi node delete 1456
-verdi node delete 1464
-verdi node delete 1465
-verdi node delete 1473
-verdi node delete 1481
-verdi node delete 1482
-verdi node delete 1490
-verdi node delete 1498
-verdi node delete 1499
-verdi node delete 1507
-verdi node delete 1515
-verdi node delete 1516
-verdi node delete 1524
-verdi node delete 1532
-verdi node delete 1533
-verdi node delete 1541
-verdi node delete 1549
-verdi node delete 1550
-verdi node delete 1553
-verdi node delete 1559
-verdi node delete 1567
-verdi node delete 1568
-verdi node delete 1571
-verdi node delete 1576
-verdi node delete 1584
-verdi node delete 1585
-verdi node delete 1588
-verdi node delete 1593
-verdi node delete 1601
-verdi node delete 1602
-verdi node delete 1610
-verdi node delete 1618
-verdi node delete 1619
-verdi node delete 1627
-verdi node delete 1635
-verdi node delete 1636
-verdi node delete 1644
-verdi node delete 1652
-verdi node delete 1653
-verdi node delete 1661
-verdi node delete 1667
-verdi node delete 1675
-verdi node delete 1676
-verdi node delete 1684
-verdi node delete 1692
-verdi node delete 1693
-verdi node delete 1701
-verdi node delete 1707
-verdi node delete 1715
-verdi node delete 1716
-verdi node delete 1724
-verdi node delete 1732
-verdi node delete 1733
-verdi node delete 1741
-verdi node delete 1749
-verdi node delete 1750
-verdi node delete 1753
-verdi node delete 1759
-verdi node delete 1767
-verdi node delete 1768
-verdi node delete 1776
-verdi node delete 1784
-verdi node delete 1785
-verdi node delete 1793
-verdi node delete 1801
-verdi node delete 1802
-verdi node delete 1807
-verdi node delete 1815
-verdi node delete 1816
-verdi node delete 1824
-verdi node delete 1832
-verdi node delete 1833
-verdi node delete 1841
-verdi node delete 1849
-verdi node delete 1850
-verdi node delete 1855
-verdi node delete 1863
-verdi node delete 1864
-verdi node delete 1872
-verdi node delete 1880
-verdi node delete 1881
-verdi node delete 1889
-verdi node delete 1897
-verdi node delete 1898
-verdi node delete 1906
-verdi node delete 1914
-verdi node delete 1915
-verdi node delete 1923
-verdi node delete 1931
-verdi node delete 1932
-verdi node delete 1940
-verdi node delete 1948
-verdi node delete 1949
-verdi node delete 1956
-verdi node delete 1961
-verdi node delete 1969
-verdi node delete 1970
-verdi node delete 1973
-verdi node delete 1978
-verdi node delete 1986
-verdi node delete 1987
-verdi node delete 1991
-verdi node delete 1999
-verdi node delete 2000
-verdi node delete 2033
-verdi node delete 2043
-verdi node delete 2050
-verdi node delete 2058
-verdi node delete 2068
-verdi node delete 2075
-verdi node delete 2082
-verdi node delete 2092
-verdi node delete 2099
-verdi node delete 2112
-verdi node delete 2122
-verdi node delete 2129
-verdi node delete 2143
-verdi node delete 2150
-verdi node delete 2163
-verdi node delete 2173
-verdi node delete 2179
-verdi node delete 2186
-verdi node delete 2196
-verdi node delete 2203
-verdi node delete 2216
-verdi node delete 2243
-verdi node delete 2249
-verdi node delete 2254
-verdi node delete 2259
-verdi node delete 2271
-verdi node delete 2279
-verdi node delete 2280
-verdi node delete 2281
-verdi node delete 2298
-verdi node delete 2306
-verdi node delete 2307
-verdi node delete 2308
-verdi node delete 2325
-verdi node delete 2326
-verdi node delete 2339
-verdi node delete 2340
-verdi node delete 2349
-verdi node delete 2354
-verdi node delete 2366
-verdi node delete 2374
-verdi node delete 2375
-verdi node delete 2376
-verdi node delete 2393
-verdi node delete 2401
-verdi node delete 2402
-verdi node delete 2403
-verdi node delete 2420
-verdi node delete 2428
-verdi node delete 2429
-verdi node delete 2430
-verdi node delete 2447
-verdi node delete 2455
-verdi node delete 2456
-verdi node delete 2457
-verdi node delete 2474
-verdi node delete 2482
-verdi node delete 2483
-verdi node delete 2484
-verdi node delete 2501
-verdi node delete 2509
-verdi node delete 2510
-verdi node delete 2511
-verdi node delete 2528
-verdi node delete 2529
-verdi node delete 2561
-verdi node delete 2566
-verdi node delete 2574
-verdi node delete 2601
-verdi node delete 2605
-verdi node delete 2610
-verdi node delete 2618
-verdi node delete 2623
-verdi node delete 2650
-verdi node delete 2655
-verdi node delete 2660
-verdi node delete 2665
-verdi node delete 2670
-verdi node delete 2678
-verdi node delete 2683
-verdi node delete 2695
-verdi node delete 2703
-verdi node delete 2704
-verdi node delete 2705
-verdi node delete 2722
-verdi node delete 2730
-verdi node delete 2731
-verdi node delete 2732
-verdi node delete 2749
-verdi node delete 2750
-verdi node delete 2763
-verdi node delete 2764
-verdi node delete 2773
-verdi node delete 2778
-verdi node delete 2790
-verdi node delete 2798
-verdi node delete 2799
-verdi node delete 2800
-verdi node delete 2817
-verdi node delete 2825
-verdi node delete 2826
-verdi node delete 2827
-verdi node delete 2844
-verdi node delete 2852
-verdi node delete 2853
-verdi node delete 2854
-verdi node delete 2871
-verdi node delete 2879
-verdi node delete 2880
-verdi node delete 2881
-verdi node delete 2898
-verdi node delete 2906
-verdi node delete 2907
-verdi node delete 2908
-verdi node delete 2925
-verdi node delete 2933
-verdi node delete 2934
-verdi node delete 2935
-verdi node delete 2952
-verdi node delete 2953
-verdi node delete 2961
-verdi node delete 2966
-verdi node delete 2978
-verdi node delete 2979
-verdi node delete 2990
-verdi node delete 2994
-verdi node delete 2999
-verdi node delete 3011
-verdi node delete 3019
-verdi node delete 3020
-verdi node delete 3021
-verdi node delete 3027
-verdi node delete 3029
-verdi node delete 3032
-verdi node delete 3034
-verdi node delete 3035
-verdi node delete 3037
-verdi node delete 3054
-verdi node delete 3062
-verdi node delete 3063
-verdi node delete 3064
-verdi node delete 3069
-verdi node delete 3071
-verdi node delete 3073
-verdi node delete 3077
-verdi node delete 3079
-verdi node delete 3081
-verdi node delete 3088
-verdi node delete 3090
-verdi node delete 3092
-verdi node delete 3105
-verdi node delete 3106
-verdi node delete 3115
-verdi node delete 3123
-verdi node delete 3124
-verdi node delete 3132
-verdi node delete 3140
-verdi node delete 3141
-verdi node delete 3149
-verdi node delete 3157
-verdi node delete 3158
-verdi node delete 3166
-verdi node delete 3174
-verdi node delete 3175
-verdi node delete 3183
-verdi node delete 3191
-verdi node delete 3192
-verdi node delete 3200
-verdi node delete 3208
-verdi node delete 3209
diff --git a/examples/workflows/benzene.xyz b/examples/workflows/benzene.xyz
new file mode 100644
index 0000000..66fd521
--- /dev/null
+++ b/examples/workflows/benzene.xyz
@@ -0,0 +1,14 @@
+12
+benzene
+C        8.39528138       6.66748365       5.00000291   
+C        8.51925793       8.06102240       5.00000509   
+C        7.12645219       6.07808103       5.00000441   
+C        5.98160072       6.88221755       5.00000464   
+C        6.10557756       8.27575669       5.00000369   
+C        7.37440623       8.86515917       5.00000734   
+H        9.50085998       8.51699940       5.00000031   
+H        9.28097172       6.04538232       5.00000427   
+H        7.03053864       5.00000000       5.00000036   
+H        5.00000000       6.42623909       5.00000439   
+H        5.21988789       8.89785938       5.00000000   
+H        7.47031658       9.94324034       5.00000253   
\ No newline at end of file
diff --git a/examples/workflows/example_gaussian_opt.py b/examples/workflows/example_gaussian_opt.py
new file mode 100644
index 0000000..f7beade
--- /dev/null
+++ b/examples/workflows/example_gaussian_opt.py
@@ -0,0 +1,55 @@
+import os
+
+import ase.io
+import numpy as np
+from aiida.engine import run_get_node
+from aiida.orm import Bool,Int, Str, StructureData, load_code
+from aiida.plugins import WorkflowFactory
+
+import aiida_nanotech_empa.utils.gaussian_wcs_postprocess as pp
+
+GaussianRelaxWorkChain = WorkflowFactory("nanotech_empa.gaussian.relax")
+
+DATA_DIR = os.path.dirname(os.path.abspath(__file__))
+OUTPUT_DIR = os.path.dirname(os.path.abspath(__file__))
+
+
+def _example_gaussian_spin(gaussian_code):#, formchk_code, cubegen_code):
+    ase_geom = ase.io.read(os.path.join(DATA_DIR, "benzene.xyz"))
+    ase_geom.cell = np.diag([10.0, 10.0, 10.0])
+
+    builder = GaussianRelaxWorkChain.get_builder()
+    builder.gaussian_code = gaussian_code
+    #builder.formchk_code = formchk_code
+    #builder.cubegen_code = cubegen_code
+    builder.structure = StructureData(ase=ase_geom)
+    builder.functional = Str("B3LYP")
+    builder.empirical_dispersion = Str("GD3")
+    builder.basis_set = Str("6-311+G(d,p)")
+    #builder.basis_set_scf = Str("STO-3G")
+    builder.multiplicity = Int(0)
+    builder.tight = Bool(True)
+    builder.options = Dict(
+        {
+            "resources": {
+                "tot_num_mpiprocs": 4,
+                "num_machines": 1,
+            },
+            "max_wallclock_seconds": 1 * 60 * 60,
+            "max_memory_kb": 4 * 1024 * 1024,
+        }
+    )    
+
+    _, wc_node = run_get_node(builder)
+
+    assert wc_node.is_finished_ok
+
+    #pp.make_report(wc_node, nb=False, save_image_loc=OUTPUT_DIR)
+
+
+if __name__ == "__main__":
+    _example_gaussian_spin(
+        load_code("gaussian@tigu"),
+        #load_code("formchk@localhost"),
+        #load_code("cubegen@localhost"),
+    )

From 71c1b563472f88c08338fd2d052722eb97d5fba6 Mon Sep 17 00:00:00 2001
From: Carlo Antonio Pignedoli <c.pignedoli@gmail.com>
Date: Wed, 11 Sep 2024 11:45:01 +0000
Subject: [PATCH 2/9] not much done

---
 .../workflows/gaussian/nics_workchain.py      | 61 ++++++++++++++-----
 examples/workflows/example_gaussian_nics.py   | 39 ++++++++++++
 2 files changed, 84 insertions(+), 16 deletions(-)
 create mode 100644 examples/workflows/example_gaussian_nics.py

diff --git a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
index 14020c8..58d5f1e 100644
--- a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
@@ -1,15 +1,40 @@
-# pylint: disable=invalid-name
-"""Run simple DFT calculation"""
-
-
-import sys
-
-import ase.io
-import click
-from aiida.common import NotExistent
-from aiida.engine import run_get_node
-from aiida.orm import Code, Dict, StructureData
-from aiida.plugins import CalculationFactory
+class GaussianRelaxWorkChain(engine.WorkChain):
+    @classmethod
+    def define(cls, spec):
+        super().define(spec)
+
+        spec.input("gaussian_code", valid_type=orm.Code)
+
+        spec.input(
+            "structure",
+            valid_type=orm.StructureData,
+            required=True,
+            help="input geometry",
+        )
+        spec.input(
+            "functional", valid_type=orm.Str, required=True, help="xc functional"
+        )
+
+        spec.input("basis_set", valid_type=orm.Str, required=True, help="basis_set")
+
+        spec.input(
+            "multiplicity",
+            valid_type=orm.Int,
+            required=False,
+            default=lambda: orm.Int(0),
+            help="spin multiplicity; 0 means RKS",
+        )
+
+        spec.input(
+            "wfn_stable_opt",
+            valid_type=orm.Bool,
+            required=False,
+            default=lambda: orm.Bool(False),
+            help="if true, perform wfn stability optimization",
+        )
+
+# Create a methane (CH4) molecule using the ASE molecular database
+ch4 = molecule('CH4')
 
 GaussianCalculation = CalculationFactory("gaussian")
 
@@ -18,7 +43,7 @@ def example_dft(gaussian_code):
     """Run a simple gaussian optimization"""
 
     # structure
-    structure = StructureData(ase=ase.io.read("./ch4.xyz"))
+    structure = StructureData(ase=ch4)
 
     num_cores = 1
     memory_mb = 300
@@ -31,17 +56,21 @@ def example_dft(gaussian_code):
                 "%mem": "%dMB" % memory_mb,
                 "%nprocshared": num_cores,
             },
-            "functional": "BLYP",
-            "basis_set": "6-31g",
+            "functional": "B3LYP",
+            "basis_set": "6-311g(d,p)",
             "charge": 0,
             "multiplicity": 1,
             "dieze_tag": "#P",
             "route_parameters": {
                 "scf": {
+                    "maxcycle":2048,
                     "cdiis": None,
+                    "conver":8
                 },
+                "int":"superfine",
+                "guess":"mix",
                 "nosymm": None,
-                "opt": None,
+                "opt": "tight",
             },
         }
     )
diff --git a/examples/workflows/example_gaussian_nics.py b/examples/workflows/example_gaussian_nics.py
new file mode 100644
index 0000000..541f9ff
--- /dev/null
+++ b/examples/workflows/example_gaussian_nics.py
@@ -0,0 +1,39 @@
+mport os
+
+import ase.io
+import numpy as np
+from aiida.engine import run_get_node
+from aiida.orm import Bool,Int, Str, StructureData, load_code
+from aiida.plugins import WorkflowFactory
+
+import aiida_nanotech_empa.utils.gaussian_wcs_postprocess as pp
+
+GaussianRelaxWorkChain = WorkflowFactory("nanotech_empa.gaussian.relax")
+
+DATA_DIR = os.path.dirname(os.path.abspath(__file__))
+OUTPUT_DIR = os.path.dirname(os.path.abspath(__file__))
+
+
+def _example_gaussian_spin(gaussian_code):#, formchk_code, cubegen_code):
+    ase_geom = ase.io.read(os.path.join(DATA_DIR, "benzene.xyz"))
+    ase_geom.cell = np.diag([10.0, 10.0, 10.0])
+
+    builder = GaussianRelaxWorkChain.get_builder()
+    
+    
+    
+    
+    _, wc_node = run_get_node(builder)
+
+    assert wc_node.is_finished_ok
+
+    #pp.make_report(wc_node, nb=False, save_image_loc=OUTPUT_DIR)
+
+
+if __name__ == "__main__":
+    _example_gaussian_spin(
+        load_code("gaussian@tigu"),
+        #load_code("formchk@localhost"),
+        #load_code("cubegen@localhost"),
+    )
+    
\ No newline at end of file

From e2311b8751e9771f8c51e680a6beeee451380122 Mon Sep 17 00:00:00 2001
From: Carlo Antonio Pignedoli <c.pignedoli@gmail.com>
Date: Thu, 12 Sep 2024 13:58:16 +0000
Subject: [PATCH 3/9] working

---
 .../workflows/gaussian/__init__.py            |   2 +
 .../workflows/gaussian/nics_workchain.py      | 219 +++++++++++-------
 .../workflows/gaussian/relax_workchain.py     |  38 ++-
 .../workflows/gaussian/scf_workchain.py       |  52 ++++-
 examples/workflows/example_gaussian_nics.py   |  59 +++--
 pyproject.toml                                |   1 +
 6 files changed, 263 insertions(+), 108 deletions(-)

diff --git a/aiida_nanotech_empa/workflows/gaussian/__init__.py b/aiida_nanotech_empa/workflows/gaussian/__init__.py
index d08b655..0b5d156 100644
--- a/aiida_nanotech_empa/workflows/gaussian/__init__.py
+++ b/aiida_nanotech_empa/workflows/gaussian/__init__.py
@@ -7,6 +7,7 @@
 from .relax_workchain import GaussianRelaxWorkChain
 from .scf_workchain import GaussianScfWorkChain
 from .spin_workchain import GaussianSpinWorkChain
+from .nics_workchain import GaussianNicsWorkChain
 
 __all__ = (
     "GaussianScfWorkChain",
@@ -18,4 +19,5 @@
     "GaussianConstrOptChainWorkChain",
     "GaussianCasscfWorkChain",
     "GaussianCasscfSeriesWorkChain",
+    "GaussianNicsWorkChain",
 )
diff --git a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
index 58d5f1e..966a42f 100644
--- a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
@@ -1,4 +1,30 @@
-class GaussianRelaxWorkChain(engine.WorkChain):
+import numpy as np
+import ase
+from aiida import engine, orm
+from .relax_workchain import GaussianRelaxWorkChain
+from .scf_workchain import GaussianScfWorkChain
+from ...utils import cycle_tools, nmr
+from ...utils import common_utils
+
+@engine.calcfunction
+def prepare_nmr_structure(structure,height):
+    ase_geo = structure.get_ase()
+    pos = ase_geo.positions
+    extents = np.array([
+        np.max(pos[:,0]) - np.min(pos[:,0]),
+        np.max(pos[:,1]) - np.min(pos[:,1]),
+        np.max(pos[:,2]) - np.min(pos[:,2]),
+    ])
+    inds = np.argsort(-extents)
+    new_pos = pos[:, inds]
+    ase_atoms = ase.Atoms(numbers=ase_geo.numbers, positions=new_pos)
+    ase_atoms_no_h = ase.Atoms([a for a in ase_atoms if a.symbol != 'H'])
+    cycles = cycle_tools.dumb_cycle_detection(ase_atoms_no_h, 8)
+    ref_p = nmr.find_ref_points(ase_atoms_no_h, cycles, height.value)
+    new_ase = ase_atoms + ref_p
+    return orm.StructureData(ase=new_ase)
+
+class GaussianNicsWorkChain(engine.WorkChain):
     @classmethod
     def define(cls, spec):
         super().define(spec)
@@ -11,6 +37,18 @@ def define(cls, spec):
             required=True,
             help="input geometry",
         )
+        spec.input(
+            "opt",
+            valid_type=orm.Bool,
+            required=False,
+            help="False do not optimize structure",
+        )
+        spec.input(
+            "height",
+            valid_type=orm.Float,
+            required=False,
+            help="Height of NICS centers",
+        )
         spec.input(
             "functional", valid_type=orm.Str, required=True, help="xc functional"
         )
@@ -24,7 +62,6 @@ def define(cls, spec):
             default=lambda: orm.Int(0),
             help="spin multiplicity; 0 means RKS",
         )
-
         spec.input(
             "wfn_stable_opt",
             valid_type=orm.Bool,
@@ -32,83 +69,105 @@ def define(cls, spec):
             default=lambda: orm.Bool(False),
             help="if true, perform wfn stability optimization",
         )
+        spec.input(
+            "empirical_dispersion",
+            valid_type=orm.Str,
+            required=False,
+            default=lambda: orm.Str(""),
+            help=("Include empirical dispersion corrections" '(e.g. "GD3", "GD3BJ")'),
+        )
+        spec.input(
+            "options",
+            valid_type=orm.Dict,
+            required=False,
+            help="Use custom metadata.options instead of the automatic ones.",
+        )        
+        
+        spec.outline(
+            cls.setup,
+            engine.if_(cls.should_submit_opt)(cls.submit_opt,cls.inspect_opt),            
+            cls.submit_nics,
+            cls.inspect_nics,
+            cls.finalize,
+        )
 
-# Create a methane (CH4) molecule using the ASE molecular database
-ch4 = molecule('CH4')
-
-GaussianCalculation = CalculationFactory("gaussian")
-
-
-def example_dft(gaussian_code):
-    """Run a simple gaussian optimization"""
-
-    # structure
-    structure = StructureData(ase=ch4)
-
-    num_cores = 1
-    memory_mb = 300
-
-    # Main parameters: geometry optimization
-    parameters = Dict(
-        {
-            "link0_parameters": {
-                "%chk": "aiida.chk",
-                "%mem": "%dMB" % memory_mb,
-                "%nprocshared": num_cores,
-            },
-            "functional": "B3LYP",
-            "basis_set": "6-311g(d,p)",
-            "charge": 0,
-            "multiplicity": 1,
-            "dieze_tag": "#P",
-            "route_parameters": {
-                "scf": {
-                    "maxcycle":2048,
-                    "cdiis": None,
-                    "conver":8
-                },
-                "int":"superfine",
-                "guess":"mix",
-                "nosymm": None,
-                "opt": "tight",
-            },
-        }
-    )
-
-    # Construct process builder
-
-    builder = GaussianCalculation.get_builder()
-
-    builder.structure = structure
-    builder.parameters = parameters
-    builder.code = gaussian_code
-
-    builder.metadata.options.resources = {
-        "num_machines": 1,
-        "tot_num_mpiprocs": num_cores,
-    }
-
-    # Should ask for extra +25% extra memory
-    builder.metadata.options.max_memory_kb = int(1.25 * memory_mb) * 1024
-    builder.metadata.options.max_wallclock_seconds = 5 * 60
-
-    print("Running calculation...")
-    res, _node = run_get_node(builder)
-
-    print("Final scf energy: %.4f" % res["output_parameters"]["scfenergies"][-1])
-
-
-@click.command("cli")
-@click.argument("codelabel", default="gaussian@localhost")
-def cli(codelabel):
-    """Click interface"""
-    try:
-        code = Code.get_from_string(codelabel)
-    except NotExistent:
-        print(f"The code '{codelabel}' does not exist")
-        sys.exit(1)
-    example_dft(code)
-
+        spec.outputs.dynamic = True
 
-if __name__ == "__main__":
-    cli()  # pylint: disable=no-value-for-parameter
+        spec.exit_code(
+            390,
+            "ERROR_TERMINATION",
+            message="One or more steps of the work chain failed.",
+        )
+    def setup(self):
+        self.report("Setting up...")
+        self.ctx.nmr_structure = self.inputs.structure
+        self.ctx_should_opt = getattr(self.inputs, "opt", True)
+        
+    def should_submit_opt(self):
+        return self.ctx_should_opt
+    
+    def submit_opt(self):
+        self.report("Submitting optimization...")
+        label = "geo_opt"
+        builder = GaussianRelaxWorkChain.get_builder()
+        builder.gaussian_code = self.inputs.gaussian_code
+        builder.structure = self.inputs.structure
+        builder.functional = self.inputs.functional
+        builder.empirical_dispersion = self.inputs.empirical_dispersion
+        builder.basis_set = self.inputs.basis_set
+        builder.multiplicity = self.inputs.multiplicity
+        builder.tight = orm.Bool(True)
+        builder.int = orm.Str("superfine")
+        builder.cdiis = orm.Bool(True)
+        builder.maxcycle = orm.Int(2048)
+        builder.conver = orm.Int(8)
+        if "options" in self.inputs:
+            builder.options = self.inputs.options
+
+        submitted_node = self.submit(builder)
+        submitted_node.description = label
+        self.to_context(**{label: submitted_node})
+        
+    def inspect_opt(self):
+        self.report("Inspecting optimization...")
+        label = "geo_opt"
+        # check if everything finished nicely
+        if not common_utils.check_if_calc_ok(self, self.ctx[label]):
+            return self.exit_codes.ERROR_TERMINATION        
+        self.ctx.nmr_structure = self.ctx[label].outputs.output_structure
+        return engine.ExitCode(0)
+    
+    def submit_nics(self):
+        self.report("Submitting NICS calculation...")
+        label = "nics"
+        builder = GaussianScfWorkChain.get_builder()
+        height = getattr(self.inputs, "height", 1.0)
+        builder.gaussian_code = self.inputs.gaussian_code
+        builder.structure = prepare_nmr_structure(self.ctx.nmr_structure,orm.Float(height))
+        builder.functional = self.inputs.functional
+        builder.basis_set = self.inputs.basis_set
+        builder.multiplicity = self.inputs.multiplicity
+        builder.int = orm.Str("superfine")
+        builder.cdiis = orm.Bool(True)
+        builder.nmr = orm.Bool(True)
+        builder.maxcycle = orm.Int(2048)
+        builder.conver = orm.Int(8)        
+        if "options" in self.inputs:
+            builder.options = self.inputs.options
+
+        submitted_node = self.submit(builder)
+        submitted_node.description = label
+        self.to_context(**{label: submitted_node})
+        
+    def inspect_nics(self):
+        self.report("Inspecting nics...")
+        label = "nics"
+        # check if everything finished nicely
+        if not common_utils.check_if_calc_ok(self, self.ctx[label]):
+            return self.exit_codes.ERROR_TERMINATION
+        
+        self.out("output_parameters", self.ctx[label].outputs.output_parameters)
+        return engine.ExitCode(0)   
+         
+    def finalize(self):
+        self.report("Finalizing...")
\ No newline at end of file
diff --git a/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py b/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py
index 9c03573..1ea2d76 100644
--- a/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py
@@ -49,7 +49,30 @@ def define(cls, spec):
             default=lambda: orm.Bool(False),
             help="Use tight optimization criteria.",
         )
-
+        spec.input(
+            "maxcycle",
+            valid_type=orm.Int,
+            required=False,
+            help="the maximum number of scf cycles",
+        )
+        spec.input(
+            "cdiis",
+            valid_type=orm.Bool,
+            required=False,
+            help="Conjugate Direct Inversion in the Iterative Subspace",
+        )        
+        spec.input(
+            "conver",
+            valid_type=orm.Int,
+            required=False,
+            help="the scf convergece threshold",
+        )        
+        spec.input(
+            "int",
+            valid_type=orm.Str,
+            required=False,
+            help="the integral grid",
+        )
         spec.input(
             "freq",
             valid_type=orm.Bool,
@@ -290,7 +313,18 @@ def optimization(self):
 
         if len(opt_dict) != 0:
             parameters["route_parameters"]["opt"] = opt_dict
-
+        conver = getattr(self.inputs,"conver",None)
+        maxcycle = getattr(self.inputs,"maxcycle",None)
+        grid = getattr(self.inputs,"int",None)
+        cdiis = getattr(self.inputs,"cdiis",None)
+        if grid:
+            parameters["route_parameters"]["int"] = grid
+        if maxcycle:
+            parameters["route_parameters"]["scf"]["maxcycle"] = maxcycle
+        if cdiis:
+            parameters["route_parameters"]["scf"]["cdiis"] = None
+        if conver:
+            parameters["route_parameters"]["scf"]["conver"] = conver
         builder.gaussian.parameters = parameters
         builder.gaussian.structure = self.inputs.structure
         builder.gaussian.code = self.inputs.gaussian_code
diff --git a/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py b/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py
index f38857a..db8644e 100644
--- a/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py
@@ -60,6 +60,36 @@ def define(cls, spec):
             required=False,
             help="the folder of a completed gaussian calculation",
         )
+        spec.input(
+            "maxcycle",
+            valid_type=orm.Int,
+            required=False,
+            help="the maximum number of scf cycles",
+        )
+        spec.input(
+            "conver",
+            valid_type=orm.Int,
+            required=False,
+            help="the scf convergece threshold",
+        )  
+        spec.input(
+            "cdiis",
+            valid_type=orm.Bool,
+            required=False,
+            help="Conjugate Direct Inversion in the Iterative Subspace",
+        )               
+        spec.input(
+            "int",
+            valid_type=orm.Str,
+            required=False,
+            help="the integral grid",
+        )
+        spec.input(
+            "nmr",
+            valid_type=orm.Bool,
+            required=False,
+            help="nmr calculation",
+        )
 
         # Inputs for cubes generation.
         spec.input("formchk_code", valid_type=orm.Code, required=False)
@@ -167,7 +197,7 @@ def setup(self):
 
         # Use default convergence criterion at the start
         # but switch to conver=7 in case of failure.
-        self.ctx.conver = None
+        self.ctx.conver = getattr(self.inputs, "conver", None)
         self.ctx.scf_label = "scf"
 
         return engine.ExitCode(0)
@@ -224,10 +254,28 @@ def scf(self):
                 "charge": 0,
                 "multiplicity": self.ctx.mult,
                 "route_parameters": {
-                    "scf": {"maxcycle": 140},
+                    "scf": {"maxcycle": 140 },
                 },
             }
         )
+        nmr = getattr(self.inputs,"nmr",None)
+        conver = getattr(self.inputs,"conver",None)
+        maxcycle = getattr(self.inputs,"maxcycle",None)
+        grid = getattr(self.inputs,"int",None)
+        cdiis = getattr(self.inputs,"cdiis",None)
+        if grid:
+            parameters["route_parameters"]["int"] = grid
+        if maxcycle:
+            parameters["route_parameters"]["scf"]["maxcycle"] = maxcycle
+        if cdiis:
+            parameters["route_parameters"]["scf"]["cdiis"] = None
+        if conver:
+            parameters["route_parameters"]["scf"]["conver"] = conver
+        if nmr:
+            if conver is None:
+                conver = 8
+            parameters["route_parameters"]["nmr"] = None
+            parameters["route_parameters"]["cphf"] = {"conver": conver}
         if self.inputs.wfn_stable_opt:
             parameters["route_parameters"]["stable"] = "opt"
         else:
diff --git a/examples/workflows/example_gaussian_nics.py b/examples/workflows/example_gaussian_nics.py
index 541f9ff..9c43f4b 100644
--- a/examples/workflows/example_gaussian_nics.py
+++ b/examples/workflows/example_gaussian_nics.py
@@ -1,39 +1,50 @@
-mport os
-
-import ase.io
+import os
+import click
+from ase.build import molecule
 import numpy as np
-from aiida.engine import run_get_node
-from aiida.orm import Bool,Int, Str, StructureData, load_code
+from aiida import engine, orm
 from aiida.plugins import WorkflowFactory
 
 import aiida_nanotech_empa.utils.gaussian_wcs_postprocess as pp
 
-GaussianRelaxWorkChain = WorkflowFactory("nanotech_empa.gaussian.relax")
-
-DATA_DIR = os.path.dirname(os.path.abspath(__file__))
-OUTPUT_DIR = os.path.dirname(os.path.abspath(__file__))
-
+GaussianNicsWorkChain = WorkflowFactory("nanotech_empa.gaussian.nics")
 
-def _example_gaussian_spin(gaussian_code):#, formchk_code, cubegen_code):
-    ase_geom = ase.io.read(os.path.join(DATA_DIR, "benzene.xyz"))
+def _example_gaussian_nics(gaussian_code,opt):
+    ase_geom = molecule('C6H6')
     ase_geom.cell = np.diag([10.0, 10.0, 10.0])
 
-    builder = GaussianRelaxWorkChain.get_builder()
+    builder = GaussianNicsWorkChain.get_builder()
+    builder.gaussian_code = gaussian_code
+    builder.structure = orm.StructureData(ase=ase_geom)
+    builder.functional = orm.Str("B3LYP")
+    builder.basis_set = orm.Str("6-311G(d,p)")
+    builder.opt = orm.Bool(opt)
+    builder.options = orm.Dict(
+        {
+            "resources": {
+                "tot_num_mpiprocs": 4,
+                "num_machines": 1,
+            },
+            "max_wallclock_seconds": 1 * 60 * 60,
+            "max_memory_kb": 8 * 1024 * 1024, #GB
+        }
+    )    
     
+    _, wc_node = engine.run_get_node(builder)
     
-    
-    
-    _, wc_node = run_get_node(builder)
-
     assert wc_node.is_finished_ok
+    return wc_node.pk
 
-    #pp.make_report(wc_node, nb=False, save_image_loc=OUTPUT_DIR)
+@click.command("cli")
+@click.argument("gaussian_code", default="gaussian@localhost")
+def run_nics(gaussian_code):
+    #print("#### Running Gaussian NICS WorkChain with geo_opt ####")
+    #uuid = _example_gaussian_nics(orm.load_code(gaussian_code),True)
+    #print(f"WorkChain completed uuid: {uuid}")
+    print("#### Running Gaussian NICS WorkChain without geo_opt ####")
+    uuid = _example_gaussian_nics(orm.load_code(gaussian_code),False)
+    print(f"WorkChain completed uuid: {uuid}")
 
 
 if __name__ == "__main__":
-    _example_gaussian_spin(
-        load_code("gaussian@tigu"),
-        #load_code("formchk@localhost"),
-        #load_code("cubegen@localhost"),
-    )
-    
\ No newline at end of file
+    run_nics()
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
index 17a63dd..19116f1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -70,6 +70,7 @@ dev = [
 "nanotech_empa.gaussian.constr_opt_chain" = "aiida_nanotech_empa.workflows.gaussian:GaussianConstrOptChainWorkChain"
 "nanotech_empa.gaussian.casscf" = "aiida_nanotech_empa.workflows.gaussian:GaussianCasscfWorkChain"
 "nanotech_empa.gaussian.casscf_series" = "aiida_nanotech_empa.workflows.gaussian:GaussianCasscfSeriesWorkChain"
+"nanotech_empa.gaussian.nics" = "aiida_nanotech_empa.workflows.gaussian:GaussianNicsWorkChain"
 "nanotech_empa.cp2k.geo_opt" = "aiida_nanotech_empa.workflows.cp2k:Cp2kGeoOptWorkChain"
 "nanotech_empa.cp2k.fragment_separation" = "aiida_nanotech_empa.workflows.cp2k:Cp2kFragmentSeparationWorkChain"
 "nanotech_empa.cp2k.ads_gw_ic" = "aiida_nanotech_empa.workflows.cp2k:Cp2kAdsorbedGwIcWorkChain"

From 2e48743f5bfb9e9db031a5619f5c0de32892f0dc Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Thu, 12 Sep 2024 14:01:52 +0000
Subject: [PATCH 4/9] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 aiida_nanotech_empa/utils/cycle_tools.py      | 60 ++++++++-------
 aiida_nanotech_empa/utils/nmr.py              | 74 ++++++++++---------
 .../workflows/gaussian/__init__.py            |  2 +-
 .../workflows/gaussian/nics_workchain.py      | 57 +++++++-------
 .../workflows/gaussian/oldnics_workchain.py   | 13 ++--
 .../workflows/gaussian/relax_workchain.py     | 12 +--
 .../workflows/gaussian/scf_workchain.py       | 16 ++--
 examples/workflows/benzene.xyz                | 24 +++---
 examples/workflows/example_gaussian_nics.py   | 27 +++----
 examples/workflows/example_gaussian_opt.py    | 18 ++---
 10 files changed, 156 insertions(+), 147 deletions(-)

diff --git a/aiida_nanotech_empa/utils/cycle_tools.py b/aiida_nanotech_empa/utils/cycle_tools.py
index f4bed8d..feaf9ca 100644
--- a/aiida_nanotech_empa/utils/cycle_tools.py
+++ b/aiida_nanotech_empa/utils/cycle_tools.py
@@ -1,38 +1,34 @@
-import numpy as np
-
-import scipy as sp
-import scipy.linalg
-
 import ase
 import ase.io
-import ase.visualize
 import ase.neighborlist
-
-import os
-import shutil
+import ase.visualize
+import numpy as np
+import scipy as sp
+import scipy.linalg
 
 
 def convert_neighbor_list(nl):
     new = {}
-    n_vert = np.max(nl)+1
-    
+    n_vert = np.max(nl) + 1
+
     for i_v in range(n_vert):
         new[i_v] = []
-    
+
     for i_v, j_v in zip(nl[0], nl[1]):
-        
+
         new[i_v].append(j_v)
-        
+
     return new
 
+
 def find_cycles(i_vert, cnl, max_length, cur_path, passed_edges):
-    
-    if len(cur_path)-1 == max_length:
+
+    if len(cur_path) - 1 == max_length:
         return []
-    
+
     acc_cycles = []
     sort_cycles = []
-    
+
     neighbs = cnl[i_vert]
 
     # if we are connected to something that is not the end
@@ -49,45 +45,47 @@ def find_cycles(i_vert, cnl, max_length, cur_path, passed_edges):
     for n in neighbs:
         edge = (np.min([i_vert, n]), np.max([i_vert, n]))
         if edge not in passed_edges:
-            
+
             if n == cur_path[0]:
                 # found cycle
                 return [cur_path]
-    
+
     # Continue in all possible directions
     for n in neighbs:
         edge = (np.min([i_vert, n]), np.max([i_vert, n]))
         if edge not in passed_edges:
 
-            cycs = find_cycles(n, cnl, max_length, cur_path + [n], passed_edges + [edge])
+            cycs = find_cycles(
+                n, cnl, max_length, cur_path + [n], passed_edges + [edge]
+            )
             for cyc in cycs:
                 sorted_cyc = tuple(sorted(cyc))
                 if sorted_cyc not in sort_cycles:
                     sort_cycles.append(sorted_cyc)
                     acc_cycles.append(cyc)
-    
+
     return acc_cycles
-    
+
 
 def dumb_cycle_detection(ase_atoms_no_h, max_length):
-    
+
     neighbor_list = ase.neighborlist.neighbor_list("ij", ase_atoms_no_h, 2.0)
-    
+
     cycles = []
     sorted_cycles = []
-    n_vert = np.max(neighbor_list)+1
-    
+    n_vert = np.max(neighbor_list) + 1
+
     cnl = convert_neighbor_list(neighbor_list)
-    
+
     for i_vert in range(n_vert):
-        
+
         cycs = find_cycles(i_vert, cnl, max_length, [i_vert], [])
         for cyc in cycs:
             sorted_cyc = tuple(sorted(cyc))
             if sorted_cyc not in sorted_cycles:
                 sorted_cycles.append(sorted_cyc)
                 cycles.append(cyc)
-    
+
     return cycles
 
 
@@ -99,7 +97,7 @@ def cycle_normal(cycle, h):
     u, s, v = np.linalg.svd(points.T)
     normal = u[:, -1]
     normal /= np.linalg.norm(normal)
-    if np.dot(normal, h*np.array([1, 1, 1])) < 0.0:
+    if np.dot(normal, h * np.array([1, 1, 1])) < 0.0:
         normal *= -1.0
     return normal
 
diff --git a/aiida_nanotech_empa/utils/nmr.py b/aiida_nanotech_empa/utils/nmr.py
index 4145659..6108d6b 100644
--- a/aiida_nanotech_empa/utils/nmr.py
+++ b/aiida_nanotech_empa/utils/nmr.py
@@ -1,16 +1,11 @@
-import numpy as np
-
-import scipy as sp
-import scipy.linalg
-
 import ase
 import ase.io
-import ase.visualize
 import ase.neighborlist
-
-import os
-import shutil
+import ase.visualize
 import cclib
+import numpy as np
+import scipy as sp
+import scipy.linalg
 
 from . import cycle_tools
 
@@ -18,67 +13,73 @@
 ### SETUP
 
 
-def find_ref_points(ase_atoms_no_h, cycles, h = 0.0):
+def find_ref_points(ase_atoms_no_h, cycles, h=0.0):
     """
     positive h means projection to z axis is positive and vice-versa
     """
 
-    centers, normals = cycle_tools.find_cycle_centers_and_normals(ase_atoms_no_h, cycles, h)
-    
+    centers, normals = cycle_tools.find_cycle_centers_and_normals(
+        ase_atoms_no_h, cycles, h
+    )
+
     ase_ref_p = ase.Atoms()
 
     for i_cyc in range(len(cycles)):
         if h == 0.0:
-            ase_ref_p.append(ase.Atom('X', centers[i_cyc]))
+            ase_ref_p.append(ase.Atom("X", centers[i_cyc]))
         else:
-            pos = centers[i_cyc] + np.abs(h)*normals[i_cyc]
-            ase_ref_p.append(ase.Atom('X', pos))
-            
+            pos = centers[i_cyc] + np.abs(h) * normals[i_cyc]
+            ase_ref_p.append(ase.Atom("X", pos))
+
     return ase_ref_p
 
+
 def dist(p1, p2):
     if isinstance(p1, ase.Atom):
         p1 = p1.position
     if isinstance(p2, ase.Atom):
         p2 = p2.position
-    return np.linalg.norm(p2-p1)
+    return np.linalg.norm(p2 - p1)
+
 
 def interp_pts(p1, p2, dx):
     vec = p2 - p1
     dist = np.sqrt(np.sum(vec**2))
-    num_p = int(np.round(dist/dx))
-    dx_real = dist/num_p
-    dvec = vec/dist*dx_real
-    
+    num_p = int(np.round(dist / dx))
+    dx_real = dist / num_p
+    dvec = vec / dist * dx_real
+
     points = np.outer(np.arange(0, num_p), dvec) + p1
     return points
 
+
 def build_path(ref_pts, dx=0.1):
 
     point_arr = None
 
-    for i_rp in range(len(ref_pts)-1):
+    for i_rp in range(len(ref_pts) - 1):
 
         pt1 = ref_pts[i_rp].position
-        pt2 = ref_pts[i_rp+1].position
+        pt2 = ref_pts[i_rp + 1].position
 
         points = interp_pts(pt1, pt2, dx)
 
-        if i_rp == len(ref_pts)-2:
+        if i_rp == len(ref_pts) - 2:
             points = np.concatenate([points, [pt2]], axis=0)
 
         if point_arr is None:
             point_arr = points
         else:
             point_arr = np.concatenate([point_arr, points], axis=0)
-    
-    ase_arr = ase.Atoms("X%d"%len(point_arr), point_arr)
+
+    ase_arr = ase.Atoms("X%d" % len(point_arr), point_arr)
     return ase_arr
 
 
 ### -----------------------------------------------------------------
 ### PROCESS
 
+
 def load_nics_gaussian(nics_path):
 
     sigma = []
@@ -110,14 +111,15 @@ def is_number(x):
         return True
     except ValueError:
         return False
-    
+
+
 def parse_nmr_cmo_matrix(log_file_str, property_dict):
-    
+
     lines = log_file_str.splitlines()
-    
+
     # build the object
-    n_atom = property_dict['natom']
-    n_occupied_mo = property_dict['homos'][0] + 1
+    n_atom = property_dict["natom"]
+    n_occupied_mo = property_dict["homos"][0] + 1
 
     nmr_cmo_matrix = np.zeros((n_atom, n_occupied_mo, 3, 3))
 
@@ -138,24 +140,24 @@ def parse_nmr_cmo_matrix(log_file_str, property_dict):
 
         if "Full Cartesian NMR shielding tensor (ppm) for atom" in lines[i_line]:
 
-            i_atom = int(lines[i_line].replace('(', ')').split(')')[-2]) - 1
+            i_atom = int(lines[i_line].replace("(", ")").split(")")[-2]) - 1
 
             i_line += 1
 
-            if 'Canonical MO contributions' in lines[i_line]:
+            if "Canonical MO contributions" in lines[i_line]:
 
                 for _i in range(2000):
                     i_line += 1
-                    if 'Total' in lines[i_line]:
+                    if "Total" in lines[i_line]:
                         break
                     split = lines[i_line].split()
 
                     if len(split) == 10 and is_number(split[-1]):
                         i_mo = int(split[0][:-1]) - 1
 
-                        arr = np.array([float(x) for x in split[1:]])                    
+                        arr = np.array([float(x) for x in split[1:]])
                         nmr_cmo_matrix[i_atom, i_mo, :, :] = arr.reshape(3, 3)
 
         i_line += 1
-        
+
     return nmr_cmo_matrix
diff --git a/aiida_nanotech_empa/workflows/gaussian/__init__.py b/aiida_nanotech_empa/workflows/gaussian/__init__.py
index 0b5d156..809a790 100644
--- a/aiida_nanotech_empa/workflows/gaussian/__init__.py
+++ b/aiida_nanotech_empa/workflows/gaussian/__init__.py
@@ -4,10 +4,10 @@
 from .delta_scf_workchain import GaussianDeltaScfWorkChain
 from .hf_mp2_workchain import GaussianHfMp2WorkChain
 from .natorb_workchain import GaussianNatOrbWorkChain
+from .nics_workchain import GaussianNicsWorkChain
 from .relax_workchain import GaussianRelaxWorkChain
 from .scf_workchain import GaussianScfWorkChain
 from .spin_workchain import GaussianSpinWorkChain
-from .nics_workchain import GaussianNicsWorkChain
 
 __all__ = (
     "GaussianScfWorkChain",
diff --git a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
index 966a42f..e667e83 100644
--- a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
@@ -1,29 +1,33 @@
-import numpy as np
 import ase
+import numpy as np
 from aiida import engine, orm
+
+from ...utils import common_utils, cycle_tools, nmr
 from .relax_workchain import GaussianRelaxWorkChain
 from .scf_workchain import GaussianScfWorkChain
-from ...utils import cycle_tools, nmr
-from ...utils import common_utils
+
 
 @engine.calcfunction
-def prepare_nmr_structure(structure,height):
+def prepare_nmr_structure(structure, height):
     ase_geo = structure.get_ase()
     pos = ase_geo.positions
-    extents = np.array([
-        np.max(pos[:,0]) - np.min(pos[:,0]),
-        np.max(pos[:,1]) - np.min(pos[:,1]),
-        np.max(pos[:,2]) - np.min(pos[:,2]),
-    ])
+    extents = np.array(
+        [
+            np.max(pos[:, 0]) - np.min(pos[:, 0]),
+            np.max(pos[:, 1]) - np.min(pos[:, 1]),
+            np.max(pos[:, 2]) - np.min(pos[:, 2]),
+        ]
+    )
     inds = np.argsort(-extents)
     new_pos = pos[:, inds]
     ase_atoms = ase.Atoms(numbers=ase_geo.numbers, positions=new_pos)
-    ase_atoms_no_h = ase.Atoms([a for a in ase_atoms if a.symbol != 'H'])
+    ase_atoms_no_h = ase.Atoms([a for a in ase_atoms if a.symbol != "H"])
     cycles = cycle_tools.dumb_cycle_detection(ase_atoms_no_h, 8)
     ref_p = nmr.find_ref_points(ase_atoms_no_h, cycles, height.value)
     new_ase = ase_atoms + ref_p
     return orm.StructureData(ase=new_ase)
 
+
 class GaussianNicsWorkChain(engine.WorkChain):
     @classmethod
     def define(cls, spec):
@@ -81,11 +85,11 @@ def define(cls, spec):
             valid_type=orm.Dict,
             required=False,
             help="Use custom metadata.options instead of the automatic ones.",
-        )        
-        
+        )
+
         spec.outline(
             cls.setup,
-            engine.if_(cls.should_submit_opt)(cls.submit_opt,cls.inspect_opt),            
+            engine.if_(cls.should_submit_opt)(cls.submit_opt, cls.inspect_opt),
             cls.submit_nics,
             cls.inspect_nics,
             cls.finalize,
@@ -98,14 +102,15 @@ def define(cls, spec):
             "ERROR_TERMINATION",
             message="One or more steps of the work chain failed.",
         )
+
     def setup(self):
         self.report("Setting up...")
         self.ctx.nmr_structure = self.inputs.structure
         self.ctx_should_opt = getattr(self.inputs, "opt", True)
-        
+
     def should_submit_opt(self):
         return self.ctx_should_opt
-    
+
     def submit_opt(self):
         self.report("Submitting optimization...")
         label = "geo_opt"
@@ -127,23 +132,25 @@ def submit_opt(self):
         submitted_node = self.submit(builder)
         submitted_node.description = label
         self.to_context(**{label: submitted_node})
-        
+
     def inspect_opt(self):
         self.report("Inspecting optimization...")
         label = "geo_opt"
         # check if everything finished nicely
         if not common_utils.check_if_calc_ok(self, self.ctx[label]):
-            return self.exit_codes.ERROR_TERMINATION        
+            return self.exit_codes.ERROR_TERMINATION
         self.ctx.nmr_structure = self.ctx[label].outputs.output_structure
         return engine.ExitCode(0)
-    
+
     def submit_nics(self):
         self.report("Submitting NICS calculation...")
         label = "nics"
         builder = GaussianScfWorkChain.get_builder()
         height = getattr(self.inputs, "height", 1.0)
         builder.gaussian_code = self.inputs.gaussian_code
-        builder.structure = prepare_nmr_structure(self.ctx.nmr_structure,orm.Float(height))
+        builder.structure = prepare_nmr_structure(
+            self.ctx.nmr_structure, orm.Float(height)
+        )
         builder.functional = self.inputs.functional
         builder.basis_set = self.inputs.basis_set
         builder.multiplicity = self.inputs.multiplicity
@@ -151,23 +158,23 @@ def submit_nics(self):
         builder.cdiis = orm.Bool(True)
         builder.nmr = orm.Bool(True)
         builder.maxcycle = orm.Int(2048)
-        builder.conver = orm.Int(8)        
+        builder.conver = orm.Int(8)
         if "options" in self.inputs:
             builder.options = self.inputs.options
 
         submitted_node = self.submit(builder)
         submitted_node.description = label
         self.to_context(**{label: submitted_node})
-        
+
     def inspect_nics(self):
         self.report("Inspecting nics...")
         label = "nics"
         # check if everything finished nicely
         if not common_utils.check_if_calc_ok(self, self.ctx[label]):
             return self.exit_codes.ERROR_TERMINATION
-        
+
         self.out("output_parameters", self.ctx[label].outputs.output_parameters)
-        return engine.ExitCode(0)   
-         
+        return engine.ExitCode(0)
+
     def finalize(self):
-        self.report("Finalizing...")
\ No newline at end of file
+        self.report("Finalizing...")
diff --git a/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py
index 06fe78e..7a7a002 100644
--- a/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py
@@ -1,7 +1,8 @@
 import numpy as np
 from aiida import engine, orm, plugins
 
-from ...utils import common_utils,nmr,cycle_tools
+from ...utils import common_utils, cycle_tools, nmr
+
 #from .delta_scf_workchain import GaussianDeltaScfWorkChain
 #from .natorb_workchain import GaussianNatOrbWorkChain
 #from .relax_workchain import GaussianRelaxWorkChain
@@ -14,7 +15,7 @@
 
 @engine.calcfunction
 def nics_structure(structure=None, h=1.0):
-    
+
     ase_geo = structure.get_ase()
     # orient the geometry
     pos = ase_geo.positions
@@ -25,10 +26,10 @@ def nics_structure(structure=None, h=1.0):
     ])
     inds = np.argsort(-extents)
     ase_atoms = ase.Atoms(numbers=ase_geo.numbers, positions=pos[:, inds])
-    ase_atoms_no_h = ase.Atoms([a for a in ase_atoms if a.symbol != 'H'])        
+    ase_atoms_no_h = ase.Atoms([a for a in ase_atoms if a.symbol != 'H'])
     cycles = cycle_tools.dumb_cycle_detection(ase_atoms_no_h, 8)
-    ase_geom = ase_atoms + nmr.find_ref_points(ase_atoms_no_h, cycles, h.value) 
-    return orm.StructureData(ase=ase_geom)   
+    ase_geom = ase_atoms + nmr.find_ref_points(ase_atoms_no_h, cycles, h.value)
+    return orm.StructureData(ase=ase_geom)
 
 class GaussianNicsWorkChain(engine.WorkChain):
     @classmethod
@@ -189,7 +190,7 @@ def submit_next_steps(self):
             builder.structure = self.ctx.gs_structure
             builder.functional = self.inputs.functional
             builder.empirical_dispersion = self.inputs.empirical_dispersion
-            builder.basis_set = 
+            builder.basis_set =
             builder.multiplicity = self.inputs.multiplicity
 
             #if "options" in self.inputs:
diff --git a/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py b/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py
index 1ea2d76..fe6b210 100644
--- a/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py
@@ -60,13 +60,13 @@ def define(cls, spec):
             valid_type=orm.Bool,
             required=False,
             help="Conjugate Direct Inversion in the Iterative Subspace",
-        )        
+        )
         spec.input(
             "conver",
             valid_type=orm.Int,
             required=False,
             help="the scf convergece threshold",
-        )        
+        )
         spec.input(
             "int",
             valid_type=orm.Str,
@@ -313,10 +313,10 @@ def optimization(self):
 
         if len(opt_dict) != 0:
             parameters["route_parameters"]["opt"] = opt_dict
-        conver = getattr(self.inputs,"conver",None)
-        maxcycle = getattr(self.inputs,"maxcycle",None)
-        grid = getattr(self.inputs,"int",None)
-        cdiis = getattr(self.inputs,"cdiis",None)
+        conver = getattr(self.inputs, "conver", None)
+        maxcycle = getattr(self.inputs, "maxcycle", None)
+        grid = getattr(self.inputs, "int", None)
+        cdiis = getattr(self.inputs, "cdiis", None)
         if grid:
             parameters["route_parameters"]["int"] = grid
         if maxcycle:
diff --git a/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py b/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py
index db8644e..e36adda 100644
--- a/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/scf_workchain.py
@@ -71,13 +71,13 @@ def define(cls, spec):
             valid_type=orm.Int,
             required=False,
             help="the scf convergece threshold",
-        )  
+        )
         spec.input(
             "cdiis",
             valid_type=orm.Bool,
             required=False,
             help="Conjugate Direct Inversion in the Iterative Subspace",
-        )               
+        )
         spec.input(
             "int",
             valid_type=orm.Str,
@@ -254,15 +254,15 @@ def scf(self):
                 "charge": 0,
                 "multiplicity": self.ctx.mult,
                 "route_parameters": {
-                    "scf": {"maxcycle": 140 },
+                    "scf": {"maxcycle": 140},
                 },
             }
         )
-        nmr = getattr(self.inputs,"nmr",None)
-        conver = getattr(self.inputs,"conver",None)
-        maxcycle = getattr(self.inputs,"maxcycle",None)
-        grid = getattr(self.inputs,"int",None)
-        cdiis = getattr(self.inputs,"cdiis",None)
+        nmr = getattr(self.inputs, "nmr", None)
+        conver = getattr(self.inputs, "conver", None)
+        maxcycle = getattr(self.inputs, "maxcycle", None)
+        grid = getattr(self.inputs, "int", None)
+        cdiis = getattr(self.inputs, "cdiis", None)
         if grid:
             parameters["route_parameters"]["int"] = grid
         if maxcycle:
diff --git a/examples/workflows/benzene.xyz b/examples/workflows/benzene.xyz
index 66fd521..5255bb3 100644
--- a/examples/workflows/benzene.xyz
+++ b/examples/workflows/benzene.xyz
@@ -1,14 +1,14 @@
 12
 benzene
-C        8.39528138       6.66748365       5.00000291   
-C        8.51925793       8.06102240       5.00000509   
-C        7.12645219       6.07808103       5.00000441   
-C        5.98160072       6.88221755       5.00000464   
-C        6.10557756       8.27575669       5.00000369   
-C        7.37440623       8.86515917       5.00000734   
-H        9.50085998       8.51699940       5.00000031   
-H        9.28097172       6.04538232       5.00000427   
-H        7.03053864       5.00000000       5.00000036   
-H        5.00000000       6.42623909       5.00000439   
-H        5.21988789       8.89785938       5.00000000   
-H        7.47031658       9.94324034       5.00000253   
\ No newline at end of file
+C        8.39528138       6.66748365       5.00000291
+C        8.51925793       8.06102240       5.00000509
+C        7.12645219       6.07808103       5.00000441
+C        5.98160072       6.88221755       5.00000464
+C        6.10557756       8.27575669       5.00000369
+C        7.37440623       8.86515917       5.00000734
+H        9.50085998       8.51699940       5.00000031
+H        9.28097172       6.04538232       5.00000427
+H        7.03053864       5.00000000       5.00000036
+H        5.00000000       6.42623909       5.00000439
+H        5.21988789       8.89785938       5.00000000
+H        7.47031658       9.94324034       5.00000253
diff --git a/examples/workflows/example_gaussian_nics.py b/examples/workflows/example_gaussian_nics.py
index 9c43f4b..12fa82f 100644
--- a/examples/workflows/example_gaussian_nics.py
+++ b/examples/workflows/example_gaussian_nics.py
@@ -1,16 +1,16 @@
-import os
 import click
-from ase.build import molecule
 import numpy as np
 from aiida import engine, orm
 from aiida.plugins import WorkflowFactory
+from ase.build import molecule
 
 import aiida_nanotech_empa.utils.gaussian_wcs_postprocess as pp
 
 GaussianNicsWorkChain = WorkflowFactory("nanotech_empa.gaussian.nics")
 
-def _example_gaussian_nics(gaussian_code,opt):
-    ase_geom = molecule('C6H6')
+
+def _example_gaussian_nics(gaussian_code, opt):
+    ase_geom = molecule("C6H6")
     ase_geom.cell = np.diag([10.0, 10.0, 10.0])
 
     builder = GaussianNicsWorkChain.get_builder()
@@ -26,25 +26,26 @@ def _example_gaussian_nics(gaussian_code,opt):
                 "num_machines": 1,
             },
             "max_wallclock_seconds": 1 * 60 * 60,
-            "max_memory_kb": 8 * 1024 * 1024, #GB
+            "max_memory_kb": 8 * 1024 * 1024,  # GB
         }
-    )    
-    
+    )
+
     _, wc_node = engine.run_get_node(builder)
-    
+
     assert wc_node.is_finished_ok
     return wc_node.pk
 
+
 @click.command("cli")
 @click.argument("gaussian_code", default="gaussian@localhost")
 def run_nics(gaussian_code):
-    #print("#### Running Gaussian NICS WorkChain with geo_opt ####")
-    #uuid = _example_gaussian_nics(orm.load_code(gaussian_code),True)
-    #print(f"WorkChain completed uuid: {uuid}")
+    # print("#### Running Gaussian NICS WorkChain with geo_opt ####")
+    # uuid = _example_gaussian_nics(orm.load_code(gaussian_code),True)
+    # print(f"WorkChain completed uuid: {uuid}")
     print("#### Running Gaussian NICS WorkChain without geo_opt ####")
-    uuid = _example_gaussian_nics(orm.load_code(gaussian_code),False)
+    uuid = _example_gaussian_nics(orm.load_code(gaussian_code), False)
     print(f"WorkChain completed uuid: {uuid}")
 
 
 if __name__ == "__main__":
-    run_nics()
\ No newline at end of file
+    run_nics()
diff --git a/examples/workflows/example_gaussian_opt.py b/examples/workflows/example_gaussian_opt.py
index f7beade..7a9bb70 100644
--- a/examples/workflows/example_gaussian_opt.py
+++ b/examples/workflows/example_gaussian_opt.py
@@ -3,7 +3,7 @@
 import ase.io
 import numpy as np
 from aiida.engine import run_get_node
-from aiida.orm import Bool,Int, Str, StructureData, load_code
+from aiida.orm import Bool, Int, Str, StructureData, load_code
 from aiida.plugins import WorkflowFactory
 
 import aiida_nanotech_empa.utils.gaussian_wcs_postprocess as pp
@@ -14,19 +14,19 @@
 OUTPUT_DIR = os.path.dirname(os.path.abspath(__file__))
 
 
-def _example_gaussian_spin(gaussian_code):#, formchk_code, cubegen_code):
+def _example_gaussian_spin(gaussian_code):  # , formchk_code, cubegen_code):
     ase_geom = ase.io.read(os.path.join(DATA_DIR, "benzene.xyz"))
     ase_geom.cell = np.diag([10.0, 10.0, 10.0])
 
     builder = GaussianRelaxWorkChain.get_builder()
     builder.gaussian_code = gaussian_code
-    #builder.formchk_code = formchk_code
-    #builder.cubegen_code = cubegen_code
+    # builder.formchk_code = formchk_code
+    # builder.cubegen_code = cubegen_code
     builder.structure = StructureData(ase=ase_geom)
     builder.functional = Str("B3LYP")
     builder.empirical_dispersion = Str("GD3")
     builder.basis_set = Str("6-311+G(d,p)")
-    #builder.basis_set_scf = Str("STO-3G")
+    # builder.basis_set_scf = Str("STO-3G")
     builder.multiplicity = Int(0)
     builder.tight = Bool(True)
     builder.options = Dict(
@@ -38,18 +38,18 @@ def _example_gaussian_spin(gaussian_code):#, formchk_code, cubegen_code):
             "max_wallclock_seconds": 1 * 60 * 60,
             "max_memory_kb": 4 * 1024 * 1024,
         }
-    )    
+    )
 
     _, wc_node = run_get_node(builder)
 
     assert wc_node.is_finished_ok
 
-    #pp.make_report(wc_node, nb=False, save_image_loc=OUTPUT_DIR)
+    # pp.make_report(wc_node, nb=False, save_image_loc=OUTPUT_DIR)
 
 
 if __name__ == "__main__":
     _example_gaussian_spin(
         load_code("gaussian@tigu"),
-        #load_code("formchk@localhost"),
-        #load_code("cubegen@localhost"),
+        # load_code("formchk@localhost"),
+        # load_code("cubegen@localhost"),
     )

From eebc0d06ce58a5cdd5ff257d82ddcee25d148c86 Mon Sep 17 00:00:00 2001
From: Carlo Antonio Pignedoli <c.pignedoli@gmail.com>
Date: Wed, 18 Sep 2024 11:45:48 +0000
Subject: [PATCH 5/9] added extras

---
 aiida_nanotech_empa/workflows/gaussian/nics_workchain.py | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
index 966a42f..05851c4 100644
--- a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
@@ -167,7 +167,12 @@ def inspect_nics(self):
             return self.exit_codes.ERROR_TERMINATION
         
         self.out("output_parameters", self.ctx[label].outputs.output_parameters)
+        self.out("output_structure", self.ctx.nics.inputs.structure)
         return engine.ExitCode(0)   
          
     def finalize(self):
-        self.report("Finalizing...")
\ No newline at end of file
+        self.report("Finalizing...")
+        
+        # Add extras.
+        struc = self.inputs.structure
+        common_utils.add_extras(struc, "surfaces", self.node.uuid)
\ No newline at end of file

From 75d2f93539cb4c9ad258e04fa77b145a2ea05da7 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 18 Sep 2024 13:52:50 +0000
Subject: [PATCH 6/9] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 aiida_nanotech_empa/workflows/gaussian/nics_workchain.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
index 9fd4c9e..2486a19 100644
--- a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
@@ -175,11 +175,11 @@ def inspect_nics(self):
 
         self.out("output_parameters", self.ctx[label].outputs.output_parameters)
         self.out("output_structure", self.ctx.nics.inputs.structure)
-        return engine.ExitCode(0)   
-         
+        return engine.ExitCode(0)
+
     def finalize(self):
         self.report("Finalizing...")
-        
+
         # Add extras.
         struc = self.inputs.structure
-        common_utils.add_extras(struc, "surfaces", self.node.uuid)
\ No newline at end of file
+        common_utils.add_extras(struc, "surfaces", self.node.uuid)

From 84256d6f16aa12ae05cfe1d4d65106ebed0ed2c3 Mon Sep 17 00:00:00 2001
From: Carlo Antonio Pignedoli <c.pignedoli@gmail.com>
Date: Wed, 25 Sep 2024 09:07:30 +0000
Subject: [PATCH 7/9] working

---
 .../workflows/gaussian/nics_workchain.py      |   1 -
 .../workflows/gaussian/oldnics_workchain.py   | 267 ------------------
 examples/workflows/example_gaussian_nics.py   |  25 +-
 examples/workflows/naphthalene.xyz            |  20 ++
 4 files changed, 40 insertions(+), 273 deletions(-)
 delete mode 100644 aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py
 create mode 100644 examples/workflows/naphthalene.xyz

diff --git a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
index 2486a19..fb46d69 100644
--- a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
@@ -6,7 +6,6 @@
 from .relax_workchain import GaussianRelaxWorkChain
 from .scf_workchain import GaussianScfWorkChain
 
-
 @engine.calcfunction
 def prepare_nmr_structure(structure, height):
     ase_geo = structure.get_ase()
diff --git a/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py
deleted file mode 100644
index 7a7a002..0000000
--- a/aiida_nanotech_empa/workflows/gaussian/oldnics_workchain.py
+++ /dev/null
@@ -1,267 +0,0 @@
-import numpy as np
-from aiida import engine, orm, plugins
-
-from ...utils import common_utils, cycle_tools, nmr
-
-#from .delta_scf_workchain import GaussianDeltaScfWorkChain
-#from .natorb_workchain import GaussianNatOrbWorkChain
-#from .relax_workchain import GaussianRelaxWorkChain
-f#rom .scf_workchain import GaussianScfWorkChain
-#from aiida.plugins import CalculationFactory
-
-#GaussianBaseWorkChain = plugins.WorkflowFactory("gaussian.base")
-#GaussianCubesWorkChain = plugins.WorkflowFactory("gaussian.cubes")
-GaussianCalculation = plugins.CalculationFactory("gaussian")
-
-@engine.calcfunction
-def nics_structure(structure=None, h=1.0):
-
-    ase_geo = structure.get_ase()
-    # orient the geometry
-    pos = ase_geo.positions
-    extents = np.array([
-        np.max(pos[:,0]) - np.min(pos[:,0]),
-        np.max(pos[:,1]) - np.min(pos[:,1]),
-        np.max(pos[:,2]) - np.min(pos[:,2]),
-    ])
-    inds = np.argsort(-extents)
-    ase_atoms = ase.Atoms(numbers=ase_geo.numbers, positions=pos[:, inds])
-    ase_atoms_no_h = ase.Atoms([a for a in ase_atoms if a.symbol != 'H'])
-    cycles = cycle_tools.dumb_cycle_detection(ase_atoms_no_h, 8)
-    ase_geom = ase_atoms + nmr.find_ref_points(ase_atoms_no_h, cycles, h.value)
-    return orm.StructureData(ase=ase_geom)
-
-class GaussianNicsWorkChain(engine.WorkChain):
-    @classmethod
-    def define(cls, spec):
-        super().define(spec)
-
-        spec.input("gaussian_code", valid_type=orm.Code)
-
-        spec.input(
-            "structure",
-            valid_type=orm.StructureData,
-            required=True,
-            help="input geometry",
-        )
-        spec.input(
-            "functional", valid_type=orm.Str, required=False,default=lambda: orm.Str("b3lyp"), help="xc functional"
-        )
-        spec.input(
-            "basis_set", valid_type=orm.Str, required=False, default=lambda: orm.Str("6-311+G(d,p)"), help="basis_set for opt"
-        )
-        spec.input(
-            "multiplicity",
-            valid_type=orm.Int,
-            required=True,
-            help="spin multiplicity",
-        )
-
-        spec.outline(
-            cls.setup
-            cls.submit_nics,
-            cls.finalize,
-        )
-
-        spec.outputs.dynamic = True
-
-        spec.exit_code(
-            390,
-            "ERROR_TERMINATION",
-            message="One or more steps of the work chain failed.",
-        )
-
-    def setup(self):
-        self.ctx.structure_wtih_centers = nics_structure(self.inputs.structure, h=1.0)
-        self.ctx.parameters = {
-                "link0_parameters": {
-                    "%chk": "aiida.chk",
-                    "%mem": "4096MB",
-                    "%nprocshared": 4,
-                },
-                "functional": "BLYP",
-                "basis_set": "6-311+G(d,p)",
-                "multiplicity": self.inputs.multiplicity.value,
-                "dieze_tag": "#P",
-                "route_parameters": {
-                    "scf": {
-                        "maxcycle":2048,
-                        "conver":8,
-                    },
-                    "int": "superfine",
-                    "chpf":{"conver":8}
-                    "nmr": None,
-                },
-            }
-
-    def submit_nics(self):
-        # Multiplicity 0 means RKS calculation.
-        label = f"m{self.inputs.multiplicity.value}_opt"
-
-        # Main parameters: geometry optimization
-        parameters = Dict(
-            {
-                "link0_parameters": {
-                    "%chk": "aiida.chk",
-                    "%mem": "4096MB",
-                    "%nprocshared": 4,
-                },
-                "functional": "BLYP",
-                "basis_set": "6-311+G(d,p)",
-                "multiplicity": self.inputs.multiplicity.value,
-                "dieze_tag": "#P",
-                "route_parameters": {
-                    "scf": {
-                        "maxcycle":2048,
-                        "conver":8,
-                    },
-                    "int": "superfine",
-                    "chpf":{"conver":8}
-                    "nmr": None,
-                },
-            }
-        )
-
-        # Construct process builder
-
-        builder = GaussianCalculation.get_builder()
-
-        builder.structure = structure
-        builder.parameters = parameters
-        builder.code = gaussian_code
-
-        builder.metadata.options.resources = {
-            "num_machines": 1,
-            "tot_num_mpiprocs": num_cores,
-        }
-
-        # Should ask for extra +25% extra memory
-        builder.metadata.options.max_memory_kb = int(1.25 * memory_mb) * 1024
-        builder.metadata.options.max_wallclock_seconds = 5 * 60
-
-        submitted_node = self.submit(builder)
-        submitted_node.description = label
-        self.to_context(**{label: submitted_node})
-
-    def inspect_opt(self):
-
-        label = f"m{self.inputs.multiplicity.value}_opt"
-
-        # check if everything finished nicely
-        if not common_utils.check_if_calc_ok(self, self.ctx[label]):
-            return self.exit_codes.ERROR_TERMINATION
-
-        opt_energy = self.ctx[label].outputs.scf_energy_ev
-        self.out("opt_energy", opt_energy)
-        self.out("opt_structure", self.ctx[label].outputs.output_structure)
-        self.out(
-            "opt_out_params", self.ctx[label].outputs.scf_output_parameters
-        )
-
-        return engine.ExitCode(0)
-
-    def submit_next_steps(self):
-
-        self.report("Submitting NMR cubes")
-
-        parameters = orm.Dict(
-                    {
-                        "link0_parameters": self.ctx.link0.copy(),
-                        "dieze_tag": "#P",
-                        "functional": self.ctx.functional,
-                        "basis_set": self.inputs.basis_set.value,
-                        "charge": 0,
-                        "multiplicity": self.ctx.mult,
-                        "route_parameters": {
-                            "scf": {"maxcycle": 140},
-                            "opt": None,
-                        },
-                    }
-                )
-
-        for mult in self.inputs.multiplicity_list:
-            label = "nmr"
-
-            if mult == self.ctx.gs_mult:
-                continue
-
-            builder = GaussianScfWorkChain.get_builder()
-            builder.gaussian_code = self.inputs.gaussian_code
-            builder.structure = self.ctx.gs_structure
-            builder.functional = self.inputs.functional
-            builder.empirical_dispersion = self.inputs.empirical_dispersion
-            builder.basis_set =
-            builder.multiplicity = self.inputs.multiplicity
-
-            #if "options" in self.inputs:
-            #    builder.options = self.inputs.options
-
-            submitted_node = self.submit(builder)
-            submitted_node.description = label
-            self.to_context(**{label: submitted_node})
-
-    def inspect_next_steps(self):
-        if not common_utils.check_if_calc_ok(self, self.ctx.gs_cubes):
-            return self.exit_codes.ERROR_TERMINATION
-
-        self.out("gs_cube_images", self.ctx.gs_cubes.outputs.cube_image_folder)
-        self.out("gs_cube_planes", self.ctx.gs_cubes.outputs.cube_planes_array)
-
-        if not common_utils.check_if_calc_ok(self, self.ctx.dscf):
-            return self.exit_codes.ERROR_TERMINATION
-
-        self.out("gs_ionization_potential", self.ctx.dscf.outputs.ionization_potential)
-        self.out("gs_electron_affinity", self.ctx.dscf.outputs.electron_affinity)
-
-        for mult in self.inputs.multiplicity_list:
-            label = f"m{mult}_vert"
-
-            if mult == self.ctx.gs_mult:
-                continue
-
-            # Check if everything finished nicely.
-            if not common_utils.check_if_calc_ok(self, self.ctx[label]):
-                return self.exit_codes.ERROR_TERMINATION
-
-            vert_energy = self.ctx[label].outputs.energy_ev
-            self.out(f"m{mult}_vert_energy", vert_energy)
-            self.out(
-                f"m{mult}_vert_out_params", self.ctx[label].outputs.output_parameters
-            )
-            self.out(
-                f"m{mult}_vert_cube_images", self.ctx[label].outputs.cube_image_folder
-            )
-            self.out(
-                f"m{mult}_vert_cube_planes", self.ctx[label].outputs.cube_planes_array
-            )
-
-        return engine.ExitCode(0)
-
-    def is_gs_oss(self):
-        """Is ground state an open-shell singlet?"""
-        return self.ctx.gs_mult == 1
-
-    def submit_nat_orb(self):
-        self.report("Submitting natural pop analysis")
-
-        builder = GaussianNatOrbWorkChain.get_builder()
-        builder.gaussian_code = self.inputs.gaussian_code
-        builder.parent_calc_folder = self.ctx.gs_scf_remote_folder
-        builder.parent_calc_params = self.ctx.gs_out_params
-        if "options" in self.inputs:
-            builder.options = self.inputs.options
-
-        submitted_node = self.submit(builder)
-        submitted_node.description = "natural orbitals pop"
-        self.to_context(natorb=submitted_node)
-
-    def inspect_nat_orb(self):
-        if not common_utils.check_if_calc_ok(self, self.ctx.natorb):
-            return self.exit_codes.ERROR_TERMINATION
-
-        self.out("gs_natorb_params", self.ctx.natorb.outputs.natorb_proc_parameters)
-
-        return engine.ExitCode(0)
-
-    def finalize(self):
-        self.report("Finalizing...")
diff --git a/examples/workflows/example_gaussian_nics.py b/examples/workflows/example_gaussian_nics.py
index 12fa82f..caddecd 100644
--- a/examples/workflows/example_gaussian_nics.py
+++ b/examples/workflows/example_gaussian_nics.py
@@ -1,21 +1,36 @@
+import pathlib
 import click
 import numpy as np
 from aiida import engine, orm
 from aiida.plugins import WorkflowFactory
 from ase.build import molecule
+from ase.io import read
 
 import aiida_nanotech_empa.utils.gaussian_wcs_postprocess as pp
 
 GaussianNicsWorkChain = WorkflowFactory("nanotech_empa.gaussian.nics")
-
+DATA_DIR = pathlib.Path(__file__).parent.absolute()
+GEO_FILE="naphthalene.xyz"
 
 def _example_gaussian_nics(gaussian_code, opt):
-    ase_geom = molecule("C6H6")
-    ase_geom.cell = np.diag([10.0, 10.0, 10.0])
+    # Check test geometry is already in database.
+    qb = orm.QueryBuilder()
+    qb.append(orm.Node, filters={"label": {"==": GEO_FILE}})
+    structure = None
+    for node_tuple in qb.iterall():
+        node = node_tuple[0]
+        structure = node
+    if structure is not None:
+        print("found existing structure: ", structure.pk)
+    else:
+        structure = orm.StructureData(ase=read(DATA_DIR / GEO_FILE))
+        structure.label = GEO_FILE
+        structure.store()
+        print("created new structure: ", structure.pk)
 
     builder = GaussianNicsWorkChain.get_builder()
     builder.gaussian_code = gaussian_code
-    builder.structure = orm.StructureData(ase=ase_geom)
+    builder.structure = structure
     builder.functional = orm.Str("B3LYP")
     builder.basis_set = orm.Str("6-311G(d,p)")
     builder.opt = orm.Bool(opt)
@@ -43,7 +58,7 @@ def run_nics(gaussian_code):
     # uuid = _example_gaussian_nics(orm.load_code(gaussian_code),True)
     # print(f"WorkChain completed uuid: {uuid}")
     print("#### Running Gaussian NICS WorkChain without geo_opt ####")
-    uuid = _example_gaussian_nics(orm.load_code(gaussian_code), False)
+    uuid = _example_gaussian_nics(orm.load_code(gaussian_code), True)
     print(f"WorkChain completed uuid: {uuid}")
 
 
diff --git a/examples/workflows/naphthalene.xyz b/examples/workflows/naphthalene.xyz
new file mode 100644
index 0000000..2f13db8
--- /dev/null
+++ b/examples/workflows/naphthalene.xyz
@@ -0,0 +1,20 @@
+18
+Properties=species:S:1:pos:R:3 napthalene=T pbc="F F F"
+C        0.00000000       1.42000000       0.00000000
+C        0.00000000       0.00000000       0.00000000
+C        1.22975600       2.13000000       0.00000000
+C        2.45951200       1.42000000       0.00000000
+C        2.45951200       0.00000000       0.00000000
+C        1.22975600      -0.71000000       0.00000000
+C        3.68926800      -0.71000000       0.00000000
+C        4.91902400       0.00000000       0.00000000
+C        4.91902400       1.42000000       0.00000000
+H       -0.98726900       1.99000000       0.00000000
+H       -0.98726900      -0.57000000       0.00000000
+H        1.22975600       3.27000000       0.00000000
+H        1.22975600      -1.85000000       0.00000000
+H        3.68926800      -1.85000000       0.00000000
+H        3.68926800       3.27000000       0.00000000
+H        5.90629300       1.99000000       0.00000000
+H        5.90629300      -0.57000000       0.00000000
+C        3.68926800       2.13000000       0.00000000

From bec152d961558ae2514c64ecd7dfaf57d36ed03e Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 25 Sep 2024 09:09:52 +0000
Subject: [PATCH 8/9] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 aiida_nanotech_empa/workflows/gaussian/nics_workchain.py | 1 +
 examples/workflows/example_gaussian_nics.py              | 4 +++-
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
index fb46d69..2486a19 100644
--- a/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
@@ -6,6 +6,7 @@
 from .relax_workchain import GaussianRelaxWorkChain
 from .scf_workchain import GaussianScfWorkChain
 
+
 @engine.calcfunction
 def prepare_nmr_structure(structure, height):
     ase_geo = structure.get_ase()
diff --git a/examples/workflows/example_gaussian_nics.py b/examples/workflows/example_gaussian_nics.py
index caddecd..295741d 100644
--- a/examples/workflows/example_gaussian_nics.py
+++ b/examples/workflows/example_gaussian_nics.py
@@ -1,4 +1,5 @@
 import pathlib
+
 import click
 import numpy as np
 from aiida import engine, orm
@@ -10,7 +11,8 @@
 
 GaussianNicsWorkChain = WorkflowFactory("nanotech_empa.gaussian.nics")
 DATA_DIR = pathlib.Path(__file__).parent.absolute()
-GEO_FILE="naphthalene.xyz"
+GEO_FILE = "naphthalene.xyz"
+
 
 def _example_gaussian_nics(gaussian_code, opt):
     # Check test geometry is already in database.

From dd31dd05c7782f24f41db3df62594a7eb945fcc9 Mon Sep 17 00:00:00 2001
From: Carlo Antonio Pignedoli <c.pignedoli@gmail.com>
Date: Thu, 3 Oct 2024 12:51:54 +0000
Subject: [PATCH 9/9] added extras to input structure

---
 .../workflows/gaussian/casscf_series_workchain.py             | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/aiida_nanotech_empa/workflows/gaussian/casscf_series_workchain.py b/aiida_nanotech_empa/workflows/gaussian/casscf_series_workchain.py
index 4aac31f..0baf0ec 100644
--- a/aiida_nanotech_empa/workflows/gaussian/casscf_series_workchain.py
+++ b/aiida_nanotech_empa/workflows/gaussian/casscf_series_workchain.py
@@ -267,5 +267,7 @@ def finalize(self):
                         f"{var}_cube_image_folder",
                         self.ctx[var].outputs.cube_image_folder,
                     )
-
+        # Add extras.
+        struc = self.inputs.structure
+        common_utils.add_extras(struc, "surfaces", self.node.uuid)
         return engine.ExitCode(0)