From 760183434ff40d3e73ca4a2cec42b19a0da43042 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Mon, 11 Dec 2023 16:10:47 +0000
Subject: [PATCH 01/47] initial commit

---
 autode/mol_graphs.py                 |   9 ++
 autode/opt/coordinates/__init__.py   |   2 +-
 autode/opt/coordinates/base.py       |  10 --
 autode/opt/coordinates/internals.py  | 138 ++++++++++++++++++++++++++-
 autode/opt/coordinates/primitives.py |   2 +-
 tests/test_opt/test_coordiantes.py   |   1 -
 6 files changed, 147 insertions(+), 15 deletions(-)

diff --git a/autode/mol_graphs.py b/autode/mol_graphs.py
index 9ef148633..c01ae1b7b 100644
--- a/autode/mol_graphs.py
+++ b/autode/mol_graphs.py
@@ -115,6 +115,15 @@ def node_matcher(self):
 
         return matcher
 
+    @property
+    def is_connected(self) -> bool:
+        """Is this graph fully connected (i.e. not separate parts)"""
+        return nx.is_connected(self)
+
+    def get_components(self):
+        """Generate the separate connected components"""
+        return nx.connected_components(self)
+
 
 def make_graph(
     species: "Species",
diff --git a/autode/opt/coordinates/__init__.py b/autode/opt/coordinates/__init__.py
index 2aadb5de9..36daee459 100644
--- a/autode/opt/coordinates/__init__.py
+++ b/autode/opt/coordinates/__init__.py
@@ -1,3 +1,3 @@
-from autode.opt.coordinates.base import OptCoordinates, CartesianComponent
+from autode.opt.coordinates.base import OptCoordinates
 from autode.opt.coordinates.cartesian import CartesianCoordinates
 from autode.opt.coordinates.dic import DIC, DICWithConstraints
diff --git a/autode/opt/coordinates/base.py b/autode/opt/coordinates/base.py
index 0f5135ced..747680e99 100644
--- a/autode/opt/coordinates/base.py
+++ b/autode/opt/coordinates/base.py
@@ -1,7 +1,6 @@
 # mypy: disable-error-code="has-type"
 import numpy as np
 from copy import deepcopy
-from enum import IntEnum, unique
 from typing import Optional, Union, Sequence, List, TYPE_CHECKING
 from abc import ABC, abstractmethod
 
@@ -15,15 +14,6 @@
     from autode.hessians import Hessian
 
 
-@unique
-class CartesianComponent(IntEnum):
-    """Cartesian component in 3D space"""
-
-    x = 0
-    y = 1
-    z = 2
-
-
 class OptCoordinates(ValueArray, ABC):
     """Coordinates used to perform optimisations"""
 
diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 9e7462e16..73522cbcb 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -7,19 +7,26 @@
 q : Primitive internal coordinates
 G : Spectroscopic G matrix
 """
-import numpy as np
+import copy
 
+import numpy as np
+from enum import Enum
+import itertools
 from typing import Any, Optional, Type, List, TYPE_CHECKING
 from abc import ABC, abstractmethod
-from autode.opt.coordinates.base import OptCoordinates, CartesianComponent
+from autode.values import Angle
+from autode.opt.coordinates.base import OptCoordinates
 from autode.opt.coordinates.primitives import (
     PrimitiveInverseDistance,
     Primitive,
     PrimitiveDistance,
+    ConstrainedPrimitiveDistance,
+    PrimitiveBondAngle,
     PrimitiveDihedralAngle,
 )
 
 if TYPE_CHECKING:
+    from autode.species import Species
     from autode.opt.coordinates.cartesian import CartesianCoordinates
     from autode.opt.coordinates.primitives import (
         ConstrainedPrimitive,
@@ -234,3 +241,130 @@ def _primitive_type(self):
 class AnyPIC(PIC):
     def _populate_all(self, x: np.ndarray) -> None:
         raise RuntimeError("Cannot populate all on an AnyPIC instance")
+
+
+class LinearAngleType(Enum):
+    linear_bend = 1
+    cos_angle = 2
+    remove = 3
+
+
+def build_pic_from_species(
+    species,
+    aux_bonds=False,
+    linear_angles=LinearAngleType.remove,
+    robust_dihedrals=False,
+):
+    pic = AnyPIC()
+    core_graph = species.graph.copy()
+    _add_distances_from_species(pic, species, core_graph, aux_bonds=aux_bonds)
+    _add_bends_from_species(
+        pic, species, core_graph, linear_angles=linear_angles
+    )
+    _add_dihedrals_from_species(
+        pic, species, core_graph, robust_dihedrals=robust_dihedrals
+    )
+
+
+def _get_connected_graph_from_species(mol: "Species"):
+    if mol.graph is None:
+        mol.reset_graph()
+
+    assert mol.graph is not None
+    core_graph = copy.deepcopy(mol.graph)
+
+    # join disconnected graphs
+    if core_graph.is_connected:
+        for components in core_graph.get_components():
+            pass
+
+    # The constraints should be counted as bonds
+    if mol.constraints.distance is not None:
+        for i, j in mol.constraints.distance:
+            if not core_graph.has_edge(i, j):
+                core_graph.add_edge(i, j, pi=False, active=False)
+
+
+def _add_distances_from_species(pic, species, core_graph, aux_bonds=True):
+    n = 0
+    for i, j in sorted(core_graph.edges):
+        if (
+            species.constraints.distance is not None
+            and (i, j) in species.constraints.distance
+        ):
+            r = species.constraints.distance[(i, j)]
+            pic.append(ConstrainedPrimitiveDistance(i, j, r))
+            n += 1
+        else:
+            pic.append(PrimitiveDistance(i, j))
+    assert n == species.constraints.n_distance
+
+    if not aux_bonds:
+        return None
+
+    for i, j in itertools.combinations(range(species.n_atoms), r=2):
+        if core_graph.has_edge(i, j):
+            continue
+        if species.distance(i, j) < 2.5 * species.eqm_bond_distance(i, j):
+            pic.append(PrimitiveDistance(i, j))
+    return None
+
+
+def _add_bends_from_species(
+    pic, species, core_graph, linear_angles=LinearAngleType.linear_bend
+):
+    for o in range(species.n_atoms):
+        for n, m in itertools.combinations(core_graph.neighbors(o), r=2):
+            if species.angle(m, o, n) < Angle(175, "deg"):
+                pic.append(PrimitiveBondAngle(o=o, m=m, n=n))
+            elif linear_angles == LinearAngleType.linear_bend:
+                pass  # todo linear bends
+            elif linear_angles == LinearAngleType.cos_angle:
+                pass
+            elif linear_angles == LinearAngleType.remove:
+                pass
+            else:
+                raise Exception(
+                    "Linear angle handling method not properly defined"
+                )
+
+    return None
+
+
+def _add_dihedrals_from_species(
+    pic, species, core_graph, robust_dihedrals=False
+):
+    # no dihedrals possible with less than 4 atoms
+    if species.n_atoms < 4:
+        return
+
+    for o, p in core_graph.edges:
+        for m in species.graph.neighbors(o):
+            if m == p:
+                continue
+
+            if species.angle(m, o, p) > Angle(175, "deg"):
+                continue
+
+            for n in species.graph.neighbors(p):
+                if n == o:
+                    continue
+
+                # avoid triangle rings like cyclopropane
+                if n == m:
+                    continue
+
+                is_linear_1 = species.angle(m, o, p) > Angle(175, "deg")
+                is_linear_2 = species.angle(o, p, n) > Angle(175, "deg")
+
+                # don't add when both angles are linear
+                if is_linear_1 and is_linear_2:
+                    continue
+
+                # if only one angle almost linear, add robust dihedral
+                if (is_linear_1 or is_linear_2) and robust_dihedrals:
+                    pass  # todo robust dihedrals
+                else:
+                    pic.append(PrimitiveDihedralAngle(m, o, p, n))
+
+    return None
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 764c322b0..5d2619730 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -11,7 +11,7 @@
 )
 
 if TYPE_CHECKING:
-    from autode.opt.coordinates import CartesianCoordinates, CartesianComponent
+    from autode.opt.coordinates import CartesianCoordinates
 
 
 def _get_3d_vecs_from_atom_idxs(
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 9fe31f72e..1900267a4 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -5,7 +5,6 @@
 from autode.species.molecule import Molecule
 from autode.values import Angle
 from autode.exceptions import CoordinateTransformFailed
-from autode.opt.coordinates.base import CartesianComponent
 from autode.opt.coordinates.internals import PrimitiveInverseDistances, PIC
 from autode.opt.coordinates.cartesian import CartesianCoordinates
 from autode.opt.coordinates.dic import DIC

From b3e4b3ac8567fe54fb444764b590fbaea2a57a2a Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Wed, 13 Dec 2023 14:54:14 +0000
Subject: [PATCH 02/47] second commit

---
 autode/opt/coordinates/internals.py | 87 ++++++++++++++++++++++-------
 1 file changed, 67 insertions(+), 20 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 73522cbcb..e845c0701 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -6,6 +6,8 @@
 B : Wilson B matrix
 q : Primitive internal coordinates
 G : Spectroscopic G matrix
+
+Set-up of redundant primitives is based on J. Chem. Phys., 117, 2002, 9160
 """
 import copy
 
@@ -27,6 +29,7 @@
 
 if TYPE_CHECKING:
     from autode.species import Species
+    from autode.mol_graphs import MolecularGraph
     from autode.opt.coordinates.cartesian import CartesianCoordinates
     from autode.opt.coordinates.primitives import (
         ConstrainedPrimitive,
@@ -250,33 +253,64 @@ class LinearAngleType(Enum):
 
 
 def build_pic_from_species(
-    species,
+    mol: "Species",
     aux_bonds=False,
     linear_angles=LinearAngleType.remove,
     robust_dihedrals=False,
 ):
     pic = AnyPIC()
-    core_graph = species.graph.copy()
-    _add_distances_from_species(pic, species, core_graph, aux_bonds=aux_bonds)
-    _add_bends_from_species(
-        pic, species, core_graph, linear_angles=linear_angles
-    )
+    core_graph = _get_connected_graph_from_species(mol)
+    _add_bonds_from_species(pic, mol, core_graph, aux_bonds=aux_bonds)
+    _add_angles_from_species(pic, mol, core_graph, linear_angles=linear_angles)
     _add_dihedrals_from_species(
-        pic, species, core_graph, robust_dihedrals=robust_dihedrals
+        pic, mol, core_graph, robust_dihedrals=robust_dihedrals
     )
 
 
-def _get_connected_graph_from_species(mol: "Species"):
+def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
+    """
+    Creates a fully connected graph from a species, by (1) joining
+    disconnected fragments by their shortest distance, (2) connecting
+    constrained bonds, (3) joining hydrogen bonds, if present.
+
+    Args:
+        mol: A species containing atoms and coordinates
+
+    Returns:
+        (MolecularGraph):
+    """
+    # if graph does not exist, create one
     if mol.graph is None:
         mol.reset_graph()
 
     assert mol.graph is not None
     core_graph = copy.deepcopy(mol.graph)
 
-    # join disconnected graphs
-    if core_graph.is_connected:
-        for components in core_graph.get_components():
-            pass
+    # join disconnected graph components
+    if not core_graph.is_connected:
+        components = core_graph.get_components()
+        for comp_i, comp_j in itertools.combinations(components, r=2):
+            min_dist = float("inf")
+            min_pair = (-1, -1)
+            for i, j in itertools.product(list(comp_i), list(comp_j)):
+                if mol.distance(i, j) < min_dist:
+                    min_dist = mol.distance(i, j)
+                    min_pair = (i, j)
+            core_graph.add_edge(*min_pair, pi=False, active=False)
+
+    # join hydrogen bonds
+    h_bond_x = ["N", "O", "F", "P", "S", "Cl"]
+    for i, j in itertools.combinations(range(mol.n_atoms), r=2):
+        if (
+            mol.atoms[i].label in h_bond_x
+            and mol.atoms[j].label == "H"
+            or mol.atoms[j].label in h_bond_x
+            and mol.atoms[i].label == "H"
+        ):
+            vdw_sum = mol.atoms[i].vdw_radius + mol.atoms[j].vdw_radius
+            if mol.distance(i, j) < 0.9 * vdw_sum:
+                if not core_graph.has_edge(i, j):
+                    core_graph.add_edge(i, j, pi=False, active=False)
 
     # The constraints should be counted as bonds
     if mol.constraints.distance is not None:
@@ -284,33 +318,46 @@ def _get_connected_graph_from_species(mol: "Species"):
             if not core_graph.has_edge(i, j):
                 core_graph.add_edge(i, j, pi=False, active=False)
 
+    return core_graph
+
+
+def _add_bonds_from_species(pic, mol, core_graph, aux_bonds=False):
+    """
+    Modify the supplied AnyPIC instance in-place by adding bonds, from the
+    connectivity graph supplied
 
-def _add_distances_from_species(pic, species, core_graph, aux_bonds=True):
+    Args:
+        pic: The AnyPIC instance (modified in-place)
+        mol: The species object
+        core_graph: The connectivity graph
+        aux_bonds: Whether to add auxiliary bonds (< 2.5 * covalent radii sum)
+    """
     n = 0
     for i, j in sorted(core_graph.edges):
         if (
-            species.constraints.distance is not None
-            and (i, j) in species.constraints.distance
+            mol.constraints.distance is not None
+            and (i, j) in mol.constraints.distance
         ):
-            r = species.constraints.distance[(i, j)]
+            r = mol.constraints.distance[(i, j)]
             pic.append(ConstrainedPrimitiveDistance(i, j, r))
             n += 1
         else:
             pic.append(PrimitiveDistance(i, j))
-    assert n == species.constraints.n_distance
+    assert n == mol.constraints.n_distance
 
     if not aux_bonds:
         return None
 
-    for i, j in itertools.combinations(range(species.n_atoms), r=2):
+    # add auxiliary bonds if specified
+    for i, j in itertools.combinations(range(mol.n_atoms), r=2):
         if core_graph.has_edge(i, j):
             continue
-        if species.distance(i, j) < 2.5 * species.eqm_bond_distance(i, j):
+        if mol.distance(i, j) < 2.5 * mol.eqm_bond_distance(i, j):
             pic.append(PrimitiveDistance(i, j))
     return None
 
 
-def _add_bends_from_species(
+def _add_angles_from_species(
     pic, species, core_graph, linear_angles=LinearAngleType.linear_bend
 ):
     for o in range(species.n_atoms):

From 31659ee92faa8229efd4d286afb935a56116e202 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Wed, 13 Dec 2023 17:47:20 +0000
Subject: [PATCH 03/47] second commit

---
 autode/opt/coordinates/internals.py | 73 +++++++++++++++++------------
 1 file changed, 44 insertions(+), 29 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index e845c0701..85e74d60d 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -246,22 +246,15 @@ def _populate_all(self, x: np.ndarray) -> None:
         raise RuntimeError("Cannot populate all on an AnyPIC instance")
 
 
-class LinearAngleType(Enum):
-    linear_bend = 1
-    cos_angle = 2
-    remove = 3
-
-
 def build_pic_from_species(
     mol: "Species",
     aux_bonds=False,
-    linear_angles=LinearAngleType.remove,
     robust_dihedrals=False,
 ):
     pic = AnyPIC()
     core_graph = _get_connected_graph_from_species(mol)
     _add_bonds_from_species(pic, mol, core_graph, aux_bonds=aux_bonds)
-    _add_angles_from_species(pic, mol, core_graph, linear_angles=linear_angles)
+    _add_angles_from_species(pic, mol, core_graph)
     _add_dihedrals_from_species(
         pic, mol, core_graph, robust_dihedrals=robust_dihedrals
     )
@@ -357,43 +350,65 @@ def _add_bonds_from_species(pic, mol, core_graph, aux_bonds=False):
     return None
 
 
+class LinearAngleType(Enum):
+    linear_bend = 1
+    cos_angle = 2
+    remove = 3
+
+
 def _add_angles_from_species(
-    pic, species, core_graph, linear_angles=LinearAngleType.linear_bend
-):
-    for o in range(species.n_atoms):
+    pic: AnyPIC, mol: "Species", core_graph: "MolecularGraph"
+) -> None:
+    """
+    Modify the set of primitives in-place by adding angles, from the
+    connectivity graph supplied
+
+    Args:
+        pic: The AnyPIC instance (modified in-place)
+        mol: The species object
+        core_graph: The connectivity graph
+    """
+    for o in range(mol.n_atoms):
         for n, m in itertools.combinations(core_graph.neighbors(o), r=2):
-            if species.angle(m, o, n) < Angle(175, "deg"):
+            # avoid almost linear angles
+            if mol.angle(m, o, n) < Angle(175, "deg"):
                 pic.append(PrimitiveBondAngle(o=o, m=m, n=n))
-            elif linear_angles == LinearAngleType.linear_bend:
-                pass  # todo linear bends
-            elif linear_angles == LinearAngleType.cos_angle:
-                pass
-            elif linear_angles == LinearAngleType.remove:
-                pass
-            else:
-                raise Exception(
-                    "Linear angle handling method not properly defined"
-                )
 
     return None
 
 
 def _add_dihedrals_from_species(
-    pic, species, core_graph, robust_dihedrals=False
+    pic: AnyPIC,
+    mol: "Species",
+    core_graph: "MolecularGraph",
+    robust_dihedrals=False,
 ):
+    """
+    Modify the set of primitives in-place by adding dihedrals (torsions),
+    from the connectivity graph supplied
+
+    Args:
+        pic: The AnyPIC instance (modified in-place)
+        mol: The species
+        core_graph: The connectivity graph
+        robust_dihedrals:
+
+    Returns:
+
+    """
     # no dihedrals possible with less than 4 atoms
-    if species.n_atoms < 4:
+    if mol.n_atoms < 4:
         return
 
     for o, p in core_graph.edges:
-        for m in species.graph.neighbors(o):
+        for m in core_graph.neighbors(o):
             if m == p:
                 continue
 
-            if species.angle(m, o, p) > Angle(175, "deg"):
+            if mol.angle(m, o, p) > Angle(175, "deg"):
                 continue
 
-            for n in species.graph.neighbors(p):
+            for n in core_graph.neighbors(p):
                 if n == o:
                     continue
 
@@ -401,8 +416,8 @@ def _add_dihedrals_from_species(
                 if n == m:
                     continue
 
-                is_linear_1 = species.angle(m, o, p) > Angle(175, "deg")
-                is_linear_2 = species.angle(o, p, n) > Angle(175, "deg")
+                is_linear_1 = mol.angle(m, o, p) > Angle(175, "deg")
+                is_linear_2 = mol.angle(o, p, n) > Angle(175, "deg")
 
                 # don't add when both angles are linear
                 if is_linear_1 and is_linear_2:

From 0a88af77139bef83f89cdbf29c4357bb39286144 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Thu, 14 Dec 2023 14:45:09 +0000
Subject: [PATCH 04/47] simplify code

---
 autode/opt/coordinates/internals.py | 43 +++++++++++---------
 autode/opt/optimisers/crfo.py       | 61 ++++-------------------------
 autode/opt/optimisers/trm.py        | 21 ----------
 3 files changed, 33 insertions(+), 92 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 85e74d60d..bb6ddec90 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -249,15 +249,26 @@ def _populate_all(self, x: np.ndarray) -> None:
 def build_pic_from_species(
     mol: "Species",
     aux_bonds=False,
-    robust_dihedrals=False,
-):
+) -> AnyPIC:
+    """
+    Build a set of primitives from the species, using the graph as
+    a starting point for the connectivity of the species. Also joins
+    any disjoint parts of the graph, and adds hydrogen bonds to
+    ensure that the primitives are redundant
+
+    Args:
+        mol:
+        aux_bonds:
+
+    Returns:
+        (AnyPIC): The set of primitive internals
+    """
     pic = AnyPIC()
     core_graph = _get_connected_graph_from_species(mol)
     _add_bonds_from_species(pic, mol, core_graph, aux_bonds=aux_bonds)
     _add_angles_from_species(pic, mol, core_graph)
-    _add_dihedrals_from_species(
-        pic, mol, core_graph, robust_dihedrals=robust_dihedrals
-    )
+    _add_dihedrals_from_species(pic, mol, core_graph)
+    return pic
 
 
 def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
@@ -314,7 +325,12 @@ def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
     return core_graph
 
 
-def _add_bonds_from_species(pic, mol, core_graph, aux_bonds=False):
+def _add_bonds_from_species(
+    pic: AnyPIC,
+    mol: "Species",
+    core_graph: "MolecularGraph",
+    aux_bonds: bool = False,
+):
     """
     Modify the supplied AnyPIC instance in-place by adding bonds, from the
     connectivity graph supplied
@@ -381,8 +397,7 @@ def _add_dihedrals_from_species(
     pic: AnyPIC,
     mol: "Species",
     core_graph: "MolecularGraph",
-    robust_dihedrals=False,
-):
+) -> None:
     """
     Modify the set of primitives in-place by adding dihedrals (torsions),
     from the connectivity graph supplied
@@ -391,10 +406,6 @@ def _add_dihedrals_from_species(
         pic: The AnyPIC instance (modified in-place)
         mol: The species
         core_graph: The connectivity graph
-        robust_dihedrals:
-
-    Returns:
-
     """
     # no dihedrals possible with less than 4 atoms
     if mol.n_atoms < 4:
@@ -419,13 +430,9 @@ def _add_dihedrals_from_species(
                 is_linear_1 = mol.angle(m, o, p) > Angle(175, "deg")
                 is_linear_2 = mol.angle(o, p, n) > Angle(175, "deg")
 
-                # don't add when both angles are linear
-                if is_linear_1 and is_linear_2:
+                # if any angle is linear, don't add dihedral
+                if is_linear_1 or is_linear_2:
                     continue
-
-                # if only one angle almost linear, add robust dihedral
-                if (is_linear_1 or is_linear_2) and robust_dihedrals:
-                    pass  # todo robust dihedrals
                 else:
                     pic.append(PrimitiveDihedralAngle(m, o, p, n))
 
diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py
index 8dceaccf5..d5cc96c2a 100644
--- a/autode/opt/optimisers/crfo.py
+++ b/autode/opt/optimisers/crfo.py
@@ -1,4 +1,5 @@
-"""Constrained rational function optimisation
+"""
+Constrained rational function optimisation
 
 Notation follows:
 [1] J. Baker, J. Comput. Chem., 18, 8 1080
@@ -11,7 +12,11 @@
 from autode.log import logger
 from autode.values import GradientRMS, Angle, Distance
 from autode.opt.coordinates import CartesianCoordinates, DICWithConstraints
-from autode.opt.coordinates.internals import PIC, AnyPIC
+from autode.opt.coordinates.internals import (
+    PIC,
+    AnyPIC,
+    build_pic_from_species,
+)
 from autode.opt.optimisers.rfo import RFOptimiser
 from autode.opt.optimisers.hessian_update import BFGSDampedUpdate, NullUpdate
 from autode.opt.coordinates.primitives import (
@@ -125,15 +130,7 @@ def _build_internal_coordinates(self):
             )
 
         cartesian_coords = CartesianCoordinates(self._species.coordinates)
-        primitives = self._primitives
-
-        if len(primitives) < cartesian_coords.expected_number_of_dof:
-            logger.info(
-                "Had an incomplete set of primitives. Adding "
-                "additional distances"
-            )
-            for i, j in combinations(range(self._species.n_atoms), 2):
-                primitives.append(PrimitiveDistance(i, j))
+        primitives = build_pic_from_species(self._species)
 
         self._coords = DICWithConstraints.from_cartesian(
             x=cartesian_coords, primitives=primitives
@@ -141,48 +138,6 @@ def _build_internal_coordinates(self):
         self._coords.zero_lagrangian_multipliers()
         return None
 
-    @property
-    def _primitives(self) -> PIC:
-        """Primitive internal coordinates in this molecule"""
-        assert self._species and self._species.graph
-        logger.info("Generating primitive internal coordinates")
-        graph = self._species.graph.copy()
-
-        # Any distance constraints should also count as 'bonds' when forming
-        # the set of primitive internal coordinates, so that there is a
-        # single molecule if those distances are approaching dissociation
-        if self._species.constraints.distance is not None:
-            logger.info("Adding distance constraints as primitives")
-            for i, j in self._species.constraints.distance:
-                graph.add_edge(i, j)
-
-        pic = AnyPIC()
-
-        for i, j in sorted(graph.edges):
-            if (
-                self._species.constraints.distance
-                and (i, j) in self._species.constraints.distance
-            ):
-                r = self._species.constraints.distance[(i, j)]
-                pic.append(ConstrainedPrimitiveDistance(i, j, value=r))
-
-            else:
-                pic.append(PrimitiveDistance(i, j))
-
-        for o in range(self._species.n_atoms):
-            for n, m in combinations(graph.neighbors(o), r=2):
-                pic.append(PrimitiveBondAngle(o=o, m=m, n=n))
-
-        if self._species.n_atoms > 2 and not self._species.is_planar():
-            for dihedral in _dihedrals(self._species):
-                pic.append(dihedral)
-
-        logger.info(
-            f"Using {pic.n_constrained} constraints in {len(pic)} "
-            f"primitive internal coordinates"
-        )
-        return pic
-
     def _lambda_p_from_eigvals_and_gradient(
         self, b: np.ndarray, f: np.ndarray
     ) -> float:
diff --git a/autode/opt/optimisers/trm.py b/autode/opt/optimisers/trm.py
index 32b674143..a44a48f73 100644
--- a/autode/opt/optimisers/trm.py
+++ b/autode/opt/optimisers/trm.py
@@ -130,27 +130,6 @@ def _initialise_species_and_method(
         ), "HybridTRMOptimiser cannot work with constraints!"
         return None
 
-    def _build_internal_coordinates(self) -> None:
-        """ "Build delocalised internal coordinates"""
-        if self._species is None:
-            raise RuntimeError(
-                "Cannot set initial coordinates. No species set"
-            )
-
-        cart_coords = CartesianCoordinates(self._species.coordinates)
-        primitives = self._primitives
-
-        if len(primitives) < cart_coords.expected_number_of_dof:
-            logger.info(
-                "Had an incomplete set of primitives. Adding "
-                "additional distances"
-            )
-            for i, j in combinations(range(self._species.n_atoms), 2):
-                primitives.append(PrimitiveDistance(i, j))
-
-        self._coords = DIC.from_cartesian(x=cart_coords, primitives=primitives)
-        return None
-
     def _step(self) -> None:
         """
         Hybrid RFO/TRM step; if the RFO step is larger than the trust

From a261052a8de5e3065b4670c26d0105bc737600fa Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Thu, 14 Dec 2023 16:20:38 +0000
Subject: [PATCH 05/47] linear angle

---
 autode/opt/coordinates/_autodiff.py  | 11 +++--
 autode/opt/coordinates/primitives.py | 67 +++++++++++++++++++++++++++-
 2 files changed, 74 insertions(+), 4 deletions(-)

diff --git a/autode/opt/coordinates/_autodiff.py b/autode/opt/coordinates/_autodiff.py
index 0fa48d437..61c967fb9 100644
--- a/autode/opt/coordinates/_autodiff.py
+++ b/autode/opt/coordinates/_autodiff.py
@@ -580,7 +580,9 @@ class DifferentiableVector3D:
     hyper-dual numbers
     """
 
-    def __init__(self, items: Sequence["VectorHyperDual"]):
+    def __init__(
+        self, items: Sequence[Union["VectorHyperDual", numeric_type]]
+    ):
         """
         Initialise the 3D vector from a list of 3 hyperdual numbers
 
@@ -590,7 +592,9 @@ def __init__(self, items: Sequence["VectorHyperDual"]):
         items = list(items)
         if len(items) != 3:
             raise ValueError("A 3D vector must have only 3 components")
-        assert all(isinstance(item, VectorHyperDual) for item in items)
+        assert all(
+            isinstance(item, (VectorHyperDual, *numeric)) for item in items
+        )
         self._data = items
 
     @staticmethod
@@ -611,7 +615,7 @@ def dot(self, other: "DifferentiableVector3D") -> "VectorHyperDual":
             (VectorHyperDual): A scalar number (with derivatives)
         """
         self._check_same_type(other)
-        dot = 0
+        dot: Union[VectorHyperDual, numeric_type] = 0
         for k in range(3):
             dot = dot + self._data[k] * other._data[k]
         assert isinstance(dot, VectorHyperDual)
@@ -702,6 +706,7 @@ def cross(
         Returns:
             (DifferentiableVector3D):
         """
+        self._check_same_type(other)
         return DifferentiableVector3D(
             [
                 self._data[1] * other._data[2]
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 5d2619730..40ff0d912 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -1,7 +1,7 @@
 import numpy as np
 
 from abc import ABC, abstractmethod
-from typing import Tuple, TYPE_CHECKING, List
+from typing import Tuple, TYPE_CHECKING, List, Optional
 from autode.opt.coordinates._autodiff import (
     get_differentiable_vars,
     DifferentiableMath,
@@ -405,3 +405,68 @@ def _evaluate(
 
     def __repr__(self):
         return f"Dihedral({self.m}-{self.o}-{self.p}-{self.n})"
+
+
+class PrimitiveLinearAngle(Primitive):
+    def __init__(self, m: int, o: int, n: int, axis: int):
+        """Linear Bend: m-o-n"""
+        super().__init__(m, o, n)
+        self.m = int(m)
+        self.o = int(o)
+        self.n = int(n)
+
+        self.axis = int(axis)
+        if self.axis not in [1, 2]:
+            raise ValueError("The axis must be 1 or 2")
+
+        self.axis_vec: Optional[DifferentiableVector3D] = None
+
+    def __eq__(self, other):
+        return isinstance(other, self.__class__) and (
+            self._ordered_idxs == other._ordered_idxs
+            and self.o == other.o
+            and self.axis == other.axis
+        )
+        # todo check the sign of the bends if m and n swapped
+
+    def _init_axis(self, x: "CartesianCoordinates") -> None:
+        _x = x.reshape(-1, 3)
+        cart_axes = [
+            np.array([1.0, 0.0, 0.0]),
+            np.array([0.0, 1.0, 0.0]),
+            np.array([0.0, 0.0, 1.0]),
+        ]
+        # choose cartesian axis with the lowest overlap with m-n vector
+        _m, _n = _x[self.m], _x[self.n]
+        w = _m - _n
+        w /= np.linalg.norm(w)
+        overlaps = []
+        for axis in cart_axes:
+            overlaps.append(np.dot(w, axis))
+        cart_ax = cart_axes[np.argmin(overlaps)]
+        # make the axis completely perpendicular
+        perp_axis = cart_ax - np.dot(cart_ax, w) * w
+        self.axis_vec = DifferentiableVector3D(list(perp_axis))
+        return None
+
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
+    ):
+        """Linear Bend angle m-o-n against a Cartesian axis"""
+        if self.axis_vec is None:
+            self._init_axis(x)
+
+        assert self.axis_vec is not None
+        _x = x.ravel()
+        vec_m, vec_o, vec_n = _get_3d_vecs_from_atom_idxs(
+            self.m, self.o, self.n, x=_x, deriv_order=deriv_order
+        )
+        w = vec_m - vec_n
+
+        cross_vec = w.cross(self.axis_vec)
+        if self.axis == 2:
+            cross_vec = w.cross(cross_vec)
+
+        u = vec_m - vec_o
+        v = vec_n - vec_o
+        return cross_vec.dot(u.cross(v)) / (u.norm() * v.norm())

From 54530827d8e4b441bf04fe770212524ac0caf582 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Thu, 14 Dec 2023 16:50:52 +0000
Subject: [PATCH 06/47] linear angle

---
 autode/opt/coordinates/internals.py  | 5 +++++
 autode/opt/coordinates/primitives.py | 3 +++
 2 files changed, 8 insertions(+)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index bb6ddec90..39cafc377 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -24,6 +24,7 @@
     PrimitiveDistance,
     ConstrainedPrimitiveDistance,
     PrimitiveBondAngle,
+    PrimitiveLinearAngle,
     PrimitiveDihedralAngle,
 )
 
@@ -389,6 +390,10 @@ def _add_angles_from_species(
             # avoid almost linear angles
             if mol.angle(m, o, n) < Angle(175, "deg"):
                 pic.append(PrimitiveBondAngle(o=o, m=m, n=n))
+            else:
+                # stabilise linear angles by two orthogonal bends
+                pic.append(PrimitiveLinearAngle(m, o, n, 1))
+                pic.append(PrimitiveLinearAngle(m, o, n, 2))
 
     return None
 
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 40ff0d912..734b509a7 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -470,3 +470,6 @@ def _evaluate(
         u = vec_m - vec_o
         v = vec_n - vec_o
         return cross_vec.dot(u.cross(v)) / (u.norm() * v.norm())
+
+    def __repr__(self):
+        return f"LinearBend{self.axis}({self.m}-{self.o}-{self.n})"

From d46415c720eb2c452c375a272864cb011dc38c9a Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sun, 17 Dec 2023 15:52:08 +0000
Subject: [PATCH 07/47] test update, more fixes needed

---
 autode/opt/coordinates/_autodiff.py  | 12 +++++++
 autode/opt/coordinates/internals.py  | 13 +++-----
 autode/opt/coordinates/primitives.py | 30 +++++++++++------
 tests/test_opt/test_coordiantes.py   | 48 +++++++++++++++++++++++++---
 4 files changed, 82 insertions(+), 21 deletions(-)

diff --git a/autode/opt/coordinates/_autodiff.py b/autode/opt/coordinates/_autodiff.py
index 61c967fb9..5f635773e 100644
--- a/autode/opt/coordinates/_autodiff.py
+++ b/autode/opt/coordinates/_autodiff.py
@@ -694,6 +694,18 @@ def __rmul__(self, other):
         """Multiplication of scalar and vector is commutative"""
         return self.__mul__(other)
 
+    def __truediv__(self, other: Union[VectorHyperDual, numeric_type]):
+        """
+        Division of a 3D vector with a scalar
+
+        Args:
+            other (VectorHyperDual|float|int):
+
+        Returns:
+            (DifferentiableVector3D):
+        """
+        return self.__mul__(1 / other)
+
     def cross(
         self, other: "DifferentiableVector3D"
     ) -> "DifferentiableVector3D":
diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 39cafc377..4c7838bd2 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -26,6 +26,7 @@
     PrimitiveBondAngle,
     PrimitiveLinearAngle,
     PrimitiveDihedralAngle,
+    LinearBendType,
 )
 
 if TYPE_CHECKING:
@@ -367,12 +368,6 @@ def _add_bonds_from_species(
     return None
 
 
-class LinearAngleType(Enum):
-    linear_bend = 1
-    cos_angle = 2
-    remove = 3
-
-
 def _add_angles_from_species(
     pic: AnyPIC, mol: "Species", core_graph: "MolecularGraph"
 ) -> None:
@@ -392,8 +387,10 @@ def _add_angles_from_species(
                 pic.append(PrimitiveBondAngle(o=o, m=m, n=n))
             else:
                 # stabilise linear angles by two orthogonal bends
-                pic.append(PrimitiveLinearAngle(m, o, n, 1))
-                pic.append(PrimitiveLinearAngle(m, o, n, 2))
+                pic.append(PrimitiveLinearAngle(m, o, n, LinearBendType.BEND))
+                pic.append(
+                    PrimitiveLinearAngle(m, o, n, LinearBendType.COMPLEMENT)
+                )
 
     return None
 
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 734b509a7..afe39aed5 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -1,6 +1,7 @@
 import numpy as np
 
 from abc import ABC, abstractmethod
+from enum import Enum
 from typing import Tuple, TYPE_CHECKING, List, Optional
 from autode.opt.coordinates._autodiff import (
     get_differentiable_vars,
@@ -407,18 +408,23 @@ def __repr__(self):
         return f"Dihedral({self.m}-{self.o}-{self.p}-{self.n})"
 
 
+class LinearBendType(Enum):
+    """For linear angles, there are two orthogonal directions"""
+
+    BEND = 0
+    COMPLEMENT = 1
+
+
 class PrimitiveLinearAngle(Primitive):
-    def __init__(self, m: int, o: int, n: int, axis: int):
+    def __init__(self, m: int, o: int, n: int, axis: LinearBendType):
         """Linear Bend: m-o-n"""
         super().__init__(m, o, n)
         self.m = int(m)
         self.o = int(o)
         self.n = int(n)
 
-        self.axis = int(axis)
-        if self.axis not in [1, 2]:
-            raise ValueError("The axis must be 1 or 2")
-
+        assert isinstance(axis, LinearBendType)
+        self.axis = axis
         self.axis_vec: Optional[DifferentiableVector3D] = None
 
     def __eq__(self, other):
@@ -427,7 +433,7 @@ def __eq__(self, other):
             and self.o == other.o
             and self.axis == other.axis
         )
-        # todo check the sign of the bends if m and n swapped
+        # TODO: check the sign of the bends if m and n swapped
 
     def _init_axis(self, x: "CartesianCoordinates") -> None:
         _x = x.reshape(-1, 3)
@@ -436,6 +442,7 @@ def _init_axis(self, x: "CartesianCoordinates") -> None:
             np.array([0.0, 1.0, 0.0]),
             np.array([0.0, 0.0, 1.0]),
         ]
+
         # choose cartesian axis with the lowest overlap with m-n vector
         _m, _n = _x[self.m], _x[self.n]
         w = _m - _n
@@ -443,9 +450,11 @@ def _init_axis(self, x: "CartesianCoordinates") -> None:
         overlaps = []
         for axis in cart_axes:
             overlaps.append(np.dot(w, axis))
-        cart_ax = cart_axes[np.argmin(overlaps)]
-        # make the axis completely perpendicular
+        cart_ax = cart_axes[np.argmin(np.abs(overlaps))]
+
+        # make the axis completely perpendicular to m-n vector
         perp_axis = cart_ax - np.dot(cart_ax, w) * w
+        perp_axis /= np.linalg.norm(perp_axis)
         self.axis_vec = DifferentiableVector3D(list(perp_axis))
         return None
 
@@ -462,9 +471,12 @@ def _evaluate(
             self.m, self.o, self.n, x=_x, deriv_order=deriv_order
         )
         w = vec_m - vec_n
+        w = w / w.norm()
+        # TODO: does w need to be normalised?
 
         cross_vec = w.cross(self.axis_vec)
-        if self.axis == 2:
+        # if complement is requested, perform another cross product
+        if self.axis == LinearBendType.COMPLEMENT:
             cross_vec = w.cross(cross_vec)
 
         u = vec_m - vec_o
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 1900267a4..1bb01b45e 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -15,6 +15,8 @@
     PrimitiveBondAngle,
     ConstrainedPrimitiveBondAngle,
     PrimitiveDihedralAngle,
+    PrimitiveLinearAngle,
+    LinearBendType,
 )
 
 
@@ -611,6 +613,35 @@ def test_dihedral_equality():
     )
 
 
+def test_linear_angle():
+    acetylene = Molecule(
+        atoms=[
+            Atom("C", 0.35540, -0.20370, -0.44810),
+            Atom("C", -0.37180, 0.21470, 0.40200),
+            Atom("H", 1.01560, -0.60550, -1.23530),
+            Atom("H", -0.99920, 0.59450, 1.15720),
+        ]
+    )
+    x = CartesianCoordinates(acetylene.coordinates)
+    angle = PrimitiveLinearAngle(0, 1, 2, LinearBendType.BEND)
+    assert angle.axis_vec is None
+    angle._init_axis(x)
+    assert angle.axis_vec is not None
+    axis_vec = np.array(angle.axis_vec._data)
+    m_n_vec = acetylene.coordinates[0] - acetylene.coordinates[2]
+    assert np.dot(axis_vec, m_n_vec) < 0.001
+
+    # the atom order can be reversed without any change in the values
+    angle2 = PrimitiveLinearAngle(2, 1, 0, LinearBendType.BEND)
+    assert np.isclose(angle(x), angle2(x))
+    assert angle == angle2
+
+    # however, linear bend complement would not be the same
+    angle3 = PrimitiveLinearAngle(0, 1, 2, LinearBendType.COMPLEMENT)
+    assert angle != angle3
+    assert not np.isclose(angle(x), angle3(x), rtol=1e-3)
+
+
 def test_primitives_consistent_with_mol_values():
     # test that the primitive values are the same as the mol.distance etc.
     h2o2 = h2o2_mol()
@@ -626,7 +657,7 @@ def test_primitives_consistent_with_mol_values():
 
 
 # fmt: off
-dihedral_mols = [
+extra_mols = [
     Molecule(
         atoms=[
             Atom("C", 0.63365, 0.11934, -0.13163),
@@ -642,18 +673,27 @@ def test_primitives_consistent_with_mol_values():
             Atom("H", 1.28230, -0.63391, -0.54779),
             Atom("H", -1.08517, -1.07984, -0.05599),
         ]
-    )  # for testing dihedral derivatives over zero
+    ),  # for testing dihedral derivatives over zero
+    Molecule(
+        atoms=[
+            Atom("C", 0.35540, -0.20370, -0.44810),
+            Atom("C", -0.37180, 0.21470, 0.40200),
+            Atom("H", 1.01560, -0.60550, -1.23530),
+            Atom("H", -0.99920, 0.59450, 1.15720),
+        ]
+    )  # for testing linear angles
 ]
 
 test_mols = [
     h2o2_mol(), h2o2_mol(), water_mol(),
-    water_mol(), water_mol(), *dihedral_mols
+    water_mol(), water_mol(), *extra_mols
 ]
 test_prims = [
     PrimitiveDihedralAngle(2, 0, 1, 3), PrimitiveBondAngle(2, 0, 1),
     PrimitiveBondAngle(0, 1, 2), PrimitiveDistance(0, 1),
     PrimitiveInverseDistance(0, 1), PrimitiveDihedralAngle(2, 0, 1, 3),
-    PrimitiveDihedralAngle(2, 0, 1, 3)
+    PrimitiveDihedralAngle(2, 0, 1, 3),
+    PrimitiveLinearAngle(0, 1, 2, LinearBendType.BEND)
 ]
 # fmt: on
 

From 994e14acb158a9baca17d176c5fba2c1b795aeff Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Mon, 18 Dec 2023 16:17:54 +0000
Subject: [PATCH 08/47] unifinished update

---
 autode/opt/coordinates/internals.py | 85 +++++++++++++++++++++++------
 1 file changed, 68 insertions(+), 17 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 4c7838bd2..42d5dc21a 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -292,18 +292,6 @@ def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
     assert mol.graph is not None
     core_graph = copy.deepcopy(mol.graph)
 
-    # join disconnected graph components
-    if not core_graph.is_connected:
-        components = core_graph.get_components()
-        for comp_i, comp_j in itertools.combinations(components, r=2):
-            min_dist = float("inf")
-            min_pair = (-1, -1)
-            for i, j in itertools.product(list(comp_i), list(comp_j)):
-                if mol.distance(i, j) < min_dist:
-                    min_dist = mol.distance(i, j)
-                    min_pair = (i, j)
-            core_graph.add_edge(*min_pair, pi=False, active=False)
-
     # join hydrogen bonds
     h_bond_x = ["N", "O", "F", "P", "S", "Cl"]
     for i, j in itertools.combinations(range(mol.n_atoms), r=2):
@@ -318,6 +306,26 @@ def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
                 if not core_graph.has_edge(i, j):
                     core_graph.add_edge(i, j, pi=False, active=False)
 
+    # join disconnected graph components
+    if not core_graph.is_connected:
+        components = core_graph.get_components()
+        for comp_i, comp_j in itertools.combinations(components, r=2):
+            min_dist = float("inf")
+            min_pair = (-1, -1)
+            for i, j in itertools.product(list(comp_i), list(comp_j)):
+                if mol.distance(i, j) < min_dist:
+                    min_dist = mol.distance(i, j)
+                    min_pair = (i, j)
+            # avoid connecting distant components
+            if min_dist < 5.0:
+                core_graph.add_edge(*min_pair, pi=False, active=False)
+
+    if not core_graph.is_connected:
+        raise RuntimeError(
+            "Unable to join all the fragments, distance between "
+            "one or more pairs of fragments is too high (>5.0 Å)"
+        )
+
     # The constraints should be counted as bonds
     if mol.constraints.distance is not None:
         for i, j in mol.constraints.distance:
@@ -369,7 +377,10 @@ def _add_bonds_from_species(
 
 
 def _add_angles_from_species(
-    pic: AnyPIC, mol: "Species", core_graph: "MolecularGraph"
+    pic: AnyPIC,
+    mol: "Species",
+    core_graph: "MolecularGraph",
+    lin_thresh=Angle(175, "deg"),
 ) -> None:
     """
     Modify the set of primitives in-place by adding angles, from the
@@ -383,15 +394,25 @@ def _add_angles_from_species(
     for o in range(mol.n_atoms):
         for n, m in itertools.combinations(core_graph.neighbors(o), r=2):
             # avoid almost linear angles
-            if mol.angle(m, o, n) < Angle(175, "deg"):
+            if mol.angle(m, o, n) < lin_thresh:
                 pic.append(PrimitiveBondAngle(o=o, m=m, n=n))
             else:
+                # if central atom is connected to another, no need to include
+                other_neighbours = list(core_graph.neighbors(o))
+                other_neighbours.remove(m)
+                other_neighbours.remove(n)
+                if any(
+                    mol.angle(m, o, x) < lin_thresh for x in other_neighbours
+                ) or any(
+                    mol.angle(n, o, x) < lin_thresh for x in other_neighbours
+                ):
+                    continue
+
                 # stabilise linear angles by two orthogonal bends
                 pic.append(PrimitiveLinearAngle(m, o, n, LinearBendType.BEND))
                 pic.append(
                     PrimitiveLinearAngle(m, o, n, LinearBendType.COMPLEMENT)
                 )
-
     return None
 
 
@@ -399,6 +420,7 @@ def _add_dihedrals_from_species(
     pic: AnyPIC,
     mol: "Species",
     core_graph: "MolecularGraph",
+    lin_thresh=Angle(175, "deg"),
 ) -> None:
     """
     Modify the set of primitives in-place by adding dihedrals (torsions),
@@ -429,8 +451,8 @@ def _add_dihedrals_from_species(
                 if n == m:
                     continue
 
-                is_linear_1 = mol.angle(m, o, p) > Angle(175, "deg")
-                is_linear_2 = mol.angle(o, p, n) > Angle(175, "deg")
+                is_linear_1 = mol.angle(m, o, p) > lin_thresh
+                is_linear_2 = mol.angle(o, p, n) > lin_thresh
 
                 # if any angle is linear, don't add dihedral
                 if is_linear_1 or is_linear_2:
@@ -438,4 +460,33 @@ def _add_dihedrals_from_species(
                 else:
                     pic.append(PrimitiveDihedralAngle(m, o, p, n))
 
+    # find all linear atom chains A--B--C--D... and add dihedrals to terminal atoms
+    linear_chains = []
+    for b in range(mol.n_atoms):
+        for a, c in itertools.combinations(core_graph.neighbors(b), r=2):
+            if mol.angle(a, b, c) > lin_thresh:
+                linear_chains.append((a, b, c))
+
+    def concatenate_adjacent_chains(chains_list):
+        for chain1, chain2 in itertools.combinations(chains_list, r=2):
+            new_chain = None
+            if chain1[0] == chain2[0]:
+                new_chain = list(reversed(chain2)) + list(chain1)
+            elif chain1[0] == chain2[-1]:
+                new_chain = list(chain2) + list(chain1)
+            elif chain1[-1] == chain2[0]:
+                new_chain = list(chain1) + list(chain2)
+            elif chain1[-1] == chain2[-1]:
+                new_chain = list(chain1) + list(reversed(chain2))
+            else:
+                continue
+
+            if new_chain is not None:
+                chains_list.remove(chain1)
+                chains_list.remove(chain2)
+                chains_list.append(tuple(new_chain))
+                return concatenate_adjacent_chains(chains_list)
+
+    concatenate_adjacent_chains(linear_chains)
+    pass
     return None

From cf7de52df1f724c423e72f996b1078f6f225768f Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 18 Dec 2023 18:47:06 +0000
Subject: [PATCH 09/47] linear angle update

---
 autode/opt/coordinates/primitives.py | 24 +++++++++++++-----------
 1 file changed, 13 insertions(+), 11 deletions(-)

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index afe39aed5..8f9e08b70 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -446,16 +446,23 @@ def _init_axis(self, x: "CartesianCoordinates") -> None:
         # choose cartesian axis with the lowest overlap with m-n vector
         _m, _n = _x[self.m], _x[self.n]
         w = _m - _n
-        w /= np.linalg.norm(w)
         overlaps = []
         for axis in cart_axes:
             overlaps.append(np.dot(w, axis))
         cart_ax = cart_axes[np.argmin(np.abs(overlaps))]
 
-        # make the axis completely perpendicular to m-n vector
-        perp_axis = cart_ax - np.dot(cart_ax, w) * w
-        perp_axis /= np.linalg.norm(perp_axis)
-        self.axis_vec = DifferentiableVector3D(list(perp_axis))
+        # generate perpendicular axis by cross product
+        first_axis = np.cross(cart_ax, w)
+        if self.axis == LinearBendType.BEND:
+            axis_vec = first_axis / np.linalg.norm(first_axis)
+        elif self.axis == LinearBendType.COMPLEMENT:
+            second_axis = np.cross(first_axis, w)
+            axis_vec = second_axis / np.linalg.norm(second_axis)
+        else:
+            raise ValueError("Unknown axis type for linear bend")
+
+        self.axis_vec = DifferentiableVector3D(list(axis_vec))
+
         return None
 
     def _evaluate(
@@ -474,14 +481,9 @@ def _evaluate(
         w = w / w.norm()
         # TODO: does w need to be normalised?
 
-        cross_vec = w.cross(self.axis_vec)
-        # if complement is requested, perform another cross product
-        if self.axis == LinearBendType.COMPLEMENT:
-            cross_vec = w.cross(cross_vec)
-
         u = vec_m - vec_o
         v = vec_n - vec_o
-        return cross_vec.dot(u.cross(v)) / (u.norm() * v.norm())
+        return self.axis_vec.dot(u.cross(v)) / (u.norm() * v.norm())
 
     def __repr__(self):
         return f"LinearBend{self.axis}({self.m}-{self.o}-{self.n})"

From 30f1fc804d041547ae267f35ab9c6ca80ba4ae78 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 18 Dec 2023 18:49:12 +0000
Subject: [PATCH 10/47] unfinished update

---
 autode/opt/coordinates/primitives.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 8f9e08b70..d1befe515 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -477,9 +477,6 @@ def _evaluate(
         vec_m, vec_o, vec_n = _get_3d_vecs_from_atom_idxs(
             self.m, self.o, self.n, x=_x, deriv_order=deriv_order
         )
-        w = vec_m - vec_n
-        w = w / w.norm()
-        # TODO: does w need to be normalised?
 
         u = vec_m - vec_o
         v = vec_n - vec_o

From d02560280cda7bd8dd61b1a50218a100204cde49 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 18 Dec 2023 20:30:46 +0000
Subject: [PATCH 11/47] unfinished update

---
 autode/opt/coordinates/_autodiff.py  | 8 ++++----
 autode/opt/coordinates/primitives.py | 8 +++++---
 tests/test_opt/test_coordiantes.py   | 8 ++++----
 3 files changed, 13 insertions(+), 11 deletions(-)

diff --git a/autode/opt/coordinates/_autodiff.py b/autode/opt/coordinates/_autodiff.py
index 5f635773e..0855a441e 100644
--- a/autode/opt/coordinates/_autodiff.py
+++ b/autode/opt/coordinates/_autodiff.py
@@ -604,7 +604,9 @@ def _check_same_type(other) -> None:
             raise ValueError("Operation must be done with another 3D vector!")
         return None
 
-    def dot(self, other: "DifferentiableVector3D") -> "VectorHyperDual":
+    def dot(
+        self, other: "DifferentiableVector3D"
+    ) -> Union["VectorHyperDual", numeric_type]:
         """
         Dot product of two 3D vectors
 
@@ -618,10 +620,9 @@ def dot(self, other: "DifferentiableVector3D") -> "VectorHyperDual":
         dot: Union[VectorHyperDual, numeric_type] = 0
         for k in range(3):
             dot = dot + self._data[k] * other._data[k]
-        assert isinstance(dot, VectorHyperDual)
         return dot
 
-    def norm(self) -> "VectorHyperDual":
+    def norm(self) -> Union["VectorHyperDual", numeric_type]:
         """
         Euclidean (l2) norm of this 3D vector
 
@@ -631,7 +632,6 @@ def norm(self) -> "VectorHyperDual":
         norm = DifferentiableMath.sqrt(
             self._data[0] ** 2 + self._data[1] ** 2 + self._data[2] ** 2
         )
-        assert isinstance(norm, VectorHyperDual)
         return norm
 
     def __add__(
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index d1befe515..9c2ce9429 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -239,7 +239,7 @@ def _evaluate(
         vec_i, vec_j = _get_3d_vecs_from_atom_idxs(
             self.i, self.j, x=x, deriv_order=deriv_order
         )
-        return 1.0 / (vec_i - vec_j).norm()
+        return 1.0 / (vec_i - vec_j).norm()  # type: ignore
 
     def __repr__(self):
         return f"InverseDistance({self.i}-{self.j})"
@@ -261,7 +261,7 @@ def _evaluate(
         vec_i, vec_j = _get_3d_vecs_from_atom_idxs(
             self.i, self.j, x=x, deriv_order=deriv_order
         )
-        return (vec_i - vec_j).norm()
+        return (vec_i - vec_j).norm()  # type: ignore
 
     def __repr__(self):
         return f"Distance({self.i}-{self.j})"
@@ -446,6 +446,7 @@ def _init_axis(self, x: "CartesianCoordinates") -> None:
         # choose cartesian axis with the lowest overlap with m-n vector
         _m, _n = _x[self.m], _x[self.n]
         w = _m - _n
+        w /= np.linalg.norm(w)
         overlaps = []
         for axis in cart_axes:
             overlaps.append(np.dot(w, axis))
@@ -483,4 +484,5 @@ def _evaluate(
         return self.axis_vec.dot(u.cross(v)) / (u.norm() * v.norm())
 
     def __repr__(self):
-        return f"LinearBend{self.axis}({self.m}-{self.o}-{self.n})"
+        axis_str = "B" if self.axis == LinearBendType.BEND else "C"
+        return f"LinearBend{axis_str}({self.m}-{self.o}-{self.n})"
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 1bb01b45e..019852a7c 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -623,7 +623,7 @@ def test_linear_angle():
         ]
     )
     x = CartesianCoordinates(acetylene.coordinates)
-    angle = PrimitiveLinearAngle(0, 1, 2, LinearBendType.BEND)
+    angle = PrimitiveLinearAngle(0, 1, 3, LinearBendType.BEND)
     assert angle.axis_vec is None
     angle._init_axis(x)
     assert angle.axis_vec is not None
@@ -632,12 +632,12 @@ def test_linear_angle():
     assert np.dot(axis_vec, m_n_vec) < 0.001
 
     # the atom order can be reversed without any change in the values
-    angle2 = PrimitiveLinearAngle(2, 1, 0, LinearBendType.BEND)
+    angle2 = PrimitiveLinearAngle(3, 1, 0, LinearBendType.BEND)
     assert np.isclose(angle(x), angle2(x))
     assert angle == angle2
 
     # however, linear bend complement would not be the same
-    angle3 = PrimitiveLinearAngle(0, 1, 2, LinearBendType.COMPLEMENT)
+    angle3 = PrimitiveLinearAngle(0, 1, 3, LinearBendType.COMPLEMENT)
     assert angle != angle3
     assert not np.isclose(angle(x), angle3(x), rtol=1e-3)
 
@@ -693,7 +693,7 @@ def test_primitives_consistent_with_mol_values():
     PrimitiveBondAngle(0, 1, 2), PrimitiveDistance(0, 1),
     PrimitiveInverseDistance(0, 1), PrimitiveDihedralAngle(2, 0, 1, 3),
     PrimitiveDihedralAngle(2, 0, 1, 3),
-    PrimitiveLinearAngle(0, 1, 2, LinearBendType.BEND)
+    PrimitiveLinearAngle(0, 1, 3, LinearBendType.BEND)
 ]
 # fmt: on
 

From 865ba19d37790b4a7bb789be3cf4115d5e1c7a36 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Tue, 19 Dec 2023 15:56:44 +0000
Subject: [PATCH 12/47] unifinished update

---
 autode/opt/coordinates/internals.py | 59 ++++++++++++++---------------
 1 file changed, 28 insertions(+), 31 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 42d5dc21a..da48c4dfa 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -435,31 +435,6 @@ def _add_dihedrals_from_species(
     if mol.n_atoms < 4:
         return
 
-    for o, p in core_graph.edges:
-        for m in core_graph.neighbors(o):
-            if m == p:
-                continue
-
-            if mol.angle(m, o, p) > Angle(175, "deg"):
-                continue
-
-            for n in core_graph.neighbors(p):
-                if n == o:
-                    continue
-
-                # avoid triangle rings like cyclopropane
-                if n == m:
-                    continue
-
-                is_linear_1 = mol.angle(m, o, p) > lin_thresh
-                is_linear_2 = mol.angle(o, p, n) > lin_thresh
-
-                # if any angle is linear, don't add dihedral
-                if is_linear_1 or is_linear_2:
-                    continue
-                else:
-                    pic.append(PrimitiveDihedralAngle(m, o, p, n))
-
     # find all linear atom chains A--B--C--D... and add dihedrals to terminal atoms
     linear_chains = []
     for b in range(mol.n_atoms):
@@ -469,7 +444,6 @@ def _add_dihedrals_from_species(
 
     def concatenate_adjacent_chains(chains_list):
         for chain1, chain2 in itertools.combinations(chains_list, r=2):
-            new_chain = None
             if chain1[0] == chain2[0]:
                 new_chain = list(reversed(chain2)) + list(chain1)
             elif chain1[0] == chain2[-1]:
@@ -481,12 +455,35 @@ def concatenate_adjacent_chains(chains_list):
             else:
                 continue
 
-            if new_chain is not None:
-                chains_list.remove(chain1)
-                chains_list.remove(chain2)
-                chains_list.append(tuple(new_chain))
-                return concatenate_adjacent_chains(chains_list)
+            chains_list.remove(chain1)
+            chains_list.remove(chain2)
+            chains_list.append(tuple(new_chain))
+            return concatenate_adjacent_chains(chains_list)
 
     concatenate_adjacent_chains(linear_chains)
+    terminal_points = [(chain[0], chain[-1]) for chain in linear_chains]
+
+    # add normal + linear chain dihedrals
+    for o, p in list(core_graph.edges) + terminal_points:
+        for m in core_graph.neighbors(o):
+            if m == p:
+                continue
+
+            for n in core_graph.neighbors(p):
+                if n == o:
+                    continue
+
+                # avoid triangle rings like cyclopropane
+                if n == m:
+                    continue
+
+                is_linear_1 = mol.angle(m, o, p) > lin_thresh
+                is_linear_2 = mol.angle(o, p, n) > lin_thresh
+
+                # if any angle is linear, don't add dihedral
+                if is_linear_1 or is_linear_2:
+                    continue
+                else:
+                    pic.append(PrimitiveDihedralAngle(m, o, p, n))
     pass
     return None

From f72d68edc33c2a02fa43923d953ae4120d66ad10 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Wed, 20 Dec 2023 19:30:42 +0000
Subject: [PATCH 13/47] unfinished update

---
 autode/opt/coordinates/internals.py  |  6 +-
 autode/opt/coordinates/primitives.py | 88 ++++++++++++++++++++++------
 tests/test_opt/test_coordiantes.py   | 12 ++--
 3 files changed, 77 insertions(+), 29 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index da48c4dfa..cc6f75649 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -407,12 +407,8 @@ def _add_angles_from_species(
                     mol.angle(n, o, x) < lin_thresh for x in other_neighbours
                 ):
                     continue
-
+                # TODO: use other atoms for linear bend instead
                 # stabilise linear angles by two orthogonal bends
-                pic.append(PrimitiveLinearAngle(m, o, n, LinearBendType.BEND))
-                pic.append(
-                    PrimitiveLinearAngle(m, o, n, LinearBendType.COMPLEMENT)
-                )
     return None
 
 
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 9c2ce9429..25c82b504 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -415,17 +415,17 @@ class LinearBendType(Enum):
     COMPLEMENT = 1
 
 
-class PrimitiveLinearAngle(Primitive):
-    def __init__(self, m: int, o: int, n: int, axis: LinearBendType):
+class AbstractLinearAngle(Primitive, ABC):
+    def __init__(self, m: int, o: int, n: int, r: int, axis: LinearBendType):
         """Linear Bend: m-o-n"""
-        super().__init__(m, o, n)
+        super().__init__(m, o, n, r)
         self.m = int(m)
         self.o = int(o)
         self.n = int(n)
+        self.r = int(r)
 
         assert isinstance(axis, LinearBendType)
         self.axis = axis
-        self.axis_vec: Optional[DifferentiableVector3D] = None
 
     def __eq__(self, other):
         return isinstance(other, self.__class__) and (
@@ -435,7 +435,71 @@ def __eq__(self, other):
         )
         # TODO: check the sign of the bends if m and n swapped
 
-    def _init_axis(self, x: "CartesianCoordinates") -> None:
+    def _calc_linear_bend(
+        self,
+        m_vec: DifferentiableVector3D,
+        o_vec: DifferentiableVector3D,
+        n_vec: DifferentiableVector3D,
+        r_vec: DifferentiableVector3D,
+    ):
+        """
+        Evaluate the linear bend from the vector positions of the
+        atoms involved in the angle m, o, n, and the reference
+        atom (or dummy atom) r
+
+        Args:
+            m_vec:
+            o_vec:
+            n_vec:
+            r_vec:
+
+        Returns:
+            (VectorHyperDual): The value, optionally containing derivatives
+        """
+        # As defined in J. Comput. Chem., 20(10), 1999, 1067
+        o_m = m_vec - o_vec
+        o_n = n_vec - o_vec
+        o_r = r_vec - o_vec
+
+        # eq.(44) p 1073
+        u = o_m.cross(o_r)
+        u = u / u.norm()
+
+        # eq. (46) and (47) p 1074
+        if self.axis == LinearBendType.BEND:
+            return u.dot(o_n) / o_n.norm()
+        elif self.axis == LinearBendType.COMPLEMENT:
+            return u.dot(o_n.cross(o_m)) / (o_n.norm() * o_m.norm())
+
+
+class PrimitiveLinearAngle(AbstractLinearAngle):
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
+    ):
+        """Linear Bend angle m-o-n against a Cartesian axis"""
+
+        _x = x.ravel()
+        vec_m, vec_o, vec_n, vec_r = _get_3d_vecs_from_atom_idxs(
+            self.m, self.o, self.n, self.r, x=_x, deriv_order=deriv_order
+        )
+
+        return self._calc_linear_bend(vec_m, vec_o, vec_n, vec_r)
+
+    def __repr__(self):
+        axis_str = "B" if self.axis == LinearBendType.BEND else "C"
+        return f"LinearBend{axis_str}({self.m}-{self.o}-{self.n}, {self.r})"
+
+
+class PrimitiveDummyLinearAngle(AbstractLinearAngle):
+    """Linear bend with a dummy atom"""
+
+    def __init__(self, m: int, o: int, n: int, axis: LinearBendType):
+        super().__init__(m, o, n, -1, axis)
+
+        self._dummy: Optional[DifferentiableVector3D] = None
+
+    def _init_dummy_atom(self, x: "CartesianCoordinates") -> None:
+        """Create the dummy atom"""
         _x = x.reshape(-1, 3)
         cart_axes = [
             np.array([1.0, 0.0, 0.0]),
@@ -469,19 +533,7 @@ def _init_axis(self, x: "CartesianCoordinates") -> None:
     def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
     ):
-        """Linear Bend angle m-o-n against a Cartesian axis"""
-        if self.axis_vec is None:
-            self._init_axis(x)
-
-        assert self.axis_vec is not None
-        _x = x.ravel()
-        vec_m, vec_o, vec_n = _get_3d_vecs_from_atom_idxs(
-            self.m, self.o, self.n, x=_x, deriv_order=deriv_order
-        )
-
-        u = vec_m - vec_o
-        v = vec_n - vec_o
-        return self.axis_vec.dot(u.cross(v)) / (u.norm() * v.norm())
+        pass
 
     def __repr__(self):
         axis_str = "B" if self.axis == LinearBendType.BEND else "C"
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 019852a7c..1afabfb11 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -15,7 +15,7 @@
     PrimitiveBondAngle,
     ConstrainedPrimitiveBondAngle,
     PrimitiveDihedralAngle,
-    PrimitiveLinearAngle,
+    PrimitiveDummyLinearAngle,
     LinearBendType,
 )
 
@@ -623,21 +623,21 @@ def test_linear_angle():
         ]
     )
     x = CartesianCoordinates(acetylene.coordinates)
-    angle = PrimitiveLinearAngle(0, 1, 3, LinearBendType.BEND)
+    angle = PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.BEND)
     assert angle.axis_vec is None
-    angle._init_axis(x)
+    # TODO: check position of dummy atom
     assert angle.axis_vec is not None
     axis_vec = np.array(angle.axis_vec._data)
     m_n_vec = acetylene.coordinates[0] - acetylene.coordinates[2]
     assert np.dot(axis_vec, m_n_vec) < 0.001
 
     # the atom order can be reversed without any change in the values
-    angle2 = PrimitiveLinearAngle(3, 1, 0, LinearBendType.BEND)
+    angle2 = PrimitiveDummyLinearAngle(3, 1, 0, LinearBendType.BEND)
     assert np.isclose(angle(x), angle2(x))
     assert angle == angle2
 
     # however, linear bend complement would not be the same
-    angle3 = PrimitiveLinearAngle(0, 1, 3, LinearBendType.COMPLEMENT)
+    angle3 = PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.COMPLEMENT)
     assert angle != angle3
     assert not np.isclose(angle(x), angle3(x), rtol=1e-3)
 
@@ -693,7 +693,7 @@ def test_primitives_consistent_with_mol_values():
     PrimitiveBondAngle(0, 1, 2), PrimitiveDistance(0, 1),
     PrimitiveInverseDistance(0, 1), PrimitiveDihedralAngle(2, 0, 1, 3),
     PrimitiveDihedralAngle(2, 0, 1, 3),
-    PrimitiveLinearAngle(0, 1, 3, LinearBendType.BEND)
+    PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.BEND)
 ]
 # fmt: on
 

From e5859b56f78f499cabd76a841afa5aba47651561 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Wed, 20 Dec 2023 21:22:47 +0000
Subject: [PATCH 14/47] better linear angles

---
 autode/opt/coordinates/primitives.py | 53 +++++++++++++++-------------
 1 file changed, 29 insertions(+), 24 deletions(-)

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 25c82b504..ef2f580e2 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -415,7 +415,7 @@ class LinearBendType(Enum):
     COMPLEMENT = 1
 
 
-class AbstractLinearAngle(Primitive, ABC):
+class LinearAngleBase(Primitive, ABC):
     def __init__(self, m: int, o: int, n: int, r: int, axis: LinearBendType):
         """Linear Bend: m-o-n"""
         super().__init__(m, o, n, r)
@@ -470,13 +470,17 @@ def _calc_linear_bend(
             return u.dot(o_n) / o_n.norm()
         elif self.axis == LinearBendType.COMPLEMENT:
             return u.dot(o_n.cross(o_m)) / (o_n.norm() * o_m.norm())
+        else:
+            raise ValueError("Unknown axis for linear bend")
+
 
+class PrimitiveLinearAngle(LinearAngleBase):
+    """Linear Angle w.r.t. a reference atom"""
 
-class PrimitiveLinearAngle(AbstractLinearAngle):
     def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
     ):
-        """Linear Bend angle m-o-n against a Cartesian axis"""
+        """Linear Bend angle m-o-n against reference atom r"""
 
         _x = x.ravel()
         vec_m, vec_o, vec_n, vec_r = _get_3d_vecs_from_atom_idxs(
@@ -490,16 +494,18 @@ def __repr__(self):
         return f"LinearBend{axis_str}({self.m}-{self.o}-{self.n}, {self.r})"
 
 
-class PrimitiveDummyLinearAngle(AbstractLinearAngle):
+class PrimitiveDummyLinearAngle(LinearAngleBase):
     """Linear bend with a dummy atom"""
 
     def __init__(self, m: int, o: int, n: int, axis: LinearBendType):
         super().__init__(m, o, n, -1, axis)
 
-        self._dummy: Optional[DifferentiableVector3D] = None
+        self._vec_r: Optional[DifferentiableVector3D] = None
 
-    def _init_dummy_atom(self, x: "CartesianCoordinates") -> None:
-        """Create the dummy atom"""
+    def _get_dummy_atom(
+        self, x: "CartesianCoordinates"
+    ) -> DifferentiableVector3D:
+        """Create the dummy atom r"""
         _x = x.reshape(-1, 3)
         cart_axes = [
             np.array([1.0, 0.0, 0.0]),
@@ -507,33 +513,32 @@ def _init_dummy_atom(self, x: "CartesianCoordinates") -> None:
             np.array([0.0, 0.0, 1.0]),
         ]
 
-        # choose cartesian axis with the lowest overlap with m-n vector
-        _m, _n = _x[self.m], _x[self.n]
-        w = _m - _n
+        # choose cartesian axis with the lowest overlap with m-o vector
+        w = _x[self.m] - _x[self.o]
         w /= np.linalg.norm(w)
         overlaps = []
         for axis in cart_axes:
             overlaps.append(np.dot(w, axis))
         cart_ax = cart_axes[np.argmin(np.abs(overlaps))]
 
-        # generate perpendicular axis by cross product
-        first_axis = np.cross(cart_ax, w)
-        if self.axis == LinearBendType.BEND:
-            axis_vec = first_axis / np.linalg.norm(first_axis)
-        elif self.axis == LinearBendType.COMPLEMENT:
-            second_axis = np.cross(first_axis, w)
-            axis_vec = second_axis / np.linalg.norm(second_axis)
-        else:
-            raise ValueError("Unknown axis type for linear bend")
-
-        self.axis_vec = DifferentiableVector3D(list(axis_vec))
-
-        return None
+        # place dummy atom perpendicular to m-o bond 1 A away
+        perp_axis = np.cross(cart_ax, w)
+        _o = _x[self.o]
+        dummy_point = _o + perp_axis / np.linalg.norm(perp_axis)
+        return DifferentiableVector3D(list(dummy_point))
 
     def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
     ):
-        pass
+        if self._vec_r is None:
+            self._vec_r = self._get_dummy_atom(x)
+
+        _x = x.ravel()
+        vec_m, vec_o, vec_n = _get_3d_vecs_from_atom_idxs(
+            self.m, self.o, self.n, x=_x, deriv_order=deriv_order
+        )
+
+        return self._calc_linear_bend(vec_m, vec_o, vec_n, self._vec_r)
 
     def __repr__(self):
         axis_str = "B" if self.axis == LinearBendType.BEND else "C"

From 0fc99bb3ea16af0d54fcd2a646c6d61376efb129 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Wed, 20 Dec 2023 22:12:59 +0000
Subject: [PATCH 15/47] better linear angles

---
 autode/opt/coordinates/internals.py  | 55 ++++++++++++++++++++++++++--
 autode/opt/coordinates/primitives.py |  2 +-
 2 files changed, 52 insertions(+), 5 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index cc6f75649..86ffca2d3 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -16,7 +16,7 @@
 import itertools
 from typing import Any, Optional, Type, List, TYPE_CHECKING
 from abc import ABC, abstractmethod
-from autode.values import Angle
+from autode.values import Angle, Distance
 from autode.opt.coordinates.base import OptCoordinates
 from autode.opt.coordinates.primitives import (
     PrimitiveInverseDistance,
@@ -24,6 +24,7 @@
     PrimitiveDistance,
     ConstrainedPrimitiveDistance,
     PrimitiveBondAngle,
+    PrimitiveDummyLinearAngle,
     PrimitiveLinearAngle,
     PrimitiveDihedralAngle,
     LinearBendType,
@@ -391,6 +392,33 @@ def _add_angles_from_species(
         mol: The species object
         core_graph: The connectivity graph
     """
+
+    def get_ref_atom(a, b, c):
+        """get a reference atom for a-b-c linear angle"""
+        # all atoms in 4 A radius
+        near_atoms = [
+            idx
+            for idx in range(mol.n_atoms)
+            if mol.distance(b, idx) < Distance(4.0, "ang")
+        ]
+        deviations_from_90 = []
+        for atom in near_atoms:
+            i_b_a = mol.angle(atom, b, a)
+            if i_b_a > lin_thresh or i_b_a < 180 - lin_thresh:
+                continue
+            i_b_c = mol.angle(atom, b, c)
+            if i_b_c > lin_thresh or i_b_c < 180 - lin_thresh:
+                continue
+            deviation_a = abs(i_b_a - Angle(90, "deg"))
+            deviation_c = abs(i_b_c - Angle(90, "deg"))
+            avg_dev = (deviation_a + deviation_c) / 2
+            deviations_from_90.append(avg_dev)
+
+        if len(deviations_from_90) == 0:
+            return None
+
+        return near_atoms[np.argmin(deviations_from_90)]
+
     for o in range(mol.n_atoms):
         for n, m in itertools.combinations(core_graph.neighbors(o), r=2):
             # avoid almost linear angles
@@ -407,8 +435,28 @@ def _add_angles_from_species(
                     mol.angle(n, o, x) < lin_thresh for x in other_neighbours
                 ):
                     continue
-                # TODO: use other atoms for linear bend instead
-                # stabilise linear angles by two orthogonal bends
+
+                # for linear bends, ideally a reference atom is needed
+                r = get_ref_atom(m, o, n)
+                if r is not None:
+                    pic.append(
+                        PrimitiveLinearAngle(m, o, n, r, LinearBendType.BEND)
+                    )
+                    pic.append(
+                        PrimitiveLinearAngle(
+                            m, o, n, r, LinearBendType.COMPLEMENT
+                        )
+                    )
+                else:  # these use dummy atom for reference
+                    pic.append(
+                        PrimitiveDummyLinearAngle(m, o, n, LinearBendType.BEND)
+                    )
+                    pic.append(
+                        PrimitiveDummyLinearAngle(
+                            m, o, n, LinearBendType.COMPLEMENT
+                        )
+                    )
+
     return None
 
 
@@ -481,5 +529,4 @@ def concatenate_adjacent_chains(chains_list):
                     continue
                 else:
                     pic.append(PrimitiveDihedralAngle(m, o, p, n))
-    pass
     return None
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index ef2f580e2..44f1cc6d8 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -542,4 +542,4 @@ def _evaluate(
 
     def __repr__(self):
         axis_str = "B" if self.axis == LinearBendType.BEND else "C"
-        return f"LinearBend{axis_str}({self.m}-{self.o}-{self.n})"
+        return f"LinearBend{axis_str}({self.m}-{self.o}-{self.n}, D)"

From d561f9b25e039c9e1b9c3b1d4a3fd1e7de9a0983 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Wed, 20 Dec 2023 22:14:19 +0000
Subject: [PATCH 16/47] minor change

---
 autode/opt/coordinates/internals.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 86ffca2d3..a9720936c 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -452,9 +452,7 @@ def get_ref_atom(a, b, c):
                         PrimitiveDummyLinearAngle(m, o, n, LinearBendType.BEND)
                     )
                     pic.append(
-                        PrimitiveDummyLinearAngle(
-                            m, o, n, LinearBendType.COMPLEMENT
-                        )
+                        PrimitiveDummyLinearAngle(m, o, n, LinearBendType.COMPLEMENT)  # fmt: skip
                     )
 
     return None

From d2a768b867f7ce1cbea0f15861ccf0a2bf240ffc Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 21 Dec 2023 11:06:49 +0000
Subject: [PATCH 17/47] fix linear bend code

---
 autode/opt/coordinates/internals.py | 12 ++++++++----
 autode/opt/optimisers/crfo.py       | 28 ----------------------------
 2 files changed, 8 insertions(+), 32 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index a9720936c..09a45cb86 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -381,7 +381,7 @@ def _add_angles_from_species(
     pic: AnyPIC,
     mol: "Species",
     core_graph: "MolecularGraph",
-    lin_thresh=Angle(175, "deg"),
+    lin_thresh=Angle(170, "deg"),
 ) -> None:
     """
     Modify the set of primitives in-place by adding angles, from the
@@ -395,12 +395,14 @@ def _add_angles_from_species(
 
     def get_ref_atom(a, b, c):
         """get a reference atom for a-b-c linear angle"""
-        # all atoms in 4 A radius
+        # all atoms in 4 A radius except a, b, c
         near_atoms = [
             idx
             for idx in range(mol.n_atoms)
             if mol.distance(b, idx) < Distance(4.0, "ang")
+            and idx not in (a, b, c)
         ]
+        # get atoms closest to perpendicular
         deviations_from_90 = []
         for atom in near_atoms:
             i_b_a = mol.angle(atom, b, a)
@@ -452,7 +454,9 @@ def get_ref_atom(a, b, c):
                         PrimitiveDummyLinearAngle(m, o, n, LinearBendType.BEND)
                     )
                     pic.append(
-                        PrimitiveDummyLinearAngle(m, o, n, LinearBendType.COMPLEMENT)  # fmt: skip
+                        PrimitiveDummyLinearAngle(
+                            m, o, n, LinearBendType.COMPLEMENT
+                        )
                     )
 
     return None
@@ -462,7 +466,7 @@ def _add_dihedrals_from_species(
     pic: AnyPIC,
     mol: "Species",
     core_graph: "MolecularGraph",
-    lin_thresh=Angle(175, "deg"),
+    lin_thresh=Angle(170, "deg"),
 ) -> None:
     """
     Modify the set of primitives in-place by adding dihedrals (torsions),
diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py
index d5cc96c2a..a7eb56c19 100644
--- a/autode/opt/optimisers/crfo.py
+++ b/autode/opt/optimisers/crfo.py
@@ -173,31 +173,3 @@ def _lambda_n_from_eigvals_and_gradient(
 
         eigenvalues = np.linalg.eigvalsh(aug_h)
         return eigenvalues[0]
-
-
-def _dihedrals(species):
-    """
-    Iterator over the dihedrals in a species. Skipping those that contain
-    bond angles close to 180 degrees (tolerance <179)
-    """
-
-    for o, p in species.graph.edges:
-        for m in species.graph.neighbors(o):
-            if m == p:
-                continue
-
-            if np.isclose(species.angle(m, o, p), Angle(np.pi), atol=0.04):
-                continue  # Don't add potentially ill-defined dihedrals
-
-            for n in species.graph.neighbors(p):
-                if n == o:
-                    continue
-
-                # avoid triangular rings like cyclopropane
-                if m == n:
-                    continue
-
-                if np.isclose(species.angle(o, p, n), Angle(np.pi), atol=0.04):
-                    continue
-
-                yield PrimitiveDihedralAngle(m, o, p, n)

From 5f97251dd409cdbf9832542cf05f9d21c207089d Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 21 Dec 2023 18:03:54 +0000
Subject: [PATCH 18/47] fix some tests

---
 autode/opt/coordinates/primitives.py |  7 +++---
 autode/opt/optimisers/crfo.py        | 15 ++---------
 tests/test_opt/molecules.py          | 11 +++++++++
 tests/test_opt/test_coordiantes.py   | 37 ++++++++++++++++------------
 4 files changed, 37 insertions(+), 33 deletions(-)

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 44f1cc6d8..5eedc7647 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -377,7 +377,7 @@ def __init__(self, m: int, o: int, p: int, n: int):
         self.n = int(n)
 
     def __eq__(self, other) -> bool:
-        """Equality of two distance functions"""
+        """Equality of two dihedral angles"""
         return isinstance(other, self.__class__) and (
             self._atom_indexes == other._atom_indexes
             or self._atom_indexes == tuple(reversed(other._atom_indexes))
@@ -428,12 +428,11 @@ def __init__(self, m: int, o: int, n: int, r: int, axis: LinearBendType):
         self.axis = axis
 
     def __eq__(self, other):
+        """Equality of two linear bend angles"""
         return isinstance(other, self.__class__) and (
-            self._ordered_idxs == other._ordered_idxs
-            and self.o == other.o
+            self._atom_indexes == other._atom_indexes
             and self.axis == other.axis
         )
-        # TODO: check the sign of the bends if m and n swapped
 
     def _calc_linear_bend(
         self,
diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py
index a7eb56c19..e1973f396 100644
--- a/autode/opt/optimisers/crfo.py
+++ b/autode/opt/optimisers/crfo.py
@@ -8,23 +8,12 @@
 import numpy as np
 from typing import Union
 
-from itertools import combinations
 from autode.log import logger
-from autode.values import GradientRMS, Angle, Distance
+from autode.values import GradientRMS, Distance
 from autode.opt.coordinates import CartesianCoordinates, DICWithConstraints
-from autode.opt.coordinates.internals import (
-    PIC,
-    AnyPIC,
-    build_pic_from_species,
-)
+from autode.opt.coordinates.internals import build_pic_from_species
 from autode.opt.optimisers.rfo import RFOptimiser
 from autode.opt.optimisers.hessian_update import BFGSDampedUpdate, NullUpdate
-from autode.opt.coordinates.primitives import (
-    PrimitiveDistance,
-    PrimitiveBondAngle,
-    PrimitiveDihedralAngle,
-    ConstrainedPrimitiveDistance,
-)
 
 
 class CRFOptimiser(RFOptimiser):
diff --git a/tests/test_opt/molecules.py b/tests/test_opt/molecules.py
index 92936dfcc..b0f29d263 100644
--- a/tests/test_opt/molecules.py
+++ b/tests/test_opt/molecules.py
@@ -6,6 +6,17 @@ def h_atom():
     return Molecule(atoms=[Atom("H")], mult=2)
 
 
+def acetylene_mol():
+    return Molecule(
+        atoms=[
+            Atom("C", 0.3800, -0.2049, -0.4861),
+            Atom("C", -0.3727, 0.1836, 0.4744),
+            Atom("H", 1.0156, -0.6055, -1.2353),
+            Atom("H", -0.9992, 0.5945, 1.1572),
+        ]
+    )
+
+
 def methane_mol():
     return Molecule(
         atoms=[
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 1afabfb11..3fdfeaa9e 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -624,22 +624,27 @@ def test_linear_angle():
     )
     x = CartesianCoordinates(acetylene.coordinates)
     angle = PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.BEND)
-    assert angle.axis_vec is None
-    # TODO: check position of dummy atom
-    assert angle.axis_vec is not None
-    axis_vec = np.array(angle.axis_vec._data)
-    m_n_vec = acetylene.coordinates[0] - acetylene.coordinates[2]
-    assert np.dot(axis_vec, m_n_vec) < 0.001
-
-    # the atom order can be reversed without any change in the values
-    angle2 = PrimitiveDummyLinearAngle(3, 1, 0, LinearBendType.BEND)
-    assert np.isclose(angle(x), angle2(x))
-    assert angle == angle2
-
-    # however, linear bend complement would not be the same
-    angle3 = PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.COMPLEMENT)
-    assert angle != angle3
-    assert not np.isclose(angle(x), angle3(x), rtol=1e-3)
+    assert angle._vec_r is None
+    _ = angle(x)
+    assert angle._vec_r is not None
+    old_r_vec = angle._vec_r
+    # the dummy atom should not change after the first call
+    _ = angle(x)
+    _ = angle(x)
+    assert angle._vec_r is old_r_vec
+
+    axis_vec = np.array(np.array(angle._vec_r._data) - x.reshape(-1, 3)[1])
+    m_n_vec = acetylene.coordinates[0] - acetylene.coordinates[1]
+    assert abs(np.dot(axis_vec, m_n_vec)) < 0.001
+
+    # check that the linear bond complement does not have the same value
+    angle2 = PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.COMPLEMENT)
+    assert angle != angle2
+    assert not np.isclose(angle(x), angle2(x), rtol=1e-3)
+
+    # for linear angle, swapping the end points changes the definition
+    angle3 = PrimitiveDummyLinearAngle(3, 1, 0, LinearBendType.BEND)
+    assert angle3 != angle
 
 
 def test_primitives_consistent_with_mol_values():

From b176e3a10cec7b9086c310e813ba9021c64c00b3 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 21 Dec 2023 18:11:19 +0000
Subject: [PATCH 19/47] add tests

---
 tests/test_opt/molecules.py        | 18 ++++++++++++++++++
 tests/test_opt/test_coordiantes.py |  9 ++++++---
 2 files changed, 24 insertions(+), 3 deletions(-)

diff --git a/tests/test_opt/molecules.py b/tests/test_opt/molecules.py
index b0f29d263..145d94066 100644
--- a/tests/test_opt/molecules.py
+++ b/tests/test_opt/molecules.py
@@ -52,3 +52,21 @@ def h2o2_mol():
             Atom("H", 0.58605, 0.91107, 0.59006),
         ]
     )
+
+
+def feco5_mol():
+    return Molecule(
+        atoms=[
+            Atom("O", -1.5139, -3.5069, -0.3015),
+            Atom("C", -0.8455, -2.5982, -0.3019),
+            Atom("Fe", 0.3545, -0.9664, -0.3020),
+            Atom("C", 1.5539, 0.6660, -0.3004),
+            Atom("O", 2.2216, 1.5751, -0.2992),
+            Atom("C", 1.2470, -1.6232, 1.3934),
+            Atom("O", 1.7445, -1.9892, 2.3373),
+            Atom("C", 1.0930, -1.5083, -2.1084),
+            Atom("O", 1.5042, -1.8102, -3.1145),
+            Atom("C", -1.2751, 0.2315, -0.1941),
+            Atom("O", -2.1828, 0.8986, -0.1345),
+        ]
+    )
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 3fdfeaa9e..f32edac94 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -1,6 +1,6 @@
 import pytest
 import numpy as np
-from .molecules import h2, methane_mol, water_mol, h2o2_mol
+from .molecules import h2, methane_mol, water_mol, h2o2_mol, feco5_mol
 from autode.atoms import Atom
 from autode.species.molecule import Molecule
 from autode.values import Angle
@@ -15,6 +15,7 @@
     PrimitiveBondAngle,
     ConstrainedPrimitiveBondAngle,
     PrimitiveDihedralAngle,
+    PrimitiveLinearAngle,
     PrimitiveDummyLinearAngle,
     LinearBendType,
 )
@@ -686,7 +687,8 @@ def test_primitives_consistent_with_mol_values():
             Atom("H", 1.01560, -0.60550, -1.23530),
             Atom("H", -0.99920, 0.59450, 1.15720),
         ]
-    )  # for testing linear angles
+    ),
+    feco5_mol(), # for testing linear angles
 ]
 
 test_mols = [
@@ -698,7 +700,8 @@ def test_primitives_consistent_with_mol_values():
     PrimitiveBondAngle(0, 1, 2), PrimitiveDistance(0, 1),
     PrimitiveInverseDistance(0, 1), PrimitiveDihedralAngle(2, 0, 1, 3),
     PrimitiveDihedralAngle(2, 0, 1, 3),
-    PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.BEND)
+    PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.BEND),
+    PrimitiveLinearAngle(2, 3, 4, 8, LinearBendType.BEND),
 ]
 # fmt: on
 

