diff --git a/aiida_nanotech_empa/utils/cycle_tools.py b/aiida_nanotech_empa/utils/cycle_tools.py
new file mode 100644
index 0000000..feaf9ca
--- /dev/null
+++ b/aiida_nanotech_empa/utils/cycle_tools.py
@@ -0,0 +1,119 @@
+import ase
+import ase.io
+import ase.neighborlist
+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
+
+    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..6108d6b
--- /dev/null
+++ b/aiida_nanotech_empa/utils/nmr.py
@@ -0,0 +1,163 @@
+import ase
+import ase.io
+import ase.neighborlist
+import ase.visualize
+import cclib
+import numpy as np
+import scipy as sp
+import scipy.linalg
+
+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/__init__.py b/aiida_nanotech_empa/workflows/gaussian/__init__.py
index d08b655..809a790 100644
--- a/aiida_nanotech_empa/workflows/gaussian/__init__.py
+++ b/aiida_nanotech_empa/workflows/gaussian/__init__.py
@@ -4,6 +4,7 @@
 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
@@ -18,4 +19,5 @@
     "GaussianConstrOptChainWorkChain",
     "GaussianCasscfWorkChain",
     "GaussianCasscfSeriesWorkChain",
+    "GaussianNicsWorkChain",
 )
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)
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..2486a19
--- /dev/null
+++ b/aiida_nanotech_empa/workflows/gaussian/nics_workchain.py
@@ -0,0 +1,185 @@
+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
+
+
+@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)
+
+        spec.input("gaussian_code", valid_type=orm.Code)
+
+        spec.input(
+            "structure",
+            valid_type=orm.StructureData,
+            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"
+        )
+
+        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",
+        )
+        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,
+        )
+
+        spec.outputs.dynamic = True
+
+        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)
+        self.out("output_structure", self.ctx.nics.inputs.structure)
+        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)
diff --git a/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py b/aiida_nanotech_empa/workflows/gaussian/relax_workchain.py
index 9c03573..fe6b210 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..e36adda 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)
@@ -228,6 +258,24 @@ def scf(self):
                 },
             }
         )
+        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/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..5255bb3
--- /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
diff --git a/examples/workflows/example_gaussian_nics.py b/examples/workflows/example_gaussian_nics.py
new file mode 100644
index 0000000..295741d
--- /dev/null
+++ b/examples/workflows/example_gaussian_nics.py
@@ -0,0 +1,68 @@
+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):
+    # 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 = structure
+    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)
+
+    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 without geo_opt ####")
+    uuid = _example_gaussian_nics(orm.load_code(gaussian_code), True)
+    print(f"WorkChain completed uuid: {uuid}")
+
+
+if __name__ == "__main__":
+    run_nics()
diff --git a/examples/workflows/example_gaussian_opt.py b/examples/workflows/example_gaussian_opt.py
new file mode 100644
index 0000000..7a9bb70
--- /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"),
+    )
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
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"