Nics (#166)

* work in progress

* not much done

* working

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* added extras

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* working

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* added extras to input structure

---------

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

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)

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

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",
 )

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)