From 0757cda5ca789401a5f0cfe9dc4500934a600d27 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 21 Dec 2023 19:50:18 +0000
Subject: [PATCH 20/47] fix failing test

---
 tests/test_opt/test_coordiantes.py | 19 ++++++++++++++++++-
 tests/test_opt/test_crfo.py        |  3 ++-
 2 files changed, 20 insertions(+), 2 deletions(-)

diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index f32edac94..1c80f6e94 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -1,3 +1,5 @@
+import itertools
+
 import pytest
 import numpy as np
 from .molecules import h2, methane_mol, water_mol, h2o2_mol, feco5_mol
@@ -5,7 +7,11 @@
 from autode.species.molecule import Molecule
 from autode.values import Angle
 from autode.exceptions import CoordinateTransformFailed
-from autode.opt.coordinates.internals import PrimitiveInverseDistances, PIC
+from autode.opt.coordinates.internals import (
+    PrimitiveInverseDistances,
+    PIC,
+    build_pic_from_species,
+)
 from autode.opt.coordinates.cartesian import CartesianCoordinates
 from autode.opt.coordinates.dic import DIC
 from autode.opt.coordinates.primitives import (
@@ -756,6 +762,9 @@ def test_repr():
         PrimitiveBondAngle(0, 1, 2),
         ConstrainedPrimitiveBondAngle(0, 1, 2, value=1.0),
         PrimitiveDihedralAngle(0, 1, 2, 3),
+        PrimitiveLinearAngle(0, 1, 2, 3, LinearBendType.BEND),
+        PrimitiveLinearAngle(0, 1, 2, 3, LinearBendType.COMPLEMENT),
+        PrimitiveDummyLinearAngle(0, 1, 2, LinearBendType.BEND),
     ]
 
     for p in prims:
