Repository: duartegroup/autodE
Pull request #324 (2024-02-26 08:24:32)

V1.4.2


Patch:
diff --git a/README.md b/README.md
index 5d402b4..3b8e650 100644
--- a/README.md
+++ b/README.md
@@ -94,3 +94,4 @@ If **autodE** is used in a publication please consider citing the [paper](https:
 - Domen Pregeljc ([@dpregeljc](https://github.com/dpregeljc))
 - Jonathon Vandezande ([@jevandezande](https://github.com/jevandezande))
 - Shoubhik Maiti ([@shoubhikraj](https://github.com/shoubhikraj))
+- Daniel Hollas ([@danielhollas](https://github.com/danielhollas))
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/config.py b/autode/config.py
index fd534e4..48f2522 100644
--- a/autode/config.py
+++ b/autode/config.py
@@ -198,6 +198,9 @@ class _ConfigClass:
         #
         # Path can be unset and will be assigned if it can be found in $PATH
         path = None
+        #
+        # File extensions to copy when a calculation completes
+        copied_output_exts = [".out", ".hess", ".xyz", ".inp", ".pc"]
 
         optts_block = (
             "\n%geom\n"
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/neb/original.py b/autode/neb/original.py
index b0fd5d1..58060dd 100644
--- a/autode/neb/original.py
+++ b/autode/neb/original.py
@@ -3,7 +3,6 @@ The theory behind this original NEB implementation is taken from
 Henkelman and H. J ́onsson, J. Chem. Phys. 113, 9978 (2000)
 """
 import numpy as np
-import matplotlib.pyplot as plt
 
 from typing import Optional, Sequence, List, Any, TYPE_CHECKING, Union, Type
 from copy import deepcopy
@@ -361,6 +360,8 @@ class Images(Path):
         self, save=False, name="None", color=None, xlabel="NEB coordinate"
     ):
         """Plot the NEB surface"""
+        import matplotlib.pyplot as plt
+
         blues = plt.get_cmap("Blues")
 
         color = (
@@ -673,6 +674,8 @@ class NEB:
             etol_per_image: Energy tolerance per image to use in the L-BFGS-B
                             minimisation
         """
+        import matplotlib.pyplot as plt
+
         self.print_geometries(name=f"{name_prefix}neb_init")
 
         # Calculate energy on the first and final points as these will not be recalc-ed
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
new file mode 100644
index 0000000..0855a44
--- /dev/null
+++ b/autode/opt/coordinates/_autodiff.py
@@ -0,0 +1,731 @@
+"""
+Automatic differentiation routines in pure Python
+
+References:
+[1] P. Rehner, G. Bauer, Front. Chem. Eng., 2021, 3, 758090
+"""
+from typing import Union, Callable, Sequence, Optional
+from enum import Enum
+from copy import deepcopy
+import numpy as np
+import math
+
+numeric = (float, int)
+numeric_type = Union[float, int]
+
+
+class DerivativeOrder(Enum):
+    """Order of derivative"""
+
+    zeroth = 0
+    first = 1
+    second = 2
+
+
+def get_differentiable_vars(
+    values: Sequence[numeric_type],
+    symbols: Sequence[str],
+    deriv_order: DerivativeOrder = DerivativeOrder.second,
+):
+    """
+    Obtain differentiable variables from a series of numbers
+
+    Args:
+        values: The values of the variables (numbers)
+        symbols: List of symbols (strings) of the numbers
+        deriv_order: Order of differentiation
+
+    Returns:
+        (list[VectorHyperDual]): A list of hyper dual numbers
+    """
+    assert all(isinstance(sym, str) for sym in symbols)
+    assert len(symbols) == len(values)
+    symbols = list(symbols)
+
+    hyperduals = []
+    for symbol, value in zip(symbols, values):
+        var = VectorHyperDual.from_variable(
+            value, symbol, all_symbols=symbols, order=deriv_order
+        )
+        hyperduals.append(var)
+
+    return hyperduals
+
+
+class VectorHyperDual:
+    """
+    Hyper-dual numbers with vector infinitesimals upto the
+    second order (i.e., upto second partial derivatives)
+    """
+
+    def __init__(
+        self,
+        value: numeric_type,
+        symbols: Sequence[str],
+        first_der: Optional[np.ndarray] = None,
+        second_der: Optional[np.ndarray] = None,
+    ):
+        """
+        Create a vector hyper dual number, i.e. a scalar function
+        with one or more variables
+
+        Args:
+            value: The scalar value of the hyper-dual
+            symbols: A list of unique strings representing the variables
+            first_der: 1D array of first derivatives against variables
+            second_der: 2D matrix of second derivatives against variables
+        """
+        assert isinstance(value, numeric)
+        self._val = float(value)
+
+        assert all(isinstance(sym, str) for sym in symbols)
+        if len(set(symbols)) != len(list(symbols)):
+            raise RuntimeError("Symbols must be unique!")
+        self._symbols = tuple(symbols)
+
+        # load the derivatives with sanity checks
+        self._first_der: Optional[np.ndarray] = None
+        self._second_der: Optional[np.ndarray] = None
+        self._order = DerivativeOrder.zeroth
+        self._init_deriv_arrays(first_der, second_der)
+
+    def _init_deriv_arrays(
+        self, first_der: Optional[np.ndarray], second_der: Optional[np.ndarray]
+    ) -> None:
+        """
+        Initialise the derivative matrices, checking they have the
+        correct shape, and set the derivative order for this
+        hyper-dual number
+        """
+        if first_der is None:
+            return None
+        assert isinstance(first_der, np.ndarray)
+        first_der = first_der.flatten()
+        if not first_der.shape == (self.n_vars,):
+            raise ValueError(
+                f"Number of symbols ({self.n_vars}) does not match with"
+                f" shape of derivative array {first_der.shape}"
+            )
+        self._first_der = first_der.astype(float)
+        self._order = DerivativeOrder.first
+
+        if second_der is None:
+            return None
+        assert isinstance(second_der, np.ndarray)
+        if not second_der.shape == (self.n_vars, self.n_vars):
+            raise ValueError(
+                f"Number of symbols ({self.n_vars}) does not match with"
+                f" shape of second derivative matrix {first_der.shape}"
+            )
+        self._second_der = second_der.astype(float)
+        self._order = DerivativeOrder.second
+
+    def __repr__(self):
+        rstring = f"HyperDual({self.value}"
+        if self._order in [DerivativeOrder.first, DerivativeOrder.second]:
+            rstring += f", f'[{self.n_vars}]"
+        if self._order == DerivativeOrder.second:
+            rstring += f', f"[{self.n_vars}, {self.n_vars}]'
+        rstring += ")"
+        return rstring
+
+    @property
+    def n_vars(self) -> int:
+        """Number of variables in this hyper-dual"""
+        return len(self._symbols)
+
+    def copy(self) -> "VectorHyperDual":
+        return deepcopy(self)
+
+    def _check_compatible(self, other: "VectorHyperDual") -> None:
+        """
+        Check the compatibility of two VectorHyperDual numbers for
+        any operation that involves the two.
+
+        Args:
+            other (VectorHyperDual):
+
+        Raises:
+            (ValueError): If they are incompatible
+        """
+        if self.n_vars != other.n_vars:
+            raise ValueError(
+                "Incompatible number of differentiable variables, "
+                "cannot perform operation"
+            )
+        if self._symbols != other._symbols:
+            raise ValueError(
+                "The differentiable variable symbols do not match, "
+                "cannot perform operation!"
+            )
+        if self._order != other._order:
+            raise ValueError("The order of derivative do not match!")
+        return None
+
+    @property
+    def value(self) -> float:
+        """Return the value of the hyper-dual number"""
+        return self._val
+
+    @value.setter
+    def value(self, value: float):
+        assert isinstance(value, numeric)
+        self._val = float(value)
+
+    @classmethod
+    def from_variable(
+        cls,
+        value: float,
+        symbol: str,
+        all_symbols: Sequence[str],
+        order: DerivativeOrder,
+    ):
+        """
+        Create a hyper-dual number from one variable, requires
+        list of symbols and the symbol of this specific variable.
+        Essentially a variable x can be considered as a scalar function
+        of a list of variables - 1 * x + 0 * y + 0 * z + ...
+
+        Args:
+            value: The value of the variable (will be converted to float)
+            symbol: The symbol of the current variable, must be in all_symbols
+            all_symbols: List of strings indicating all required variables
+            order: The order of differentiation to consider
+
+        Returns:
+            (VectorHyperDual): The hyper-dual representing the variable
+        """
+        assert all(isinstance(sym, str) for sym in all_symbols)
+        assert isinstance(symbol, str)
+        assert symbol in all_symbols
+
+        val = float(value)
+        first_der = None
+        second_der = None
+        n = len(all_symbols)
+        idx = list(all_symbols).index(symbol)
+        order = DerivativeOrder(order)
+
+        if order == DerivativeOrder.first or order == DerivativeOrder.second:
+            first_der = np.zeros(shape=n, dtype=float)
+            first_der[idx] = 1.0
+        if order == DerivativeOrder.second:
+            second_der = np.zeros(shape=(n, n), dtype=float)
+
+        return VectorHyperDual(val, all_symbols, first_der, second_der)
+
+    def differentiate_wrt(
+        self,
+        symbol1: str,
+        symbol2: Union[str, None] = None,
+    ) -> Optional[float]:
+        """
+        Derivative of this hyper-dual number (scalar function) against one
+        or two variable(s) identified by their string(s).
+
+        Args:
+            symbol1 (str):
+            symbol2 (str|None):
+
+        Returns:
+            (float|None): The derivative value, or None if not available
+        """
+        assert isinstance(symbol1, str)
+        if symbol1 not in self._symbols:
+            return None
+
+        if self._order == DerivativeOrder.zeroth:
+            return None
+
+        idx_1 = self._symbols.index(symbol1)
+        assert self._first_der is not None
+        if symbol2 is None:
+            return self._first_der[idx_1]
+
+        assert isinstance(symbol2, str)
+        if symbol2 not in self._symbols:
+            return None
+        idx_2 = self._symbols.index(symbol2)
+        # check if second derivs are available
+        if self._order == DerivativeOrder.first:
+            return None
+        assert self._second_der is not None
+        return self._second_der[idx_1, idx_2]
+
+    def __add__(
+        self, other: Union["VectorHyperDual", numeric_type]
+    ) -> "VectorHyperDual":
+        """Adding a hyper dual number"""
+
+        if isinstance(other, numeric):
+            new = self.copy()
+            new._val += float(other)
+            return new
+
+        # add to another dual number
+        elif isinstance(other, VectorHyperDual):
+            self._check_compatible(other)
+
+            val = self._val + other._val
+            if self._order == DerivativeOrder.zeroth:
+                return VectorHyperDual(val, self._symbols)
+
+            assert self._first_der is not None
+            assert other._first_der is not None
+            first_der = self._first_der + other._first_der
+            if self._order == DerivativeOrder.first:
+                return VectorHyperDual(val, self._symbols, first_der)
+
+            assert self._second_der is not None
+            assert other._second_der is not None
+            second_der = self._second_der + other._second_der
+            return VectorHyperDual(val, self._symbols, first_der, second_der)
+
+        else:
+            raise TypeError("Unknown type for addition")
+
+    def __radd__(self, other):
+        """Addition is commutative"""
+        return self.__add__(other)
+
+    def __neg__(self) -> "VectorHyperDual":
+        """Unary negative operation"""
+        new = self.copy()
+        new._val = -new._val
+        if self._order == DerivativeOrder.first:
+            assert new._first_der is not None
+            new._first_der = -new._first_der
+        elif self._order == DerivativeOrder.second:
+            assert new._first_der is not None
+            assert new._second_der is not None
+            new._first_der = -new._first_der
+            new._second_der = -new._second_der
+        return new
+
+    def __sub__(self, other):
+        """Subtraction of hyper dual numbers"""
+        return self.__add__(-other)
+
+    def __rsub__(self, other):
+        """Reverse subtraction"""
+        return other + (-self)
+
+    def __mul__(self, other) -> "VectorHyperDual":
+        """Multiply a hyper dual number with float or another hyper dual"""
+        if isinstance(other, numeric):
+            new = self.copy()
+            new._val *= float(other)
+            return new
+
+        # Product rule for derivatives, Eqn (24) in ref. [1]
+        elif isinstance(other, VectorHyperDual):
+            self._check_compatible(other)
+
+            val = self._val * other._val
+            if self._order == DerivativeOrder.zeroth:
+                return VectorHyperDual(val, self._symbols)
+
+            assert self._first_der is not None
+            assert other._first_der is not None
+            first_der = (
+                self._val * other._first_der + other._val * self._first_der
+            )
+            if self._order == DerivativeOrder.first:
+                return VectorHyperDual(val, self._symbols, first_der)
+
+            assert self._second_der is not None
+            assert other._second_der is not None
+            second_der = (
+                self._val * other._second_der
+                + np.outer(self._first_der, other._first_der)
+                + np.outer(other._first_der, self._first_der)
+                + other._val * self._second_der
+            )
+            return VectorHyperDual(val, self._symbols, first_der, second_der)
+        else:
+            raise TypeError("Unknown type for multiplication")
+
+    def __rmul__(self, other):
+        return self.__mul__(other)
+
+    def __truediv__(self, other):
+        """True division, defined by multiplicative inverse"""
+        return self.__mul__(DifferentiableMath.pow(other, -1))
+
+    def __rtruediv__(self, other):
+        """Reverse true division"""
+        return DifferentiableMath.pow(self, -1).__mul__(other)
+
+    def __pow__(self, power, modulo=None) -> "VectorHyperDual":
+        if modulo is not None:
+            raise NotImplementedError("Modulo inverse is not implemented")
+
+        result = DifferentiableMath.pow(self, power)
+        assert isinstance(result, VectorHyperDual)
+        return result
+
+    def __rpow__(self, other):
+        return DifferentiableMath.pow(other, self)
+
+    @staticmethod
+    def apply_operation(
+        num: Union["VectorHyperDual", numeric_type],
+        operator: Callable[[float], float],
+        operator_first_deriv: Callable[[float], float],
+        operator_second_deriv: Callable[[float], float],
+    ) -> Union["VectorHyperDual", numeric_type]:
+        """
+        Perform an operation on the hyperdual (i.e. apply a scalar function),
+        also compatible with Python numeric (float/int) types.
+
+        Args:
+            num: Number that is hyper-dual (or float/int)
+            operator: Function that returns the value (result) of operation
+            operator_first_deriv: Should return first derivative of operation
+            operator_second_deriv: Should return second derivative of operation
+
+        Returns:
+            (VectorHyperDual|float): The result
+        """
+        # pass through numeric types
+        if isinstance(num, numeric):
+            return operator(float(num))
+
+        assert isinstance(num, VectorHyperDual)
+
+        val = operator(num._val)
+
+        if num._order == DerivativeOrder.zeroth:
+            return VectorHyperDual(val, num._symbols)
+
+        # Eqn (25) in reference [1]
+        assert num._first_der is not None
+        f_dash_x0 = operator_first_deriv(num._val)
+        first_der = num._first_der * f_dash_x0
+
+        if num._order == DerivativeOrder.first:
+            return VectorHyperDual(val, num._symbols, first_der)
+
+        assert num._second_der is not None
+        second_der = np.outer(num._first_der, num._first_der) * (
+            operator_second_deriv(num._val)
+        )
+        second_der += num._second_der * f_dash_x0
+        return VectorHyperDual(val, num._symbols, first_der, second_der)
+
+
+class DifferentiableMath:
+    """
+    Class defining math functions that can be used on
+    hyper dual numbers (i.e. differentiable functions),
+    as well as standard numeric types (float and int)
+    """
+
+    @staticmethod
+    def sqrt(
+        num: Union[VectorHyperDual, numeric_type]
+    ) -> Union[VectorHyperDual, numeric_type]:
+        """Calculate the square root of a hyperdual number"""
+
+        if isinstance(num, numeric):
+            assert num > 0
+        else:
+            assert num.value > 0
+
+        return VectorHyperDual.apply_operation(
+            num,
+            operator=lambda x0: math.sqrt(x0),
+            operator_first_deriv=lambda x0: 1 / (2 * math.sqrt(x0)),
+            operator_second_deriv=lambda x0: -1 / (4 * math.pow(x0, 3 / 2)),
+        )
+
+    @staticmethod
+    def exp(
+        num: Union[VectorHyperDual, numeric_type]
+    ) -> Union[VectorHyperDual, numeric_type]:
+        """Raise e to the power of num"""
+
+        return VectorHyperDual.apply_operation(
+            num,
+            operator=lambda x0: math.exp(x0),
+            operator_first_deriv=lambda x0: math.exp(x0),
+            operator_second_deriv=lambda x0: math.exp(x0),
+        )
+
+    @staticmethod
+    def pow(
+        num: Union[VectorHyperDual, numeric_type],
+        power: Union[VectorHyperDual, numeric_type],
+    ) -> Union[VectorHyperDual, numeric_type]:
+        """Exponentiation of one hyperdual to another"""
+
+        if isinstance(num, numeric) and isinstance(power, numeric):
+            return math.pow(num, power)
+
+        elif isinstance(num, VectorHyperDual) and isinstance(power, numeric):
+            if num.value < 0 and isinstance(power, float):
+                raise AssertionError(
+                    "Math error, can't raise negative number to fractional power"
+                )
+            return VectorHyperDual.apply_operation(
+                num,
+                operator=lambda x0: math.pow(x0, power),  # type: ignore
+                operator_first_deriv=lambda x0: power  # type: ignore
+                * math.pow(x0, power - 1),
+                operator_second_deriv=lambda x0: power  # type: ignore
+                * (power - 1)
+                * math.pow(x0, power - 2),
+            )
+
+        elif isinstance(power, VectorHyperDual) and isinstance(
+            num, (numeric, VectorHyperDual)
+        ):
+            if (isinstance(num, numeric) and num < 0) or (
+                isinstance(num, VectorHyperDual) and num.value < 0
+            ):
+                raise AssertionError(
+                    "Only positive numbers can be used with"
+                    " differentiable exponent"
+                )
+            # use identity x^y = e^(y log_x) for x > 0
+            return DifferentiableMath.exp(power * DifferentiableMath.log(num))
+
+        else:
+            raise TypeError("Unknown type for exponentiation")
+
+    @staticmethod
+    def log(
+        num: Union[VectorHyperDual, numeric_type]
+    ) -> Union[VectorHyperDual, numeric_type]:
+        """Natural logarithm"""
+
+        if isinstance(num, numeric):
+            assert num > 0
+        else:
+            assert num.value > 0
+
+        return VectorHyperDual.apply_operation(
+            num,
+            operator=lambda x0: math.log(x0),
+            operator_first_deriv=lambda x0: 1.0 / x0,
+            operator_second_deriv=lambda x0: -1.0 / (x0**2),
+        )
+
+    @staticmethod
+    def acos(
+        num: Union[VectorHyperDual, numeric_type]
+    ) -> Union[VectorHyperDual, numeric_type]:
+        """Calculate the arccosine of a hyperdual number"""
+
+        if isinstance(num, VectorHyperDual):
+            assert -1 < num.value < 1
+        else:
+            assert -1 < num < 1
+
+        return VectorHyperDual.apply_operation(
+            num,
+            operator=lambda x0: math.acos(x0),
+            operator_first_deriv=lambda x0: -1 / math.sqrt(1 - x0**2),
+            operator_second_deriv=lambda x0: -x0
+            / math.pow(1 - x0**2, 3 / 2),
+        )
+
+    @staticmethod
+    def atan(
+        num: Union[VectorHyperDual, numeric_type]
+    ) -> Union[VectorHyperDual, numeric_type]:
+        """Calculate the arctangent of a hyperdual number"""
+
+        return VectorHyperDual.apply_operation(
+            num,
+            operator=lambda x0: math.atan(x0),
+            operator_first_deriv=lambda x0: 1 / (1 + x0**2),
+            operator_second_deriv=lambda x0: (-2 * x0) / (x0**2 + 1) ** 2,
+        )
+
+    @staticmethod
+    def atan2(
+        num_y: Union[VectorHyperDual, numeric_type],
+        num_x: Union[VectorHyperDual, numeric_type],
+    ) -> Union[VectorHyperDual, numeric_type]:
+        """Calculate the arctan2 of two hyper dual numbers"""
+        if isinstance(num_y, numeric) and isinstance(num_x, numeric):
+            return math.atan2(num_y, num_x)
+
+        # https://en.wikipedia.org/wiki/Atan2 four overlapping half-planes
+        def atan2_derivs_x_not_0(y, x):
+            return DifferentiableMath.atan(y / x)
+
+        def atan2_derivs_x_close_0(y, x):
+            return -DifferentiableMath.atan(x / y)
+
+        x_val = float(num_x) if isinstance(num_x, numeric) else num_x.value
+        y_val = float(num_y) if isinstance(num_y, numeric) else num_y.value
+        res_val = math.atan2(y_val, x_val)
+
+        # when atan2(y,x)->pi/2, x->0 or y/x->inf, use other formula for derivs
+        if math.isclose(abs(res_val), math.pi / 2, abs_tol=0.1):
+            res = atan2_derivs_x_close_0(num_y, num_x)
+            res.value = res_val
+            return res
+        else:
+            res = atan2_derivs_x_not_0(num_y, num_x)
+            res.value = res_val
+            return res
+
+
+class DifferentiableVector3D:
+    """
+    Convenience class to represent a 3D vector of differentiable
+    hyper-dual numbers
+    """
+
+    def __init__(
+        self, items: Sequence[Union["VectorHyperDual", numeric_type]]
+    ):
+        """
+        Initialise the 3D vector from a list of 3 hyperdual numbers
+
+        Args:
+            items: A list of 3 hyper-dual numbers
+        """
+        items = list(items)
+        if len(items) != 3:
+            raise ValueError("A 3D vector must have only 3 components")
+        assert all(
+            isinstance(item, (VectorHyperDual, *numeric)) for item in items
+        )
+        self._data = items
+
+    @staticmethod
+    def _check_same_type(other) -> None:
+        """Check that another object is also a 3D differentiable vector"""
+        if not isinstance(other, DifferentiableVector3D):
+            raise ValueError("Operation must be done with another 3D vector!")
+        return None
+
+    def dot(
+        self, other: "DifferentiableVector3D"
+    ) -> Union["VectorHyperDual", numeric_type]:
+        """
+        Dot product of two 3D vectors
+
+        Args:
+            other (DifferentiableVector3D):
+
+        Returns:
+            (VectorHyperDual): A scalar number (with derivatives)
+        """
+        self._check_same_type(other)
+        dot: Union[VectorHyperDual, numeric_type] = 0
+        for k in range(3):
+            dot = dot + self._data[k] * other._data[k]
+        return dot
+
+    def norm(self) -> Union["VectorHyperDual", numeric_type]:
+        """
+        Euclidean (l2) norm of this 3D vector
+
+        Returns:
+            (VectorHyperDual): A scalar number (with derivatives)
+        """
+        norm = DifferentiableMath.sqrt(
+            self._data[0] ** 2 + self._data[1] ** 2 + self._data[2] ** 2
+        )
+        return norm
+
+    def __add__(
+        self, other: "DifferentiableVector3D"
+    ) -> "DifferentiableVector3D":
+        """
+        Vector addition in 3D, returns a vector
+
+        Args:
+            other (DifferentiableVector3D):
+
+        Returns:
+            (DifferentiableVector3D):
+        """
+        self._check_same_type(other)
+        return DifferentiableVector3D(
+            [self._data[k] + other._data[k] for k in range(3)]
+        )
+
+    def __neg__(self) -> "DifferentiableVector3D":
+        """
+        Unary negation of a vector, returns another vector
+
+        Returns:
+            (DifferentiableVector3D):
+        """
+        return DifferentiableVector3D([-self._data[k] for k in range(3)])
+
+    def __sub__(self, other) -> "DifferentiableVector3D":
+        """
+        Vector subtraction in 3D, defined in terms of addition
+        and negation
+
+        Args:
+            other (DifferentiableVector3D):
+
+        Returns:
+            (DifferentiableVector3D):
+        """
+        return self.__add__(-other)
+
+    def __mul__(
+        self, other: Union[VectorHyperDual, numeric_type]
+    ) -> "DifferentiableVector3D":
+        """
+        Multiplication of a 3D vector with a scalar
+
+        Args:
+            other (VectorHyperDual|float|int):
+
+        Returns:
+            (DifferentiableVector3D):
+        """
+        assert isinstance(other, numeric) or isinstance(other, VectorHyperDual)
+        return DifferentiableVector3D(
+            [self._data[k] * other for k in range(3)]
+        )
+
+    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":
+        """
+        Cross-product of two 3D vectors, produces another vector
+
+        Args:
+            other (DifferentiableVector3D):
+
+        Returns:
+            (DifferentiableVector3D):
+        """
+        self._check_same_type(other)
+        return DifferentiableVector3D(
+            [
+                self._data[1] * other._data[2]
+                - self._data[2] * other._data[1],
+                self._data[2] * other._data[0]
+                - self._data[0] * other._data[2],
+                self._data[0] * other._data[1]
+                - self._data[1] * other._data[0],
+            ]
+        )
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/dic.py b/autode/opt/coordinates/dic.py
index b857f2a..2d52778 100644
--- a/autode/opt/coordinates/dic.py
+++ b/autode/opt/coordinates/dic.py
@@ -216,34 +216,42 @@ class DIC(InternalCoordinates):  # lgtm [py/missing-equals]
         q_init = self.primitives(x_k)
         x_1 = self.to("cartesian") + np.matmul(self.B_T_inv, value)
 
+        success = False
         for i in range(1, _max_back_transform_iterations + 1):
-            x_k = x_k + np.matmul(self.B_T_inv, (s_new - s_k))
+            try:
+                x_k = x_k + np.matmul(self.B_T_inv, (s_new - s_k))
 
-            # Rebuild the primitives & DIC from the back-transformed Cartesians
-            q_k = self.primitives.close_to(x_k, q_init)
-            s_k = np.matmul(self.U.T, q_k)
-            self.B = np.matmul(self.U.T, self.primitives.B)
-            self.B_T_inv = np.linalg.pinv(self.B)
+                # Rebuild the primitives & DIC from the back-transformed Cartesians
+                q_k = self.primitives.close_to(x_k, q_init)
+                s_k = np.matmul(self.U.T, q_k)
+                self.B = np.matmul(self.U.T, self.primitives.B)
+                self.B_T_inv = np.linalg.pinv(self.B)
 
-            rms_s = np.sqrt(np.mean(np.square(s_k - s_new)))
+                rms_s = np.sqrt(np.mean(np.square(s_k - s_new)))
+
+            # for ill-conditioned primitives, there might be math error
+            except ArithmeticError:
+                break
 
             if rms_s < 1e-10:
-                logger.info(
-                    f"DIC transformation converged in {i} cycle(s) "
-                    f"in {time() - start_time:.4f} s"
-                )
+                success = True
                 break
 
-            if i == _max_back_transform_iterations:
-                logger.warning(
-                    f"Failed to transform in {i} cycles. "
-                    f"Final RMS(s) = {rms_s:.8f}"
+        if success:
+            logger.info(
+                f"DIC transformation converged in {i} cycle(s) "
+                f"in {time() - start_time:.4f} s"
+            )
+        else:
+            logger.warning(
+                f"Failed to transform in {i} cycles. "
+                f"Final RMS(s) = {rms_s:.8f}"
+            )
+            x_k = x_1
+            if not self.allow_unconverged_back_transform:
+                raise CoordinateTransformFailed(
+                    "DIC->Cart iterative back-transform did not converge"
                 )
-                x_k = x_1
-                if not self.allow_unconverged_back_transform:
-                    raise CoordinateTransformFailed(
-                        "DIC->Cart iterative back-transform did not converge"
-                    )
 
         s_k = np.matmul(self.U.T, self.primitives(x_k))
         self.B = np.matmul(self.U.T, self.primitives.B)
diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 675ca1f..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,6 +98,19 @@ 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"""
+        raise NotImplementedError(
+            "Please use PIC.add() to add new primitives to the set"
+        )
+
     @property
     def B(self) -> np.ndarray:
         """Wilson B matrix"""
@@ -173,22 +202,12 @@ class PIC(list, ABC):
                 "primitive internal coordinates"
             )
 
-        cart_coords = x.reshape((-1, 3))
+        cart_coords = x.ravel()
 
-        n_atoms, _ = cart_coords.shape
-        B = np.zeros(shape=(len(self), 3 * n_atoms))
+        B = np.zeros(shape=(len(self), len(cart_coords)))
 
         for i, primitive in enumerate(self):
-            for j in range(n_atoms):
-                B[i, 3 * j + 0] = primitive.derivative(
-                    j, CartesianComponent.x, x=cart_coords
-                )
-                B[i, 3 * j + 1] = primitive.derivative(
-                    j, CartesianComponent.y, x=cart_coords
-                )
-                B[i, 3 * j + 2] = primitive.derivative(
-                    j, CartesianComponent.z, x=cart_coords
-                )
+            B[i] = primitive.derivative(x=cart_coords)
 
         self._B = B
         return None
@@ -215,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
 
@@ -239,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 22bfabb..52efca8 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -1,10 +1,56 @@
 import numpy as np
 
 from abc import ABC, abstractmethod
-from typing import Tuple, TYPE_CHECKING
+from enum import Enum
+from typing import Tuple, TYPE_CHECKING, List, Optional
+from autode.opt.coordinates._autodiff import (
+    get_differentiable_vars,
+    DifferentiableMath,
+    DifferentiableVector3D,
+    DerivativeOrder,
+    VectorHyperDual,
+)
 
 if TYPE_CHECKING:
-    from autode.opt.coordinates import CartesianCoordinates, CartesianComponent
+    from autode.opt.coordinates import CartesianCoordinates
+
+
+def _get_3d_vecs_from_atom_idxs(
+    *args: int,
+    x: "CartesianCoordinates",
+    deriv_order: DerivativeOrder,
+) -> List[DifferentiableVector3D]:
+    """
+    Obtain differentiable 3D vectors from the Cartesian components
+    of each atom, given by atomic indices in order. The symbols are
+    strings denoting their position in the flat Cartesian coordinate.
+
+    Args:
+        *args: Integers denoting the atom positions
+        x: Cartesian coordinate array
+        deriv_order: Order of derivatives for initialising variables
+
+    Returns:
+        (list[VectorHyperDual]): A list of differentiable variables
+    """
+    assert all(isinstance(idx, int) and idx >= 0 for idx in args)
+    # get positions in the flat Cartesian array
+    _x = x.ravel()
+    cart_idxs = []
+    for atom_idx in args:
+        for k in range(3):
+            cart_idxs.append(3 * atom_idx + k)
+    variables = get_differentiable_vars(
+        values=[_x[idx] for idx in cart_idxs],
+        symbols=[str(idx) for idx in cart_idxs],
+        deriv_order=deriv_order,
+    )
+    atom_vecs = []
+    for pos_idx in range(len(args)):
+        atom_vecs.append(
+            DifferentiableVector3D(variables[pos_idx * 3 : pos_idx * 3 + 3])
+        )
+    return atom_vecs
 
 
 class Primitive(ABC):
@@ -17,18 +63,35 @@ class Primitive(ABC):
         self._atom_indexes = atom_indexes
 
     @abstractmethod
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
+    ):
+        """
+        The function that performs the main evaluation of the PIC,
+        and optionally returns derivative or second derivatives.
+        The returned hyper-dual must have the proper cartesian idxs
+        set.
+
+        Args:
+            x: Cartesian coordinates
+            deriv_order: The order of derivatives requested - 0, 1 or 2
+
+        Returns:
+            (VectorHyperDual): The result, optionally containing derivatives
+        """
+
     def __call__(self, x: "CartesianCoordinates") -> float:
         """Return the value of this PIC given a set of cartesian coordinates"""
+        _x = x.ravel()
+        res = self._evaluate(_x, deriv_order=DerivativeOrder.zeroth)
+        return res.value
 
-    @abstractmethod
     def derivative(
         self,
-        i: int,
-        component: "CartesianComponent",
         x: "CartesianCoordinates",
-    ) -> float:
+    ) -> np.ndarray:
         r"""
-        Calculate the derivative with respect to a cartesian coordinate
+        Calculate the derivatives with respect to cartesian coordinates
 
         .. math::
 
@@ -40,17 +103,56 @@ class Primitive(ABC):
 
         -----------------------------------------------------------------------
         Arguments:
-            i: Cartesian index to take the derivative with respect to. [0-N),
-               for N atoms
 
-            component: Cartesian component (x, y, z) to take the derivative
-                       with respect to
+            x: Cartesian coordinate array of shape (N, )
 
-            x: Cartesian coordinates
+        Returns:
+            (np.ndarray): Derivative array of shape (N, )
+        """
+        _x = x.ravel()
+        res = self._evaluate(_x, deriv_order=DerivativeOrder.first)
+        derivs = np.zeros_like(_x, dtype=float)
+        for i in range(_x.shape[0]):
+            dqdx_i = res.differentiate_wrt(str(i))
+            if dqdx_i is not None:
+                derivs[i] = dqdx_i
+
+        return derivs
+
+    def second_derivative(
+        self,
+        x: "CartesianCoordinates",
+    ) -> np.ndarray:
+        r"""
+        Calculate the second derivatives with respect to cartesian coordinates
+
+        .. math::
+
+            \frac{d^2 q}
+                  {d\boldsymbol{X}_{i, k}^2} {\Bigg\rvert}_{X=X0}
+
+        where :math:`q` is the primitive coordinate and :math:`\boldsymbol{X}`
+        are the cartesian coordinates.
+
+        -----------------------------------------------------------------------
+        Arguments:
+
+            x: Cartesian coordinate array of shape (N, )
 
         Returns:
-            (float): Derivative
+            (np.ndarray): Second derivative matrix of shape (N, N)
         """
+        _x = x.ravel()
+        x_n = _x.shape[0]
+        res = self._evaluate(_x, deriv_order=DerivativeOrder.second)
+        derivs = np.zeros(shape=(x_n, x_n), dtype=float)
+        for i in range(x_n):
+            for j in range(x_n):
+                d2q_dx2_ij = res.differentiate_wrt(str(i), str(j))
+                if d2q_dx2_ij is not None:
+                    derivs[i, j] = d2q_dx2_ij
+
+        return derivs
 
     @abstractmethod
     def __eq__(self, other):
@@ -121,8 +223,8 @@ class _DistanceFunction(Primitive, ABC):
 
 
 class PrimitiveInverseDistance(_DistanceFunction):
-    """
-    Inverse distance between to atoms:
+    r"""
+    Inverse distance between two atoms:
 
     .. math::
 
@@ -130,40 +232,21 @@ class PrimitiveInverseDistance(_DistanceFunction):
                 {|\boldsymbol{X}_i - \boldsymbol{X}_j|}
     """
 
-    def derivative(
-        self,
-        i: int,
-        component: "CartesianComponent",
-        x: "CartesianCoordinates",
-    ) -> float:
-        """
-        Derivative with respect to Cartesian displacement
-
-        -----------------------------------------------------------------------
-        See Also:
-            :py:meth:`Primitive.derivative <Primitive.derivative>`
-        """
-
-        _x = x.reshape((-1, 3))
-        k = int(component)
-
-        if i != self.i and i != self.j:
-            return 0  # Atom does not form part of this distance
-
-        elif i == self.i:
-            return -(_x[i, k] - _x[self.j, k]) * self(x) ** 3
-
-        else:  # i == self.idx_j:
-            return (_x[self.i, k] - _x[self.j, k]) * self(x) ** 3
-
-    def __call__(self, x: "CartesianCoordinates") -> float:
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
+    ) -> "VectorHyperDual":
         """1 / |x_i - x_j|"""
-        _x = x.reshape((-1, 3))
-        return 1.0 / np.linalg.norm(_x[self.i] - _x[self.j])
+        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()  # type: ignore
+
+    def __repr__(self):
+        return f"InverseDistance({self.i}-{self.j})"
 
 
 class PrimitiveDistance(_DistanceFunction):
-    """
+    r"""
     Distance between two atoms:
 
     .. math::
@@ -171,33 +254,14 @@ class PrimitiveDistance(_DistanceFunction):
         q = |\boldsymbol{X}_i - \boldsymbol{X}_j|
     """
 
-    def derivative(
-        self,
-        i: int,
-        component: "CartesianComponent",
-        x: "CartesianCoordinates",
-    ) -> float:
-        """
-        Derivative with respect to Cartesian displacement
-
-        -----------------------------------------------------------------------
-        See Also:
-            :py:meth:`Primitive.derivative <Primitive.derivative>`
-        """
-        _x = x.reshape((-1, 3))
-        k = int(component)
-
-        if i != self.i and i != self.j:
-            return 0  # Atom does not form part of this distance
-
-        val = (_x[self.i, k] - _x[self.j, k]) / self(x)
-
-        return val if i == self.i else -val
-
-    def __call__(self, x: "CartesianCoordinates") -> float:
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
+    ) -> "VectorHyperDual":
         """|x_i - x_j|"""
-        _x = x.reshape((-1, 3))
-        return np.linalg.norm(_x[self.i] - _x[self.j])
+        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()  # type: ignore
 
     def __repr__(self):
         return f"Distance({self.i}-{self.j})"
@@ -230,13 +294,18 @@ class ConstrainedPrimitiveDistance(ConstrainedPrimitive, PrimitiveDistance):
 
 
 class PrimitiveBondAngle(Primitive):
-    def __init__(self, o: int, m: int, n: int):
+    """
+    Bond angle between three atoms, calculated with the
+    arccosine of the normalised dot product
+    """
+
+    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 = o
-        self.m = m
-        self.n = n
+        self.m = int(m)
+        self.o = int(o)
+        self.n = int(n)
 
     def __eq__(self, other) -> bool:
         """Equality of two distance functions"""
@@ -247,80 +316,38 @@ class PrimitiveBondAngle(Primitive):
             and other._ordered_idxs == self._ordered_idxs
         )
 
-    def __call__(self, x: "CartesianCoordinates") -> float:
-        _x = x.reshape((-1, 3))
-        u = _x[self.m, :] - _x[self.o, :]
-        v = _x[self.n, :] - _x[self.o, :]
-
-        theta = np.arccos(u.dot(v) / (np.linalg.norm(u) * np.linalg.norm(v)))
-        return theta
-
-    def derivative(
-        self,
-        i: int,
-        component: "CartesianComponent",
-        x: "CartesianCoordinates",
-    ) -> float:
-        if i not in (self.o, self.m, self.n):
-            return 0.0
-
-        k = int(component)
-
-        _x = x.reshape((-1, 3))
-        u = _x[self.m, :] - _x[self.o, :]
-        lambda_u = np.linalg.norm(u)
-        u /= lambda_u
-
-        v = _x[self.n, :] - _x[self.o, :]
-        lambda_v = np.linalg.norm(v)
-        v /= lambda_v
-
-        t0, t1 = np.array([1.0, -1.0, 1.0]), np.array([-1.0, 1.0, 1.0])
-
-        if not np.isclose(np.abs(np.arccos(u.dot(v))), 1.0):
-            w = np.cross(u, v)
-        elif not np.isclose(
-            np.abs(np.arccos(u.dot(t0))), 1.0
-        ) and not np.isclose(np.abs(np.arccos(v.dot(t0))), 1.0):
-            w = np.cross(u, t0)
-        else:
-            w = np.cross(u, t1)
-
-        w /= np.linalg.norm(w)
-
-        dqdx = 0.0
-
-        if i in (self.m, self.o):
-            sign = 1 if i == self.m else -1
-            dqdx += sign * np.cross(u, w)[k] / lambda_u
-
-        if i in (self.n, self.o):
-            sign = 1 if i == self.n else -1
-            dqdx += sign * np.cross(w, v)[k] / lambda_v
-
-        return dqdx
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
+    ):
+        """m - o - n angle"""
+        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 DifferentiableMath.acos(u.dot(v) / (u.norm() * v.norm()))
 
     def __repr__(self):
         return f"Angle({self.m}-{self.o}-{self.n})"
 
 
 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
 
@@ -344,88 +371,174 @@ class PrimitiveDihedralAngle(Primitive):
         """Dihedral angle: m-o-p-n"""
         super().__init__(m, o, p, n)
 
-        self.m = m
-        self.o = o
-        self.p = p
-        self.n = n
-
-    def __call__(self, x: "CartesianCoordinates") -> float:
-        """Value of the dihedral"""
-        return self._value(x, return_derivative=False)
-
-    def derivative(
-        self,
-        i: int,
-        component: "CartesianComponent",
-        x: "CartesianCoordinates",
-    ) -> float:
-        return self._value(x, i=i, component=component, return_derivative=True)
+        self.m = int(m)
+        self.o = int(o)
+        self.p = int(p)
+        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))
         )
 
-    def _value(self, x, i=None, component=None, return_derivative=False):
-        """Evaluate either the value or the derivative. Shared function
-        to reuse local variables"""
-
-        _x = x.reshape((-1, 3))
-        u = _x[self.m, :] - _x[self.o, :]
-        lambda_u = np.linalg.norm(u)
-        u /= lambda_u
-
-        v = _x[self.n, :] - _x[self.p, :]
-        lambda_v = np.linalg.norm(v)
-        v /= lambda_v
-
-        w = _x[self.p, :] - _x[self.o, :]
-        lambda_w = np.linalg.norm(w)
-        w /= lambda_w
-
-        phi_u = np.arccos(u.dot(w))
-        phi_v = np.arccos(w.dot(v))
-
-        if not return_derivative:
-            v1 = np.cross(u, w)
-            v2 = np.cross(-w, v)
-            return -np.arctan2(np.cross(v1, w).dot(v2), (v1.dot(v2)))
-
-        # are now computing the derivative..
-        if i not in self._atom_indexes:
-            return 0.0
-
-        k = int(component)
-        dqdx = 0.0
-
-        if i in (self.m, self.o):
-            sign = 1 if i == self.m else -1
-            dqdx += sign * (
-                np.cross(u, w)[k] / (lambda_u * np.sin(phi_u) ** 2)
-            )
-
-        if i in (self.p, self.n):
-            sign = 1 if i == self.p else -1
-            dqdx += sign * (
-                np.cross(v, w)[k] / (lambda_v * np.sin(phi_v) ** 2)
-            )
-
-        if i in (self.o, self.p):
-            sign = 1 if i == self.o else -1
-            dqdx += sign * (
-                (
-                    (np.cross(u, w)[k] * np.cos(phi_u))
-                    / (lambda_w * np.sin(phi_u) ** 2)
-                )
-                - (
-                    (np.cross(v, w)[k] * np.cos(phi_v))
-                    / (lambda_w * np.sin(phi_v) ** 2)
-                )
-            )
-
-        return dqdx
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
+    ) -> "VectorHyperDual":
+        """Dihedral m-o-p-n"""
+        # https://en.wikipedia.org/wiki/Dihedral_angle#In_polymer_physics
+        _x = x.ravel()
+        vec_m, vec_o, vec_p, vec_n = _get_3d_vecs_from_atom_idxs(
+            self.m, self.o, self.p, self.n, x=_x, deriv_order=deriv_order
+        )
+        u_1 = vec_o - vec_m
+        u_2 = vec_p - vec_o
+        u_3 = vec_n - vec_p
+
+        norm_u2 = u_2.norm()
+        v1 = u_2.cross(u_3)
+        v2 = u_1.cross(u_2)
+        v3 = u_1 * norm_u2
+        dihedral = DifferentiableMath.atan2(v3.dot(v1), v2.dot(v1))
+        assert isinstance(dihedral, VectorHyperDual)
+        return dihedral
 
     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 f6056aa..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,27 +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
-
-                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/path/path.py b/autode/path/path.py
index 81ab756..64939cc 100644
--- a/autode/path/path.py
+++ b/autode/path/path.py
@@ -1,5 +1,4 @@
 import numpy as np
-import matplotlib.pyplot as plt
 
 from autode.species import Species
 from autode.input_output import atoms_to_xyz_file
@@ -147,6 +146,8 @@ class Path(list):
         self, save: bool, name: str, color: str, xlabel: str
     ) -> None:
         """Plot this path"""
+        import matplotlib.pyplot as plt
+
         if len(self) == 0 or any(item.energy is None for item in self):
             logger.error("Could not plot a surface, an energy was None")
             return
diff --git a/autode/pes/pes_nd.py b/autode/pes/pes_nd.py
index edc5933..7741c59 100644
--- a/autode/pes/pes_nd.py
+++ b/autode/pes/pes_nd.py
@@ -5,7 +5,6 @@ surface and connecting minima and saddle points
 """
 import numpy as np
 import itertools as it
-import matplotlib.pyplot as plt
 
 from abc import ABC, abstractmethod
 from typing import (
@@ -196,6 +195,7 @@ class PESnD(ABC):
 
             units: Units of the surface. One of {'Ha', 'eV', 'kcal', 'kJ'}
         """
+        import matplotlib.pyplot as plt
 
         if interp_factor < 0:
             raise ValueError(
@@ -559,6 +559,8 @@ class PESnD(ABC):
             interp_factor:
             units_name:
         """
+        import matplotlib.pyplot as plt
+
         r_x = self._rs[0]
         energies, units = self._energies, energy_unit_from_name(units_name)
         energies = units.times * (energies - np.min(energies))
@@ -605,6 +607,7 @@ class PESnD(ABC):
             interp_factor:
             units_name:
         """
+        import matplotlib.pyplot as plt
         from mpl_toolkits.mplot3d import Axes3D
         from matplotlib.ticker import FormatStrFormatter
 
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/autode/wrappers/ORCA.py b/autode/wrappers/ORCA.py
index 9336192..1732de6 100644
--- a/autode/wrappers/ORCA.py
+++ b/autode/wrappers/ORCA.py
@@ -220,7 +220,7 @@ class ORCA(autode.wrappers.methods.ExternalMethodOEGH):
     def execute(self, calc):
         @work_in_tmp_dir(
             filenames_to_copy=calc.input.filenames,
-            kept_file_exts=(".out", ".hess", ".xyz", ".inp", ".pc"),
+            kept_file_exts=Config.ORCA.copied_output_exts,
         )
         def execute_orca():
             run_external(
diff --git a/doc/changelog.rst b/doc/changelog.rst
index 7d4c49f..924a896 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -1,6 +1,24 @@
 Changelog
 =========
 
+1.4.2
+------
+-------
+
+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
+*********
+- Fixes triangular rings being incorrectly treated as dihedral angles
+
+Usability improvements/Changes
+******************************
+- Faster import of autode package by lazily loading matplotlib
+- ORCA output files copied after a calculation are configurable
 
 1.4.1
 ------
diff --git a/setup.py b/setup.py
index 0e243e7..300adb9 100644
--- a/setup.py
+++ b/setup.py
@@ -32,7 +32,7 @@ extensions = [
 
 setup(
     name="autode",
-    version="1.4.1",
+    version="1.4.2",
     python_requires=">3.7",
     packages=[
         "autode",


Test patch:
diff --git a/tests/test_import.py b/tests/test_import.py
new file mode 100644
index 0000000..db90f0f
--- /dev/null
+++ b/tests/test_import.py
@@ -0,0 +1,32 @@
+"""Tests for the import speed of autode."""
+import sys
+
+import pytest
+
+SLOW_IMPORTS = ["matplotlib"]
+
+
+@pytest.fixture
+def unimport_slow_imports():
+    """Remove modules in ``SLOW_IMPORTS`` from ``sys.modules``."""
+    for module in SLOW_IMPORTS:
+        if module in sys.modules:
+            del sys.modules[module]
+
+
+@pytest.mark.usefixtures("unimport_slow_imports")
+def test_slow_imports_during_tab_completion():
+    """Check that importing autode does not import certain python modules that would make import slow."""
+
+    # Let's double check that the undesired imports are not already loaded
+    for modulename in SLOW_IMPORTS:
+        assert (
+            modulename not in sys.modules
+        ), f"Module `{modulename}` was not properly unloaded"
+
+    import autode
+
+    for modulename in SLOW_IMPORTS:
+        assert (
+            modulename not in sys.modules
+        ), f"Detected loaded module {modulename} after autode import"
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_autodiff.py b/tests/test_opt/test_autodiff.py
new file mode 100644
index 0000000..eec7065
--- /dev/null
+++ b/tests/test_opt/test_autodiff.py
@@ -0,0 +1,190 @@
+import math
+import numpy as np
+import pytest
+from autode.opt.coordinates._autodiff import (
+    get_differentiable_vars,
+    DifferentiableMath,
+    VectorHyperDual,
+    DerivativeOrder,
+    DifferentiableVector3D,
+)
+
+
+def test_hyperdual_sanity_checks():
+    a, b = get_differentiable_vars(
+        values=[1, 2], symbols=["a", "b"], deriv_order=DerivativeOrder.first
+    )
+    assert repr(a) is not None
+
+    c = get_differentiable_vars(
+        values=[1], symbols=["c"], deriv_order=DerivativeOrder.first
+    )[0]
+    with pytest.raises(ValueError, match="Incompatible number"):
+        _ = a + c
+
+    c, d = get_differentiable_vars(
+        values=[1, 2], symbols=["c", "d"], deriv_order=DerivativeOrder.first
+    )
+    with pytest.raises(ValueError, match="symbols do not match"):
+        _ = a * c
+
+    a2, b2 = get_differentiable_vars(
+        values=[1, 2], symbols=["a", "b"], deriv_order=DerivativeOrder.second
+    )
+    with pytest.raises(ValueError, match="order of derivative do not match!"):
+        _ = DifferentiableMath.atan2(a2, b)
+
+    with pytest.raises(RuntimeError, match="Symbols must be unique"):
+        _ = get_differentiable_vars(
+            values=[1, 2],
+            symbols=["a", "a"],
+            deriv_order=DerivativeOrder.first,
+        )
+
+
+def x_to_the_y_derivs(x, y):
+    d_dx = y * x ** (y - 1)
+    d_dy = x**y * math.log(x)
+    return d_dx, d_dy
+
+
+def x_to_the_y_second_derivs(x, y):
+    d2_dx2 = y * (y - 1) * x ** (y - 2)
+    d2_dxdy = math.log(x) * y * x ** (y - 1) + x**y * (1 / x)
+    d2_dydx = x ** (y - 1) + y * x ** (y - 1) * math.log(x)
+    d2_dy2 = math.log(x) * math.log(x) * x**y
+    return d2_dx2, d2_dxdy, d2_dydx, d2_dy2
+
+
+def test_autodiff_exponential_func():
+    x_val, y_val = 1.0, 2.0
+    x, y = get_differentiable_vars(
+        values=[x_val, y_val],
+        symbols=["x", "y"],
+    )
+    result = x**y
+    d_dx, d_dy = x_to_the_y_derivs(x_val, y_val)
+    assert math.isclose(result.differentiate_wrt("x"), d_dx)
+    assert math.isclose(result.differentiate_wrt("y"), d_dy)
+
+    d2_dx2, d2_dxdy, d2_dydx, d2_dy2 = x_to_the_y_second_derivs(x_val, y_val)
+    assert math.isclose(result.differentiate_wrt("x", "x"), d2_dx2)
+    assert math.isclose(result.differentiate_wrt("x", "y"), d2_dxdy)
+    assert math.isclose(result.differentiate_wrt("y", "x"), d2_dydx)
+    assert math.isclose(result.differentiate_wrt("y", "y"), d2_dy2)
+    # check the second derivatives are symmetric
+    assert math.isclose(d2_dydx, d2_dxdy)
+
+
+def test_exponential_math_sanity_checks():
+    x, y = get_differentiable_vars([-0.1, -0.2], ["x", "y"])
+    assert x**2 is not None
+    assert x**-1 is not None
+    # negative number raised to fractional power is complex
+    with pytest.raises(AssertionError):
+        _ = x**0.2
+    # negative number cannot be raised to differentiable power
+    with pytest.raises(AssertionError):
+        _ = (-2) ** y
+    with pytest.raises(AssertionError):
+        _ = x**y
+
+
+def test_math_funcs_work_with_native_types():
+    # python float or int types should be passed through
+    assert math.isclose(DifferentiableMath.acos(0), math.pi / 2)
+    (y,) = get_differentiable_vars([1], ["y"])
+    res = DifferentiableMath.atan2(y, 0)
+    assert isinstance(res, VectorHyperDual)
+    assert math.isclose(res.value, math.pi / 2)
+    res = DifferentiableMath.atan2(1, 0)
+    assert isinstance(res, float)
+    assert math.isclose(DifferentiableMath.pow(0.9, 1.3), math.pow(0.9, 1.3))
+
+    # however, python's complex type is not supported
+    with pytest.raises(TypeError, match="Unknown type for addition"):
+        _ = y + (4 + 1j)
+    with pytest.raises(TypeError, match="Unknown type for multiplication"):
+        _ = y * (1 + 2j)
+    with pytest.raises(TypeError, match="Unknown type for exponentiation"):
+        _ = y ** (1 + 2j)
+
+
+def test_hyperdual_init_checks():
+    with pytest.raises(ValueError):
+        _ = VectorHyperDual(
+            0.1, symbols=["x", "y"], first_der=np.array([0.1, 0.2, 0.3])
+        )
+
+    with pytest.raises(ValueError):
+        _ = VectorHyperDual(
+            0.1,
+            symbols=["x", "y"],
+            first_der=np.array([0.1, 0.2]),
+            second_der=np.zeros(shape=(2, 3)),
+        )
+
+
+def test_hyperdual_order():
+    symbols = ["x", "y"]
+    x = VectorHyperDual(
+        0.1,
+        symbols=symbols,
+    )
+    assert x._order == DerivativeOrder.zeroth
+    x = VectorHyperDual(0.1, symbols=symbols, first_der=np.zeros(2))
+    assert x._order == DerivativeOrder.first
+    # will ignore second derivatives if first is not present
+    x = VectorHyperDual(0.1, symbols=symbols, second_der=np.zeros((3, 2)))
+    assert x._order == DerivativeOrder.zeroth
+    x = VectorHyperDual(
+        0.1,
+        symbols=symbols,
+        first_der=np.zeros(2),
+        second_der=np.zeros((2, 2)),
+    )
+    assert x._order == DerivativeOrder.second
+
+
+def test_derivative_not_available():
+    x, y = get_differentiable_vars(
+        [1.0, 2.0], symbols=["x", "y"], deriv_order=DerivativeOrder.first
+    )
+    res = 1 - x**2 + y
+    assert isinstance(res, VectorHyperDual)
+    assert math.isclose(res.value, 2)
+    assert math.isclose(res.differentiate_wrt("x"), -2)
+    # higher order derivatives are not available
+    assert res.differentiate_wrt("x", "y") is None
+    # unknown variables
+    assert res.differentiate_wrt("z") is None
+    assert res.differentiate_wrt("y", "z") is None
+
+    # with zero order, only value is present
+    (x,) = get_differentiable_vars(
+        [1.0], symbols=["x"], deriv_order=DerivativeOrder.zeroth
+    )
+    res = 1 + x**2
+    assert isinstance(res, VectorHyperDual)
+    assert res.differentiate_wrt("x") is None
+
+
+def test_hyperdual_3d_vector_sanity_checks():
+    x, y, z = get_differentiable_vars(
+        [0.1, 0.2, 0.3],
+        ["x", "y", "z"],
+    )
+
+    with pytest.raises(ValueError):
+        _ = DifferentiableVector3D([x, y])
+
+    vec1 = DifferentiableVector3D([x, y, z])
+    vec2 = DifferentiableVector3D([x, y, z])
+
+    with pytest.raises(ValueError):
+        _ = vec1.dot(2)
+
+    with pytest.raises(AssertionError):
+        _ = vec1 * vec2
+
+    assert 2 * vec1 is not None
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 2f775e8..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,
 )
 
 
@@ -28,18 +45,12 @@ def test_inv_dist_primitives():
     assert np.isclose(inv_dist(x), 0.5)  # 1/2.0 = 0.5 Å-1
 
     # Check a couple of derivatives by hand
-    assert np.isclose(
-        inv_dist.derivative(0, CartesianComponent.x, x=x), 2 * inv_dist(x) ** 3
-    )
-    assert np.isclose(
-        inv_dist.derivative(1, CartesianComponent.x, x=x),
-        -2 * inv_dist(x) ** 3,
-    )
+    derivs = inv_dist.derivative(x=x)
+    assert np.isclose(derivs[3 * 0 + 0], 2 * inv_dist(x) ** 3)
+    assert np.isclose(derivs[3 * 1 + 0], -2 * inv_dist(x) ** 3)
 
     # Derivatives with respect to zero components
-    assert np.isclose(inv_dist.derivative(0, CartesianComponent.y, x=x), 0)
-    # or those that are not present in the system should be zero
-    assert np.isclose(inv_dist.derivative(2, CartesianComponent.x, x=x), 0)
+    assert np.isclose(derivs[3 * 0 + 1], 0)
 
 
 def test_dist_primitives():
@@ -50,18 +61,12 @@ def test_dist_primitives():
     inv_dist = PrimitiveDistance(0, 1)
     assert np.isclose(inv_dist(x), 2.0)
 
-    assert np.isclose(
-        inv_dist.derivative(0, CartesianComponent.x, x=x), -2 / 2
-    )
+    derivs = inv_dist.derivative(x)
+    assert np.isclose(derivs[3 * 0 + 0], -2 / 2)
+    assert np.isclose(derivs[3 * 1 + 0], +2 / 2)
 
-    assert np.isclose(
-        inv_dist.derivative(1, CartesianComponent.x, x=x), +2 / 2
-    )
-
-    for component in (CartesianComponent.y, CartesianComponent.z):
-        assert np.isclose(inv_dist.derivative(1, component, x=x), 0)
-
-    assert np.isclose(inv_dist.derivative(2, CartesianComponent.x, x=x), 0)
+    for k in (1, 2):
+        assert np.isclose(derivs[3 * 1 + k], 0)
 
 
 def test_primitive_equality():
@@ -539,6 +544,23 @@ def test_pic_b_no_primitives():
         c._calc_B(np.arange(6, dtype=float).reshape(2, 3))
 
 
+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.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():
     d = ConstrainedPrimitiveDistance(0, 1, value=1.0)
 
@@ -562,25 +584,20 @@ 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 CartesianComponent:
-            analytic = angle.derivative(atom_idx, component, init_coords)
+        for component in (0, 1, 2):
+            analytic = derivs[3 * atom_idx + component]
 
             assert np.isclose(
                 analytic, numerical_derivative(atom_idx, component), atol=1e-6
             )
 
-    # Derivative should be zero for an atom not present the bond angle
-    assert np.isclose(
-        angle.derivative(3, CartesianComponent.x, init_coords), 0.0
-    )
-
 
 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():
@@ -605,17 +622,11 @@ def test_dihedral_primitive_derivative():
     init_coords = m.coordinates.copy()
 
     dihedral = PrimitiveDihedralAngle(2, 0, 1, 3)
-
+    analytic = dihedral.derivative(init_coords)
     for atom_idx in (0, 1, 2, 3):
-        for component in CartesianComponent:
-            analytic = dihedral.derivative(atom_idx, component, init_coords)
-            numerical = numerical_derivative(atom_idx, component)
-
-            assert np.isclose(analytic, numerical, atol=1e-6)
-
-    assert np.isclose(
-        dihedral.derivative(5, CartesianComponent.x, init_coords), 0.0
-    )
+        for k in [0, 1, 2]:
+            numerical = numerical_derivative(atom_idx, k)
+            assert np.isclose(analytic[3 * atom_idx + k], numerical, atol=1e-6)
 
 
 def test_dihedral_equality():
@@ -627,40 +638,136 @@ def test_dihedral_equality():
     )
 
 
-@pytest.mark.parametrize(
-    "h_coord",
-    [
-        np.array([1.08517, 1.07993, 0.05600]),
-        np.array([1.28230, -0.63391, -0.54779]),
-    ],
-)
-def test_dihedral_primitive_derivative_over_zero(h_coord):
-    def numerical_derivative(a, b, h=1e-6):
-        coords = init_coords.copy()
-        coords[a, int(b)] += h
-        y_plus = dihedral(coords)
-        coords[a, int(b)] -= 2 * h
-        y_minus = dihedral(coords)
-        return (y_plus - y_minus) / (2 * h)
+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()
+    coords = h2o2.coordinates
+    dist = PrimitiveDistance(0, 1)
+    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(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)
+
 
-    m = Molecule(
+# fmt: off
+extra_mols = [
+    Molecule(
         atoms=[
             Atom("C", 0.63365, 0.11934, -0.13163),
             Atom("C", -0.63367, -0.11938, 0.13153),
-            Atom("H", 0.0, 0.0, 0.0),
+            Atom("H", 1.08517, 1.07993, 0.05600),
             Atom("H", -1.08517, -1.07984, -0.05599),
         ]
-    )
-    m.atoms[2].coord = h_coord
-    init_coords = m.coordinates
-
-    dihedral = PrimitiveDihedralAngle(2, 0, 1, 3)
-
-    for atom_idx in (0, 1, 2, 3):
-        for component in CartesianComponent:
-            analytic = dihedral.derivative(atom_idx, component, init_coords)
-            numerical = numerical_derivative(atom_idx, component)
-            assert np.isclose(analytic, numerical, atol=1e-6)
+    ),
+    Molecule(
+        atoms=[
+            Atom("C", 0.63365, 0.11934, -0.13163),
+            Atom("C", -0.63367, -0.11938, 0.13153),
+            Atom("H", 1.28230, -0.63391, -0.54779),
+            Atom("H", -1.08517, -1.07984, -0.05599),
+        ]
+    ),  # 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(), *extra_mols
+]
+test_prims = [
+    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),
+    PrimitiveLinearAngle(2, 3, 4, 8, LinearBendType.BEND),
+]
+# fmt: on
+
+
+@pytest.mark.parametrize("mol,prim", list(zip(test_mols, test_prims)))
+def test_primitive_first_derivs(mol, prim):
+    init_coords = CartesianCoordinates(mol.coordinates)
+    init_prim = prim(init_coords)
+
+    def numerical_first_deriv(coords, h=1e-8):
+        coords = coords.flatten()
+        derivs = np.zeros_like(coords)
+        for i in range(coords.shape[0]):
+            coords[i] += h
+            derivs[i] = (prim(coords) - init_prim) / h
+            coords[i] -= h
+        return derivs
+
+    analytic = prim.derivative(init_coords)
+    numeric = numerical_first_deriv(init_coords)
+    assert np.allclose(analytic, numeric, atol=1e-6)
+
+
+@pytest.mark.parametrize("mol,prim", list(zip(test_mols, test_prims)))
+def test_primitve_second_deriv(mol, prim):
+    init_coords = CartesianCoordinates(mol.coordinates)
+    init_first_der = prim.derivative(init_coords)
+
+    def numerical_second_deriv(coords, h=1e-8):
+        coords = coords.flatten()
+        derivs = np.zeros((coords.shape[0], coords.shape[0]))
+        for i in range(coords.shape[0]):
+            coords[i] += h
+            derivs[i] = (prim.derivative(coords) - init_first_der) / h
+            coords[i] -= h
+        return derivs
+
+    analytic = prim.second_derivative(init_coords)
+    # second derivative matrix should be symmetric
+    assert np.allclose(analytic, analytic.T)
+    numeric = numerical_second_deriv(init_coords)
+    assert np.allclose(analytic, numeric, atol=1e-6)
 
 
 def test_repr():
@@ -670,23 +777,18 @@ 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:
         assert repr(p) is not None
 
 
-@pytest.mark.parametrize("sign", [1, -1])
-def test_angle_normal(sign):
-    angle = PrimitiveBondAngle(0, 1, 2)
-    x = np.array([[0.0, 0.0, 0.0], [1.0, sign * 1.0, 1.0], [1.0, 0.0, 1.0]])
-
-    assert not np.isinf(angle.derivative(0, 1, x))
-
-
 def test_dic_large_step_allowed_unconverged_back_transform():
     x = CartesianCoordinates(water_mol().coordinates)
     dic = DIC.from_cartesian(x)
@@ -705,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)
@@ -714,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
 
@@ -729,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")