Repository: duartegroup/autodE
Pull request #314 (2023-12-21 18:47:17)

More comprehensive primitive generation
Partial fix for #312, #261 and #313.


***

## Checklist

* [ ] The changes include an associated explanation of how/why
* [ ] Test pass
* [ ] Documentation has been updated
* [ ] Changelog has been updated


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) Report
Attention: `39 lines` in your changes are missing coverage. Please review.
> Comparison is base [(`1075b84`)](https://app.codecov.io/gh/duartegroup/autodE/commit/1075b84d22e1bcc0ceba69152305ce155aed4494?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) 97.37% compared to head [(`0757cda`)](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) 97.31%.

| [Files](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) | Patch % | Lines |
|---|---|---|
| [autode/opt/coordinates/internals.py](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-YXV0b2RlL29wdC9jb29yZGluYXRlcy9pbnRlcm5hbHMucHk=) | 77.85% | [31 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [autode/opt/coordinates/\_autodiff.py](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-YXV0b2RlL29wdC9jb29yZGluYXRlcy9fYXV0b2RpZmYucHk=) | 98.29% | [5 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [autode/mol\_graphs.py](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-YXV0b2RlL21vbF9ncmFwaHMucHk=) | 80.00% | [1 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [autode/opt/coordinates/primitives.py](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-YXV0b2RlL29wdC9jb29yZGluYXRlcy9wcmltaXRpdmVzLnB5) | 99.23% | [1 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |
| [tests/test\_opt/molecules.py](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#diff-dGVzdHMvdGVzdF9vcHQvbW9sZWN1bGVzLnB5) | 75.00% | [1 Missing :warning: ](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##           master     #314      +/-   ##
==========================================
- Coverage   97.37%   97.31%   -0.07%     
==========================================
  Files         206      208       +2     
  Lines       23048    23616     +568     
==========================================
+ Hits        22443    22981     +538     
- Misses        605      635      +30     
```

| [Flag](https://app.codecov.io/gh/duartegroup/autodE/pull/314/flags?src=pr&el=flags&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) | Coverage Δ | |
|---|---|---|
| [unittests](https://app.codecov.io/gh/duartegroup/autodE/pull/314/flags?src=pr&el=flag&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup) | `97.31% <95.08%> (-0.07%)` | :arrow_down: |

Flags with carried forward coverage won't be shown. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup#carryforward-flags-in-the-pull-request-comment) to find out more.


</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/duartegroup/autodE/pull/314?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=duartegroup).


Patch:
diff --git a/autode/atoms.py b/autode/atoms.py
index 699e7f6..6e88d8f 100644
--- a/autode/atoms.py
+++ b/autode/atoms.py
@@ -1041,9 +1041,10 @@ class AtomCollection:
                 f"least one zero vector"
             )
 
-        value = np.arccos(np.dot(vec1, vec2) / norms)
+        # 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 9ef1486..521036a 100644
--- a/autode/mol_graphs.py
+++ b/autode/mol_graphs.py
@@ -115,6 +115,15 @@ class MolecularGraph(nx.Graph):
 
         return matcher
 
+    @property
+    def is_connected(self) -> bool:
+        """Is this graph fully connected (i.e. not separate parts)"""
+        return nx.is_connected(self)
+
+    def connected_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 2aadb5d..36daee4 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/_autodiff.py b/autode/opt/coordinates/_autodiff.py
index 0fa48d4..0855a44 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 @@ class DifferentiableVector3D:
         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
@@ -600,7 +604,9 @@ class DifferentiableVector3D:
             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
 
@@ -611,13 +617,12 @@ class DifferentiableVector3D:
             (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)
         return dot
 
-    def norm(self) -> "VectorHyperDual":
+    def norm(self) -> Union["VectorHyperDual", numeric_type]:
         """
         Euclidean (l2) norm of this 3D vector
 
@@ -627,7 +632,6 @@ class DifferentiableVector3D:
         norm = DifferentiableMath.sqrt(
             self._data[0] ** 2 + self._data[1] ** 2 + self._data[2] ** 2
         )
-        assert isinstance(norm, VectorHyperDual)
         return norm
 
     def __add__(
@@ -690,6 +694,18 @@ class DifferentiableVector3D:
         """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":
@@ -702,6 +718,7 @@ class DifferentiableVector3D:
         Returns:
             (DifferentiableVector3D):
         """
+        self._check_same_type(other)
         return DifferentiableVector3D(
             [
                 self._data[1] * other._data[2]
diff --git a/autode/opt/coordinates/base.py b/autode/opt/coordinates/base.py
index 0f5135c..747680e 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 @@ if TYPE_CHECKING:
     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 9e7462e..3de919b 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -6,20 +6,29 @@ x : Cartesian coordinates
 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 numpy as np
-
+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, Distance
+from autode.opt.coordinates.base import OptCoordinates
 from autode.opt.coordinates.primitives import (
     PrimitiveInverseDistance,
     Primitive,
     PrimitiveDistance,
+    ConstrainedPrimitiveDistance,
+    PrimitiveBondAngle,
+    PrimitiveDummyLinearAngle,
+    PrimitiveLinearAngle,
     PrimitiveDihedralAngle,
+    LinearBendType,
 )
 
 if TYPE_CHECKING:
+    from autode.species import Species
     from autode.opt.coordinates.cartesian import CartesianCoordinates
     from autode.opt.coordinates.primitives import (
         ConstrainedPrimitive,
@@ -27,12 +36,19 @@ if TYPE_CHECKING:
     )
 
 
+# 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"""
 
         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))
 
@@ -82,10 +98,18 @@ class PIC(list, ABC):
                 f"from {args}. Must be primitive internals"
             )
 
+    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"""
-        assert isinstance(item, Primitive), "Must be a Primitive type!"
-        super().append(item)
+        raise NotImplementedError(
+            "Please use PIC.add() to add new primitives to the set"
+        )
 
     @property
     def B(self) -> np.ndarray:
@@ -210,7 +234,7 @@ class _FunctionOfDistances(PIC):
         # 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
 
@@ -234,3 +258,371 @@ class PrimitiveDistances(_FunctionOfDistances):
 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.
+
+        Args:
+            mol: The species object
+
+        Returns:
+            (AnyPIC): The set of primitive internals
+        """
+        pic = cls()
+        # 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)
+        pic._add_chain_dihedrals_from_species(mol)
+        return pic
+
+    def _add_bonds_from_species(
+        self,
+        mol: "Species",
+    ):
+        """
+        Add bonds to the current set of primitives, from the
+        connectivity graph of the species
+
+        Args:
+            mol: The species object
+        """
+        assert mol.graph is not None
+
+        n = 0
+        for i, j in sorted(mol.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,
+        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:
+            a:
+            b:
+            c:
+            bonded: Whether to look for only atoms bonded to the central
+                    atom (b) for reference
+            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
+        """
+
+        # only check bonded atoms if requested
+        if bonded:
+            near_atoms = list(mol.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_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
+
+        return min(deviations_from_90, key=deviations_from_90.get)  # type: ignore
+
+    def _add_angles_from_species(
+        self,
+        mol: "Species",
+    ) -> None:
+        """
+        Modify the set of primitives in-place by adding angles, from the
+        connectivity graph supplied
+
+        Args:
+            mol (Species): The species object
+        """
+        assert mol.graph is not None
+
+        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:
+                    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, 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, 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 _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]]:
+        """
+        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_chain_dihedrals_from_species(self, mol: "Species"):
+        """
+        Add extra dihedrals for chain molecules like allene, which
+        are required to cover all degrees of freedom
+
+        Args:
+            mol:
+
+        Returns:
+
+        """
+        assert mol.graph is not None
+        linear_chains = self._get_linear_chains(mol)
+
+        for chain in linear_chains:
+            o, p = chain[0], chain[-1]
+            for m in mol.graph.neighbors(o):
+                if m == p:
+                    continue
+
+                if m in chain:
+                    continue
+
+                for n in mol.graph.neighbors(p):
+                    if n == o:
+                        continue
+
+                    if n == m:
+                        continue
+
+                    if n in chain:
+                        continue
+
+                    if _is_dihedral_well_defined(mol, 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)
+            mol.graph.add_edge(*min_pair, pi=False, active=False)
+
+    assert mol.graph.is_connected, "Unknown error in connecting graph"
+
+    # 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
+
+
+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/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 764c322..52efca8 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -1,7 +1,8 @@
 import numpy as np
 
 from abc import ABC, abstractmethod
-from typing import Tuple, TYPE_CHECKING, List
+from enum import Enum
+from typing import Tuple, TYPE_CHECKING, List, Optional
 from autode.opt.coordinates._autodiff import (
     get_differentiable_vars,
     DifferentiableMath,
@@ -11,7 +12,7 @@ from autode.opt.coordinates._autodiff import (
 )
 
 if TYPE_CHECKING:
-    from autode.opt.coordinates import CartesianCoordinates, CartesianComponent
+    from autode.opt.coordinates import CartesianCoordinates
 
 
 def _get_3d_vecs_from_atom_idxs(
@@ -238,7 +239,7 @@ class PrimitiveInverseDistance(_DistanceFunction):
         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})"
@@ -260,7 +261,7 @@ class PrimitiveDistance(_DistanceFunction):
         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})"
@@ -298,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:
@@ -331,22 +332,22 @@ class PrimitiveBondAngle(Primitive):
 
 
 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
 
@@ -376,7 +377,7 @@ class PrimitiveDihedralAngle(Primitive):
         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))
@@ -405,3 +406,139 @@ class PrimitiveDihedralAngle(Primitive):
 
     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 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)
+        self.m = int(m)
+        self.o = int(o)
+        self.n = int(n)
+        self.r = int(r)
+
+        assert isinstance(axis, LinearBendType)
+        self.axis = axis
+
+    def __eq__(self, other):
+        """Equality of two linear bend angles"""
+        return isinstance(other, self.__class__) and (
+            self._atom_indexes == other._atom_indexes
+            and self.axis == other.axis
+        )
+
+    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())
+        else:
+            raise ValueError("Unknown axis for linear bend")
+
+
+class PrimitiveLinearAngle(LinearAngleBase):
+    """Linear Angle w.r.t. a reference atom"""
+
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
+    ):
+        """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(
+            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(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._vec_r: Optional[DifferentiableVector3D] = None
+
+    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]),
+            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-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))]
+
+        # 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
+    ):
+        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"
+        return f"LinearBend{axis_str}({self.m}-{self.o}-{self.n}, D)"
diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py
index 8dceacc..038f996 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
@@ -7,18 +8,14 @@ Notation follows:
 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
+from autode.opt.coordinates.internals import AnyPIC
 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,
+from autode.opt.optimisers.hessian_update import (
+    BFGSDampedUpdate,
+    BFGSSR1Update,
 )
 
 
@@ -41,7 +38,7 @@ class CRFOptimiser(RFOptimiser):
 
         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"""
@@ -125,15 +122,7 @@ class CRFOptimiser(RFOptimiser):
             )
 
         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 = AnyPIC.from_species(self._species)
 
         self._coords = DICWithConstraints.from_cartesian(
             x=cartesian_coords, primitives=primitives
@@ -141,48 +130,6 @@ class CRFOptimiser(RFOptimiser):
         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:
@@ -218,31 +165,3 @@ class CRFOptimiser(RFOptimiser):
 
         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)
diff --git a/autode/opt/optimisers/trm.py b/autode/opt/optimisers/trm.py
index 32b6741..a44a48f 100644
--- a/autode/opt/optimisers/trm.py
+++ b/autode/opt/optimisers/trm.py
@@ -130,27 +130,6 @@ class HybridTRMOptimiser(CRFOptimiser):
         ), "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
diff --git a/autode/values.py b/autode/values.py
index f361962..ea771c6 100644
--- a/autode/values.py
+++ b/autode/values.py
@@ -240,6 +240,10 @@ class Value(ABC, float):
     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/doc/changelog.rst b/doc/changelog.rst
index 1d2b977..3113ebb 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
 *********


Test patch:
diff --git a/tests/test_opt/molecules.py b/tests/test_opt/molecules.py
index 92936df..43956fc 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=[
@@ -41,3 +52,37 @@ 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.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.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),
+        ]
+    )
+
+
+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 9fe31f7..eb8e210 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -1,12 +1,26 @@
+import itertools
+
 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,
+    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
 from autode.exceptions import CoordinateTransformFailed
-from autode.opt.coordinates.base import CartesianComponent
-from autode.opt.coordinates.internals import PrimitiveInverseDistances, PIC
+from autode.opt.coordinates.internals import (
+    PrimitiveInverseDistances,
+    PIC,
+    AnyPIC,
+)
 from autode.opt.coordinates.cartesian import CartesianCoordinates
 from autode.opt.coordinates.dic import DIC
 from autode.opt.coordinates.primitives import (
@@ -16,6 +30,9 @@ from autode.opt.coordinates.primitives import (
     PrimitiveBondAngle,
     ConstrainedPrimitiveBondAngle,
     PrimitiveDihedralAngle,
+    PrimitiveLinearAngle,
+    PrimitiveDummyLinearAngle,
+    LinearBendType,
 )
 
 
@@ -527,12 +544,21 @@ def test_pic_b_no_primitives():
         c._calc_B(np.arange(6, dtype=float).reshape(2, 3))
 
 
-def test_pic_append_type_checking():
-    c = PIC()
-    # append should check for primitive type
-    c.append(PrimitiveDistance(0, 1))
+def test_pic_add_sanity_checking():
+    c = AnyPIC()
+    # pic add should check for primitive type
+    c.add(PrimitiveDistance(0, 1))
     with pytest.raises(AssertionError):
-        c.append(3)
+        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))
+    assert len(c) == 1
 
 
 def test_constrained_distance_satisfied():
@@ -558,7 +584,7 @@ def test_angle_primitive_derivative():
     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):
@@ -570,8 +596,8 @@ def test_angle_primitive_derivative():
 
 
 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():
@@ -612,6 +638,40 @@ 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 = PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.BEND)
+    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():
     # test that the primitive values are the same as the mol.distance etc.
     h2o2 = h2o2_mol()
@@ -620,14 +680,14 @@ 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)
 
 
 # fmt: off
-dihedral_mols = [
+extra_mols = [
     Molecule(
         atoms=[
             Atom("C", 0.63365, 0.11934, -0.13163),
@@ -643,18 +703,29 @@ dihedral_mols = [
             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),
+        ]
+    ),
+    feco5_mol(), # 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),
+    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)
+    PrimitiveDihedralAngle(2, 0, 1, 3),
+    PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.BEND),
+    PrimitiveLinearAngle(2, 3, 4, 8, LinearBendType.BEND),
 ]
 # fmt: on
 
@@ -706,9 +777,12 @@ 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),
+        PrimitiveDummyLinearAngle(0, 1, 2, LinearBendType.BEND),
     ]
 
     for p in prims:
@@ -733,7 +807,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)
@@ -742,8 +816,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
 
@@ -757,3 +831,121 @@ def test_dics_cannot_be_built_with_incomplete_primitives():
 
     with pytest.raises(RuntimeError):
         _ = DIC.from_cartesian(x=x, primitives=primitives)
+
+
+def test_pic_generation_linear_angle_ref():
+    # Fe(CO)5 with linear Fe-C-O bonds
+    m = feco5_mol()
+    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))
+    # check that linear bends use reference atoms, not dummy
+    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, 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())
+    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 = AnyPIC.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 = (
+        "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 = AnyPIC.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
+    # 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 = AnyPIC.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():
+    # extra dihedrals are needed for ends of linear chains like allene
+    cumulene = cumulene_mol()
+    pic = AnyPIC.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
+
+
+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 = AnyPIC.from_species(ptcl4)
+    _ = pic(ptcl4.coordinates.flatten())
+    assert np.linalg.matrix_rank(pic.B) == 3 * ptcl4.n_atoms - 6
diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py
index 50854b6..f040278 100644
--- a/tests/test_opt/test_crfo.py
+++ b/tests/test_opt/test_crfo.py
@@ -11,7 +11,7 @@ from autode.opt.optimisers.crfo import CRFOptimiser
 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, feco5_mol, cumulene_mol
 from ..testutils import requires_working_xtb_install
 
 
@@ -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
@@ -160,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)
@@ -199,25 +200,25 @@ 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)
     for pair in ((2, 0), (3, 0), (4, 1), (5, 1)):
-        pic.append(prim.PrimitiveDistance(*pair))
+        pic.add(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))
+        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
@@ -339,3 +340,41 @@ 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")
+
+
+@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=15, gtol=1e-5, etol=1e-6)
+    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")