@@ -804,3 +813,11 @@ def test_dics_cannot_be_built_with_incomplete_primitives():
 
     with pytest.raises(RuntimeError):
         _ = DIC.from_cartesian(x=x, primitives=primitives)
+
+
+def test_pic_generation():
+    m = feco5_mol()
+    pic = build_pic_from_species(m)
+
+    # check that there are no duplicates
+    assert not any(ic1 == ic2 for ic1, ic2 in itertools.combinations(pic, r=2))
diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index 50854b689..599815f56 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -116,7 +116,8 @@ def test_primitive_projection_discard():
     r_initial = optimiser._species.distance(0, 1)
 
     x = CartesianCoordinates(optimiser._species.coordinates)
-    s = DICWithConstraints.from_cartesian(x, optimiser._primitives)
+    optimiser._build_internal_coordinates()
+    s = optimiser._coords
     assert len(s) == 3
 
     # Shift on the first couple of DIC but nothing on the final one

From de52e9e722ca50114e71c008f7bdd3ee82c1f3d6 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Fri, 22 Dec 2023 10:14:08 +0000
Subject: [PATCH 21/47] add tests

---
 autode/opt/coordinates/internals.py | 11 ++++++-----
 tests/test_opt/molecules.py         |  4 ++--
 tests/test_opt/test_coordiantes.py  |  5 ++++-
 tests/test_opt/test_crfo.py         | 16 +++++++++++++++-
 4 files changed, 27 insertions(+), 9 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 09a45cb86..3ef16c1f6 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -388,9 +388,10 @@ def _add_angles_from_species(
     connectivity graph supplied
 
     Args:
-        pic: The AnyPIC instance (modified in-place)
-        mol: The species object
-        core_graph: The connectivity graph
+        pic (AnyPIC): The AnyPIC instance (modified in-place)
+        mol (Species): The species object
+        core_graph (MolecularGraph): The connectivity graph
+        lin_thresh (Angle): The angle threshold for linearity
     """
 
     def get_ref_atom(a, b, c):
@@ -406,10 +407,10 @@ def get_ref_atom(a, b, c):
         deviations_from_90 = []
         for atom in near_atoms:
             i_b_a = mol.angle(atom, b, a)
-            if i_b_a > lin_thresh or i_b_a < 180 - lin_thresh:
+            if i_b_a > lin_thresh or i_b_a < (Angle(180, "deg") - lin_thresh):
                 continue
             i_b_c = mol.angle(atom, b, c)
-            if i_b_c > lin_thresh or i_b_c < 180 - lin_thresh:
+            if i_b_c > lin_thresh or i_b_c < (Angle(180, "deg") - lin_thresh):
                 continue
             deviation_a = abs(i_b_a - Angle(90, "deg"))
             deviation_c = abs(i_b_c - Angle(90, "deg"))
diff --git a/tests/test_opt/molecules.py b/tests/test_opt/molecules.py
index 145d94066..fdf0e6eb9 100644
--- a/tests/test_opt/molecules.py
+++ b/tests/test_opt/molecules.py
@@ -60,11 +60,11 @@ def feco5_mol():
             Atom("O", -1.5139, -3.5069, -0.3015),
             Atom("C", -0.8455, -2.5982, -0.3019),
             Atom("Fe", 0.3545, -0.9664, -0.3020),
-            Atom("C", 1.5539, 0.6660, -0.3004),
+            Atom("C", 1.4039, 0.5711, -0.2183),
             Atom("O", 2.2216, 1.5751, -0.2992),
             Atom("C", 1.2470, -1.6232, 1.3934),
             Atom("O", 1.7445, -1.9892, 2.3373),
-            Atom("C", 1.0930, -1.5083, -2.1084),
+            Atom("C", 1.0993, -1.3993, -2.0287),
             Atom("O", 1.5042, -1.8102, -3.1145),
             Atom("C", -1.2751, 0.2315, -0.1941),
             Atom("O", -2.1828, 0.8986, -0.1345),
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 1c80f6e94..9e1fbf6fe 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -815,9 +815,12 @@ def test_dics_cannot_be_built_with_incomplete_primitives():
         _ = DIC.from_cartesian(x=x, primitives=primitives)
 
 
-def test_pic_generation():
+def test_pic_generation_linear_angle_ref():
+    # Fe(CO)5 with linear Fe-C-O bonds
     m = feco5_mol()
     pic = build_pic_from_species(m)
 
     # check that there are no duplicates
     assert not any(ic1 == ic2 for ic1, ic2 in itertools.combinations(pic, r=2))
+    # check that linear bends use reference atoms, not dummy
+    assert not any(isinstance(ic, PrimitiveDummyLinearAngle) for ic in pic)
diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index 599815f56..3f52770b7 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -11,7 +11,7 @@
 from autode.opt.coordinates import CartesianCoordinates, DICWithConstraints
 from autode.opt.coordinates.primitives import PrimitiveDihedralAngle
 from autode.utils import work_in_tmp_dir
-from .molecules import h2o2_mol
+from .molecules import h2o2_mol, acetylene_mol
 from ..testutils import requires_working_xtb_install
 
 
@@ -340,3 +340,17 @@ def test_linear_dihedrals_are_removed():
     assert not any(
         isinstance(q, PrimitiveDihedralAngle) for q in dic.primitives
     )
+
+
+@requires_working_xtb_install
+@work_in_tmp_dir()
+def test_optimise_linear_molecule():
+    mol = acetylene_mol()
+    # the two H-C-C angles are almost linear
+    assert val.Angle(170, "deg") < mol.angle(0, 1, 3) < val.Angle(176, "deg")
+    assert val.Angle(170, "deg") < mol.angle(2, 0, 1) < val.Angle(176, "deg")
+    opt = CRFOptimiser(maxiter=10, gtol=1e-4, etol=1e-5)
+    opt.run(mol, XTB())
+    assert opt.converged
+    assert mol.angle(0, 1, 3) > val.Angle(179, "deg")
+    assert mol.angle(2, 0, 1) > val.Angle(179, "deg")

From 40abfd7cc01afbe9a8e09c1751ad6f78745f4693 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Fri, 22 Dec 2023 11:14:11 +0000
Subject: [PATCH 22/47] minor dihedral fix

---
 autode/opt/coordinates/internals.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 3ef16c1f6..d60ab95a4 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -523,9 +523,15 @@ def concatenate_adjacent_chains(chains_list):
                 # avoid triangle rings like cyclopropane
                 if n == m:
                     continue
-
-                is_linear_1 = mol.angle(m, o, p) > lin_thresh
-                is_linear_2 = mol.angle(o, p, n) > lin_thresh
+                zero_angle_thresh = Angle(180, "deg") - lin_thresh
+                is_linear_1 = (
+                    mol.angle(m, o, p) > lin_thresh
+                    or mol.angle(m, o, p) < zero_angle_thresh
+                )
+                is_linear_2 = (
+                    mol.angle(o, p, n) > lin_thresh
+                    or mol.angle(o, p, n) < zero_angle_thresh
+                )
 
                 # if any angle is linear, don't add dihedral
                 if is_linear_1 or is_linear_2:

From f600b3cede566617c84cf11b6a218467c99521f5 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Fri, 22 Dec 2023 20:32:42 +0000
Subject: [PATCH 23/47] dihedral code fix

---
 autode/opt/coordinates/internals.py | 104 +++++++++++++++++-----------
 autode/opt/optimisers/crfo.py       |   7 +-
 2 files changed, 70 insertions(+), 41 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index d60ab95a4..1085a212c 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -482,60 +482,86 @@ def _add_dihedrals_from_species(
     if mol.n_atoms < 4:
         return
 
+    zero_angle_thresh = Angle(180, "deg") - lin_thresh
+
+    def is_dihedral_well_defined(w, x, y, z):
+        """A dihedral is well-defined if any angle is not linear"""
+        is_linear_1 = (
+            mol.angle(w, x, y) > lin_thresh
+            or mol.angle(w, x, y) < zero_angle_thresh
+        )
+        is_linear_2 = (
+            mol.angle(x, y, z) > lin_thresh
+            or mol.angle(x, y, z) < zero_angle_thresh
+        )
+        return not (is_linear_1 or is_linear_2)
+
+    # add normal dihedrals
+    for o, p in list(core_graph.edges):
+        for m in core_graph.neighbors(o):
+            if m == p:
+                continue
+
+            for n in core_graph.neighbors(p):
+                if n == o:
+                    continue
+
+                # avoid triangle rings like cyclopropane
+                if n == m:
+                    continue
+
+                if is_dihedral_well_defined(m, o, p, n):
+                    pic.append(PrimitiveDihedralAngle(m, o, p, n))
+
     # find all linear atom chains A--B--C--D... and add dihedrals to terminal atoms
-    linear_chains = []
+
+    def extend_chain(chain: List[int]):
+        for idx in range(mol.n_atoms):
+            if idx in chain:
+                continue
+            # if idx -- 0 -- 1 > 170 degrees
+            if mol.angle(chain[1], chain[0], idx) > lin_thresh:
+                chain.insert(0, idx)
+                continue
+            # if (-2) -- (-1) -- idx > 170 degrees
+            if mol.angle(chain[-2], chain[-1], idx) > lin_thresh:
+                chain.append(idx)
+                continue
+
+    linear_chains: List[list] = []
     for b in range(mol.n_atoms):
+        if any(b in chain for chain in linear_chains):
+            continue
         for a, c in itertools.combinations(core_graph.neighbors(b), r=2):
-            if mol.angle(a, b, c) > lin_thresh:
-                linear_chains.append((a, b, c))
-
-    def concatenate_adjacent_chains(chains_list):
-        for chain1, chain2 in itertools.combinations(chains_list, r=2):
-            if chain1[0] == chain2[0]:
-                new_chain = list(reversed(chain2)) + list(chain1)
-            elif chain1[0] == chain2[-1]:
-                new_chain = list(chain2) + list(chain1)
-            elif chain1[-1] == chain2[0]:
-                new_chain = list(chain1) + list(chain2)
-            elif chain1[-1] == chain2[-1]:
-                new_chain = list(chain1) + list(reversed(chain2))
-            else:
+            if any(a in chain for chain in linear_chains) or any(
+                c in chain for chain in linear_chains
+            ):
                 continue
+            if mol.angle(a, b, c) > lin_thresh:
+                chain = [a, b, c]
+                extend_chain(chain)
+                linear_chains.append(chain)
 
-            chains_list.remove(chain1)
-            chains_list.remove(chain2)
-            chains_list.append(tuple(new_chain))
-            return concatenate_adjacent_chains(chains_list)
-
-    concatenate_adjacent_chains(linear_chains)
-    terminal_points = [(chain[0], chain[-1]) for chain in linear_chains]
-
-    # add normal + linear chain dihedrals
-    for o, p in list(core_graph.edges) + terminal_points:
+    # add linear chain dihedrals
+    for chain in linear_chains:
+        o, p = chain[0], chain[-1]
         for m in core_graph.neighbors(o):
             if m == p:
                 continue
 
+            if m in chain:
+                continue
+
             for n in core_graph.neighbors(p):
                 if n == o:
                     continue
 
-                # avoid triangle rings like cyclopropane
                 if n == m:
                     continue
-                zero_angle_thresh = Angle(180, "deg") - lin_thresh
-                is_linear_1 = (
-                    mol.angle(m, o, p) > lin_thresh
-                    or mol.angle(m, o, p) < zero_angle_thresh
-                )
-                is_linear_2 = (
-                    mol.angle(o, p, n) > lin_thresh
-                    or mol.angle(o, p, n) < zero_angle_thresh
-                )
-
-                # if any angle is linear, don't add dihedral
-                if is_linear_1 or is_linear_2:
+
+                if n in chain:
                     continue
-                else:
+
+                if is_dihedral_well_defined(m, o, p, n):
                     pic.append(PrimitiveDihedralAngle(m, o, p, n))
     return None
diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py
index e1973f396..5c86640db 100644
--- a/autode/opt/optimisers/crfo.py
+++ b/autode/opt/optimisers/crfo.py
@@ -13,7 +13,10 @@
 from autode.opt.coordinates import CartesianCoordinates, DICWithConstraints
 from autode.opt.coordinates.internals import build_pic_from_species
 from autode.opt.optimisers.rfo import RFOptimiser
-from autode.opt.optimisers.hessian_update import BFGSDampedUpdate, NullUpdate
+from autode.opt.optimisers.hessian_update import (
+    BFGSDampedUpdate,
+    BFGSSR1Update,
+)
 
 
 class CRFOptimiser(RFOptimiser):
@@ -35,7 +38,7 @@ def __init__(
 
         self.alpha = Distance(init_alpha, units="ang")
         assert self.alpha > 0
-        self._hessian_update_types = [BFGSDampedUpdate, NullUpdate]
+        self._hessian_update_types = [BFGSDampedUpdate, BFGSSR1Update]
 
     def _step(self) -> None:
         """Partitioned rational function step"""

From 113f621ac2e38ade6b10d228f4e7b651987f9bb3 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Fri, 22 Dec 2023 20:41:43 +0000
Subject: [PATCH 24/47] add test for linear bend with ref

---
 tests/test_opt/test_coordiantes.py |  2 ++
 tests/test_opt/test_crfo.py        | 14 +++++++++++++-
 2 files changed, 15 insertions(+), 1 deletion(-)

diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 9e1fbf6fe..5ca47ae1b 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -824,3 +824,5 @@ def test_pic_generation_linear_angle_ref():
     assert not any(ic1 == ic2 for ic1, ic2 in itertools.combinations(pic, r=2))
     # check that linear bends use reference atoms, not dummy
     assert not any(isinstance(ic, PrimitiveDummyLinearAngle) for ic in pic)
+    # there should not be any dihedral for this geometry
+    assert not any(isinstance(ic, PrimitiveDihedralAngle) for ic in pic)
diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index 3f52770b7..275f9e4d4 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -11,7 +11,7 @@
 from autode.opt.coordinates import CartesianCoordinates, DICWithConstraints
 from autode.opt.coordinates.primitives import PrimitiveDihedralAngle
 from autode.utils import work_in_tmp_dir
-from .molecules import h2o2_mol, acetylene_mol
+from .molecules import h2o2_mol, acetylene_mol, feco5_mol
 from ..testutils import requires_working_xtb_install
 
 
@@ -354,3 +354,15 @@ def test_optimise_linear_molecule():
     assert opt.converged
     assert mol.angle(0, 1, 3) > val.Angle(179, "deg")
     assert mol.angle(2, 0, 1) > val.Angle(179, "deg")
+
+
+@requires_working_xtb_install
+@work_in_tmp_dir()
+def test_optimise_linear_bend_with_ref():
+    mol = feco5_mol()
+    # the Fe-C-O angle are manually deviated
+    assert val.Angle(170, "deg") < mol.angle(2, 3, 4) < val.Angle(176, "deg")
+    # large molecule so allow few iters, no need to converge fully
+    opt = CRFOptimiser(maxiter=10, gtol=1e-4, etol=1e-5)
+    opt.run(mol, XTB())
+    assert mol.angle(2, 3, 4) > val.Angle(178, "deg")

From 902aa749307033263e64b82e07adad941fa92e4b Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Fri, 22 Dec 2023 20:44:50 +0000
Subject: [PATCH 25/47] add test for linear bend with ref

---
 tests/test_opt/test_crfo.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index 275f9e4d4..d613d8bc7 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -161,7 +161,7 @@ def test_xtb_opt_with_distance_constraint():
 
     assert np.isclose(water.distance(0, 1), 0.99, atol=0.01)
 
-    CRFOptimiser.optimise(species=water, method=XTB())
+    CRFOptimiser.optimise(species=water, method=XTB(), etol=1e-6)
 
     # Optimisation should generate an O-H distance *very* close to 1.1 Å
     assert np.isclose(water.distance(0, 1).to("Å"), 1.1, atol=1e-4)

From de46ccc032a7e4f6827bda9da5527b409b1e2f1a Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 23 Dec 2023 12:25:12 +0000
Subject: [PATCH 26/47] add tests

---
 autode/opt/coordinates/internals.py |  4 +--
 tests/test_opt/test_coordiantes.py  | 39 +++++++++++++++++++++++++++++
 2 files changed, 41 insertions(+), 2 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 1085a212c..5f262b45c 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -318,13 +318,13 @@ def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
                     min_dist = mol.distance(i, j)
                     min_pair = (i, j)
             # avoid connecting distant components
-            if min_dist < 5.0:
+            if min_dist < Distance(4.0, "ang"):
                 core_graph.add_edge(*min_pair, pi=False, active=False)
 
     if not core_graph.is_connected:
         raise RuntimeError(
             "Unable to join all the fragments, distance between "
-            "one or more pairs of fragments is too high (>5.0 Å)"
+            "one or more pairs of fragments is too high (>4.0 Å)"
         )
 
     # The constraints should be counted as bonds
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 5ca47ae1b..6576a2972 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -826,3 +826,42 @@ def test_pic_generation_linear_angle_ref():
     assert not any(isinstance(ic, PrimitiveDummyLinearAngle) for ic in pic)
     # there should not be any dihedral for this geometry
     assert not any(isinstance(ic, PrimitiveDihedralAngle) for ic in pic)
+
+
+def test_pic_generation_disjoin_graph():
+    # the algorithm should fully connect the graph
+    xyz_string = (
+        "16\n\n"
+        "C      -0.00247        1.65108        0.05872\n"
+        "C       1.19010        1.11169        0.27709\n"
+        "C       1.58519       -0.30014        0.31049\n"
+        "C       0.05831       -1.54292       -0.45110\n"
+        "C      -1.18798       -1.04262        0.13551\n"
+        "C      -1.28206        0.99883       -0.23631\n"
+        "H      -0.07432        2.73634        0.08639\n"
+        "H       2.01755        1.78921        0.47735\n"
+        "H       1.70503       -0.70916        1.30550\n"
+        "H       2.40398       -0.55376       -0.34855\n"
+        "H       0.44229       -2.48695       -0.08638\n"
+        "H       0.15289       -1.41865       -1.51944\n"
+        "H      -1.25410       -1.13318        1.21833\n"
+        "H      -2.09996       -1.35918       -0.36715\n"
+        "H      -2.09462        1.29055        0.41495\n"
+        "H      -1.56001        1.00183       -1.28217\n"
+    )
+    with open("diels_alder_complex.xyz", "w") as fh:
+        fh.write(xyz_string)
+
+    mol = Molecule("diels_alder_complex.xyz")
+    assert not mol.graph.is_connected
+    pic = build_pic_from_species(mol)
+
+    # shortest bond is between 4, 5 which should also generate angle, torsion
+    assert PrimitiveDistance(4, 5) in pic
+    assert PrimitiveBondAngle(3, 4, 5) in pic
+    assert PrimitiveBondAngle(4, 5, 0) in pic
+    assert PrimitiveDihedralAngle(3, 4, 5, 0) in pic
+
+    # the other distance between fragments is 2, 3 which should not be connected
+    assert PrimitiveDistance(2, 3) not in pic
+    assert PrimitiveBondAngle(1, 2, 3) not in pic

From bbe1f1458a7c4cbc58a1972b043872ea4b96b470 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 23 Dec 2023 12:33:41 +0000
Subject: [PATCH 27/47] make bond angle consistent

---
 autode/opt/coordinates/internals.py  |  2 +-
 autode/opt/coordinates/primitives.py | 14 +++++++-------
 tests/test_opt/test_coordiantes.py   | 22 +++++++++++-----------
 tests/test_opt/test_crfo.py          | 14 +++++++-------
 4 files changed, 26 insertions(+), 26 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 5f262b45c..3fd94d71d 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -426,7 +426,7 @@ def get_ref_atom(a, b, c):
         for n, m in itertools.combinations(core_graph.neighbors(o), r=2):
             # avoid almost linear angles
             if mol.angle(m, o, n) < lin_thresh:
-                pic.append(PrimitiveBondAngle(o=o, m=m, n=n))
+                pic.append(PrimitiveBondAngle(m=m, o=o, n=n))
             else:
                 # if central atom is connected to another, no need to include
                 other_neighbours = list(core_graph.neighbors(o))
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 5eedc7647..52efca8a1 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -299,12 +299,12 @@ class PrimitiveBondAngle(Primitive):
     arccosine of the normalised dot product
     """
 
-    def __init__(self, o: int, m: int, n: int):
+    def __init__(self, m: int, o: int, n: int):
         """Bond angle m-o-n"""
-        super().__init__(o, m, n)
+        super().__init__(m, o, n)
 
-        self.o = int(o)
         self.m = int(m)
+        self.o = int(o)
         self.n = int(n)
 
     def __eq__(self, other) -> bool:
@@ -332,22 +332,22 @@ def __repr__(self):
 
 
 class ConstrainedPrimitiveBondAngle(ConstrainedPrimitive, PrimitiveBondAngle):
-    def __init__(self, o: int, m: int, n: int, value: float):
+    def __init__(self, m: int, o: int, n: int, value: float):
         """
         Angle (m-o-n) constrained to a value (in radians)
 
         -----------------------------------------------------------------------
         Arguments:
 
-            o: Atom index
-
             m: Atom index
 
+            o: Atom index
+
             n: Atom index
 
             value: Required value of the constrained angle
         """
-        super().__init__(o=o, m=m, n=n)
+        super().__init__(m=m, o=o, n=n)
 
         self._theta0 = value
 
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 6576a2972..4796dbd4a 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -566,7 +566,7 @@ def numerical_derivative(a, b, h=1e-8):
     m = water_mol()
     init_coords = m.coordinates.copy()
 
-    angle = PrimitiveBondAngle(0, 1, 2)
+    angle = PrimitiveBondAngle(1, 0, 2)
     derivs = angle.derivative(init_coords)
     for atom_idx in (0, 1, 2):
         for component in (0, 1, 2):
@@ -578,8 +578,8 @@ def numerical_derivative(a, b, h=1e-8):
 
 
 def test_angle_primitive_equality():
-    assert PrimitiveBondAngle(0, 1, 2) == PrimitiveBondAngle(0, 2, 1)
-    assert PrimitiveBondAngle(0, 1, 2) != PrimitiveBondAngle(2, 1, 0)
+    assert PrimitiveBondAngle(1, 0, 2) == PrimitiveBondAngle(2, 0, 1)
+    assert PrimitiveBondAngle(1, 0, 2) != PrimitiveBondAngle(1, 2, 0)
 
 
 def test_dihedral_value():
@@ -662,7 +662,7 @@ def test_primitives_consistent_with_mol_values():
     assert np.isclose(dist(coords), h2o2.distance(0, 1), rtol=1e-8)
     invdist = PrimitiveInverseDistance(1, 2)
     assert np.isclose(invdist(coords), 1 / h2o2.distance(1, 2), rtol=1e-8)
-    ang = PrimitiveBondAngle(2, 0, 1)  # bond is defined in a different way
+    ang = PrimitiveBondAngle(0, 2, 1)
     assert np.isclose(ang(coords), h2o2.angle(0, 2, 1), rtol=1e-8)
     dihedral = PrimitiveDihedralAngle(2, 0, 1, 3)
     assert np.isclose(dihedral(coords), h2o2.dihedral(2, 0, 1, 3), rtol=1e-8)
@@ -702,8 +702,8 @@ def test_primitives_consistent_with_mol_values():
     water_mol(), water_mol(), *extra_mols
 ]
 test_prims = [
-    PrimitiveDihedralAngle(2, 0, 1, 3), PrimitiveBondAngle(2, 0, 1),
-    PrimitiveBondAngle(0, 1, 2), PrimitiveDistance(0, 1),
+    PrimitiveDihedralAngle(2, 0, 1, 3), PrimitiveBondAngle(0, 2, 1),
+    PrimitiveBondAngle(1, 0, 2), PrimitiveDistance(0, 1),
     PrimitiveInverseDistance(0, 1), PrimitiveDihedralAngle(2, 0, 1, 3),
     PrimitiveDihedralAngle(2, 0, 1, 3),
     PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.BEND),
