diff --git a/autode/opt/coordinates/_autodiff.py b/autode/opt/coordinates/_autodiff.py
new file mode 100644
index 000000000..0fa48d437
--- /dev/null
+++ b/autode/opt/coordinates/_autodiff.py
@@ -0,0 +1,714 @@
+"""
+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["VectorHyperDual"]):
+        """
+        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) 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") -> "VectorHyperDual":
+        """
+        Dot product of two 3D vectors
+
+        Args:
+            other (DifferentiableVector3D):
+
+        Returns:
+            (VectorHyperDual): A scalar number (with derivatives)
+        """
+        self._check_same_type(other)
+        dot = 0
+        for k in range(3):
+            dot = dot + self._data[k] * other._data[k]
+        assert isinstance(dot, VectorHyperDual)
+        return dot
+
+    def norm(self) -> "VectorHyperDual":
+        """
+        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
+        )
+        assert isinstance(norm, VectorHyperDual)
+        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 cross(
+        self, other: "DifferentiableVector3D"
+    ) -> "DifferentiableVector3D":
+        """
+        Cross-product of two 3D vectors, produces another vector
+
+        Args:
+            other (DifferentiableVector3D):
+
+        Returns:
+            (DifferentiableVector3D):
+        """
+        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/dic.py b/autode/opt/coordinates/dic.py
index b857f2adb..2d5277854 100644
--- a/autode/opt/coordinates/dic.py
+++ b/autode/opt/coordinates/dic.py
@@ -216,34 +216,42 @@ def iadd(self, value: np.ndarray) -> "OptCoordinates":
         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 675ca1f1c..9e7462e16 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -82,6 +82,11 @@ def __init__(self, *args: Any):
                 f"from {args}. Must be primitive internals"
             )
 
+    def append(self, item: Primitive) -> None:
+        """Append an item to this set of primitives"""
+        assert isinstance(item, Primitive), "Must be a Primitive type!"
+        super().append(item)
+
     @property
     def B(self) -> np.ndarray:
         """Wilson B matrix"""
@@ -173,22 +178,12 @@ def _calc_B(self, x: np.ndarray) -> None:
                 "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
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 22bfabb86..764c322b0 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -1,12 +1,57 @@
 import numpy as np
 
 from abc import ABC, abstractmethod
-from typing import Tuple, TYPE_CHECKING
+from typing import Tuple, TYPE_CHECKING, List
+from autode.opt.coordinates._autodiff import (
+    get_differentiable_vars,
+    DifferentiableMath,
+    DifferentiableVector3D,
+    DerivativeOrder,
+    VectorHyperDual,
+)
 
 if TYPE_CHECKING:
     from autode.opt.coordinates import CartesianCoordinates, CartesianComponent
 
 
+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):
     """Primitive internal coordinate"""
 
@@ -17,18 +62,35 @@ def __init__(self, *atom_indexes: int):
         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 +102,56 @@ def derivative(
 
         -----------------------------------------------------------------------
         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 +222,8 @@ def __eq__(self, other) -> bool:
 
 
 class PrimitiveInverseDistance(_DistanceFunction):
-    """
-    Inverse distance between to atoms:
+    r"""
+    Inverse distance between two atoms:
 
     .. math::
 
@@ -130,40 +231,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()
+
+    def __repr__(self):
+        return f"InverseDistance({self.i}-{self.j})"
 
 
 class PrimitiveDistance(_DistanceFunction):
-    """
+    r"""
     Distance between two atoms:
 
     .. math::
@@ -171,33 +253,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()
 
     def __repr__(self):
         return f"Distance({self.i}-{self.j})"
@@ -230,13 +293,18 @@ def __repr__(self):
 
 
 class PrimitiveBondAngle(Primitive):
+    """
+    Bond angle between three atoms, calculated with the
+    arccosine of the normalised dot product
+    """
+
     def __init__(self, o: int, m: int, n: int):
         """Bond angle m-o-n"""
         super().__init__(o, m, n)
 
-        self.o = o
-        self.m = m
-        self.n = n
+        self.o = int(o)
+        self.m = int(m)
+        self.n = int(n)
 
     def __eq__(self, other) -> bool:
         """Equality of two distance functions"""
@@ -247,58 +315,16 @@ def __eq__(self, other) -> bool:
             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})"
@@ -344,22 +370,10 @@ def __init__(self, m: int, o: int, p: int, n: int):
         """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"""
@@ -368,64 +382,26 @@ def __eq__(self, other) -> bool:
             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})"
diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py
index f6056aa68..8dceaccf5 100644
--- a/autode/opt/optimisers/crfo.py
+++ b/autode/opt/optimisers/crfo.py
@@ -238,6 +238,10 @@ def _dihedrals(species):
                 if n == o:
                     continue
 
+                # avoid triangular rings like cyclopropane
+                if m == n:
+                    continue
+
                 if np.isclose(species.angle(o, p, n), Angle(np.pi), atol=0.04):
                     continue
 
diff --git a/doc/changelog.rst b/doc/changelog.rst
index 460ed3494..1d2b97732 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -7,11 +7,11 @@ Changelog
 
 Functionality improvements
 **************************
-- ...
+- Replaces hard-coded derivatives for primitive internal coordinates with automatic differentiation
 
 Bug Fixes
 *********
-- ...
+- Fixes triangular rings being incorrectly treated as dihedral angles
 
 
 1.4.1
diff --git a/tests/test_opt/test_autodiff.py b/tests/test_opt/test_autodiff.py
new file mode 100644
index 000000000..eec7065d1
--- /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 2f775e81b..9fe31f72e 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -28,18 +28,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 +44,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 +527,14 @@ def test_pic_b_no_primitives():
         c._calc_B(np.arange(6, dtype=float).reshape(2, 3))
 
 
+def test_pic_append_type_checking():
+    c = PIC()
+    # append should check for primitive type
+    c.append(PrimitiveDistance(0, 1))
+    with pytest.raises(AssertionError):
+        c.append(3)
+
+
 def test_constrained_distance_satisfied():
     d = ConstrainedPrimitiveDistance(0, 1, value=1.0)
 
@@ -563,20 +559,15 @@ def numerical_derivative(a, b, h=1e-8):
     init_coords = m.coordinates.copy()
 
     angle = PrimitiveBondAngle(0, 1, 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)
@@ -605,17 +596,11 @@ def numerical_derivative(a, b, h=1e-8):
     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 +612,91 @@ 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_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(2, 0, 1)  # bond is defined in a different way
+    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
+dihedral_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
+]
+
+test_mols = [
+    h2o2_mol(), h2o2_mol(), water_mol(),
+    water_mol(), water_mol(), *dihedral_mols
+]
+test_prims = [
+    PrimitiveDihedralAngle(2, 0, 1, 3), PrimitiveBondAngle(2, 0, 1),
+    PrimitiveBondAngle(0, 1, 2), PrimitiveDistance(0, 1),
+    PrimitiveInverseDistance(0, 1), PrimitiveDihedralAngle(2, 0, 1, 3),
+    PrimitiveDihedralAngle(2, 0, 1, 3)
+]
+# 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():
@@ -679,14 +715,6 @@ def test_repr():
         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)