@@ -759,8 +759,8 @@ def test_repr():
         PrimitiveInverseDistance(0, 1),
         PrimitiveDistance(0, 1),
         ConstrainedPrimitiveDistance(0, 1, value=1e-3),
-        PrimitiveBondAngle(0, 1, 2),
-        ConstrainedPrimitiveBondAngle(0, 1, 2, value=1.0),
+        PrimitiveBondAngle(1, 0, 2),
+        ConstrainedPrimitiveBondAngle(1, 0, 2, value=1.0),
         PrimitiveDihedralAngle(0, 1, 2, 3),
         PrimitiveLinearAngle(0, 1, 2, 3, LinearBendType.BEND),
         PrimitiveLinearAngle(0, 1, 2, 3, LinearBendType.COMPLEMENT),
@@ -789,7 +789,7 @@ def test_dic_large_step_allowed_unconverged_back_transform():
 
 
 def test_constrained_angle_delta():
-    q = ConstrainedPrimitiveBondAngle(0, 1, 2, value=np.pi)
+    q = ConstrainedPrimitiveBondAngle(1, 0, 2, value=np.pi)
     mol = water_mol()
     theta = mol.angle(1, 0, 2)
     x = CartesianCoordinates(mol.coordinates)
@@ -798,8 +798,8 @@ def test_constrained_angle_delta():
 
 
 def test_constrained_angle_equality():
-    a = ConstrainedPrimitiveBondAngle(0, 1, 2, value=np.pi)
-    b = ConstrainedPrimitiveBondAngle(0, 2, 1, value=np.pi)
+    a = ConstrainedPrimitiveBondAngle(1, 0, 2, value=np.pi)
+    b = ConstrainedPrimitiveBondAngle(2, 0, 1, value=np.pi)
 
     assert a == b
 
diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index d613d8bc7..794e0f297 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -200,7 +200,7 @@ def test_baker1997_example():
     r1 = prim.ConstrainedPrimitiveDistance(0, 1, value=1.5)
     r2 = prim.ConstrainedPrimitiveDistance(3, 4, value=2.5)
     theta = prim.ConstrainedPrimitiveBondAngle(
-        1, 0, 5, value=val.Angle(123.0, "º").to("rad")
+        0, 1, 5, value=val.Angle(123.0, "º").to("rad")
     )
 
     pic = PIC(r1, r2, theta)
@@ -208,12 +208,12 @@ def test_baker1997_example():
         pic.append(prim.PrimitiveDistance(*pair))
 
     for triple in (
-        (1, 0, 3),
-        (1, 0, 3),
-        (2, 0, 3),
-        (4, 1, 0),
-        (5, 1, 0),
-        (4, 1, 5),
+        (0, 1, 3),
+        (0, 1, 3),
+        (0, 2, 3),
+        (1, 4, 0),
+        (1, 5, 0),
+        (1, 4, 5),
     ):
         pic.append(prim.PrimitiveBondAngle(*triple))
 

From c1a90feea8e6897bf0b3f0498032986cb2196afe Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 25 Dec 2023 10:58:37 +0000
Subject: [PATCH 28/47] add tests and fix value bug

---
 autode/opt/coordinates/internals.py |  6 +++---
 autode/values.py                    |  4 ++++
 tests/test_opt/molecules.py         | 16 ++++++++++++++++
 tests/test_opt/test_coordiantes.py  | 26 ++++++++++++++++++++++++--
 4 files changed, 47 insertions(+), 5 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 3fd94d71d..6dde80bd9 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -404,7 +404,7 @@ def get_ref_atom(a, b, c):
             and idx not in (a, b, c)
         ]
         # get atoms closest to perpendicular
-        deviations_from_90 = []
+        deviations_from_90 = {}
         for atom in near_atoms:
             i_b_a = mol.angle(atom, b, a)
             if i_b_a > lin_thresh or i_b_a < (Angle(180, "deg") - lin_thresh):
@@ -415,12 +415,12 @@ def get_ref_atom(a, b, c):
             deviation_a = abs(i_b_a - Angle(90, "deg"))
             deviation_c = abs(i_b_c - Angle(90, "deg"))
             avg_dev = (deviation_a + deviation_c) / 2
-            deviations_from_90.append(avg_dev)
+            deviations_from_90[atom] = avg_dev
 
         if len(deviations_from_90) == 0:
             return None
 
-        return near_atoms[np.argmin(deviations_from_90)]
+        return min(deviations_from_90, key=deviations_from_90.get)
 
     for o in range(mol.n_atoms):
         for n, m in itertools.combinations(core_graph.neighbors(o), r=2):
diff --git a/autode/values.py b/autode/values.py
index f36196276..ea771c631 100644
--- a/autode/values.py
+++ b/autode/values.py
@@ -240,6 +240,10 @@ def __radd__(self, other) -> TypeValue:
     def __sub__(self, other) -> TypeValue:
         return self.__add__(-other)
 
+    def __neg__(self) -> TypeValue:
+        """Unary negation operation"""
+        return self._like_self_from_float(-float(self))
+
     def __floordiv__(self, other) -> Union[float, TypeValue]:
         x = float(self) // self._other_same_units(other)
         return x if isinstance(other, Value) else self._like_self_from_float(x)
diff --git a/tests/test_opt/molecules.py b/tests/test_opt/molecules.py
index fdf0e6eb9..43956fc06 100644
--- a/tests/test_opt/molecules.py
+++ b/tests/test_opt/molecules.py
@@ -70,3 +70,19 @@ def feco5_mol():
             Atom("O", -2.1828, 0.8986, -0.1345),
         ]
     )
+
+
+def cumulene_mol():
+    return Molecule(
+        atoms=[
+            Atom("C", -1.3968, 3.7829, 0.9582),
+            Atom("C", -0.8779, 2.7339, 0.6902),
+            Atom("C", -1.9158, 4.8319, 1.2262),
+            Atom("C", -2.4770, 5.9661, 1.5160),
+            Atom("C", -0.3167, 1.5996, 0.4004),
+            Atom("H", -2.1197, 6.8879, 1.0720),
+            Atom("H", -3.3113, 6.0085, 2.2063),
+            Atom("H", 0.1364, 1.4405, -0.5711),
+            Atom("H", -0.2928, 0.7947, 1.1256),
+        ]
+    )
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 4796dbd4a..54f0d411f 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -2,7 +2,14 @@
 
 import pytest
 import numpy as np
-from .molecules import h2, methane_mol, water_mol, h2o2_mol, feco5_mol
+from .molecules import (
+    h2,
+    methane_mol,
+    water_mol,
+    h2o2_mol,
+    feco5_mol,
+    cumulene_mol,
+)
 from autode.atoms import Atom
 from autode.species.molecule import Molecule
 from autode.values import Angle
@@ -828,7 +835,7 @@ def test_pic_generation_linear_angle_ref():
     assert not any(isinstance(ic, PrimitiveDihedralAngle) for ic in pic)
 
 
-def test_pic_generation_disjoin_graph():
+def test_pic_generation_disjoint_graph():
     # the algorithm should fully connect the graph
     xyz_string = (
         "16\n\n"
@@ -865,3 +872,18 @@ def test_pic_generation_disjoin_graph():
     # the other distance between fragments is 2, 3 which should not be connected
     assert PrimitiveDistance(2, 3) not in pic
     assert PrimitiveBondAngle(1, 2, 3) not in pic
+
+
+def test_pic_generation_chain_dihedrals():
+    # extra dihedrals are needed for ends of linear chains like allene
+    cumulene = cumulene_mol()
+    pic = build_pic_from_species(cumulene)
+
+    assert PrimitiveDihedralAngle(5, 3, 4, 8) in pic
+    assert PrimitiveDihedralAngle(6, 3, 4, 7) in pic
+    assert PrimitiveDihedralAngle(8, 4, 3, 6) in pic
+    assert PrimitiveDihedralAngle(7, 4, 3, 6) in pic
+
+    # check that the 3N-6 degrees of freedom are maintained
+    _ = pic(cumulene.coordinates.flatten())
+    assert np.linalg.matrix_rank(pic.B) == 3 * cumulene.n_atoms - 6

From f2e5f65565762c09451e4a838ff56e1e6e47244f Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 25 Dec 2023 15:40:36 +0000
Subject: [PATCH 29/47] add tests

---
 autode/opt/coordinates/internals.py | 93 ++++++++++++++---------------
 tests/test_opt/test_coordiantes.py  | 36 ++++++++++-
 2 files changed, 80 insertions(+), 49 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 6dde80bd9..f9e1496fe 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -12,7 +12,6 @@
 import copy
 
 import numpy as np
-from enum import Enum
 import itertools
 from typing import Any, Optional, Type, List, TYPE_CHECKING
 from abc import ABC, abstractmethod
@@ -98,7 +97,9 @@ def __init__(self, *args: Any):
     def append(self, item: Primitive) -> None:
         """Append an item to this set of primitives"""
         assert isinstance(item, Primitive), "Must be a Primitive type!"
-        super().append(item)
+        # prevent duplicate primitives
+        if item not in self:
+            super().append(item)
 
     @property
     def B(self) -> np.ndarray:
@@ -251,7 +252,6 @@ def _populate_all(self, x: np.ndarray) -> None:
 
 def build_pic_from_species(
     mol: "Species",
-    aux_bonds=False,
 ) -> AnyPIC:
     """
     Build a set of primitives from the species, using the graph as
@@ -261,14 +261,13 @@ def build_pic_from_species(
 
     Args:
         mol:
-        aux_bonds:
 
     Returns:
         (AnyPIC): The set of primitive internals
     """
     pic = AnyPIC()
     core_graph = _get_connected_graph_from_species(mol)
-    _add_bonds_from_species(pic, mol, core_graph, aux_bonds=aux_bonds)
+    _add_bonds_from_species(pic, mol, core_graph)
     _add_angles_from_species(pic, mol, core_graph)
     _add_dihedrals_from_species(pic, mol, core_graph)
     return pic
@@ -340,7 +339,6 @@ def _add_bonds_from_species(
     pic: AnyPIC,
     mol: "Species",
     core_graph: "MolecularGraph",
-    aux_bonds: bool = False,
 ):
     """
     Modify the supplied AnyPIC instance in-place by adding bonds, from the
@@ -350,7 +348,6 @@ def _add_bonds_from_species(
         pic: The AnyPIC instance (modified in-place)
         mol: The species object
         core_graph: The connectivity graph
-        aux_bonds: Whether to add auxiliary bonds (< 2.5 * covalent radii sum)
     """
     n = 0
     for i, j in sorted(core_graph.edges):
@@ -365,15 +362,6 @@ def _add_bonds_from_species(
             pic.append(PrimitiveDistance(i, j))
     assert n == mol.constraints.n_distance
 
-    if not aux_bonds:
-        return None
-
-    # add auxiliary bonds if specified
-    for i, j in itertools.combinations(range(mol.n_atoms), r=2):
-        if core_graph.has_edge(i, j):
-            continue
-        if mol.distance(i, j) < 2.5 * mol.eqm_bond_distance(i, j):
-            pic.append(PrimitiveDistance(i, j))
     return None
 
 
@@ -381,7 +369,7 @@ def _add_angles_from_species(
     pic: AnyPIC,
     mol: "Species",
     core_graph: "MolecularGraph",
-    lin_thresh=Angle(170, "deg"),
+    lin_thresh: Angle = Angle(170, "deg"),
 ) -> None:
     """
     Modify the set of primitives in-place by adding angles, from the
@@ -393,27 +381,36 @@ def _add_angles_from_species(
         core_graph (MolecularGraph): The connectivity graph
         lin_thresh (Angle): The angle threshold for linearity
     """
+    lin_thresh = lin_thresh.to("rad")
 
-    def get_ref_atom(a, b, c):
+    def get_ref_atom(a, b, c, bonded=False):
         """get a reference atom for a-b-c linear angle"""
-        # all atoms in 4 A radius except a, b, c
-        near_atoms = [
-            idx
-            for idx in range(mol.n_atoms)
-            if mol.distance(b, idx) < Distance(4.0, "ang")
-            and idx not in (a, b, c)
-        ]
+        # only check bonded atoms if requested
+        if bonded:
+            near_atoms = list(core_graph.neighbors(b))
+            near_atoms.remove(a)
+            near_atoms.remove(c)
+
+        # otherwise get all atoms in 4 A radius except a, b, c
+        else:
+            near_atoms = [
+                idx
+                for idx in range(mol.n_atoms)
+                if mol.distance(b, idx) < Distance(4.0, "ang")
+                and idx not in (a, b, c)
+            ]
+
         # get atoms closest to perpendicular
         deviations_from_90 = {}
         for atom in near_atoms:
             i_b_a = mol.angle(atom, b, a)
-            if i_b_a > lin_thresh or i_b_a < (Angle(180, "deg") - lin_thresh):
+            if i_b_a > lin_thresh or i_b_a < (np.pi - lin_thresh):
                 continue
             i_b_c = mol.angle(atom, b, c)
-            if i_b_c > lin_thresh or i_b_c < (Angle(180, "deg") - lin_thresh):
+            if i_b_c > lin_thresh or i_b_c < (np.pi - lin_thresh):
                 continue
-            deviation_a = abs(i_b_a - Angle(90, "deg"))
-            deviation_c = abs(i_b_c - Angle(90, "deg"))
+            deviation_a = abs(i_b_a - np.pi / 2)
+            deviation_c = abs(i_b_c - np.pi / 2)
             avg_dev = (deviation_a + deviation_c) / 2
             deviations_from_90[atom] = avg_dev
 
@@ -424,23 +421,20 @@ def get_ref_atom(a, b, c):
 
     for o in range(mol.n_atoms):
         for n, m in itertools.combinations(core_graph.neighbors(o), r=2):
-            # avoid almost linear angles
             if mol.angle(m, o, n) < lin_thresh:
                 pic.append(PrimitiveBondAngle(m=m, o=o, n=n))
             else:
-                # if central atom is connected to another, no need to include
-                other_neighbours = list(core_graph.neighbors(o))
-                other_neighbours.remove(m)
-                other_neighbours.remove(n)
-                if any(
-                    mol.angle(m, o, x) < lin_thresh for x in other_neighbours
-                ) or any(
-                    mol.angle(n, o, x) < lin_thresh for x in other_neighbours
-                ):
+                # If central atom is connected to another atom, then the
+                # linear angle is skipped and instead an out-of-plane
+                # (improper dihedral) coordinate is used
+                r = get_ref_atom(m, o, n, bonded=True)
+                if r is not None:
+                    pic.append(PrimitiveDihedralAngle(m, r, o, n))
                     continue
 
-                # for linear bends, ideally a reference atom is needed
-                r = get_ref_atom(m, o, n)
+                # Otherwise, we use a nearby (< 4.0 A) reference atom to
+                # define two orthogonal linear bends
+                r = get_ref_atom(m, o, n, bonded=False)
                 if r is not None:
                     pic.append(
                         PrimitiveLinearAngle(m, o, n, r, LinearBendType.BEND)
@@ -450,7 +444,10 @@ def get_ref_atom(a, b, c):
                             m, o, n, r, LinearBendType.COMPLEMENT
                         )
                     )
-                else:  # these use dummy atom for reference
+
+                # For completely linear molecules (CO2), there will be no such
+                # reference atoms, so use dummy atoms instead
+                else:
                     pic.append(
                         PrimitiveDummyLinearAngle(m, o, n, LinearBendType.BEND)
                     )
@@ -467,22 +464,24 @@ def _add_dihedrals_from_species(
     pic: AnyPIC,
     mol: "Species",
     core_graph: "MolecularGraph",
-    lin_thresh=Angle(170, "deg"),
+    lin_thresh: Angle = Angle(170, "deg"),
 ) -> None:
     """
     Modify the set of primitives in-place by adding dihedrals (torsions),
     from the connectivity graph supplied
 
     Args:
-        pic: The AnyPIC instance (modified in-place)
-        mol: The species
-        core_graph: The connectivity graph
+        pic (AnyPIC): The AnyPIC instance (modified in-place)
+        mol (Species): The species
+        core_graph (MolecularGraph): The connectivity graph
+        lin_thresh (Angle): The threshold for linearity
     """
     # no dihedrals possible with less than 4 atoms
     if mol.n_atoms < 4:
         return
 
-    zero_angle_thresh = Angle(180, "deg") - lin_thresh
+    lin_thresh = lin_thresh.to("rad")
+    zero_angle_thresh = np.pi - lin_thresh
 
     def is_dihedral_well_defined(w, x, y, z):
         """A dihedral is well-defined if any angle is not linear"""
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 54f0d411f..c538d0151 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -9,7 +9,9 @@
     h2o2_mol,
     feco5_mol,
     cumulene_mol,
+    acetylene_mol,
 )
+from autode.utils import work_in_tmp_dir
 from autode.atoms import Atom
 from autode.species.molecule import Molecule
 from autode.values import Angle
@@ -831,10 +833,31 @@ def test_pic_generation_linear_angle_ref():
     assert not any(ic1 == ic2 for ic1, ic2 in itertools.combinations(pic, r=2))
     # check that linear bends use reference atoms, not dummy
     assert not any(isinstance(ic, PrimitiveDummyLinearAngle) for ic in pic)
-    # there should not be any dihedral for this geometry
-    assert not any(isinstance(ic, PrimitiveDihedralAngle) for ic in pic)
+    # for C-Fe-C, one out-of-plane dihedral should be present
+    assert PrimitiveDihedralAngle(3, 5, 2, 1) in pic
+    # check degrees of freedom = 3N - 6
+    _ = pic(m.coordinates.flatten())
+    assert np.linalg.matrix_rank(pic.B) == 3 * m.n_atoms - 6
 
 
+def test_pic_generation_linear_angle_dummy():
+    # acetylene molecule
+    mol = acetylene_mol()
+    pic = build_pic_from_species(mol)
+
+    # there should not be any usual bond angles
+    assert not any(isinstance(ic, PrimitiveBondAngle) for ic in pic)
+    # there should not be any linear angles with reference atom
+    assert not any(isinstance(ic, PrimitiveLinearAngle) for ic in pic)
+    # there should be linear angles with dummy
+    assert any(isinstance(ic, PrimitiveDummyLinearAngle) for ic in pic)
+
+    # degrees of freedom = 3N - 5 for linear molecules
+    _ = pic(mol.coordinates.flatten())
+    assert np.linalg.matrix_rank(pic.B) == 3 * mol.n_atoms - 5
+
+
+@work_in_tmp_dir()
 def test_pic_generation_disjoint_graph():
     # the algorithm should fully connect the graph
     xyz_string = (
@@ -872,6 +895,15 @@ def test_pic_generation_disjoint_graph():
     # the other distance between fragments is 2, 3 which should not be connected
     assert PrimitiveDistance(2, 3) not in pic
     assert PrimitiveBondAngle(1, 2, 3) not in pic
+    # check degrees of freedom = 3N - 6
+    _ = pic(mol.coordinates.flatten())
+    assert np.linalg.matrix_rank(pic.B) == 3 * mol.n_atoms - 6
+
+    # if the bond between 2, 3 is made into a constraint, it will generate angles
+    mol.constraints.distance = {(2, 3): mol.distance(2, 3)}
+    pic = build_pic_from_species(mol)
+    assert ConstrainedPrimitiveDistance(2, 3, mol.distance(2, 3)) in pic
+    assert PrimitiveBondAngle(1, 2, 3) in pic
 
 
 def test_pic_generation_chain_dihedrals():

From 3e826c91abf353db7be359576fa221807e023a6b Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 25 Dec 2023 15:51:33 +0000
Subject: [PATCH 30/47] fix test

---
 tests/test_opt/test_crfo.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index 794e0f297..c9a4d95dd 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -363,6 +363,6 @@ def test_optimise_linear_bend_with_ref():
     # the Fe-C-O angle are manually deviated
     assert val.Angle(170, "deg") < mol.angle(2, 3, 4) < val.Angle(176, "deg")
     # large molecule so allow few iters, no need to converge fully
-    opt = CRFOptimiser(maxiter=10, gtol=1e-4, etol=1e-5)
+    opt = CRFOptimiser(maxiter=15, gtol=1e-4, etol=1e-5)
     opt.run(mol, XTB())
     assert mol.angle(2, 3, 4) > val.Angle(178, "deg")

From 483fbe6d1adaf08b5884930fe61e6164ddcfdf84 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Wed, 27 Dec 2023 14:58:05 +0000
Subject: [PATCH 31/47] change backup hessian update scheme for crfo

---
 autode/opt/optimisers/crfo.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py
index 5c86640db..05f3f9886 100644
--- a/autode/opt/optimisers/crfo.py
+++ b/autode/opt/optimisers/crfo.py
@@ -15,7 +15,7 @@
 from autode.opt.optimisers.rfo import RFOptimiser
 from autode.opt.optimisers.hessian_update import (
     BFGSDampedUpdate,
-    BFGSSR1Update,
+    BofillUpdate,
 )
 
 
@@ -38,7 +38,7 @@ def __init__(
 
         self.alpha = Distance(init_alpha, units="ang")
         assert self.alpha > 0
-        self._hessian_update_types = [BFGSDampedUpdate, BFGSSR1Update]
+        self._hessian_update_types = [BFGSDampedUpdate, BofillUpdate]
 
     def _step(self) -> None:
         """Partitioned rational function step"""

From 46ce93a9deef29c6e72bd238325645d2abf9922c Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Wed, 27 Dec 2023 19:59:24 +0000
Subject: [PATCH 32/47] revert change to hessian update scheme

---
 autode/opt/optimisers/crfo.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py
index 05f3f9886..5c86640db 100644
--- a/autode/opt/optimisers/crfo.py
+++ b/autode/opt/optimisers/crfo.py
@@ -15,7 +15,7 @@
 from autode.opt.optimisers.rfo import RFOptimiser
 from autode.opt.optimisers.hessian_update import (
     BFGSDampedUpdate,
-    BofillUpdate,
+    BFGSSR1Update,
 )
 
 
@@ -38,7 +38,7 @@ def __init__(
 
         self.alpha = Distance(init_alpha, units="ang")
         assert self.alpha > 0
-        self._hessian_update_types = [BFGSDampedUpdate, BofillUpdate]
+        self._hessian_update_types = [BFGSDampedUpdate, BFGSSR1Update]
 
     def _step(self) -> None:
         """Partitioned rational function step"""

From acddba4efeee8c94c2d8e07dd1c97b1cba8265a0 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 28 Dec 2023 14:47:54 +0000
Subject: [PATCH 33/47] correction to angle code

---
 autode/atoms.py                    | 13 ++++++++++++-
 tests/test_opt/test_coordiantes.py | 19 +++++++++++++++++++
 2 files changed, 31 insertions(+), 1 deletion(-)

diff --git a/autode/atoms.py b/autode/atoms.py
index 699e7f67f..916ce8e19 100644
--- a/autode/atoms.py
+++ b/autode/atoms.py
@@ -1041,7 +1041,18 @@ def angle(self, i: int, j: int, k: int) -> Angle:
                 f"least one zero vector"
             )
 
-        value = np.arccos(np.dot(vec1, vec2) / norms)
+        # Catch errors due to incomplete float precision
+        cos_value = np.dot(vec1, vec2) / norms
+        if -1 < cos_value < 1:
+            pass
+        elif cos_value > 1 and np.isclose(cos_value, 1, rtol=1e-8):
+            cos_value = 1
+        elif cos_value < -1 and np.isclose(cos_value, -1, rtol=1e-8):
+            cos_value = -1
+        else:
+            raise ValueError("Cos(angle) must be in [-1, 1] range!")
+
+        value = np.arccos(cos_value)
 
         return Angle(value)
 
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index c538d0151..eb4ce675f 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -919,3 +919,22 @@ def test_pic_generation_chain_dihedrals():
     # check that the 3N-6 degrees of freedom are maintained
     _ = pic(cumulene.coordinates.flatten())
     assert np.linalg.matrix_rank(pic.B) == 3 * cumulene.n_atoms - 6
+
+
+def test_pic_generation_square_planar():
+    ptcl4 = Molecule(
+        atoms=[
+            Atom("Pt", -0.1467, -0.2594, -0.0294),
+            Atom("Cl", -0.4597, -2.5963, -0.0523),
+            Atom("Cl", 2.1804, -0.5689, -0.2496),
+            Atom("Cl", -2.4738, 0.0501, 0.1908),
+            Atom("Cl", 0.1663, 2.0776, -0.0066),
+        ],
+        charge=-2,
+    )
+
+    # for sq planar, out-of-plane dihedrals are needed to have
+    # all degrees of freedom
+    pic = build_pic_from_species(ptcl4)
+    _ = pic(ptcl4.coordinates.flatten())
+    assert np.linalg.matrix_rank(pic.B) == 3 * 5 - 6

From 05305b03093afe3729d53ecdeb12c357ad58d75d Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 28 Dec 2023 14:59:16 +0000
Subject: [PATCH 34/47] assert only for graph

---
 autode/opt/coordinates/internals.py | 15 ++++++---------
 tests/test_opt/test_crfo.py         | 14 +++++++++++++-
 2 files changed, 19 insertions(+), 10 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index f9e1496fe..71cdec06a 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -275,21 +275,18 @@ def build_pic_from_species(
 
 def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
     """
-    Creates a fully connected graph from a species, by (1) joining
-    disconnected fragments by their shortest distance, (2) connecting
-    constrained bonds, (3) joining hydrogen bonds, if present.
+    Creates a fully connected graph from the graph of a species, by
+    (1) joining disconnected fragments by their shortest distance,
+    (2) connecting constrained bonds, (3) joining hydrogen bonds,
+    if present.
 
     Args:
-        mol: A species containing atoms and coordinates
+        mol: A species (must have atoms and graph)
 
     Returns:
         (MolecularGraph):
     """
-    # if graph does not exist, create one
-    if mol.graph is None:
-        mol.reset_graph()
-
-    assert mol.graph is not None
+    assert mol.graph is not None, "Species must have graph!"
     core_graph = copy.deepcopy(mol.graph)
 
     # join hydrogen bonds
diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index c9a4d95dd..ab9aa0ab8 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -11,7 +11,7 @@
 from autode.opt.coordinates import CartesianCoordinates, DICWithConstraints
 from autode.opt.coordinates.primitives import PrimitiveDihedralAngle
 from autode.utils import work_in_tmp_dir
-from .molecules import h2o2_mol, acetylene_mol, feco5_mol
+from .molecules import h2o2_mol, acetylene_mol, feco5_mol, cumulene_mol
 from ..testutils import requires_working_xtb_install
 
 
@@ -366,3 +366,15 @@ def test_optimise_linear_bend_with_ref():
     opt = CRFOptimiser(maxiter=15, gtol=1e-4, etol=1e-5)
     opt.run(mol, XTB())
     assert mol.angle(2, 3, 4) > val.Angle(178, "deg")
+
+
+@requires_working_xtb_install
+@work_in_tmp_dir()
+def test_optimise_chain_dihedrals():
+    mol = cumulene_mol()
+    assert abs(mol.dihedral(6, 3, 4, 8)) < val.Angle(40, "deg")
+    opt = CRFOptimiser(maxiter=15, gtol=1e-4, etol=1e-4)
+    opt.run(mol, XTB())
+    # 5-C chain, should be close to 90 degrees
+    assert abs(mol.dihedral(6, 3, 4, 8)) > val.Angle(85, "deg")
+    assert abs(mol.dihedral(6, 3, 4, 8)) < val.Angle(95, "deg")

From a18f9d26ab35fea206bd089eec92bd4fb22a0623 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 28 Dec 2023 15:16:26 +0000
Subject: [PATCH 35/47] fix tests

---
 autode/atoms.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/autode/atoms.py b/autode/atoms.py
index 916ce8e19..6f4034dc3 100644
--- a/autode/atoms.py
+++ b/autode/atoms.py
@@ -1043,7 +1043,7 @@ def angle(self, i: int, j: int, k: int) -> Angle:
 
         # Catch errors due to incomplete float precision
         cos_value = np.dot(vec1, vec2) / norms
-        if -1 < cos_value < 1:
+        if -1 <= cos_value <= 1:
             pass
         elif cos_value > 1 and np.isclose(cos_value, 1, rtol=1e-8):
             cos_value = 1

From 800089b16015e741ba90bd67008673e1317c8456 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Fri, 29 Dec 2023 20:10:12 +0000
Subject: [PATCH 36/47] pr suggestions 1

---
 autode/atoms.py                     | 16 +++---------
 autode/mol_graphs.py                |  2 +-
 autode/opt/coordinates/internals.py | 38 +++++++++++++++++------------
 tests/test_opt/test_coordiantes.py  | 13 +++++++---
 tests/test_opt/test_crfo.py         |  6 ++---
 5 files changed, 38 insertions(+), 37 deletions(-)

diff --git a/autode/atoms.py b/autode/atoms.py
index 6f4034dc3..6e88d8f9f 100644
--- a/autode/atoms.py
+++ b/autode/atoms.py
@@ -1041,20 +1041,10 @@ def angle(self, i: int, j: int, k: int) -> Angle:
                 f"least one zero vector"
             )
 
-        # Catch errors due to incomplete float precision
-        cos_value = np.dot(vec1, vec2) / norms
-        if -1 <= cos_value <= 1:
-            pass
-        elif cos_value > 1 and np.isclose(cos_value, 1, rtol=1e-8):
-            cos_value = 1
-        elif cos_value < -1 and np.isclose(cos_value, -1, rtol=1e-8):
-            cos_value = -1
-        else:
-            raise ValueError("Cos(angle) must be in [-1, 1] range!")
-
-        value = np.arccos(cos_value)
+        # Cos(theta) must lie within [-1, 1]
+        cos_value = np.clip(np.dot(vec1, vec2) / norms, a_min=-1, a_max=1)
 
-        return Angle(value)
+        return Angle(np.arccos(cos_value))
 
     def dihedral(self, w: int, x: int, y: int, z: int) -> Angle:
         r"""
diff --git a/autode/mol_graphs.py b/autode/mol_graphs.py
index c01ae1b7b..521036a90 100644
--- a/autode/mol_graphs.py
+++ b/autode/mol_graphs.py
@@ -120,7 +120,7 @@ def is_connected(self) -> bool:
         """Is this graph fully connected (i.e. not separate parts)"""
         return nx.is_connected(self)
 
-    def get_components(self):
+    def connected_components(self):
         """Generate the separate connected components"""
         return nx.connected_components(self)
 
diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 71cdec06a..8b772f8c1 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -94,13 +94,19 @@ def __init__(self, *args: Any):
                 f"from {args}. Must be primitive internals"
             )
 
-    def append(self, item: Primitive) -> None:
-        """Append an item to this set of primitives"""
-        assert isinstance(item, Primitive), "Must be a Primitive type!"
-        # prevent duplicate primitives
+    def add(self, item: Primitive) -> None:
+        """Add a primitive to this set of primitive coordinates"""
+        assert isinstance(item, Primitive), "Must be a primitive"
+        # prevent duplication of primitives
         if item not in self:
             super().append(item)
 
+    def append(self, item: Primitive) -> None:
+        """Append an item to this set of primitives"""
+        raise NotImplementedError(
+            "Please use PIC.add() to add new primitives to the set"
+        )
+
     @property
     def B(self) -> np.ndarray:
         """Wilson B matrix"""
@@ -224,7 +230,7 @@ def _populate_all(self, x: np.ndarray):
         # Add all the unique inverse distances (i < j)
         for i in range(n_atoms):
             for j in range(i + 1, n_atoms):
-                self.append(self._primitive_type(i, j))
+                self.add(self._primitive_type(i, j))
 
         return None
 
@@ -305,7 +311,7 @@ def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
 
     # join disconnected graph components
     if not core_graph.is_connected:
-        components = core_graph.get_components()
+        components = core_graph.connected_components()
         for comp_i, comp_j in itertools.combinations(components, r=2):
             min_dist = float("inf")
             min_pair = (-1, -1)
@@ -353,10 +359,10 @@ def _add_bonds_from_species(
             and (i, j) in mol.constraints.distance
         ):
             r = mol.constraints.distance[(i, j)]
-            pic.append(ConstrainedPrimitiveDistance(i, j, r))
+            pic.add(ConstrainedPrimitiveDistance(i, j, r))
             n += 1
         else:
-            pic.append(PrimitiveDistance(i, j))
+            pic.add(PrimitiveDistance(i, j))
     assert n == mol.constraints.n_distance
 
     return None
@@ -419,24 +425,24 @@ def get_ref_atom(a, b, c, bonded=False):
     for o in range(mol.n_atoms):
         for n, m in itertools.combinations(core_graph.neighbors(o), r=2):
             if mol.angle(m, o, n) < lin_thresh:
-                pic.append(PrimitiveBondAngle(m=m, o=o, n=n))
+                pic.add(PrimitiveBondAngle(m=m, o=o, n=n))
             else:
                 # If central atom is connected to another atom, then the
                 # linear angle is skipped and instead an out-of-plane
                 # (improper dihedral) coordinate is used
                 r = get_ref_atom(m, o, n, bonded=True)
                 if r is not None:
-                    pic.append(PrimitiveDihedralAngle(m, r, o, n))
+                    pic.add(PrimitiveDihedralAngle(m, r, o, n))
                     continue
 
                 # Otherwise, we use a nearby (< 4.0 A) reference atom to
                 # define two orthogonal linear bends
                 r = get_ref_atom(m, o, n, bonded=False)
                 if r is not None:
-                    pic.append(
+                    pic.add(
                         PrimitiveLinearAngle(m, o, n, r, LinearBendType.BEND)
                     )
-                    pic.append(
+                    pic.add(
                         PrimitiveLinearAngle(
                             m, o, n, r, LinearBendType.COMPLEMENT
                         )
@@ -445,10 +451,10 @@ def get_ref_atom(a, b, c, bonded=False):
                 # For completely linear molecules (CO2), there will be no such
                 # reference atoms, so use dummy atoms instead
                 else:
-                    pic.append(
+                    pic.add(
                         PrimitiveDummyLinearAngle(m, o, n, LinearBendType.BEND)
                     )
-                    pic.append(
+                    pic.add(
                         PrimitiveDummyLinearAngle(
                             m, o, n, LinearBendType.COMPLEMENT
                         )
@@ -507,7 +513,7 @@ def is_dihedral_well_defined(w, x, y, z):
                     continue
 
                 if is_dihedral_well_defined(m, o, p, n):
-                    pic.append(PrimitiveDihedralAngle(m, o, p, n))
+                    pic.add(PrimitiveDihedralAngle(m, o, p, n))
 
     # find all linear atom chains A--B--C--D... and add dihedrals to terminal atoms
 
@@ -559,5 +565,5 @@ def extend_chain(chain: List[int]):
                     continue
 
                 if is_dihedral_well_defined(m, o, p, n):
-                    pic.append(PrimitiveDihedralAngle(m, o, p, n))
+                    pic.add(PrimitiveDihedralAngle(m, o, p, n))
     return None
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index eb4ce675f..42e56673e 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -544,12 +544,17 @@ def test_pic_b_no_primitives():
         c._calc_B(np.arange(6, dtype=float).reshape(2, 3))
 
 
-def test_pic_append_type_checking():
+def test_pic_add_sanity_checking():
     c = PIC()
-    # append should check for primitive type
-    c.append(PrimitiveDistance(0, 1))
+    # pic add should check for primitive type
+    c.add(PrimitiveDistance(0, 1))
     with pytest.raises(AssertionError):
-        c.append(3)
+        c.add(3)
+
+    # pic should not allow duplicate coordinates to be added
+    assert len(c) == 1
+    c.add(PrimitiveDistance(1, 0))
+    assert len(c) == 1
 
 
 def test_constrained_distance_satisfied():
diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index ab9aa0ab8..e1060592c 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -205,7 +205,7 @@ def test_baker1997_example():
 
     pic = PIC(r1, r2, theta)
     for pair in ((2, 0), (3, 0), (4, 1), (5, 1)):
-        pic.append(prim.PrimitiveDistance(*pair))
+        pic.add(prim.PrimitiveDistance(*pair))
 
     for triple in (
         (0, 1, 3),
@@ -215,10 +215,10 @@ def test_baker1997_example():
         (1, 5, 0),
         (1, 4, 5),
     ):
-        pic.append(prim.PrimitiveBondAngle(*triple))
+        pic.add(prim.PrimitiveBondAngle(*triple))
 
     for quadruple in ((4, 1, 0, 2), (4, 1, 0, 3), (5, 1, 0, 2), (5, 1, 0, 3)):
-        pic.append(prim.PrimitiveDihedralAngle(*quadruple))
+        pic.add(prim.PrimitiveDihedralAngle(*quadruple))
 
     dic = DICWithConstraints.from_cartesian(
         x=CartesianCoordinates(c2h3f.coordinates), primitives=pic

From e51ca4d76723588659238a4128a64e88aa7d8533 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 30 Dec 2023 11:40:13 +0000
Subject: [PATCH 37/47] pr suggestions 2 partial

---
 autode/opt/coordinates/internals.py | 316 +++++++++++++++-------------
 autode/opt/optimisers/crfo.py       |   4 +-
 2 files changed, 174 insertions(+), 146 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 8b772f8c1..964b18cdd 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -255,28 +255,181 @@ class AnyPIC(PIC):
     def _populate_all(self, x: np.ndarray) -> None:
         raise RuntimeError("Cannot populate all on an AnyPIC instance")
 
+    @classmethod
+    def from_species(cls, mol: "Species") -> "AnyPIC":
+        """
+        Build a set of primitives from the species, using the graph as
+        a starting point for the connectivity of the species. Also joins
+        any disjoint parts of the graph, and adds hydrogen bonds to
+        ensure that the primitives are redundant.
 
-def build_pic_from_species(
-    mol: "Species",
-) -> AnyPIC:
-    """
-    Build a set of primitives from the species, using the graph as
-    a starting point for the connectivity of the species. Also joins
-    any disjoint parts of the graph, and adds hydrogen bonds to
-    ensure that the primitives are redundant
+        Args:
+            mol: The species object
 
-    Args:
-        mol:
+        Returns:
+            (AnyPIC): The set of primitive internals
+        """
+        pic = cls()
+        core_graph = _get_connected_graph_from_species(mol)
+        pic._add_bonds_from_species(mol, core_graph)
+        pic._add_angles_from_species(mol, core_graph)
+        _add_dihedrals_from_species(pic, mol, core_graph)
+        return pic
 
-    Returns:
-        (AnyPIC): The set of primitive internals
-    """
-    pic = AnyPIC()
-    core_graph = _get_connected_graph_from_species(mol)
-    _add_bonds_from_species(pic, mol, core_graph)
-    _add_angles_from_species(pic, mol, core_graph)
-    _add_dihedrals_from_species(pic, mol, core_graph)
-    return pic
+    def _add_bonds_from_species(
+        self,
+        mol: "Species",
+        core_graph: "MolecularGraph",
+    ):
+        """
+        Add bonds to the current set of primitives, from the
+        connectivity graph supplied
+
+        Args:
+            mol: The species object
+            core_graph: The connectivity graph
+        """
+        n = 0
+        for i, j in sorted(core_graph.edges):
+            if (
+                mol.constraints.distance is not None
+                and (i, j) in mol.constraints.distance
+            ):
+                r = mol.constraints.distance[(i, j)]
+                self.add(ConstrainedPrimitiveDistance(i, j, r))
+                n += 1
+            else:
+                self.add(PrimitiveDistance(i, j))
+        assert n == mol.constraints.n_distance
+
+        return None
+
+    @staticmethod
+    def _get_ref_for_linear_angle(
+        mol,
+        core_graph,
+        lin_thresh,
+        a,
+        b,
+        c,
+        bonded: bool,
+        dist_thresh=Distance(4, "ang"),
+    ) -> Optional[int]:
+        """
+        Get a reference atom for describing a linear angle, which
+        must not itself be linear to the atoms in the angle in
+        any combination. The linear angle is a--b--c here.
+
+        Args:
+            mol:
+            core_graph:
+            lin_thresh:
+            a:
+            b:
+            c:
+            bonded: Whether to look for only atoms bonded to the central
+                    atom (b) for reference
+
+        Returns:
+            (int|None): The index of the ref. atom if found, else None
+        """
+        # only check bonded atoms if requested
+        if bonded:
+            near_atoms = list(core_graph.neighbors(b))
+            near_atoms.remove(a)
+            near_atoms.remove(c)
+
+        # otherwise get all atoms in 4 A radius except a, b, c
+        else:
+            near_atoms = [
+                idx
+                for idx in range(mol.n_atoms)
+                if mol.distance(b, idx) < dist_thresh and idx not in (a, b, c)
+            ]
+
+        # get atoms closest to perpendicular
+        deviations_from_90 = {}
+        for atom in near_atoms:
+            i_b_a = mol.angle(atom, b, a)
+            if i_b_a > lin_thresh or i_b_a < (np.pi - lin_thresh):
+                continue
+            i_b_c = mol.angle(atom, b, c)
+            if i_b_c > lin_thresh or i_b_c < (np.pi - lin_thresh):
+                continue
+            deviation_a = abs(i_b_a - np.pi / 2)
+            deviation_c = abs(i_b_c - np.pi / 2)
+            avg_dev = (deviation_a + deviation_c) / 2
+            deviations_from_90[atom] = avg_dev
+
+        if len(deviations_from_90) == 0:
+            return None
+
+        return min(deviations_from_90, key=deviations_from_90.get)  # type: ignore
+
+    def _add_angles_from_species(
+        self,
+        mol: "Species",
+        core_graph: "MolecularGraph",
+        lin_thresh: Angle = Angle(170, "deg"),
+    ) -> None:
+        """
+        Modify the set of primitives in-place by adding angles, from the
+        connectivity graph supplied
+
+        Args:
+            mol (Species): The species object
+            core_graph (MolecularGraph): The connectivity graph
+            lin_thresh (Angle): The angle threshold for linearity
+        """
+        lin_thresh = lin_thresh.to("rad")
+
+        for o in range(mol.n_atoms):
+            for n, m in itertools.combinations(core_graph.neighbors(o), r=2):
+                if mol.angle(m, o, n) < lin_thresh:
+                    self.add(PrimitiveBondAngle(m=m, o=o, n=n))
+                else:
+                    # If central atom is connected to another atom, then the
+                    # linear angle is skipped and instead an out-of-plane
+                    # (improper dihedral) coordinate is used
+                    r = self._get_ref_for_linear_angle(
+                        mol, core_graph, lin_thresh, m, o, n, bonded=True
+                    )
+                    if r is not None:
+                        self.add(PrimitiveDihedralAngle(m, r, o, n))
+                        continue
+
+                    # Otherwise, we use a nearby (< 4.0 A) reference atom to
+                    # define two orthogonal linear bends
+                    r = self._get_ref_for_linear_angle(
+                        mol, core_graph, lin_thresh, m, o, n, bonded=False
+                    )
+                    if r is not None:
+                        self.add(
+                            PrimitiveLinearAngle(
+                                m, o, n, r, LinearBendType.BEND
+                            )
+                        )
+                        self.add(
+                            PrimitiveLinearAngle(
+                                m, o, n, r, LinearBendType.COMPLEMENT
+                            )
+                        )
+
+                    # For completely linear molecules (CO2), there will be no such
+                    # reference atoms, so use dummy atoms instead
+                    else:
+                        self.add(
+                            PrimitiveDummyLinearAngle(
+                                m, o, n, LinearBendType.BEND
+                            )
+                        )
+                        self.add(
+                            PrimitiveDummyLinearAngle(
+                                m, o, n, LinearBendType.COMPLEMENT
+                            )
+                        )
+
+        return None
 
 
 def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
@@ -338,131 +491,6 @@ def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
     return core_graph
 
 
-def _add_bonds_from_species(
-    pic: AnyPIC,
-    mol: "Species",
-    core_graph: "MolecularGraph",
-):
-    """
-    Modify the supplied AnyPIC instance in-place by adding bonds, from the
-    connectivity graph supplied
-
-    Args:
-        pic: The AnyPIC instance (modified in-place)
-        mol: The species object
-        core_graph: The connectivity graph
-    """
-    n = 0
-    for i, j in sorted(core_graph.edges):
-        if (
-            mol.constraints.distance is not None
-            and (i, j) in mol.constraints.distance
-        ):
-            r = mol.constraints.distance[(i, j)]
-            pic.add(ConstrainedPrimitiveDistance(i, j, r))
-            n += 1
-        else:
-            pic.add(PrimitiveDistance(i, j))
-    assert n == mol.constraints.n_distance
-
-    return None
-
-
-def _add_angles_from_species(
-    pic: AnyPIC,
-    mol: "Species",
-    core_graph: "MolecularGraph",
-    lin_thresh: Angle = Angle(170, "deg"),
-) -> None:
-    """
-    Modify the set of primitives in-place by adding angles, from the
-    connectivity graph supplied
-
-    Args:
-        pic (AnyPIC): The AnyPIC instance (modified in-place)
-        mol (Species): The species object
-        core_graph (MolecularGraph): The connectivity graph
-        lin_thresh (Angle): The angle threshold for linearity
-    """
-    lin_thresh = lin_thresh.to("rad")
-
-    def get_ref_atom(a, b, c, bonded=False):
-        """get a reference atom for a-b-c linear angle"""
-        # only check bonded atoms if requested
-        if bonded:
-            near_atoms = list(core_graph.neighbors(b))
-            near_atoms.remove(a)
-            near_atoms.remove(c)
-
-        # otherwise get all atoms in 4 A radius except a, b, c
-        else:
-            near_atoms = [
-                idx
-                for idx in range(mol.n_atoms)
-                if mol.distance(b, idx) < Distance(4.0, "ang")
-                and idx not in (a, b, c)
-            ]
-
-        # get atoms closest to perpendicular
-        deviations_from_90 = {}
-        for atom in near_atoms:
-            i_b_a = mol.angle(atom, b, a)
-            if i_b_a > lin_thresh or i_b_a < (np.pi - lin_thresh):
-                continue
-            i_b_c = mol.angle(atom, b, c)
-            if i_b_c > lin_thresh or i_b_c < (np.pi - lin_thresh):
-                continue
-            deviation_a = abs(i_b_a - np.pi / 2)
-            deviation_c = abs(i_b_c - np.pi / 2)
-            avg_dev = (deviation_a + deviation_c) / 2
-            deviations_from_90[atom] = avg_dev
-
-        if len(deviations_from_90) == 0:
-            return None
-
-        return min(deviations_from_90, key=deviations_from_90.get)
-
-    for o in range(mol.n_atoms):
-        for n, m in itertools.combinations(core_graph.neighbors(o), r=2):
-            if mol.angle(m, o, n) < lin_thresh:
-                pic.add(PrimitiveBondAngle(m=m, o=o, n=n))
-            else:
-                # If central atom is connected to another atom, then the
-                # linear angle is skipped and instead an out-of-plane
-                # (improper dihedral) coordinate is used
-                r = get_ref_atom(m, o, n, bonded=True)
-                if r is not None:
-                    pic.add(PrimitiveDihedralAngle(m, r, o, n))
-                    continue
-
-                # Otherwise, we use a nearby (< 4.0 A) reference atom to
-                # define two orthogonal linear bends
-                r = get_ref_atom(m, o, n, bonded=False)
-                if r is not None:
-                    pic.add(
-                        PrimitiveLinearAngle(m, o, n, r, LinearBendType.BEND)
-                    )
-                    pic.add(
-                        PrimitiveLinearAngle(
-                            m, o, n, r, LinearBendType.COMPLEMENT
-                        )
-                    )
-
-                # For completely linear molecules (CO2), there will be no such
-                # reference atoms, so use dummy atoms instead
-                else:
-                    pic.add(
-                        PrimitiveDummyLinearAngle(m, o, n, LinearBendType.BEND)
-                    )
-                    pic.add(
-                        PrimitiveDummyLinearAngle(
-                            m, o, n, LinearBendType.COMPLEMENT
-                        )
-                    )
-
-    return None
-
-
 def _add_dihedrals_from_species(
     pic: AnyPIC,
     mol: "Species",
diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py
index 5c86640db..038f99666 100644
--- a/autode/opt/optimisers/crfo.py
+++ b/autode/opt/optimisers/crfo.py
@@ -11,7 +11,7 @@
 from autode.log import logger
 from autode.values import GradientRMS, Distance
 from autode.opt.coordinates import CartesianCoordinates, DICWithConstraints
-from autode.opt.coordinates.internals import build_pic_from_species
+from autode.opt.coordinates.internals import AnyPIC
 from autode.opt.optimisers.rfo import RFOptimiser
 from autode.opt.optimisers.hessian_update import (
     BFGSDampedUpdate,
@@ -122,7 +122,7 @@ def _build_internal_coordinates(self):
             )
 
         cartesian_coords = CartesianCoordinates(self._species.coordinates)
-        primitives = build_pic_from_species(self._species)
+        primitives = AnyPIC.from_species(self._species)
 
         self._coords = DICWithConstraints.from_cartesian(
             x=cartesian_coords, primitives=primitives

From 7fa4454fd603220a6d63f58804f96d4fa0d4ec43 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 30 Dec 2023 11:52:38 +0000
Subject: [PATCH 38/47] test fix

---
 tests/test_opt/test_coordiantes.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 42e56673e..8e4b88741 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -19,7 +19,7 @@
 from autode.opt.coordinates.internals import (
     PrimitiveInverseDistances,
     PIC,
-    build_pic_from_species,
+    AnyPIC,
 )
 from autode.opt.coordinates.cartesian import CartesianCoordinates
 from autode.opt.coordinates.dic import DIC
@@ -545,7 +545,7 @@ def test_pic_b_no_primitives():
 
 
 def test_pic_add_sanity_checking():
-    c = PIC()
+    c = AnyPIC()
     # pic add should check for primitive type
     c.add(PrimitiveDistance(0, 1))
     with pytest.raises(AssertionError):
@@ -832,7 +832,7 @@ def test_dics_cannot_be_built_with_incomplete_primitives():
 def test_pic_generation_linear_angle_ref():
     # Fe(CO)5 with linear Fe-C-O bonds
     m = feco5_mol()
-    pic = build_pic_from_species(m)
+    pic = AnyPIC.from_species(m)
 
     # check that there are no duplicates
     assert not any(ic1 == ic2 for ic1, ic2 in itertools.combinations(pic, r=2))
@@ -848,7 +848,7 @@ def test_pic_generation_linear_angle_ref():
 def test_pic_generation_linear_angle_dummy():
     # acetylene molecule
     mol = acetylene_mol()
-    pic = build_pic_from_species(mol)
+    pic = AnyPIC.from_species(mol)
 
     # there should not be any usual bond angles
     assert not any(isinstance(ic, PrimitiveBondAngle) for ic in pic)
@@ -889,7 +889,7 @@ def test_pic_generation_disjoint_graph():
 
     mol = Molecule("diels_alder_complex.xyz")
     assert not mol.graph.is_connected
-    pic = build_pic_from_species(mol)
+    pic = AnyPIC.from_species(mol)
 
     # shortest bond is between 4, 5 which should also generate angle, torsion
     assert PrimitiveDistance(4, 5) in pic
@@ -906,7 +906,7 @@ def test_pic_generation_disjoint_graph():
 
     # if the bond between 2, 3 is made into a constraint, it will generate angles
     mol.constraints.distance = {(2, 3): mol.distance(2, 3)}
-    pic = build_pic_from_species(mol)
+    pic = AnyPIC.from_species(mol)
     assert ConstrainedPrimitiveDistance(2, 3, mol.distance(2, 3)) in pic
     assert PrimitiveBondAngle(1, 2, 3) in pic
 
@@ -914,7 +914,7 @@ def test_pic_generation_disjoint_graph():
 def test_pic_generation_chain_dihedrals():
     # extra dihedrals are needed for ends of linear chains like allene
     cumulene = cumulene_mol()
-    pic = build_pic_from_species(cumulene)
+    pic = AnyPIC.from_species(cumulene)
 
     assert PrimitiveDihedralAngle(5, 3, 4, 8) in pic
     assert PrimitiveDihedralAngle(6, 3, 4, 7) in pic
@@ -940,6 +940,6 @@ def test_pic_generation_square_planar():
 
     # for sq planar, out-of-plane dihedrals are needed to have
     # all degrees of freedom
-    pic = build_pic_from_species(ptcl4)
+    pic = AnyPIC.from_species(ptcl4)
     _ = pic(ptcl4.coordinates.flatten())
     assert np.linalg.matrix_rank(pic.B) == 3 * 5 - 6

From 6887a875744d23c58290993f4b8423dd635309df Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 30 Dec 2023 20:03:27 +0000
Subject: [PATCH 39/47] put methods into anypic class

---
 autode/opt/coordinates/internals.py | 297 ++++++++++++++--------------
 1 file changed, 148 insertions(+), 149 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 964b18cdd..1325d5755 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -270,12 +270,71 @@ def from_species(cls, mol: "Species") -> "AnyPIC":
             (AnyPIC): The set of primitive internals
         """
         pic = cls()
-        core_graph = _get_connected_graph_from_species(mol)
+        core_graph = pic._get_connected_graph_from_species(mol)
         pic._add_bonds_from_species(mol, core_graph)
         pic._add_angles_from_species(mol, core_graph)
-        _add_dihedrals_from_species(pic, mol, core_graph)
+        pic._add_dihedrals_from_species(mol, core_graph)
         return pic
 
+    @staticmethod
+    def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
+        """
+        Creates a fully connected graph from the graph of a species, by
+        (1) joining disconnected fragments by their shortest distance,
+        (2) connecting constrained bonds, (3) joining hydrogen bonds,
+        if present.
+
+        Args:
+            mol: A species (must have atoms and graph)
+
+        Returns:
+            (MolecularGraph):
+        """
+        assert mol.graph is not None, "Species must have graph!"
+        core_graph = copy.deepcopy(mol.graph)
+
+        # join hydrogen bonds
+        h_bond_x = ["N", "O", "F", "P", "S", "Cl"]
+        for i, j in itertools.combinations(range(mol.n_atoms), r=2):
+            if (
+                mol.atoms[i].label in h_bond_x
+                and mol.atoms[j].label == "H"
+                or mol.atoms[j].label in h_bond_x
+                and mol.atoms[i].label == "H"
+            ):
+                vdw_sum = mol.atoms[i].vdw_radius + mol.atoms[j].vdw_radius
+                if mol.distance(i, j) < 0.9 * vdw_sum:
+                    if not core_graph.has_edge(i, j):
+                        core_graph.add_edge(i, j, pi=False, active=False)
+
+        # join disconnected graph components
+        if not core_graph.is_connected:
+            components = core_graph.connected_components()
+            for comp_i, comp_j in itertools.combinations(components, r=2):
+                min_dist = float("inf")
+                min_pair = (-1, -1)
+                for i, j in itertools.product(list(comp_i), list(comp_j)):
+                    if mol.distance(i, j) < min_dist:
+                        min_dist = mol.distance(i, j)
+                        min_pair = (i, j)
+                # avoid connecting distant components
+                if min_dist < Distance(4.0, "ang"):
+                    core_graph.add_edge(*min_pair, pi=False, active=False)
+
+        if not core_graph.is_connected:
+            raise RuntimeError(
+                "Unable to join all the fragments, distance between "
+                "one or more pairs of fragments is too high (>4.0 Å)"
+            )
+
+        # The constraints should be counted as bonds
+        if mol.constraints.distance is not None:
+            for i, j in mol.constraints.distance:
+                if not core_graph.has_edge(i, j):
+                    core_graph.add_edge(i, j, pi=False, active=False)
+
+        return core_graph
+
     def _add_bonds_from_species(
         self,
         mol: "Species",
@@ -329,6 +388,7 @@ def _get_ref_for_linear_angle(
             c:
             bonded: Whether to look for only atoms bonded to the central
                     atom (b) for reference
+            dist_thresh: The distance threshold to connect
 
         Returns:
             (int|None): The index of the ref. atom if found, else None
@@ -431,167 +491,106 @@ def _add_angles_from_species(
 
         return None
 
+    def _add_dihedrals_from_species(
+        self,
+        mol: "Species",
+        core_graph: "MolecularGraph",
+        lin_thresh: Angle = Angle(170, "deg"),
+    ) -> None:
+        """
+        Modify the set of primitives in-place by adding dihedrals (torsions),
+        from the connectivity graph supplied
 
-def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
-    """
-    Creates a fully connected graph from the graph of a species, by
-    (1) joining disconnected fragments by their shortest distance,
-    (2) connecting constrained bonds, (3) joining hydrogen bonds,
-    if present.
-
-    Args:
-        mol: A species (must have atoms and graph)
-
-    Returns:
-        (MolecularGraph):
-    """
-    assert mol.graph is not None, "Species must have graph!"
-    core_graph = copy.deepcopy(mol.graph)
-
-    # join hydrogen bonds
-    h_bond_x = ["N", "O", "F", "P", "S", "Cl"]
-    for i, j in itertools.combinations(range(mol.n_atoms), r=2):
-        if (
-            mol.atoms[i].label in h_bond_x
-            and mol.atoms[j].label == "H"
-            or mol.atoms[j].label in h_bond_x
-            and mol.atoms[i].label == "H"
-        ):
-            vdw_sum = mol.atoms[i].vdw_radius + mol.atoms[j].vdw_radius
-            if mol.distance(i, j) < 0.9 * vdw_sum:
-                if not core_graph.has_edge(i, j):
-                    core_graph.add_edge(i, j, pi=False, active=False)
-
-    # join disconnected graph components
-    if not core_graph.is_connected:
-        components = core_graph.connected_components()
-        for comp_i, comp_j in itertools.combinations(components, r=2):
-            min_dist = float("inf")
-            min_pair = (-1, -1)
-            for i, j in itertools.product(list(comp_i), list(comp_j)):
-                if mol.distance(i, j) < min_dist:
-                    min_dist = mol.distance(i, j)
-                    min_pair = (i, j)
-            # avoid connecting distant components
-            if min_dist < Distance(4.0, "ang"):
-                core_graph.add_edge(*min_pair, pi=False, active=False)
-
-    if not core_graph.is_connected:
-        raise RuntimeError(
-            "Unable to join all the fragments, distance between "
-            "one or more pairs of fragments is too high (>4.0 Å)"
-        )
-
-    # The constraints should be counted as bonds
-    if mol.constraints.distance is not None:
-        for i, j in mol.constraints.distance:
-            if not core_graph.has_edge(i, j):
-                core_graph.add_edge(i, j, pi=False, active=False)
-
-    return core_graph
-
-
-def _add_dihedrals_from_species(
-    pic: AnyPIC,
-    mol: "Species",
-    core_graph: "MolecularGraph",
-    lin_thresh: Angle = Angle(170, "deg"),
-) -> None:
-    """
-    Modify the set of primitives in-place by adding dihedrals (torsions),
-    from the connectivity graph supplied
-
-    Args:
-        pic (AnyPIC): The AnyPIC instance (modified in-place)
-        mol (Species): The species
-        core_graph (MolecularGraph): The connectivity graph
-        lin_thresh (Angle): The threshold for linearity
-    """
-    # no dihedrals possible with less than 4 atoms
-    if mol.n_atoms < 4:
-        return
-
-    lin_thresh = lin_thresh.to("rad")
-    zero_angle_thresh = np.pi - lin_thresh
+        Args:
+            mol (Species): The species
+            core_graph (MolecularGraph): The connectivity graph
+            lin_thresh (Angle): The threshold for linearity
+        """
+        # no dihedrals possible with less than 4 atoms
+        if mol.n_atoms < 4:
+            return
 
-    def is_dihedral_well_defined(w, x, y, z):
-        """A dihedral is well-defined if any angle is not linear"""
-        is_linear_1 = (
-            mol.angle(w, x, y) > lin_thresh
-            or mol.angle(w, x, y) < zero_angle_thresh
-        )
-        is_linear_2 = (
-            mol.angle(x, y, z) > lin_thresh
-            or mol.angle(x, y, z) < zero_angle_thresh
-        )
-        return not (is_linear_1 or is_linear_2)
+        lin_thresh = lin_thresh.to("rad")
+        zero_angle_thresh = np.pi - lin_thresh
 
-    # add normal dihedrals
-    for o, p in list(core_graph.edges):
-        for m in core_graph.neighbors(o):
-            if m == p:
-                continue
+        def is_dihedral_well_defined(w, x, y, z):
+            """A dihedral is well-defined if any angle is not linear"""
+            is_linear_1 = (
+                mol.angle(w, x, y) > lin_thresh
+                or mol.angle(w, x, y) < zero_angle_thresh
+            )
+            is_linear_2 = (
+                mol.angle(x, y, z) > lin_thresh
+                or mol.angle(x, y, z) < zero_angle_thresh
+            )
+            return not (is_linear_1 or is_linear_2)
 
-            for n in core_graph.neighbors(p):
-                if n == o:
+        # add normal dihedrals
+        for o, p in list(core_graph.edges):
+            for m in core_graph.neighbors(o):
+                if m == p:
                     continue
 
-                # avoid triangle rings like cyclopropane
-                if n == m:
-                    continue
+                for n in core_graph.neighbors(p):
+                    if n == o:
+                        continue
 
-                if is_dihedral_well_defined(m, o, p, n):
-                    pic.add(PrimitiveDihedralAngle(m, o, p, n))
+                    # avoid triangle rings like cyclopropane
+                    if n == m:
+                        continue
 
-    # find all linear atom chains A--B--C--D... and add dihedrals to terminal atoms
+                    if is_dihedral_well_defined(m, o, p, n):
+                        self.add(PrimitiveDihedralAngle(m, o, p, n))
 
-    def extend_chain(chain: List[int]):
-        for idx in range(mol.n_atoms):
-            if idx in chain:
-                continue
-            # if idx -- 0 -- 1 > 170 degrees
-            if mol.angle(chain[1], chain[0], idx) > lin_thresh:
-                chain.insert(0, idx)
-                continue
-            # if (-2) -- (-1) -- idx > 170 degrees
-            if mol.angle(chain[-2], chain[-1], idx) > lin_thresh:
-                chain.append(idx)
-                continue
+        # find all linear atom chains A--B--C--D... and add dihedrals to terminal atoms
 
-    linear_chains: List[list] = []
-    for b in range(mol.n_atoms):
-        if any(b in chain for chain in linear_chains):
-            continue
-        for a, c in itertools.combinations(core_graph.neighbors(b), r=2):
-            if any(a in chain for chain in linear_chains) or any(
-                c in chain for chain in linear_chains
-            ):
-                continue
-            if mol.angle(a, b, c) > lin_thresh:
-                chain = [a, b, c]
-                extend_chain(chain)
-                linear_chains.append(chain)
-
-    # add linear chain dihedrals
-    for chain in linear_chains:
-        o, p = chain[0], chain[-1]
-        for m in core_graph.neighbors(o):
-            if m == p:
-                continue
+        def extend_chain(chain: List[int]):
+            for idx in range(mol.n_atoms):
+                if idx in chain:
+                    continue
+                # if idx -- 0 -- 1 > 170 degrees
+                if mol.angle(chain[1], chain[0], idx) > lin_thresh:
+                    chain.insert(0, idx)
+                    continue
+                # if (-2) -- (-1) -- idx > 170 degrees
+                if mol.angle(chain[-2], chain[-1], idx) > lin_thresh:
+                    chain.append(idx)
+                    continue
 
-            if m in chain:
+        linear_chains: List[list] = []
+        for b in range(mol.n_atoms):
+            if any(b in chain for chain in linear_chains):
                 continue
-
-            for n in core_graph.neighbors(p):
-                if n == o:
+            for a, c in itertools.combinations(core_graph.neighbors(b), r=2):
+                if any(a in chain for chain in linear_chains) or any(
+                    c in chain for chain in linear_chains
+                ):
                     continue
-
-                if n == m:
+                if mol.angle(a, b, c) > lin_thresh:
+                    chain = [a, b, c]
+                    extend_chain(chain)
+                    linear_chains.append(chain)
+
+        # add linear chain dihedrals
+        for chain in linear_chains:
+            o, p = chain[0], chain[-1]
+            for m in core_graph.neighbors(o):
+                if m == p:
                     continue
 
-                if n in chain:
+                if m in chain:
                     continue
 
-                if is_dihedral_well_defined(m, o, p, n):
-                    pic.add(PrimitiveDihedralAngle(m, o, p, n))
-    return None
+                for n in core_graph.neighbors(p):
+                    if n == o:
+                        continue
+
+                    if n == m:
+                        continue
+
+                    if n in chain:
+                        continue
+
+                    if is_dihedral_well_defined(m, o, p, n):
+                        self.add(PrimitiveDihedralAngle(m, o, p, n))
+        return None

From cff5079f6834a4081b31e6e9129bd1e1b406080d Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 30 Dec 2023 20:19:47 +0000
Subject: [PATCH 40/47] pr update 3

---
 autode/opt/coordinates/internals.py | 166 ++++++++++++++--------------
 1 file changed, 80 insertions(+), 86 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 1325d5755..13012a139 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -9,8 +9,6 @@
 
 Set-up of redundant primitives is based on J. Chem. Phys., 117, 2002, 9160
 """
-import copy
-
 import numpy as np
 import itertools
 from typing import Any, Optional, Type, List, TYPE_CHECKING
@@ -31,7 +29,6 @@
 
 if TYPE_CHECKING:
     from autode.species import Species
-    from autode.mol_graphs import MolecularGraph
     from autode.opt.coordinates.cartesian import CartesianCoordinates
     from autode.opt.coordinates.primitives import (
         ConstrainedPrimitive,
@@ -45,6 +42,9 @@ def __new__(cls, input_array) -> "InternalCoordinates":
 
         arr = super().__new__(cls, input_array, units="Å")
 
+        arr._x = None
+        arr.primitives = None
+
         for attr in ("_x", "primitives"):
             setattr(arr, attr, getattr(input_array, attr, None))
 
@@ -270,86 +270,29 @@ def from_species(cls, mol: "Species") -> "AnyPIC":
             (AnyPIC): The set of primitive internals
         """
         pic = cls()
-        core_graph = pic._get_connected_graph_from_species(mol)
-        pic._add_bonds_from_species(mol, core_graph)
-        pic._add_angles_from_species(mol, core_graph)
-        pic._add_dihedrals_from_species(mol, core_graph)
+        # take a copy of mol as mol.graph might be changed
+        mol = mol.copy()
+        _connect_graph_for_species(mol)
+        pic._add_bonds_from_species(mol)
+        pic._add_angles_from_species(mol)
+        pic._add_dihedrals_from_species(mol)
         return pic
 
-    @staticmethod
-    def _get_connected_graph_from_species(mol: "Species") -> "MolecularGraph":
-        """
-        Creates a fully connected graph from the graph of a species, by
-        (1) joining disconnected fragments by their shortest distance,
-        (2) connecting constrained bonds, (3) joining hydrogen bonds,
-        if present.
-
-        Args:
-            mol: A species (must have atoms and graph)
-
-        Returns:
-            (MolecularGraph):
-        """
-        assert mol.graph is not None, "Species must have graph!"
-        core_graph = copy.deepcopy(mol.graph)
-
-        # join hydrogen bonds
-        h_bond_x = ["N", "O", "F", "P", "S", "Cl"]
-        for i, j in itertools.combinations(range(mol.n_atoms), r=2):
-            if (
-                mol.atoms[i].label in h_bond_x
-                and mol.atoms[j].label == "H"
-                or mol.atoms[j].label in h_bond_x
-                and mol.atoms[i].label == "H"
-            ):
-                vdw_sum = mol.atoms[i].vdw_radius + mol.atoms[j].vdw_radius
-                if mol.distance(i, j) < 0.9 * vdw_sum:
-                    if not core_graph.has_edge(i, j):
-                        core_graph.add_edge(i, j, pi=False, active=False)
-
-        # join disconnected graph components
-        if not core_graph.is_connected:
-            components = core_graph.connected_components()
-            for comp_i, comp_j in itertools.combinations(components, r=2):
-                min_dist = float("inf")
-                min_pair = (-1, -1)
-                for i, j in itertools.product(list(comp_i), list(comp_j)):
-                    if mol.distance(i, j) < min_dist:
-                        min_dist = mol.distance(i, j)
-                        min_pair = (i, j)
-                # avoid connecting distant components
-                if min_dist < Distance(4.0, "ang"):
-                    core_graph.add_edge(*min_pair, pi=False, active=False)
-
-        if not core_graph.is_connected:
-            raise RuntimeError(
-                "Unable to join all the fragments, distance between "
-                "one or more pairs of fragments is too high (>4.0 Å)"
-            )
-
-        # The constraints should be counted as bonds
-        if mol.constraints.distance is not None:
-            for i, j in mol.constraints.distance:
-                if not core_graph.has_edge(i, j):
-                    core_graph.add_edge(i, j, pi=False, active=False)
-
-        return core_graph
-
     def _add_bonds_from_species(
         self,
         mol: "Species",
-        core_graph: "MolecularGraph",
     ):
         """
         Add bonds to the current set of primitives, from the
-        connectivity graph supplied
+        connectivity graph of the species
 
         Args:
             mol: The species object
-            core_graph: The connectivity graph
         """
+        assert mol.graph is not None
+
         n = 0
-        for i, j in sorted(core_graph.edges):
+        for i, j in sorted(mol.graph.edges):
             if (
                 mol.constraints.distance is not None
                 and (i, j) in mol.constraints.distance
@@ -366,7 +309,6 @@ def _add_bonds_from_species(
     @staticmethod
     def _get_ref_for_linear_angle(
         mol,
-        core_graph,
         lin_thresh,
         a,
         b,
@@ -381,7 +323,6 @@ def _get_ref_for_linear_angle(
 
         Args:
             mol:
-            core_graph:
             lin_thresh:
             a:
             b:
@@ -395,7 +336,7 @@ def _get_ref_for_linear_angle(
         """
         # only check bonded atoms if requested
         if bonded:
-            near_atoms = list(core_graph.neighbors(b))
+            near_atoms = list(mol.graph.neighbors(b))
             near_atoms.remove(a)
             near_atoms.remove(c)
 
@@ -429,7 +370,6 @@ def _get_ref_for_linear_angle(
     def _add_angles_from_species(
         self,
         mol: "Species",
-        core_graph: "MolecularGraph",
         lin_thresh: Angle = Angle(170, "deg"),
     ) -> None:
         """
@@ -438,13 +378,13 @@ def _add_angles_from_species(
 
         Args:
             mol (Species): The species object
-            core_graph (MolecularGraph): The connectivity graph
             lin_thresh (Angle): The angle threshold for linearity
         """
+        assert mol.graph is not None
         lin_thresh = lin_thresh.to("rad")
 
         for o in range(mol.n_atoms):
-            for n, m in itertools.combinations(core_graph.neighbors(o), r=2):
+            for n, m in itertools.combinations(mol.graph.neighbors(o), r=2):
                 if mol.angle(m, o, n) < lin_thresh:
                     self.add(PrimitiveBondAngle(m=m, o=o, n=n))
                 else:
@@ -452,7 +392,7 @@ def _add_angles_from_species(
                     # linear angle is skipped and instead an out-of-plane
                     # (improper dihedral) coordinate is used
                     r = self._get_ref_for_linear_angle(
-                        mol, core_graph, lin_thresh, m, o, n, bonded=True
+                        mol, lin_thresh, m, o, n, bonded=True
                     )
                     if r is not None:
                         self.add(PrimitiveDihedralAngle(m, r, o, n))
@@ -461,7 +401,7 @@ def _add_angles_from_species(
                     # Otherwise, we use a nearby (< 4.0 A) reference atom to
                     # define two orthogonal linear bends
                     r = self._get_ref_for_linear_angle(
-                        mol, core_graph, lin_thresh, m, o, n, bonded=False
+                        mol, lin_thresh, m, o, n, bonded=False
                     )
                     if r is not None:
                         self.add(
@@ -494,7 +434,6 @@ def _add_angles_from_species(
     def _add_dihedrals_from_species(
         self,
         mol: "Species",
-        core_graph: "MolecularGraph",
         lin_thresh: Angle = Angle(170, "deg"),
     ) -> None:
         """
@@ -503,13 +442,13 @@ def _add_dihedrals_from_species(
 
         Args:
             mol (Species): The species
-            core_graph (MolecularGraph): The connectivity graph
             lin_thresh (Angle): The threshold for linearity
         """
         # no dihedrals possible with less than 4 atoms
         if mol.n_atoms < 4:
             return
 
+        assert mol.graph is not None
         lin_thresh = lin_thresh.to("rad")
         zero_angle_thresh = np.pi - lin_thresh
 
@@ -526,12 +465,12 @@ def is_dihedral_well_defined(w, x, y, z):
             return not (is_linear_1 or is_linear_2)
 
         # add normal dihedrals
-        for o, p in list(core_graph.edges):
-            for m in core_graph.neighbors(o):
+        for o, p in list(mol.graph.edges):
+            for m in mol.graph.neighbors(o):
                 if m == p:
                     continue
 
-                for n in core_graph.neighbors(p):
+                for n in mol.graph.neighbors(p):
                     if n == o:
                         continue
 
@@ -561,7 +500,7 @@ def extend_chain(chain: List[int]):
         for b in range(mol.n_atoms):
             if any(b in chain for chain in linear_chains):
                 continue
-            for a, c in itertools.combinations(core_graph.neighbors(b), r=2):
+            for a, c in itertools.combinations(mol.graph.neighbors(b), r=2):
                 if any(a in chain for chain in linear_chains) or any(
                     c in chain for chain in linear_chains
                 ):
@@ -574,14 +513,14 @@ def extend_chain(chain: List[int]):
         # add linear chain dihedrals
         for chain in linear_chains:
             o, p = chain[0], chain[-1]
-            for m in core_graph.neighbors(o):
+            for m in mol.graph.neighbors(o):
                 if m == p:
                     continue
 
                 if m in chain:
                     continue
 
-                for n in core_graph.neighbors(p):
+                for n in mol.graph.neighbors(p):
                     if n == o:
                         continue
 
@@ -594,3 +533,58 @@ def extend_chain(chain: List[int]):
                     if is_dihedral_well_defined(m, o, p, n):
                         self.add(PrimitiveDihedralAngle(m, o, p, n))
         return None
+
+
+def _connect_graph_for_species(mol: "Species") -> None:
+    """
+    Creates a fully connected graph from the graph of a species, by
+    (1) joining disconnected fragments by their shortest distance,
+    (2) connecting constrained bonds, (3) joining hydrogen bonds,
+    if present. The molecular graph is modified in-place.
+
+    Args:
+        mol: A species (must have atoms and graph)
+    """
+    assert mol.graph is not None, "Species must have graph!"
+
+    # join hydrogen bonds
+    h_bond_x = ["N", "O", "F", "P", "S", "Cl"]
+    for i, j in itertools.combinations(range(mol.n_atoms), r=2):
+        if (
+            mol.atoms[i].label in h_bond_x
+            and mol.atoms[j].label == "H"
+            or mol.atoms[j].label in h_bond_x
+            and mol.atoms[i].label == "H"
+        ):
+            vdw_sum = mol.atoms[i].vdw_radius + mol.atoms[j].vdw_radius
+            if mol.distance(i, j) < 0.9 * vdw_sum:
+                if not mol.graph.has_edge(i, j):
+                    mol.graph.add_edge(i, j, pi=False, active=False)
+
+    # join disconnected graph components
+    if not mol.graph.is_connected:
+        components = mol.graph.connected_components()
+        for comp_i, comp_j in itertools.combinations(components, r=2):
+            min_dist = float("inf")
+            min_pair = (-1, -1)
+            for i, j in itertools.product(list(comp_i), list(comp_j)):
+                if mol.distance(i, j) < min_dist:
+                    min_dist = mol.distance(i, j)
+                    min_pair = (i, j)
+            # avoid connecting distant components
+            if min_dist < Distance(4.0, "ang"):
+                mol.graph.add_edge(*min_pair, pi=False, active=False)
+
+    if not mol.graph.is_connected:
+        raise RuntimeError(
+            "Unable to join all the fragments, distance between "
+            "one or more pairs of fragments is too high (>4.0 Å)"
+        )
+
+    # The constraints should be counted as bonds
+    if mol.constraints.distance is not None:
+        for i, j in mol.constraints.distance:
+            if not mol.graph.has_edge(i, j):
+                mol.graph.add_edge(i, j, pi=False, active=False)
+
+    return None

From 32ddb137579c07df59744de6738a12d8448b5180 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 30 Dec 2023 20:26:59 +0000
Subject: [PATCH 41/47] pr update 4

---
 autode/opt/coordinates/internals.py | 35 +++++++++++++----------------
 1 file changed, 16 insertions(+), 19 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 13012a139..68a1dc4ea 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -36,6 +36,10 @@
     )
 
 
+# Angle threshold for linearity (should be in radians)
+_lin_thresh = Angle(170, "deg").to("rad")
+
+
 class InternalCoordinates(OptCoordinates, ABC):  # lgtm [py/missing-equals]
     def __new__(cls, input_array) -> "InternalCoordinates":
         """New instance of these internal coordinates"""
@@ -309,7 +313,6 @@ def _add_bonds_from_species(
     @staticmethod
     def _get_ref_for_linear_angle(
         mol,
-        lin_thresh,
         a,
         b,
         c,
@@ -323,7 +326,6 @@ def _get_ref_for_linear_angle(
 
         Args:
             mol:
-            lin_thresh:
             a:
             b:
             c:
@@ -334,6 +336,7 @@ def _get_ref_for_linear_angle(
         Returns:
             (int|None): The index of the ref. atom if found, else None
         """
+
         # only check bonded atoms if requested
         if bonded:
             near_atoms = list(mol.graph.neighbors(b))
@@ -352,10 +355,10 @@ def _get_ref_for_linear_angle(
         deviations_from_90 = {}
         for atom in near_atoms:
             i_b_a = mol.angle(atom, b, a)
-            if i_b_a > lin_thresh or i_b_a < (np.pi - lin_thresh):
+            if i_b_a > _lin_thresh or i_b_a < (np.pi - _lin_thresh):
                 continue
             i_b_c = mol.angle(atom, b, c)
-            if i_b_c > lin_thresh or i_b_c < (np.pi - lin_thresh):
+            if i_b_c > _lin_thresh or i_b_c < (np.pi - _lin_thresh):
                 continue
             deviation_a = abs(i_b_a - np.pi / 2)
             deviation_c = abs(i_b_c - np.pi / 2)
@@ -370,7 +373,6 @@ def _get_ref_for_linear_angle(
     def _add_angles_from_species(
         self,
         mol: "Species",
-        lin_thresh: Angle = Angle(170, "deg"),
     ) -> None:
         """
         Modify the set of primitives in-place by adding angles, from the
@@ -378,21 +380,19 @@ def _add_angles_from_species(
 
         Args:
             mol (Species): The species object
-            lin_thresh (Angle): The angle threshold for linearity
         """
         assert mol.graph is not None
-        lin_thresh = lin_thresh.to("rad")
 
         for o in range(mol.n_atoms):
             for n, m in itertools.combinations(mol.graph.neighbors(o), r=2):
-                if mol.angle(m, o, n) < lin_thresh:
+                if mol.angle(m, o, n) < _lin_thresh:
                     self.add(PrimitiveBondAngle(m=m, o=o, n=n))
                 else:
                     # If central atom is connected to another atom, then the
                     # linear angle is skipped and instead an out-of-plane
                     # (improper dihedral) coordinate is used
                     r = self._get_ref_for_linear_angle(
-                        mol, lin_thresh, m, o, n, bonded=True
+                        mol, m, o, n, bonded=True
                     )
                     if r is not None:
                         self.add(PrimitiveDihedralAngle(m, r, o, n))
@@ -401,7 +401,7 @@ def _add_angles_from_species(
                     # Otherwise, we use a nearby (< 4.0 A) reference atom to
                     # define two orthogonal linear bends
                     r = self._get_ref_for_linear_angle(
-                        mol, lin_thresh, m, o, n, bonded=False
+                        mol, m, o, n, bonded=False
                     )
                     if r is not None:
                         self.add(
@@ -434,7 +434,6 @@ def _add_angles_from_species(
     def _add_dihedrals_from_species(
         self,
         mol: "Species",
-        lin_thresh: Angle = Angle(170, "deg"),
     ) -> None:
         """
         Modify the set of primitives in-place by adding dihedrals (torsions),
@@ -442,24 +441,22 @@ def _add_dihedrals_from_species(
 
         Args:
             mol (Species): The species
-            lin_thresh (Angle): The threshold for linearity
         """
         # no dihedrals possible with less than 4 atoms
         if mol.n_atoms < 4:
             return
 
         assert mol.graph is not None
-        lin_thresh = lin_thresh.to("rad")
-        zero_angle_thresh = np.pi - lin_thresh
+        zero_angle_thresh = np.pi - _lin_thresh
 
         def is_dihedral_well_defined(w, x, y, z):
             """A dihedral is well-defined if any angle is not linear"""
             is_linear_1 = (
-                mol.angle(w, x, y) > lin_thresh
+                mol.angle(w, x, y) > _lin_thresh
                 or mol.angle(w, x, y) < zero_angle_thresh
             )
             is_linear_2 = (
-                mol.angle(x, y, z) > lin_thresh
+                mol.angle(x, y, z) > _lin_thresh
                 or mol.angle(x, y, z) < zero_angle_thresh
             )
             return not (is_linear_1 or is_linear_2)
@@ -488,11 +485,11 @@ def extend_chain(chain: List[int]):
                 if idx in chain:
                     continue
                 # if idx -- 0 -- 1 > 170 degrees
-                if mol.angle(chain[1], chain[0], idx) > lin_thresh:
+                if mol.angle(chain[1], chain[0], idx) > _lin_thresh:
                     chain.insert(0, idx)
                     continue
                 # if (-2) -- (-1) -- idx > 170 degrees
-                if mol.angle(chain[-2], chain[-1], idx) > lin_thresh:
+                if mol.angle(chain[-2], chain[-1], idx) > _lin_thresh:
                     chain.append(idx)
                     continue
 
@@ -505,7 +502,7 @@ def extend_chain(chain: List[int]):
                     c in chain for chain in linear_chains
                 ):
                     continue
-                if mol.angle(a, b, c) > lin_thresh:
+                if mol.angle(a, b, c) > _lin_thresh:
                     chain = [a, b, c]
                     extend_chain(chain)
                     linear_chains.append(chain)

From d6fabb435ec53d44e3f077199893b97fc841d306 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sun, 31 Dec 2023 21:41:23 +0000
Subject: [PATCH 42/47] pr update 5

---
 autode/opt/coordinates/internals.py | 81 ++++++++++++++++++-----------
 1 file changed, 50 insertions(+), 31 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 68a1dc4ea..8e6fb8a2f 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -331,7 +331,8 @@ def _get_ref_for_linear_angle(
             c:
             bonded: Whether to look for only atoms bonded to the central
                     atom (b) for reference
-            dist_thresh: The distance threshold to connect
+            dist_thresh: The distance threshold to check for atoms if
+                    bonded = False
 
         Returns:
             (int|None): The index of the ref. atom if found, else None
@@ -431,6 +432,51 @@ def _add_angles_from_species(
 
         return None
 
+    @staticmethod
+    def _get_linear_chains(mol: "Species") -> List[List[int]]:
+        """
+        Obtain a list of all the continuous chains of linear atoms
+        present in the species i.e. A--B--C--D--E...
+
+        Args:
+            mol: The species
+
+        Returns:
+            (list[list[int]): A list of lists, each containing indices
+                            of the atoms of linear chains in order
+        """
+        assert mol.graph is not None
+
+        def extend_chain(chain: List[int]):
+            """Extend a chain in-place"""
+            for idx in range(mol.n_atoms):
+                if idx in chain:
+                    continue
+
+                if mol.angle(chain[1], chain[0], idx) > _lin_thresh:
+                    chain.insert(0, idx)
+                    continue
+
+                if mol.angle(chain[-2], chain[-1], idx) > _lin_thresh:
+                    chain.append(idx)
+                    continue
+            return None
+
+        linear_chains: List[list] = []
+        for b in range(mol.n_atoms):
+            for a, c in itertools.combinations(mol.graph.neighbors(b), r=2):
+                if any(
+                    a in chain and b in chain and c in chain
+                    for chain in linear_chains
+                ):
+                    continue
+                if mol.angle(a, b, c) > _lin_thresh:
+                    chain = [a, b, c]
+                    extend_chain(chain)
+                    linear_chains.append(chain)
+
+        return linear_chains
+
     def _add_dihedrals_from_species(
         self,
         mol: "Species",
@@ -450,7 +496,7 @@ def _add_dihedrals_from_species(
         zero_angle_thresh = np.pi - _lin_thresh
 
         def is_dihedral_well_defined(w, x, y, z):
-            """A dihedral is well-defined if any angle is not linear"""
+            """A dihedral is not well-defined if any angle is linear"""
             is_linear_1 = (
                 mol.angle(w, x, y) > _lin_thresh
                 or mol.angle(w, x, y) < zero_angle_thresh
@@ -478,36 +524,9 @@ def is_dihedral_well_defined(w, x, y, z):
                     if is_dihedral_well_defined(m, o, p, n):
                         self.add(PrimitiveDihedralAngle(m, o, p, n))
 
-        # find all linear atom chains A--B--C--D... and add dihedrals to terminal atoms
-
-        def extend_chain(chain: List[int]):
-            for idx in range(mol.n_atoms):
-                if idx in chain:
-                    continue
-                # if idx -- 0 -- 1 > 170 degrees
-                if mol.angle(chain[1], chain[0], idx) > _lin_thresh:
-                    chain.insert(0, idx)
-                    continue
-                # if (-2) -- (-1) -- idx > 170 degrees
-                if mol.angle(chain[-2], chain[-1], idx) > _lin_thresh:
-                    chain.append(idx)
-                    continue
-
-        linear_chains: List[list] = []
-        for b in range(mol.n_atoms):
-            if any(b in chain for chain in linear_chains):
-                continue
-            for a, c in itertools.combinations(mol.graph.neighbors(b), r=2):
-                if any(a in chain for chain in linear_chains) or any(
-                    c in chain for chain in linear_chains
-                ):
-                    continue
-                if mol.angle(a, b, c) > _lin_thresh:
-                    chain = [a, b, c]
-                    extend_chain(chain)
-                    linear_chains.append(chain)
+        # add linear chain terminal dihedrals
+        linear_chains = self._get_linear_chains(mol)
 
-        # add linear chain dihedrals
         for chain in linear_chains:
             o, p = chain[0], chain[-1]
             for m in mol.graph.neighbors(o):

From 239c08a4a603c16bd0782d4210dd6b0c1248862d Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sun, 31 Dec 2023 21:43:47 +0000
Subject: [PATCH 43/47] pr update 6

---
 autode/opt/coordinates/internals.py | 10 ++--------
 1 file changed, 2 insertions(+), 8 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 8e6fb8a2f..a996f6c11 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -587,15 +587,9 @@ def _connect_graph_for_species(mol: "Species") -> None:
                 if mol.distance(i, j) < min_dist:
                     min_dist = mol.distance(i, j)
                     min_pair = (i, j)
-            # avoid connecting distant components
-            if min_dist < Distance(4.0, "ang"):
-                mol.graph.add_edge(*min_pair, pi=False, active=False)
+            mol.graph.add_edge(*min_pair, pi=False, active=False)
 
-    if not mol.graph.is_connected:
-        raise RuntimeError(
-            "Unable to join all the fragments, distance between "
-            "one or more pairs of fragments is too high (>4.0 Å)"
-        )
+    assert mol.graph.is_connected, "Unknown error in connecting graph"
 
     # The constraints should be counted as bonds
     if mol.constraints.distance is not None:

From 16282cbcabc8b44225b0bcd341831ab40dde2c31 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sun, 31 Dec 2023 22:13:15 +0000
Subject: [PATCH 44/47] pr update 7

---
 autode/opt/coordinates/internals.py | 122 +++++++++++++++++-----------
 tests/test_opt/test_coordiantes.py  |   6 +-
 2 files changed, 78 insertions(+), 50 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index a996f6c11..757505447 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -280,6 +280,7 @@ def from_species(cls, mol: "Species") -> "AnyPIC":
         pic._add_bonds_from_species(mol)
         pic._add_angles_from_species(mol)
         pic._add_dihedrals_from_species(mol)
+        pic._add_chain_dihedrals_from_species(mol)
         return pic
 
     def _add_bonds_from_species(
@@ -361,10 +362,8 @@ def _get_ref_for_linear_angle(
             i_b_c = mol.angle(atom, b, c)
             if i_b_c > _lin_thresh or i_b_c < (np.pi - _lin_thresh):
                 continue
-            deviation_a = abs(i_b_a - np.pi / 2)
-            deviation_c = abs(i_b_c - np.pi / 2)
-            avg_dev = (deviation_a + deviation_c) / 2
-            deviations_from_90[atom] = avg_dev
+            deviation = abs(i_b_a - np.pi / 2)
+            deviations_from_90[atom] = deviation
 
         if len(deviations_from_90) == 0:
             return None
@@ -432,6 +431,41 @@ def _add_angles_from_species(
 
         return None
 
+    def _add_dihedrals_from_species(
+        self,
+        mol: "Species",
+    ) -> None:
+        """
+        Modify the set of primitives in-place by adding dihedrals (torsions),
+        from the connectivity graph supplied
+
+        Args:
+            mol (Species): The species
+        """
+        # no dihedrals possible with less than 4 atoms
+        if mol.n_atoms < 4:
+            return
+
+        assert mol.graph is not None
+
+        for o, p in list(mol.graph.edges):
+            for m in mol.graph.neighbors(o):
+                if m == p:
+                    continue
+
+                for n in mol.graph.neighbors(p):
+                    if n == o:
+                        continue
+
+                    # avoid triangle rings like cyclopropane
+                    if n == m:
+                        continue
+
+                    if _is_dihedral_well_defined(mol, m, o, p, n):
+                        self.add(PrimitiveDihedralAngle(m, o, p, n))
+
+        return None
+
     @staticmethod
     def _get_linear_chains(mol: "Species") -> List[List[int]]:
         """
@@ -477,54 +511,18 @@ def extend_chain(chain: List[int]):
 
         return linear_chains
 
-    def _add_dihedrals_from_species(
-        self,
-        mol: "Species",
-    ) -> None:
+    def _add_chain_dihedrals_from_species(self, mol: "Species"):
         """
-        Modify the set of primitives in-place by adding dihedrals (torsions),
-        from the connectivity graph supplied
+        Add extra dihedrals for chain molecules like allene, which
+        are required to cover all degrees of freedom
 
         Args:
-            mol (Species): The species
-        """
-        # no dihedrals possible with less than 4 atoms
-        if mol.n_atoms < 4:
-            return
-
-        assert mol.graph is not None
-        zero_angle_thresh = np.pi - _lin_thresh
-
-        def is_dihedral_well_defined(w, x, y, z):
-            """A dihedral is not well-defined if any angle is linear"""
-            is_linear_1 = (
-                mol.angle(w, x, y) > _lin_thresh
-                or mol.angle(w, x, y) < zero_angle_thresh
-            )
-            is_linear_2 = (
-                mol.angle(x, y, z) > _lin_thresh
-                or mol.angle(x, y, z) < zero_angle_thresh
-            )
-            return not (is_linear_1 or is_linear_2)
-
-        # add normal dihedrals
-        for o, p in list(mol.graph.edges):
-            for m in mol.graph.neighbors(o):
-                if m == p:
-                    continue
-
-                for n in mol.graph.neighbors(p):
-                    if n == o:
-                        continue
-
-                    # avoid triangle rings like cyclopropane
-                    if n == m:
-                        continue
+            mol:
 
-                    if is_dihedral_well_defined(m, o, p, n):
-                        self.add(PrimitiveDihedralAngle(m, o, p, n))
+        Returns:
 
-        # add linear chain terminal dihedrals
+        """
+        assert mol.graph is not None
         linear_chains = self._get_linear_chains(mol)
 
         for chain in linear_chains:
@@ -546,8 +544,9 @@ def is_dihedral_well_defined(w, x, y, z):
                     if n in chain:
                         continue
 
-                    if is_dihedral_well_defined(m, o, p, n):
+                    if _is_dihedral_well_defined(mol, m, o, p, n):
                         self.add(PrimitiveDihedralAngle(m, o, p, n))
+
         return None
 
 
@@ -598,3 +597,30 @@ def _connect_graph_for_species(mol: "Species") -> None:
                 mol.graph.add_edge(i, j, pi=False, active=False)
 
     return None
+
+
+def _is_dihedral_well_defined(mol, a, b, c, d) -> bool:
+    """
+    A dihedral a--b--c--d is only well-defined when the
+    two constituent angles are not linear
+
+    Args:
+        mol:
+        a:
+        b:
+        c:
+        d:
+
+    Returns:
+        (bool): True if well-defined otherwise False
+    """
+    zero_angle_thresh = np.pi - _lin_thresh
+    is_linear_1 = (
+        mol.angle(a, b, c) > _lin_thresh
+        or mol.angle(a, b, c) < zero_angle_thresh
+    )
+    is_linear_2 = (
+        mol.angle(b, c, d) > _lin_thresh
+        or mol.angle(b, c, d) < zero_angle_thresh
+    )
+    return not (is_linear_1 or is_linear_2)
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 8e4b88741..8c315f107 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -838,8 +838,10 @@ def test_pic_generation_linear_angle_ref():
     assert not any(ic1 == ic2 for ic1, ic2 in itertools.combinations(pic, r=2))
     # check that linear bends use reference atoms, not dummy
     assert not any(isinstance(ic, PrimitiveDummyLinearAngle) for ic in pic)
-    # for C-Fe-C, one out-of-plane dihedral should be present
-    assert PrimitiveDihedralAngle(3, 5, 2, 1) in pic
+    assert PrimitiveLinearAngle(4, 3, 2, 8, LinearBendType.BEND) in pic
+    # for C-Fe-C, only one out-of-plane dihedral should be present
+    assert PrimitiveDihedralAngle(3, 7, 2, 1) in pic
+    assert sum(isinstance(ic, PrimitiveDihedralAngle) for ic in pic) == 1
     # check degrees of freedom = 3N - 6
     _ = pic(m.coordinates.flatten())
     assert np.linalg.matrix_rank(pic.B) == 3 * m.n_atoms - 6

From 3943421811c4b9bb3a514a054f3590fdb1af7962 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 1 Jan 2024 10:43:52 +0000
Subject: [PATCH 45/47] pr update 8

---
 autode/opt/coordinates/internals.py | 6 ++++--
 tests/test_opt/test_coordiantes.py  | 2 +-
 tests/test_opt/test_crfo.py         | 2 +-
 3 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 757505447..3de919b13 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -362,8 +362,10 @@ def _get_ref_for_linear_angle(
             i_b_c = mol.angle(atom, b, c)
             if i_b_c > _lin_thresh or i_b_c < (np.pi - _lin_thresh):
                 continue
-            deviation = abs(i_b_a - np.pi / 2)
-            deviations_from_90[atom] = deviation
+            deviation_a = abs(i_b_a - np.pi / 2)
+            deviation_b = abs(i_b_c - np.pi / 2)
+            avg_dev = (deviation_a + deviation_b) / 2
+            deviations_from_90[atom] = avg_dev
 
         if len(deviations_from_90) == 0:
             return None
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 8c315f107..c58600231 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -840,7 +840,7 @@ def test_pic_generation_linear_angle_ref():
     assert not any(isinstance(ic, PrimitiveDummyLinearAngle) for ic in pic)
     assert PrimitiveLinearAngle(4, 3, 2, 8, LinearBendType.BEND) in pic
     # for C-Fe-C, only one out-of-plane dihedral should be present
-    assert PrimitiveDihedralAngle(3, 7, 2, 1) in pic
+    assert PrimitiveDihedralAngle(3, 5, 2, 1) in pic
     assert sum(isinstance(ic, PrimitiveDihedralAngle) for ic in pic) == 1
     # check degrees of freedom = 3N - 6
     _ = pic(m.coordinates.flatten())
diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index e1060592c..f040278a8 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -363,7 +363,7 @@ def test_optimise_linear_bend_with_ref():
     # the Fe-C-O angle are manually deviated
     assert val.Angle(170, "deg") < mol.angle(2, 3, 4) < val.Angle(176, "deg")
     # large molecule so allow few iters, no need to converge fully
-    opt = CRFOptimiser(maxiter=15, gtol=1e-4, etol=1e-5)
+    opt = CRFOptimiser(maxiter=15, gtol=1e-5, etol=1e-6)
     opt.run(mol, XTB())
     assert mol.angle(2, 3, 4) > val.Angle(178, "deg")
 

From 2547c03cd0aef03a97eec37b2ac3fdf8b09ddc09 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 1 Jan 2024 11:07:35 +0000
Subject: [PATCH 46/47] pr update

---
 tests/test_opt/test_coordiantes.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index c58600231..dae461932 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -944,4 +944,4 @@ def test_pic_generation_square_planar():
     # all degrees of freedom
     pic = AnyPIC.from_species(ptcl4)
     _ = pic(ptcl4.coordinates.flatten())
-    assert np.linalg.matrix_rank(pic.B) == 3 * 5 - 6
+    assert np.linalg.matrix_rank(pic.B) == 3 * ptcl4.n_atoms - 6

From f8187e8a4c9ad2da08b8d11cd2b1d7d27273e0ba Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Wed, 3 Jan 2024 14:55:16 +0000
Subject: [PATCH 47/47] pr update, changelog update

---
 doc/changelog.rst                  | 2 ++
 tests/test_opt/test_coordiantes.py | 4 ++++
 2 files changed, 6 insertions(+)

diff --git a/doc/changelog.rst b/doc/changelog.rst
index 1d2b97732..3113ebb9d 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -8,6 +8,8 @@ Changelog
 Functionality improvements
 **************************
 - Replaces hard-coded derivatives for primitive internal coordinates with automatic differentiation
+- More comprehensive primitive internal coordinate generation
+- Adds the capacity to handle linear molecules in internal coordinate system
 
 Bug Fixes
 *********
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index dae461932..eb8e21040 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -551,6 +551,10 @@ def test_pic_add_sanity_checking():
     with pytest.raises(AssertionError):
         c.add(3)
 
+    # pic append is disallowed
+    with pytest.raises(NotImplementedError, match="Please use PIC.add()"):
+        c.append(PrimitiveDistance(0, 1))
+
     # pic should not allow duplicate coordinates to be added
     assert len(c) == 1
     c.add(PrimitiveDistance(1, 0))