From a25db2e37947aca7d03cbc752cfa3a63ea7e05d8 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 16 Nov 2023 11:41:05 +0000
Subject: [PATCH 01/37] initial commit autodiff

---
 autode/opt/coordinates/autodiff.py | 523 +++++++++++++++++++++++++++++
 1 file changed, 523 insertions(+)
 create mode 100644 autode/opt/coordinates/autodiff.py

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
new file mode 100644
index 000000000..1da2f2787
--- /dev/null
+++ b/autode/opt/coordinates/autodiff.py
@@ -0,0 +1,523 @@
+"""
+Automatic differentiation routines in pure Python
+"""
+from typing import Union, Callable, Sequence, Optional
+from copy import deepcopy
+import numpy as np
+import math
+
+numeric = (float, int)
+numeric_type = Union[float, int]
+
+
+def get_differentiable_vars(*args, deriv_order: int = 2, symbols=None):
+    """
+    Obtain differentiable variables from a series of numbers
+
+    Args:
+        *args: The values of the variables (numbers)
+        deriv_order: Order of differentiation
+        symbols: Optional list of symbols for the numbers,
+                 if not given, numeric strings ('0', '1', '2', ...)
+                 are generated in order
+
+    Returns:
+        (list[VectorHyperDual]): A list of hyper dual numbers
+    """
+    if symbols is None:
+        symbols = [str(idx) for idx in range(len(args))]
+    assert all(isinstance(sym, str) for sym in symbols)
+    assert len(symbols) == len(args)
+    symbols = list(symbols)
+
+    hyperduals = []
+    for symbol, value in zip(symbols, args):
+        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 = 0
+        if first_der is None:
+            return
+        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 = 1
+
+        if second_der is None:
+            return
+        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 = 2
+
+    @property
+    def n_vars(self) -> int:
+        """Number of variables in this hyper-dual"""
+        return len(self._symbols)
+
+    def copy(self):
+        return deepcopy(self)
+
+    def _check_compatible(self, other: "VectorHyperDual"):
+        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):
+        return self._val
+
+    @value.setter
+    def value(self, value):
+        assert isinstance(value, numeric)
+        self._val = float(value)
+
+    @classmethod
+    def from_variable(
+        cls, value: float, symbol: str, all_symbols: Sequence[str], order: int
+    ):
+        """
+        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)
+
+        if order >= 1:
+            first_der = np.zeros(shape=n, dtype=float)
+            first_der[idx] = 1.0
+        if order == 2:
+            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 == 0:
+            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 == 1:
+            return None
+        assert self._second_der is not None
+        return self._second_der[idx_1, idx_2]
+
+    def __add__(self, other):
+        """Adding a hyper dual number"""
+
+        # add to a float
+        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 == 0:
+                return VectorHyperDual(val, self._symbols)
+
+            first_der = self._first_der + other._first_der
+            if self._order == 1:
+                return VectorHyperDual(val, self._symbols, first_der)
+
+            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):
+        """Unary negative operation"""
+        new = self.copy()
+        new._val = -new._val
+        if self._order >= 1:
+            new._first_der = -new._first_der
+        if self._order == 2:
+            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):
+        """Multiply a hyper dual number with float or another hyper dual"""
+        if isinstance(other, numeric):
+            new = self.copy()
+            new._val *= float(other)
+            return new
+
+        # multiply with another dual
+        elif isinstance(other, VectorHyperDual):
+            self._check_compatible(other)
+
+            val = self._val * other._val
+            if self._order == 0:
+                return VectorHyperDual(val, self._symbols)
+
+            first_der = (
+                self._val * other._first_der + other._val * self._first_der
+            )
+            if self._order == 1:
+                return VectorHyperDual(val, self._symbols, first_der)
+
+            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):
+        if modulo is not None:
+            raise NotImplementedError("Modulo inverse is not implemented")
+
+        return DifferentiableMath.pow(self, power)
+
+    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],
+    ):
+        """
+        Perform an operation on the hyperdual (i.e. apply a scalar function)
+
+        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))
+
+        val = operator(num._val)
+
+        if num._order == 0:
+            return VectorHyperDual(val, num._symbols)
+
+        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 == 1:
+            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]):
+        """Calculate the square root of a hyperdual number"""
+
+        if isinstance(num, numeric):
+            assert num > 0
+        else:
+            assert num.value > 0
+
+        def value(x0: float):
+            return math.sqrt(x0)
+
+        def derivative(x0: float):
+            return 1 / (2 * math.sqrt(x0))
+
+        def second_derivative(x0: float):
+            return -1 / (4 * math.pow(x0, 3 / 2))
+
+        return VectorHyperDual.apply_operation(
+            num, value, derivative, second_derivative
+        )
+
+    @staticmethod
+    def exp(num: Union[VectorHyperDual, numeric_type]):
+        """Raise e to the power of num"""
+
+        def value(x0: float):
+            return math.exp(x0)
+
+        def derivative(x0: float):
+            return math.exp(x0)
+
+        def second_derivative(x0: float):
+            return math.exp(x0)
+
+        return VectorHyperDual.apply_operation(
+            num, value, derivative, second_derivative
+        )
+
+    @staticmethod
+    def pow(
+        num: Union[VectorHyperDual, numeric_type],
+        power: Union[VectorHyperDual, numeric_type],
+    ):
+        """Exponentiation"""
+
+        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 ValueError(
+                    "Math error, can't raise negative number to fractional power"
+                )
+
+            def value(x0: float):
+                return math.pow(x0, power)
+
+            def derivative(x0: float):
+                return power * math.pow(x0, power - 1)
+
+            def second_derivative(x0: float):
+                return power * (power - 1) * math.pow(x0, power - 2)
+
+            return VectorHyperDual.apply_operation(
+                num, value, derivative, second_derivative
+            )
+
+        elif isinstance(power, VectorHyperDual):
+            if (isinstance(num, numeric) and num) < 0 or (
+                isinstance(num, VectorHyperDual) and num.value < 0
+            ):
+                raise ValueError(
+                    "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))
+
+    @staticmethod
+    def log(num: Union[VectorHyperDual, numeric_type]):
+        """Natural logarithm"""
+
+        if isinstance(num, numeric):
+            assert num > 0
+        else:
+            assert num.value > 0
+
+        def value(x0: float):
+            return math.log(x0)
+
+        def derivative(x0: float):
+            return 1.0 / x0
+
+        def second_derivative(x0: float):
+            return -1.0 / (x0**2)
+
+        return VectorHyperDual.apply_operation(
+            num, value, derivative, second_derivative
+        )
+
+    @staticmethod
+    def acos(num: 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
+
+        def value(x0: float):
+            return math.acos(x0)
+
+        def derivative(x0: float):
+            return -1 / math.sqrt(1 - x0**2)
+
+        def second_derivative(x0: float):
+            return -x0 / math.pow(1 - x0**2, 3 / 2)
+
+        return VectorHyperDual.apply_operation(
+            num, value, derivative, second_derivative
+        )
+
+    @staticmethod
+    def atan(num: Union[VectorHyperDual, numeric_type]):
+        """Calculate the arctangent of a hyperdual number"""
+
+        def value(x0: float):
+            return math.atan(x0)
+
+        def derivative(x0: float):
+            return 1 / (1 + x0**2)
+
+        def second_derivative(x0: float):
+            return (-2 * x0) / (x0**2 + 1) ** 2
+
+        return VectorHyperDual.apply_operation(
+            num, value, derivative, second_derivative
+        )
+
+    @staticmethod
+    def atan2(
+        num_y: Union[VectorHyperDual, numeric_type],
+        num_x: Union[VectorHyperDual, numeric_type],
+    ):
+        """Calculate the arctan2 of two hyper dual numbers"""
+
+        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

From 1bbb7e489268bf4424330f55f718785513805ccc Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 16 Nov 2023 11:56:29 +0000
Subject: [PATCH 02/37] second commit autodiff

---
 autode/opt/coordinates/autodiff.py   | 14 +++++---
 autode/opt/coordinates/primitives.py | 50 +++++++++++++++++++++++++++-
 2 files changed, 58 insertions(+), 6 deletions(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index 1da2f2787..ab4a2348b 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -10,12 +10,16 @@
 numeric_type = Union[float, int]
 
 
-def get_differentiable_vars(*args, deriv_order: int = 2, symbols=None):
+def get_differentiable_vars(
+    values: Sequence[numeric_type],
+    deriv_order: int = 2,
+    symbols: Optional[Sequence[str]] = None,
+):
     """
     Obtain differentiable variables from a series of numbers
 
     Args:
-        *args: The values of the variables (numbers)
+        values: The values of the variables (numbers)
         deriv_order: Order of differentiation
         symbols: Optional list of symbols for the numbers,
                  if not given, numeric strings ('0', '1', '2', ...)
@@ -25,13 +29,13 @@ def get_differentiable_vars(*args, deriv_order: int = 2, symbols=None):
         (list[VectorHyperDual]): A list of hyper dual numbers
     """
     if symbols is None:
-        symbols = [str(idx) for idx in range(len(args))]
+        symbols = [str(idx) for idx in range(len(values))]
     assert all(isinstance(sym, str) for sym in symbols)
-    assert len(symbols) == len(args)
+    assert len(symbols) == len(values)
     symbols = list(symbols)
 
     hyperduals = []
-    for symbol, value in zip(symbols, args):
+    for symbol, value in zip(symbols, values):
         var = VectorHyperDual.from_variable(
             value, symbol, all_symbols=symbols, order=deriv_order
         )
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 22bfabb86..7cda7f430 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -1,10 +1,27 @@
 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,
+)
 
 if TYPE_CHECKING:
     from autode.opt.coordinates import CartesianCoordinates, CartesianComponent
+    from autode.opt.coordinates.autodiff import VectorHyperDual
+
+
+def _norm_vec3(vec: List["VectorHyperDual"]):
+    """Norm of a 3D vector"""
+    assert len(vec) == 3
+    return DifferentiableMath.sqrt(vec[0] ** 2 + vec[1] ** 2 + vec[2] ** 2)
+
+
+def _sub_vec3(vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]):
+    """Evaluate vec1 - vec2 for 3D vectors"""
+    assert len(vec1) == len(vec2) == 3
+    return [vec1[i] - vec2[i] for i in range(3)]
 
 
 class Primitive(ABC):
@@ -16,6 +33,21 @@ def __init__(self, *atom_indexes: int):
         """A primitive internal coordinate that involves a number of atoms"""
         self._atom_indexes = atom_indexes
 
+    def _evaluate(self, x: "CartesianCoordinates", deriv_order: int):
+        """
+        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
+        """
+
     @abstractmethod
     def __call__(self, x: "CartesianCoordinates") -> float:
         """Return the value of this PIC given a set of cartesian coordinates"""
@@ -130,6 +162,22 @@ class PrimitiveInverseDistance(_DistanceFunction):
                 {|\boldsymbol{X}_i - \boldsymbol{X}_j|}
     """
 
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: int
+    ) -> "VectorHyperDual":
+        """Inverse distance"""
+        _x = x.ravel()
+        idxs = []
+        for atom in [self.i, self.j]:
+            for k in range(3):
+                idxs.append(3 * atom + k)
+        i_0, i_1, i_2, j_0, j_1, j_2 = get_differentiable_vars(
+            values=[_x[idx] for idx in idxs],
+            deriv_order=deriv_order,
+            symbols=[str(idx) for idx in idxs],
+        )
+        return 1.0 / _norm_vec3(_sub_vec3([i_0, i_1, i_2], [j_0, j_1, j_2]))
+
     def derivative(
         self,
         i: int,

From 00966afc5ca5c632d6b5a299e10cf02c54c53b39 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 16 Nov 2023 18:26:45 +0000
Subject: [PATCH 03/37] third commit autodiff

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

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 7cda7f430..41c8a8a8f 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -1,7 +1,7 @@
 import numpy as np
 
 from abc import ABC, abstractmethod
-from typing import Tuple, TYPE_CHECKING, List
+from typing import Tuple, TYPE_CHECKING, List, Optional
 from autode.opt.coordinates.autodiff import (
     get_differentiable_vars,
     DifferentiableMath,
@@ -24,6 +24,44 @@ def _sub_vec3(vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]):
     return [vec1[i] - vec2[i] for i in range(3)]
 
 
+def _dot_vec3(vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]):
+    """Evaluate vec1.vec2 for 3D vectors"""
+    assert len(vec1) == len(vec2) == 3
+    return vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2]
+
+
+def get_vars_from_atom_idxs(
+    *args,
+    x: "CartesianCoordinates",
+    deriv_order: int,
+) -> List["VectorHyperDual"]:
+    """
+    Obtain differentiable variables 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(atom, int) and atom > 0 for atom in args)
+    # get positions in the flat Cartesian array
+    _x = x.ravel()
+    cart_idxs = []
+    for atom in args:
+        for k in range(3):
+            cart_idxs.append(3 * atom + k)
+    return get_differentiable_vars(
+        values=[_x[idx] for idx in cart_idxs],
+        deriv_order=deriv_order,
+        symbols=[str(idx) for idx in cart_idxs],
+    )
+
+
 class Primitive(ABC):
     """Primitive internal coordinate"""
 
@@ -165,16 +203,9 @@ class PrimitiveInverseDistance(_DistanceFunction):
     def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: int
     ) -> "VectorHyperDual":
-        """Inverse distance"""
-        _x = x.ravel()
-        idxs = []
-        for atom in [self.i, self.j]:
-            for k in range(3):
-                idxs.append(3 * atom + k)
-        i_0, i_1, i_2, j_0, j_1, j_2 = get_differentiable_vars(
-            values=[_x[idx] for idx in idxs],
-            deriv_order=deriv_order,
-            symbols=[str(idx) for idx in idxs],
+        """1 / |x_i - x_j|"""
+        i_0, i_1, i_2, j_0, j_1, j_2 = get_vars_from_atom_idxs(
+            self.i, self.j, x=x, deriv_order=deriv_order
         )
         return 1.0 / _norm_vec3(_sub_vec3([i_0, i_1, i_2], [j_0, j_1, j_2]))
 
@@ -242,6 +273,15 @@ def derivative(
 
         return val if i == self.i else -val
 
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: int
+    ) -> "VectorHyperDual":
+        """|x_i - x_j|"""
+        i_0, i_1, i_2, j_0, j_1, j_2 = get_vars_from_atom_idxs(
+            self.i, self.j, x=x, deriv_order=deriv_order
+        )
+        return _norm_vec3(_sub_vec3([i_0, i_1, i_2], [j_0, j_1, j_2]))
+
     def __call__(self, x: "CartesianCoordinates") -> float:
         """|x_i - x_j|"""
         _x = x.reshape((-1, 3))
@@ -295,6 +335,18 @@ def __eq__(self, other) -> bool:
             and other._ordered_idxs == self._ordered_idxs
         )
 
+    def _evaluate(self, x: "CartesianCoordinates", deriv_order: int):
+        """m - o - n angle"""
+        m_0, m_1, m_2, o_0, o_1, o_2, n_0, n_1, n_2 = get_vars_from_atom_idxs(
+            self.m, self.o, self.n, x=x, deriv_order=deriv_order
+        )
+        u = _sub_vec3([m_0, m_1, m_2], [o_0, o_1, o_2])
+        v = _sub_vec3([n_0, n_1, n_2], [o_0, o_1, o_2])
+        theta = DifferentiableMath.acos(
+            _dot_vec3(u, v) / (_norm_vec3(u) * _norm_vec3(v))
+        )
+        return theta
+
     def __call__(self, x: "CartesianCoordinates") -> float:
         _x = x.reshape((-1, 3))
         u = _x[self.m, :] - _x[self.o, :]

From b5a602feabd4def08348fc1ece913350a6671323 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 16 Nov 2023 19:25:04 +0000
Subject: [PATCH 04/37] autodiff update

---
 autode/opt/coordinates/primitives.py | 34 ++++++++++++++++------------
 1 file changed, 19 insertions(+), 15 deletions(-)

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 41c8a8a8f..465b66fce 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -89,14 +89,14 @@ def _evaluate(self, x: "CartesianCoordinates", deriv_order: int):
     @abstractmethod
     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=0)
+        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
 
@@ -110,17 +110,21 @@ 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 coordinates
+            x: Cartesian coordinate array of shape (N, )
 
         Returns:
-            (float): Derivative
+            (np.ndarray): Derivative array of shape (N, )
         """
+        _x = x.ravel()
+        res = self._evaluate(_x, deriv_order=1)
+        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
 
     @abstractmethod
     def __eq__(self, other):
@@ -209,7 +213,7 @@ def _evaluate(
         )
         return 1.0 / _norm_vec3(_sub_vec3([i_0, i_1, i_2], [j_0, j_1, j_2]))
 
-    def derivative(
+    def legacy_derivative(
         self,
         i: int,
         component: "CartesianComponent",
@@ -250,7 +254,7 @@ class PrimitiveDistance(_DistanceFunction):
         q = |\boldsymbol{X}_i - \boldsymbol{X}_j|
     """
 
-    def derivative(
+    def legacy_derivative(
         self,
         i: int,
         component: "CartesianComponent",
@@ -355,7 +359,7 @@ def __call__(self, x: "CartesianCoordinates") -> float:
         theta = np.arccos(u.dot(v) / (np.linalg.norm(u) * np.linalg.norm(v)))
         return theta
 
-    def derivative(
+    def legacy_derivative(
         self,
         i: int,
         component: "CartesianComponent",
@@ -453,7 +457,7 @@ def __call__(self, x: "CartesianCoordinates") -> float:
         """Value of the dihedral"""
         return self._value(x, return_derivative=False)
 
-    def derivative(
+    def legacy_derivative(
         self,
         i: int,
         component: "CartesianComponent",

From 850c22990c99e0b2b7f9654fcb44976b6cf7000b Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 16 Nov 2023 21:16:38 +0000
Subject: [PATCH 05/37] autodiff update

---
 autode/opt/coordinates/primitives.py | 44 ++++++++++++++++++++++++++++
 1 file changed, 44 insertions(+)

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 465b66fce..d76018ce1 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -30,6 +30,17 @@ def _dot_vec3(vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]):
     return vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2]
 
 
+def _cross_vec3(vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]):
+    """Evaluate vec1 x vec2 for 3D vectors"""
+    assert len(vec1) == len(vec2) == 3
+    # a cross b = a1 b2 - a2 b1, a2 b0 - a0 b2, a0 b1 - a1 b0
+    return [
+        vec1[1] * vec2[2] - vec1[2] * vec2[1],
+        vec1[2] * vec2[0] - vec1[0] * vec2[2],
+        vec1[0] * vec2[1] - vec1[1] * vec2[0],
+    ]
+
+
 def get_vars_from_atom_idxs(
     *args,
     x: "CartesianCoordinates",
@@ -472,6 +483,39 @@ def __eq__(self, other) -> bool:
             or self._atom_indexes == tuple(reversed(other._atom_indexes))
         )
 
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: int
+    ) -> "VectorHyperDual":
+        """Dihedral m-o-p-n"""
+        _x = x.ravel()
+        variables = get_vars_from_atom_idxs(
+            self.m, self.o, self.p, self.n, x=_x, deriv_order=deriv_order
+        )
+        _m = variables[:3]
+        _o = variables[3:6]
+        _p = variables[6:9]
+        _n = variables[9:12]
+
+        # todo do we need to normalise these vectors?
+        u = _sub_vec3(_m, _o)
+        lambda_u = _norm_vec3(u)
+        u = [k / lambda_u for k in u]
+
+        v = _sub_vec3(_n, _p)
+        lambda_v = _norm_vec3(v)
+        v = [k / lambda_v for k in v]
+
+        w = _sub_vec3(_p, _o)
+        lambda_w = _norm_vec3(w)
+        w = [k / lambda_w for k in w]
+
+        v1 = _cross_vec3(u, w)
+        v2 = _cross_vec3([-k for k in w], v)
+
+        return DifferentiableMath.atan2(
+            _dot_vec3(_cross_vec3(v1, w), v2), _dot_vec3(v1, v2)
+        )
+
     def _value(self, x, i=None, component=None, return_derivative=False):
         """Evaluate either the value or the derivative. Shared function
         to reuse local variables"""

From 60ebc13cc7ac6f5f073ebcd93fb42d5119923519 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Mon, 20 Nov 2023 14:07:57 +0000
Subject: [PATCH 06/37] doc update

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

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index d76018ce1..9f677470c 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -97,7 +97,6 @@ def _evaluate(self, x: "CartesianCoordinates", deriv_order: int):
             (VectorHyperDual): The result, optionally containing derivatives
         """
 
-    @abstractmethod
     def __call__(self, x: "CartesianCoordinates") -> float:
         """Return the value of this PIC given a set of cartesian coordinates"""
         _x = x.ravel()
@@ -109,7 +108,7 @@ def derivative(
         x: "CartesianCoordinates",
     ) -> np.ndarray:
         r"""
-        Calculate the derivative with respect to a cartesian coordinate
+        Calculate the derivatives with respect to cartesian coordinates
 
         .. math::
 
@@ -137,6 +136,41 @@ def derivative(
 
         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:
+            (np.ndarray): Second derivative matrix of shape (N, N)
+        """
+        _x = x.ravel()
+        x_n, _ = _x.shape
+        res = self._evaluate(_x, deriv_order=2)
+        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):
         """Comparison of two primitive coordinates"""
@@ -207,7 +241,7 @@ def __eq__(self, other) -> bool:
 
 class PrimitiveInverseDistance(_DistanceFunction):
     """
-    Inverse distance between to atoms:
+    Inverse distance between two atoms:
 
     .. math::
 

From 5748900f285d52376f50b59f335bae59398f46d1 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 20 Nov 2023 22:08:07 +0000
Subject: [PATCH 07/37] autodiff update

---
 autode/opt/coordinates/primitives.py | 27 +++++++++------------------
 1 file changed, 9 insertions(+), 18 deletions(-)

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 9f677470c..9e6458031 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -531,24 +531,15 @@ def _evaluate(
         _n = variables[9:12]
 
         # todo do we need to normalise these vectors?
-        u = _sub_vec3(_m, _o)
-        lambda_u = _norm_vec3(u)
-        u = [k / lambda_u for k in u]
-
-        v = _sub_vec3(_n, _p)
-        lambda_v = _norm_vec3(v)
-        v = [k / lambda_v for k in v]
-
-        w = _sub_vec3(_p, _o)
-        lambda_w = _norm_vec3(w)
-        w = [k / lambda_w for k in w]
-
-        v1 = _cross_vec3(u, w)
-        v2 = _cross_vec3([-k for k in w], v)
-
-        return DifferentiableMath.atan2(
-            _dot_vec3(_cross_vec3(v1, w), v2), _dot_vec3(v1, v2)
-        )
+        u_1 = _sub_vec3(_o, _m)
+        u_2 = _sub_vec3(_p, _o)
+        u_3 = _sub_vec3(_n, _p)
+
+        norm_u2 = _norm_vec3(u_2)
+        v1 = _cross_vec3(u_2, u_3)
+        v2 = _cross_vec3(u_1, u_2)
+        v3 = [k * norm_u2 for k in u_1]
+        return DifferentiableMath.atan2(_dot_vec3(v3, v1), _dot_vec3(v2, v1))
 
     def _value(self, x, i=None, component=None, return_derivative=False):
         """Evaluate either the value or the derivative. Shared function

From d8bbd3629cf9191936887718d3e32ec9c36c7458 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 20 Nov 2023 22:17:56 +0000
Subject: [PATCH 08/37] test update

---
 autode/opt/coordinates/primitives.py |  2 +-
 tests/test_opt/test_coordiantes.py   | 14 ++++----------
 2 files changed, 5 insertions(+), 11 deletions(-)

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 9e6458031..e98f28311 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -59,7 +59,7 @@ def get_vars_from_atom_idxs(
     Returns:
         (list[VectorHyperDual]): A list of differentiable variables
     """
-    assert all(isinstance(atom, int) and atom > 0 for atom in args)
+    assert all(isinstance(atom, int) and atom >= 0 for atom in args)
     # get positions in the flat Cartesian array
     _x = x.ravel()
     cart_idxs = []
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 2f775e81b..c75450fac 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():

From a76f310ae3ee9e16bf04fafa85848f0a89b94626 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 23 Nov 2023 20:50:22 +0000
Subject: [PATCH 09/37] type checking and B matrix

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

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 675ca1f1c..8b83b6d2f 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -82,6 +82,16 @@ def __init__(self, *args: Any):
                 f"from {args}. Must be primitive internals"
             )
 
+    def append(self, item) -> None:
+        """Append an item to this set of primitives"""
+        if isinstance(item, Primitive):
+            super().append(item)
+        else:
+            raise TypeError(
+                f"Can only append Primitive type but"
+                f" {type(item)} was provided"
+            )
+
     @property
     def B(self) -> np.ndarray:
         """Wilson B matrix"""
@@ -179,16 +189,7 @@ def _calc_B(self, x: np.ndarray) -> None:
         B = np.zeros(shape=(len(self), 3 * n_atoms))
 
         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

From 3b9d3b27de2d844fee5bd53a16633317a5c9ef13 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 23 Nov 2023 20:58:20 +0000
Subject: [PATCH 10/37] fix some tests

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

diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index c75450fac..1ba92aaac 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -44,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
-    )
-
-    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)
+    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(2, CartesianComponent.x, x=x), 0)
+    for k in (1, 2):
+        assert np.isclose(derivs[3 * 1 + k], 0)
 
 
 def test_primitive_equality():
@@ -557,20 +551,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)

From 2fef780eec107e2266c289821501c83412583497 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Thu, 23 Nov 2023 21:32:43 +0000
Subject: [PATCH 11/37] fix some tests

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

diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 1ba92aaac..1857b6390 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -646,6 +646,22 @@ def numerical_derivative(a, b, h=1e-6):
             assert np.isclose(analytic, numerical, atol=1e-6)
 
 
+def test_primitive_first_derivs(init_coords, prim):
+    def numerical_first_deriv(coords, h=1e-6):
+        coords = coords.flatten()
+        init_prim = prim(coords)
+        derivs = np.zeros_like(coords)
+        for i in range(coords.shape[0]):
+            coords[i] += h
+            derivs[i] = prim(coords) - init_prim
+            coords[i] -= h
+        return derivs
+
+    analytic = prim.derivative(init_coords)
+    numeric = numerical_first_deriv(init_coords)
+    assert np.allclose(analytic, numeric, atol=1e-6)
+
+
 def test_repr():
     """Test that each primitive has a representation"""
 

From 78214a7fc37fc598b089123655eed663998d0734 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Fri, 24 Nov 2023 14:55:12 +0000
Subject: [PATCH 12/37] add tests

---
 autode/opt/coordinates/primitives.py |   2 +-
 tests/test_opt/test_coordiantes.py   | 111 ++++++++++++++-------------
 2 files changed, 60 insertions(+), 53 deletions(-)

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index e98f28311..c44ead5b9 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -160,7 +160,7 @@ def second_derivative(
             (np.ndarray): Second derivative matrix of shape (N, N)
         """
         _x = x.ravel()
-        x_n, _ = _x.shape
+        x_n = _x.shape[0]
         res = self._evaluate(_x, deriv_order=2)
         derivs = np.zeros(shape=(x_n, x_n), dtype=float)
         for i in range(x_n):
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 1857b6390..54554302f 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -588,17 +588,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():
@@ -610,55 +604,76 @@ 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)
-
-    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
+    ),
+    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
 
-    dihedral = PrimitiveDihedralAngle(2, 0, 1, 3)
+    analytic = prim.derivative(init_coords)
+    numeric = numerical_first_deriv(init_coords)
+    assert np.allclose(analytic, numeric, atol=1e-6)
 
-    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)
 
+@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 test_primitive_first_derivs(init_coords, prim):
-    def numerical_first_deriv(coords, h=1e-6):
+    def numerical_second_deriv(coords, h=1e-8):
         coords = coords.flatten()
-        init_prim = prim(coords)
-        derivs = np.zeros_like(coords)
+        derivs = np.zeros((coords.shape[0], coords.shape[0]))
         for i in range(coords.shape[0]):
             coords[i] += h
-            derivs[i] = prim(coords) - init_prim
+            derivs[i] = (prim.derivative(coords) - init_first_der) / h
             coords[i] -= h
         return derivs
 
-    analytic = prim.derivative(init_coords)
-    numeric = numerical_first_deriv(init_coords)
+    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)
 
 
@@ -678,14 +693,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)

From da883e31709337b1e4fcae7e962aadac28a81776 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Fri, 24 Nov 2023 15:13:49 +0000
Subject: [PATCH 13/37] fix test

---
 autode/opt/coordinates/primitives.py | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index c44ead5b9..aef33f5a9 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -371,9 +371,9 @@ 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"""
@@ -493,10 +493,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
+        self.m = int(m)
+        self.o = int(o)
+        self.p = int(p)
+        self.n = int(n)
 
     def __call__(self, x: "CartesianCoordinates") -> float:
         """Value of the dihedral"""

From 74b3419bb0ff1b40a334321655036901d4bc69dc Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Fri, 24 Nov 2023 20:24:59 +0000
Subject: [PATCH 14/37] fix tests 2

---
 autode/opt/coordinates/dic.py | 48 ++++++++++++++++++++---------------
 1 file changed, 28 insertions(+), 20 deletions(-)

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)

From b282539d327583387ecbd93aeb0a321e47a616a1 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Fri, 24 Nov 2023 22:22:58 +0000
Subject: [PATCH 15/37] fix tests 2

---
 autode/opt/coordinates/primitives.py | 27 ++++-----------------------
 1 file changed, 4 insertions(+), 23 deletions(-)

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index aef33f5a9..22f246437 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -258,6 +258,9 @@ def _evaluate(
         )
         return 1.0 / _norm_vec3(_sub_vec3([i_0, i_1, i_2], [j_0, j_1, j_2]))
 
+    def __repr__(self):
+        return f"InverseDistance({self.i}-{self.j})"
+
     def legacy_derivative(
         self,
         i: int,
@@ -284,11 +287,6 @@ def legacy_derivative(
         else:  # i == self.idx_j:
             return (_x[self.i, k] - _x[self.j, k]) * self(x) ** 3
 
-    def __call__(self, x: "CartesianCoordinates") -> float:
-        """1 / |x_i - x_j|"""
-        _x = x.reshape((-1, 3))
-        return 1.0 / np.linalg.norm(_x[self.i] - _x[self.j])
-
 
 class PrimitiveDistance(_DistanceFunction):
     """
@@ -331,11 +329,6 @@ def _evaluate(
         )
         return _norm_vec3(_sub_vec3([i_0, i_1, i_2], [j_0, j_1, j_2]))
 
-    def __call__(self, x: "CartesianCoordinates") -> float:
-        """|x_i - x_j|"""
-        _x = x.reshape((-1, 3))
-        return np.linalg.norm(_x[self.i] - _x[self.j])
-
     def __repr__(self):
         return f"Distance({self.i}-{self.j})"
 
@@ -396,14 +389,6 @@ def _evaluate(self, x: "CartesianCoordinates", deriv_order: int):
         )
         return theta
 
-    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 legacy_derivative(
         self,
         i: int,
@@ -497,10 +482,7 @@ def __init__(self, m: int, o: int, p: int, n: int):
         self.o = int(o)
         self.p = int(p)
         self.n = int(n)
-
-    def __call__(self, x: "CartesianCoordinates") -> float:
-        """Value of the dihedral"""
-        return self._value(x, return_derivative=False)
+        # TODO: check the sign matches with mol.dihedral()
 
     def legacy_derivative(
         self,
@@ -530,7 +512,6 @@ def _evaluate(
         _p = variables[6:9]
         _n = variables[9:12]
 
-        # todo do we need to normalise these vectors?
         u_1 = _sub_vec3(_o, _m)
         u_2 = _sub_vec3(_p, _o)
         u_3 = _sub_vec3(_n, _p)

From 946d4863ebd2e8ad47ff7722c5bdb8ea0fe6ddd0 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Fri, 24 Nov 2023 22:28:25 +0000
Subject: [PATCH 16/37] doc update

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

diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 22f246437..2f72bedde 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -12,25 +12,31 @@
     from autode.opt.coordinates.autodiff import VectorHyperDual
 
 
-def _norm_vec3(vec: List["VectorHyperDual"]):
+def _norm_vec3(vec: List["VectorHyperDual"]) -> "VectorHyperDual":
     """Norm of a 3D vector"""
     assert len(vec) == 3
     return DifferentiableMath.sqrt(vec[0] ** 2 + vec[1] ** 2 + vec[2] ** 2)
 
 
-def _sub_vec3(vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]):
+def _sub_vec3(
+    vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]
+) -> List["VectorHyperDual"]:
     """Evaluate vec1 - vec2 for 3D vectors"""
     assert len(vec1) == len(vec2) == 3
     return [vec1[i] - vec2[i] for i in range(3)]
 
 
-def _dot_vec3(vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]):
+def _dot_vec3(
+    vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]
+) -> "VectorHyperDual":
     """Evaluate vec1.vec2 for 3D vectors"""
     assert len(vec1) == len(vec2) == 3
     return vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2]
 
 
-def _cross_vec3(vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]):
+def _cross_vec3(
+    vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]
+) -> List["VectorHyperDual"]:
     """Evaluate vec1 x vec2 for 3D vectors"""
     assert len(vec1) == len(vec2) == 3
     # a cross b = a1 b2 - a2 b1, a2 b0 - a0 b2, a0 b1 - a1 b0
@@ -41,7 +47,7 @@ def _cross_vec3(vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]):
     ]
 
 
-def get_vars_from_atom_idxs(
+def _get_vars_from_atom_idxs(
     *args,
     x: "CartesianCoordinates",
     deriv_order: int,
@@ -253,7 +259,7 @@ def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: int
     ) -> "VectorHyperDual":
         """1 / |x_i - x_j|"""
-        i_0, i_1, i_2, j_0, j_1, j_2 = get_vars_from_atom_idxs(
+        i_0, i_1, i_2, j_0, j_1, j_2 = _get_vars_from_atom_idxs(
             self.i, self.j, x=x, deriv_order=deriv_order
         )
         return 1.0 / _norm_vec3(_sub_vec3([i_0, i_1, i_2], [j_0, j_1, j_2]))
@@ -324,7 +330,7 @@ def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: int
     ) -> "VectorHyperDual":
         """|x_i - x_j|"""
-        i_0, i_1, i_2, j_0, j_1, j_2 = get_vars_from_atom_idxs(
+        i_0, i_1, i_2, j_0, j_1, j_2 = _get_vars_from_atom_idxs(
             self.i, self.j, x=x, deriv_order=deriv_order
         )
         return _norm_vec3(_sub_vec3([i_0, i_1, i_2], [j_0, j_1, j_2]))
@@ -379,11 +385,14 @@ def __eq__(self, other) -> bool:
 
     def _evaluate(self, x: "CartesianCoordinates", deriv_order: int):
         """m - o - n angle"""
-        m_0, m_1, m_2, o_0, o_1, o_2, n_0, n_1, n_2 = get_vars_from_atom_idxs(
+        variables = _get_vars_from_atom_idxs(
             self.m, self.o, self.n, x=x, deriv_order=deriv_order
         )
-        u = _sub_vec3([m_0, m_1, m_2], [o_0, o_1, o_2])
-        v = _sub_vec3([n_0, n_1, n_2], [o_0, o_1, o_2])
+        _m = variables[0:3]
+        _o = variables[3:6]
+        _n = variables[6:9]
+        u = _sub_vec3(_m, _o)
+        v = _sub_vec3(_n, _o)
         theta = DifferentiableMath.acos(
             _dot_vec3(u, v) / (_norm_vec3(u) * _norm_vec3(v))
         )
@@ -504,7 +513,7 @@ def _evaluate(
     ) -> "VectorHyperDual":
         """Dihedral m-o-p-n"""
         _x = x.ravel()
-        variables = get_vars_from_atom_idxs(
+        variables = _get_vars_from_atom_idxs(
             self.m, self.o, self.p, self.n, x=_x, deriv_order=deriv_order
         )
         _m = variables[:3]

From 19f8147cd2f249e3d6d36552ea92caa120834b70 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 25 Nov 2023 10:33:47 +0000
Subject: [PATCH 17/37] add test

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

diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 54554302f..853c27cea 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -604,6 +604,20 @@ def test_dihedral_equality():
     )
 
 
+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)
+
+
 # fmt: off
 dihedral_mols = [
     Molecule(

From eb8da4cc7a6ff5cf2aacdd0b3c5af90f72e46e05 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 25 Nov 2023 16:12:44 +0000
Subject: [PATCH 18/37] minor refactor autodiff

---
 autode/opt/coordinates/autodiff.py   |  8 ++-----
 autode/opt/coordinates/primitives.py |  2 +-
 tests/test_opt/test_autodiff.py      | 32 ++++++++++++++++++++++++++++
 3 files changed, 35 insertions(+), 7 deletions(-)
 create mode 100644 tests/test_opt/test_autodiff.py

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index ab4a2348b..a4795aff8 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -12,24 +12,20 @@
 
 def get_differentiable_vars(
     values: Sequence[numeric_type],
+    symbols: Sequence[str],
     deriv_order: int = 2,
-    symbols: Optional[Sequence[str]] = None,
 ):
     """
     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
-        symbols: Optional list of symbols for the numbers,
-                 if not given, numeric strings ('0', '1', '2', ...)
-                 are generated in order
 
     Returns:
         (list[VectorHyperDual]): A list of hyper dual numbers
     """
-    if symbols is None:
-        symbols = [str(idx) for idx in range(len(values))]
     assert all(isinstance(sym, str) for sym in symbols)
     assert len(symbols) == len(values)
     symbols = list(symbols)
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 2f72bedde..ddd0535e6 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -74,8 +74,8 @@ def _get_vars_from_atom_idxs(
             cart_idxs.append(3 * atom + k)
     return get_differentiable_vars(
         values=[_x[idx] for idx in cart_idxs],
-        deriv_order=deriv_order,
         symbols=[str(idx) for idx in cart_idxs],
+        deriv_order=deriv_order,
     )
 
 
diff --git a/tests/test_opt/test_autodiff.py b/tests/test_opt/test_autodiff.py
new file mode 100644
index 000000000..dc7c49968
--- /dev/null
+++ b/tests/test_opt/test_autodiff.py
@@ -0,0 +1,32 @@
+import pytest
+from autode.opt.coordinates.autodiff import (
+    get_differentiable_vars,
+    DifferentiableMath,
+)
+
+
+def test_hyperdual_sanity_checks():
+    a, b = get_differentiable_vars(
+        values=[1, 2], symbols=["a", "b"], deriv_order=1
+    )
+
+    c = get_differentiable_vars(values=[1], symbols=["c"], deriv_order=1)[0]
+    with pytest.raises(ValueError, match="Incompatible number"):
+        _ = a + c
+
+    c, d = get_differentiable_vars(
+        values=[1, 2], symbols=["c", "d"], deriv_order=1
+    )
+    with pytest.raises(ValueError, match="symbols do not match"):
+        _ = a * c
+
+    a2, b2 = get_differentiable_vars(
+        values=[1, 2], symbols=["a", "b"], deriv_order=2
+    )
+    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=1
+        )

From a2e9761859962b43c7dfa20bf972407493cc874b Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 25 Nov 2023 17:57:17 +0000
Subject: [PATCH 19/37] add more tests

---
 autode/opt/coordinates/autodiff.py | 20 ++++++++---
 tests/test_opt/test_autodiff.py    | 53 ++++++++++++++++++++++++++++++
 2 files changed, 69 insertions(+), 4 deletions(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index a4795aff8..ef18c27f2 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -205,7 +205,9 @@ def differentiate_wrt(
         assert self._second_der is not None
         return self._second_der[idx_1, idx_2]
 
-    def __add__(self, other):
+    def __add__(
+        self, other: Union["VectorHyperDual", numeric_type]
+    ) -> "VectorHyperDual":
         """Adding a hyper dual number"""
 
         # add to a float
@@ -222,10 +224,14 @@ def __add__(self, other):
             if self._order == 0:
                 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 == 1:
                 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)
 
@@ -236,7 +242,7 @@ def __radd__(self, other):
         """Addition is commutative"""
         return self.__add__(other)
 
-    def __neg__(self):
+    def __neg__(self) -> "VectorHyperDual":
         """Unary negative operation"""
         new = self.copy()
         new._val = -new._val
@@ -254,7 +260,7 @@ def __rsub__(self, other):
         """Reverse subtraction"""
         return other + (-self)
 
-    def __mul__(self, other):
+    def __mul__(self, other) -> "VectorHyperDual":
         """Multiply a hyper dual number with float or another hyper dual"""
         if isinstance(other, numeric):
             new = self.copy()
@@ -269,12 +275,16 @@ def __mul__(self, other):
             if self._order == 0:
                 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 == 1:
                 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)
@@ -296,7 +306,7 @@ def __rtruediv__(self, other):
         """Reverse true division"""
         return DifferentiableMath.pow(self, -1).__mul__(other)
 
-    def __pow__(self, power, modulo=None):
+    def __pow__(self, power, modulo=None) -> "VectorHyperDual":
         if modulo is not None:
             raise NotImplementedError("Modulo inverse is not implemented")
 
@@ -501,6 +511,8 @@ def atan2(
         num_x: 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)
 
         def atan2_derivs_x_not_0(y, x):
             return DifferentiableMath.atan(y / x)
diff --git a/tests/test_opt/test_autodiff.py b/tests/test_opt/test_autodiff.py
index dc7c49968..6193e921e 100644
--- a/tests/test_opt/test_autodiff.py
+++ b/tests/test_opt/test_autodiff.py
@@ -1,7 +1,10 @@
+import math
+
 import pytest
 from autode.opt.coordinates.autodiff import (
     get_differentiable_vars,
     DifferentiableMath,
+    VectorHyperDual,
 )
 
 
@@ -30,3 +33,53 @@ def test_hyperdual_sanity_checks():
         _ = get_differentiable_vars(
             values=[1, 2], symbols=["a", "a"], deriv_order=1
         )
+
+
+def x_to_the_y(x, y):
+    return x**y
+
+
+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_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))

From 35d5e4ddec8544ca8e07d438a5bb730a3365d57f Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 25 Nov 2023 20:49:52 +0000
Subject: [PATCH 20/37] remove legacy derivatives and add test

---
 autode/opt/coordinates/internals.py  |   5 +-
 autode/opt/coordinates/primitives.py | 172 ++-------------------------
 tests/test_opt/test_coordiantes.py   |   8 ++
 3 files changed, 18 insertions(+), 167 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index 8b83b6d2f..fb181dacc 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -183,10 +183,9 @@ 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):
             B[i] = primitive.derivative(x=cart_coords)
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index ddd0535e6..4b89d9b0a 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -246,7 +246,7 @@ def __eq__(self, other) -> bool:
 
 
 class PrimitiveInverseDistance(_DistanceFunction):
-    """
+    r"""
     Inverse distance between two atoms:
 
     .. math::
@@ -267,35 +267,9 @@ def _evaluate(
     def __repr__(self):
         return f"InverseDistance({self.i}-{self.j})"
 
-    def legacy_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
-
 
 class PrimitiveDistance(_DistanceFunction):
-    """
+    r"""
     Distance between two atoms:
 
     .. math::
@@ -303,29 +277,6 @@ class PrimitiveDistance(_DistanceFunction):
         q = |\boldsymbol{X}_i - \boldsymbol{X}_j|
     """
 
-    def legacy_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 _evaluate(
         self, x: "CartesianCoordinates", deriv_order: int
     ) -> "VectorHyperDual":
@@ -366,6 +317,11 @@ 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)
@@ -398,51 +354,6 @@ def _evaluate(self, x: "CartesianCoordinates", deriv_order: int):
         )
         return theta
 
-    def legacy_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 __repr__(self):
         return f"Angle({self.m}-{self.o}-{self.n})"
 
@@ -491,15 +402,6 @@ def __init__(self, m: int, o: int, p: int, n: int):
         self.o = int(o)
         self.p = int(p)
         self.n = int(n)
-        # TODO: check the sign matches with mol.dihedral()
-
-    def legacy_derivative(
-        self,
-        i: int,
-        component: "CartesianComponent",
-        x: "CartesianCoordinates",
-    ) -> float:
-        return self._value(x, i=i, component=component, return_derivative=True)
 
     def __eq__(self, other) -> bool:
         """Equality of two distance functions"""
@@ -512,6 +414,7 @@ def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: int
     ) -> "VectorHyperDual":
         """Dihedral m-o-p-n"""
+        # https://en.wikipedia.org/wiki/Dihedral_angle#In_polymer_physics
         _x = x.ravel()
         variables = _get_vars_from_atom_idxs(
             self.m, self.o, self.p, self.n, x=_x, deriv_order=deriv_order
@@ -531,64 +434,5 @@ def _evaluate(
         v3 = [k * norm_u2 for k in u_1]
         return DifferentiableMath.atan2(_dot_vec3(v3, v1), _dot_vec3(v2, v1))
 
-    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 __repr__(self):
         return f"Dihedral({self.m}-{self.o}-{self.p}-{self.n})"
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 853c27cea..9f07c4699 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -527,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(TypeError):
+        c.append(3)
+
+
 def test_constrained_distance_satisfied():
     d = ConstrainedPrimitiveDistance(0, 1, value=1.0)
 

From 396bdc8ffc57afa752b72f523c183a3b2bdf5a61 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 27 Nov 2023 20:32:03 +0000
Subject: [PATCH 21/37] pr suggestions 1

---
 autode/opt/coordinates/internals.py  | 11 +++--------
 autode/opt/coordinates/primitives.py |  8 ++++----
 2 files changed, 7 insertions(+), 12 deletions(-)

diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py
index fb181dacc..9e7462e16 100644
--- a/autode/opt/coordinates/internals.py
+++ b/autode/opt/coordinates/internals.py
@@ -82,15 +82,10 @@ def __init__(self, *args: Any):
                 f"from {args}. Must be primitive internals"
             )
 
-    def append(self, item) -> None:
+    def append(self, item: Primitive) -> None:
         """Append an item to this set of primitives"""
-        if isinstance(item, Primitive):
-            super().append(item)
-        else:
-            raise TypeError(
-                f"Can only append Primitive type but"
-                f" {type(item)} was provided"
-            )
+        assert isinstance(item, Primitive), "Must be a Primitive type!"
+        super().append(item)
 
     @property
     def B(self) -> np.ndarray:
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 4b89d9b0a..3c3c69867 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -48,7 +48,7 @@ def _cross_vec3(
 
 
 def _get_vars_from_atom_idxs(
-    *args,
+    *args: int,
     x: "CartesianCoordinates",
     deriv_order: int,
 ) -> List["VectorHyperDual"]:
@@ -65,13 +65,13 @@ def _get_vars_from_atom_idxs(
     Returns:
         (list[VectorHyperDual]): A list of differentiable variables
     """
-    assert all(isinstance(atom, int) and atom >= 0 for atom in args)
+    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 in args:
+    for atom_idx in args:
         for k in range(3):
-            cart_idxs.append(3 * atom + k)
+            cart_idxs.append(3 * atom_idx + k)
     return get_differentiable_vars(
         values=[_x[idx] for idx in cart_idxs],
         symbols=[str(idx) for idx in cart_idxs],

From a9eeb52c0471f1fb657acff3320990d0610a5cb9 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Mon, 27 Nov 2023 20:57:30 +0000
Subject: [PATCH 22/37] Use enum for derivative order

---
 autode/opt/coordinates/autodiff.py   | 48 ++++++++++++++++++----------
 autode/opt/coordinates/primitives.py | 26 +++++++++------
 2 files changed, 47 insertions(+), 27 deletions(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index ef18c27f2..dc7d51ca3 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -2,6 +2,7 @@
 Automatic differentiation routines in pure Python
 """
 from typing import Union, Callable, Sequence, Optional
+from enum import Enum
 from copy import deepcopy
 import numpy as np
 import math
@@ -10,10 +11,18 @@
 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: int = 2,
+    deriv_order: DerivativeOrder = DerivativeOrder.second,
 ):
     """
     Obtain differentiable variables from a series of numbers
@@ -74,7 +83,7 @@ def __init__(
         # load the derivatives with sanity checks
         self._first_der: Optional[np.ndarray] = None
         self._second_der: Optional[np.ndarray] = None
-        self._order = 0
+        self._order = DerivativeOrder.zeroth
         if first_der is None:
             return
         assert isinstance(first_der, np.ndarray)
@@ -85,7 +94,7 @@ def __init__(
                 f" shape of derivative array {first_der.shape}"
             )
         self._first_der = first_der.astype(float)
-        self._order = 1
+        self._order = DerivativeOrder.first
 
         if second_der is None:
             return
@@ -96,7 +105,7 @@ def __init__(
                 f" shape of second derivative matrix {first_der.shape}"
             )
         self._second_der = second_der.astype(float)
-        self._order = 2
+        self._order = DerivativeOrder.second
 
     @property
     def n_vars(self) -> int:
@@ -132,7 +141,11 @@ def value(self, value):
 
     @classmethod
     def from_variable(
-        cls, value: float, symbol: str, all_symbols: Sequence[str], order: int
+        cls,
+        value: float,
+        symbol: str,
+        all_symbols: Sequence[str],
+        order: DerivativeOrder,
     ):
         """
         Create a hyper-dual number from one variable, requires
@@ -159,10 +172,10 @@ def from_variable(
         n = len(all_symbols)
         idx = list(all_symbols).index(symbol)
 
-        if order >= 1:
+        if order == DerivativeOrder.first or order == DerivativeOrder.second:
             first_der = np.zeros(shape=n, dtype=float)
             first_der[idx] = 1.0
-        if order == 2:
+        if order == DerivativeOrder.second:
             second_der = np.zeros(shape=(n, n), dtype=float)
 
         return VectorHyperDual(val, all_symbols, first_der, second_der)
@@ -187,7 +200,7 @@ def differentiate_wrt(
         if symbol1 not in self._symbols:
             return None
 
-        if self._order == 0:
+        if self._order == DerivativeOrder.zeroth:
             return None
 
         idx_1 = self._symbols.index(symbol1)
@@ -200,7 +213,7 @@ def differentiate_wrt(
             return None
         idx_2 = self._symbols.index(symbol2)
         # check if second derivs are available
-        if self._order == 1:
+        if self._order == DerivativeOrder.first:
             return None
         assert self._second_der is not None
         return self._second_der[idx_1, idx_2]
@@ -221,13 +234,13 @@ def __add__(
             self._check_compatible(other)
 
             val = self._val + other._val
-            if self._order == 0:
+            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 == 1:
+            if self._order == DerivativeOrder.first:
                 return VectorHyperDual(val, self._symbols, first_der)
 
             assert self._second_der is not None
@@ -246,9 +259,10 @@ def __neg__(self) -> "VectorHyperDual":
         """Unary negative operation"""
         new = self.copy()
         new._val = -new._val
-        if self._order >= 1:
+        if self._order == DerivativeOrder.first:
+            new._first_der = -new._first_der
+        elif self._order == DerivativeOrder.second:
             new._first_der = -new._first_der
-        if self._order == 2:
             new._second_der = -new._second_der
         return new
 
@@ -272,7 +286,7 @@ def __mul__(self, other) -> "VectorHyperDual":
             self._check_compatible(other)
 
             val = self._val * other._val
-            if self._order == 0:
+            if self._order == DerivativeOrder.zeroth:
                 return VectorHyperDual(val, self._symbols)
 
             assert self._first_der is not None
@@ -280,7 +294,7 @@ def __mul__(self, other) -> "VectorHyperDual":
             first_der = (
                 self._val * other._first_der + other._val * self._first_der
             )
-            if self._order == 1:
+            if self._order == DerivativeOrder.first:
                 return VectorHyperDual(val, self._symbols, first_der)
 
             assert self._second_der is not None
@@ -340,14 +354,14 @@ def apply_operation(
 
         val = operator(num._val)
 
-        if num._order == 0:
+        if num._order == DerivativeOrder.zeroth:
             return VectorHyperDual(val, num._symbols)
 
         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 == 1:
+        if num._order == DerivativeOrder.first:
             return VectorHyperDual(val, num._symbols, first_der)
 
         assert num._second_der is not None
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 3c3c69867..6aea394c8 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -1,10 +1,11 @@
 import numpy as np
 
 from abc import ABC, abstractmethod
-from typing import Tuple, TYPE_CHECKING, List, Optional
+from typing import Tuple, TYPE_CHECKING, List
 from autode.opt.coordinates.autodiff import (
     get_differentiable_vars,
     DifferentiableMath,
+    DerivativeOrder,
 )
 
 if TYPE_CHECKING:
@@ -50,7 +51,7 @@ def _cross_vec3(
 def _get_vars_from_atom_idxs(
     *args: int,
     x: "CartesianCoordinates",
-    deriv_order: int,
+    deriv_order: DerivativeOrder,
 ) -> List["VectorHyperDual"]:
     """
     Obtain differentiable variables from the Cartesian components
@@ -88,7 +89,10 @@ def __init__(self, *atom_indexes: int):
         """A primitive internal coordinate that involves a number of atoms"""
         self._atom_indexes = atom_indexes
 
-    def _evaluate(self, x: "CartesianCoordinates", deriv_order: int):
+    @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.
@@ -106,7 +110,7 @@ def _evaluate(self, x: "CartesianCoordinates", deriv_order: int):
     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=0)
+        res = self._evaluate(_x, deriv_order=DerivativeOrder.zeroth)
         return res.value
 
     def derivative(
@@ -133,7 +137,7 @@ def derivative(
             (np.ndarray): Derivative array of shape (N, )
         """
         _x = x.ravel()
-        res = self._evaluate(_x, deriv_order=1)
+        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))
@@ -167,7 +171,7 @@ def second_derivative(
         """
         _x = x.ravel()
         x_n = _x.shape[0]
-        res = self._evaluate(_x, deriv_order=2)
+        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):
@@ -256,7 +260,7 @@ class PrimitiveInverseDistance(_DistanceFunction):
     """
 
     def _evaluate(
-        self, x: "CartesianCoordinates", deriv_order: int
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
     ) -> "VectorHyperDual":
         """1 / |x_i - x_j|"""
         i_0, i_1, i_2, j_0, j_1, j_2 = _get_vars_from_atom_idxs(
@@ -278,7 +282,7 @@ class PrimitiveDistance(_DistanceFunction):
     """
 
     def _evaluate(
-        self, x: "CartesianCoordinates", deriv_order: int
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
     ) -> "VectorHyperDual":
         """|x_i - x_j|"""
         i_0, i_1, i_2, j_0, j_1, j_2 = _get_vars_from_atom_idxs(
@@ -339,7 +343,9 @@ def __eq__(self, other) -> bool:
             and other._ordered_idxs == self._ordered_idxs
         )
 
-    def _evaluate(self, x: "CartesianCoordinates", deriv_order: int):
+    def _evaluate(
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
+    ):
         """m - o - n angle"""
         variables = _get_vars_from_atom_idxs(
             self.m, self.o, self.n, x=x, deriv_order=deriv_order
@@ -411,7 +417,7 @@ def __eq__(self, other) -> bool:
         )
 
     def _evaluate(
-        self, x: "CartesianCoordinates", deriv_order: int
+        self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
     ) -> "VectorHyperDual":
         """Dihedral m-o-p-n"""
         # https://en.wikipedia.org/wiki/Dihedral_angle#In_polymer_physics

From 07e671f663a160dd7116224386c330f9d5410609 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Wed, 29 Nov 2023 16:05:39 +0000
Subject: [PATCH 23/37] fix tests

---
 autode/opt/coordinates/autodiff.py |  1 +
 tests/test_opt/test_autodiff.py    | 15 ++++++++++-----
 tests/test_opt/test_coordiantes.py |  2 +-
 3 files changed, 12 insertions(+), 6 deletions(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index dc7d51ca3..0cd4e2013 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -171,6 +171,7 @@ def from_variable(
         second_der = None
         n = len(all_symbols)
         idx = list(all_symbols).index(symbol)
+        assert isinstance(order, DerivativeOrder)
 
         if order == DerivativeOrder.first or order == DerivativeOrder.second:
             first_der = np.zeros(shape=n, dtype=float)
diff --git a/tests/test_opt/test_autodiff.py b/tests/test_opt/test_autodiff.py
index 6193e921e..a44c04de8 100644
--- a/tests/test_opt/test_autodiff.py
+++ b/tests/test_opt/test_autodiff.py
@@ -5,33 +5,38 @@
     get_differentiable_vars,
     DifferentiableMath,
     VectorHyperDual,
+    DerivativeOrder,
 )
 
 
 def test_hyperdual_sanity_checks():
     a, b = get_differentiable_vars(
-        values=[1, 2], symbols=["a", "b"], deriv_order=1
+        values=[1, 2], symbols=["a", "b"], deriv_order=DerivativeOrder.first
     )
 
-    c = get_differentiable_vars(values=[1], symbols=["c"], deriv_order=1)[0]
+    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=1
+        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=2
+        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=1
+            values=[1, 2],
+            symbols=["a", "a"],
+            deriv_order=DerivativeOrder.first,
         )
 
 
diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py
index 9f07c4699..9fe31f72e 100644
--- a/tests/test_opt/test_coordiantes.py
+++ b/tests/test_opt/test_coordiantes.py
@@ -531,7 +531,7 @@ def test_pic_append_type_checking():
     c = PIC()
     # append should check for primitive type
     c.append(PrimitiveDistance(0, 1))
-    with pytest.raises(TypeError):
+    with pytest.raises(AssertionError):
         c.append(3)
 
 

From 6a778771cd29e2dbc9c7e1768f4daad90c214dc9 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Wed, 29 Nov 2023 16:39:59 +0000
Subject: [PATCH 24/37] doc update and use lambda func

---
 autode/opt/coordinates/autodiff.py | 99 +++++++++++-------------------
 1 file changed, 35 insertions(+), 64 deletions(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index 0cd4e2013..6ffa3c253 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -1,5 +1,8 @@
 """
 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
@@ -224,7 +227,6 @@ def __add__(
     ) -> "VectorHyperDual":
         """Adding a hyper dual number"""
 
-        # add to a float
         if isinstance(other, numeric):
             new = self.copy()
             new._val += float(other)
@@ -282,7 +284,7 @@ def __mul__(self, other) -> "VectorHyperDual":
             new._val *= float(other)
             return new
 
-        # multiply with another dual
+        # Product rule for derivatives, Eqn (24) in ref. [1]
         elif isinstance(other, VectorHyperDual):
             self._check_compatible(other)
 
@@ -358,6 +360,7 @@ def apply_operation(
         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
@@ -389,34 +392,22 @@ def sqrt(num: Union[VectorHyperDual, numeric_type]):
         else:
             assert num.value > 0
 
-        def value(x0: float):
-            return math.sqrt(x0)
-
-        def derivative(x0: float):
-            return 1 / (2 * math.sqrt(x0))
-
-        def second_derivative(x0: float):
-            return -1 / (4 * math.pow(x0, 3 / 2))
-
         return VectorHyperDual.apply_operation(
-            num, value, derivative, second_derivative
+            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]):
         """Raise e to the power of num"""
 
-        def value(x0: float):
-            return math.exp(x0)
-
-        def derivative(x0: float):
-            return math.exp(x0)
-
-        def second_derivative(x0: float):
-            return math.exp(x0)
-
         return VectorHyperDual.apply_operation(
-            num, value, derivative, second_derivative
+            num,
+            operator=lambda x0: math.exp(x0),
+            operator_first_deriv=lambda x0: math.exp(x0),
+            operator_second_deriv=lambda x0: math.exp(x0),
         )
 
     @staticmethod
@@ -424,7 +415,7 @@ def pow(
         num: Union[VectorHyperDual, numeric_type],
         power: Union[VectorHyperDual, numeric_type],
     ):
-        """Exponentiation"""
+        """Exponentiation of one hyperdual to another"""
 
         if isinstance(num, numeric) and isinstance(power, numeric):
             return math.pow(num, power)
@@ -434,18 +425,14 @@ def pow(
                 raise ValueError(
                     "Math error, can't raise negative number to fractional power"
                 )
-
-            def value(x0: float):
-                return math.pow(x0, power)
-
-            def derivative(x0: float):
-                return power * math.pow(x0, power - 1)
-
-            def second_derivative(x0: float):
-                return power * (power - 1) * math.pow(x0, power - 2)
-
             return VectorHyperDual.apply_operation(
-                num, value, derivative, second_derivative
+                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):
@@ -468,17 +455,11 @@ def log(num: Union[VectorHyperDual, numeric_type]):
         else:
             assert num.value > 0
 
-        def value(x0: float):
-            return math.log(x0)
-
-        def derivative(x0: float):
-            return 1.0 / x0
-
-        def second_derivative(x0: float):
-            return -1.0 / (x0**2)
-
         return VectorHyperDual.apply_operation(
-            num, value, derivative, second_derivative
+            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
@@ -490,34 +471,23 @@ def acos(num: Union[VectorHyperDual, numeric_type]):
         else:
             assert -1 < num < 1
 
-        def value(x0: float):
-            return math.acos(x0)
-
-        def derivative(x0: float):
-            return -1 / math.sqrt(1 - x0**2)
-
-        def second_derivative(x0: float):
-            return -x0 / math.pow(1 - x0**2, 3 / 2)
-
         return VectorHyperDual.apply_operation(
-            num, value, derivative, second_derivative
+            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]):
         """Calculate the arctangent of a hyperdual number"""
 
-        def value(x0: float):
-            return math.atan(x0)
-
-        def derivative(x0: float):
-            return 1 / (1 + x0**2)
-
-        def second_derivative(x0: float):
-            return (-2 * x0) / (x0**2 + 1) ** 2
-
         return VectorHyperDual.apply_operation(
-            num, value, derivative, second_derivative
+            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
@@ -529,6 +499,7 @@ def atan2(
         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)
 

From b4bd9eb964618a905ac6a25c67cc9added5b7e27 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Wed, 29 Nov 2023 20:06:12 +0000
Subject: [PATCH 25/37] differentiable vector class

---
 autode/opt/coordinates/autodiff.py   | 71 +++++++++++++++++++-
 autode/opt/coordinates/primitives.py | 97 +++++++++-------------------
 2 files changed, 98 insertions(+), 70 deletions(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index 6ffa3c253..1fd125598 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -338,9 +338,10 @@ def apply_operation(
         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)
+        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)
@@ -519,3 +520,69 @@ def atan2_derivs_x_close_0(y, x):
             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"]):
+        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 = list(items)
+
+    @staticmethod
+    def _check_same_type(other):
+        if not isinstance(other, DifferentiableVector3D):
+            raise ValueError("Operation must be done with another 3D vector!")
+
+    def dot(self, other: "DifferentiableVector3D") -> "VectorHyperDual":
+        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":
+        return DifferentiableMath.sqrt(
+            self._data[0] ** 2 + self._data[1] ** 2 + self._data[2] ** 2
+        )
+
+    def __add__(
+        self, other: "DifferentiableVector3D"
+    ) -> "DifferentiableVector3D":
+        self._check_same_type(other)
+        return DifferentiableVector3D(
+            [self._data[k] + other._data[k] for k in range(3)]
+        )
+
+    def __neg__(self) -> "DifferentiableVector3D":
+        return DifferentiableVector3D([-self._data[k] for k in range(3)])
+
+    def __sub__(self, other) -> "DifferentiableVector3D":
+        return self.__add__(-other)
+
+    def __mul__(self, other: Union[VectorHyperDual, numeric_type]):
+        assert isinstance(other, numeric) or isinstance(other, VectorHyperDual)
+        return DifferentiableVector3D(
+            [self._data[k] * other for k in range(3)]
+        )
+
+    def __rmul__(self, other):
+        return self.__mul__(other)
+
+    def cross(self, other) -> "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/primitives.py b/autode/opt/coordinates/primitives.py
index 6aea394c8..8f14f914d 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -5,6 +5,7 @@
 from autode.opt.coordinates.autodiff import (
     get_differentiable_vars,
     DifferentiableMath,
+    DifferentiableVector3D,
     DerivativeOrder,
 )
 
@@ -13,48 +14,13 @@
     from autode.opt.coordinates.autodiff import VectorHyperDual
 
 
-def _norm_vec3(vec: List["VectorHyperDual"]) -> "VectorHyperDual":
-    """Norm of a 3D vector"""
-    assert len(vec) == 3
-    return DifferentiableMath.sqrt(vec[0] ** 2 + vec[1] ** 2 + vec[2] ** 2)
-
-
-def _sub_vec3(
-    vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]
-) -> List["VectorHyperDual"]:
-    """Evaluate vec1 - vec2 for 3D vectors"""
-    assert len(vec1) == len(vec2) == 3
-    return [vec1[i] - vec2[i] for i in range(3)]
-
-
-def _dot_vec3(
-    vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]
-) -> "VectorHyperDual":
-    """Evaluate vec1.vec2 for 3D vectors"""
-    assert len(vec1) == len(vec2) == 3
-    return vec1[0] * vec2[0] + vec1[1] * vec2[1] + vec1[2] * vec2[2]
-
-
-def _cross_vec3(
-    vec1: List["VectorHyperDual"], vec2: List["VectorHyperDual"]
-) -> List["VectorHyperDual"]:
-    """Evaluate vec1 x vec2 for 3D vectors"""
-    assert len(vec1) == len(vec2) == 3
-    # a cross b = a1 b2 - a2 b1, a2 b0 - a0 b2, a0 b1 - a1 b0
-    return [
-        vec1[1] * vec2[2] - vec1[2] * vec2[1],
-        vec1[2] * vec2[0] - vec1[0] * vec2[2],
-        vec1[0] * vec2[1] - vec1[1] * vec2[0],
-    ]
-
-
-def _get_vars_from_atom_idxs(
+def _get_3d_vecs_from_atom_idxs(
     *args: int,
     x: "CartesianCoordinates",
     deriv_order: DerivativeOrder,
-) -> List["VectorHyperDual"]:
+) -> List[DifferentiableVector3D]:
     """
-    Obtain differentiable variables from the Cartesian components
+    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.
 
@@ -73,11 +39,17 @@ def _get_vars_from_atom_idxs(
     for atom_idx in args:
         for k in range(3):
             cart_idxs.append(3 * atom_idx + k)
-    return get_differentiable_vars(
+    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):
@@ -263,10 +235,10 @@ def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
     ) -> "VectorHyperDual":
         """1 / |x_i - x_j|"""
-        i_0, i_1, i_2, j_0, j_1, j_2 = _get_vars_from_atom_idxs(
+        vec_i, vec_j = _get_3d_vecs_from_atom_idxs(
             self.i, self.j, x=x, deriv_order=deriv_order
         )
-        return 1.0 / _norm_vec3(_sub_vec3([i_0, i_1, i_2], [j_0, j_1, j_2]))
+        return 1.0 / (vec_i - vec_j).norm()
 
     def __repr__(self):
         return f"InverseDistance({self.i}-{self.j})"
@@ -285,10 +257,10 @@ def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
     ) -> "VectorHyperDual":
         """|x_i - x_j|"""
-        i_0, i_1, i_2, j_0, j_1, j_2 = _get_vars_from_atom_idxs(
+        vec_i, vec_j = _get_3d_vecs_from_atom_idxs(
             self.i, self.j, x=x, deriv_order=deriv_order
         )
-        return _norm_vec3(_sub_vec3([i_0, i_1, i_2], [j_0, j_1, j_2]))
+        return (vec_i - vec_j).norm()
 
     def __repr__(self):
         return f"Distance({self.i}-{self.j})"
@@ -347,18 +319,12 @@ def _evaluate(
         self, x: "CartesianCoordinates", deriv_order: DerivativeOrder
     ):
         """m - o - n angle"""
-        variables = _get_vars_from_atom_idxs(
+        vec_m, vec_o, vec_n = _get_3d_vecs_from_atom_idxs(
             self.m, self.o, self.n, x=x, deriv_order=deriv_order
         )
-        _m = variables[0:3]
-        _o = variables[3:6]
-        _n = variables[6:9]
-        u = _sub_vec3(_m, _o)
-        v = _sub_vec3(_n, _o)
-        theta = DifferentiableMath.acos(
-            _dot_vec3(u, v) / (_norm_vec3(u) * _norm_vec3(v))
-        )
-        return theta
+        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})"
@@ -422,23 +388,18 @@ def _evaluate(
         """Dihedral m-o-p-n"""
         # https://en.wikipedia.org/wiki/Dihedral_angle#In_polymer_physics
         _x = x.ravel()
-        variables = _get_vars_from_atom_idxs(
+        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
         )
-        _m = variables[:3]
-        _o = variables[3:6]
-        _p = variables[6:9]
-        _n = variables[9:12]
-
-        u_1 = _sub_vec3(_o, _m)
-        u_2 = _sub_vec3(_p, _o)
-        u_3 = _sub_vec3(_n, _p)
-
-        norm_u2 = _norm_vec3(u_2)
-        v1 = _cross_vec3(u_2, u_3)
-        v2 = _cross_vec3(u_1, u_2)
-        v3 = [k * norm_u2 for k in u_1]
-        return DifferentiableMath.atan2(_dot_vec3(v3, v1), _dot_vec3(v2, v1))
+        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
+        return DifferentiableMath.atan2(v3.dot(v1), v2.dot(v1))
 
     def __repr__(self):
         return f"Dihedral({self.m}-{self.o}-{self.p}-{self.n})"

From 15b0f6dbfc61b6946fc9d0268288cf08d3e94d22 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Thu, 30 Nov 2023 14:24:44 +0000
Subject: [PATCH 26/37] doc update

---
 autode/opt/coordinates/autodiff.py | 77 ++++++++++++++++++++++++++++--
 1 file changed, 74 insertions(+), 3 deletions(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index 1fd125598..d2e1ae1d6 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -529,17 +529,34 @@ class DifferentiableVector3D:
     """
 
     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
+        """
         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 = list(items)
 
     @staticmethod
-    def _check_same_type(other):
+    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):
@@ -548,6 +565,12 @@ def dot(self, other: "DifferentiableVector3D") -> "VectorHyperDual":
         return dot
 
     def norm(self) -> "VectorHyperDual":
+        """
+        Euclidean (l2) norm of this 3D vector
+
+        Returns:
+            (VectorHyperDual): A scalar number (with derivatives)
+        """
         return DifferentiableMath.sqrt(
             self._data[0] ** 2 + self._data[1] ** 2 + self._data[2] ** 2
         )
@@ -555,27 +578,75 @@ def norm(self) -> "VectorHyperDual":
     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]):
+    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":
+    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]

From 4186393639b00d3d11c6a4c599b7a05b731c2bdc Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com>
Date: Thu, 30 Nov 2023 15:04:35 +0000
Subject: [PATCH 27/37] doc update

---
 autode/opt/coordinates/autodiff.py | 40 ++++++++++++++++++++++++------
 1 file changed, 32 insertions(+), 8 deletions(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index d2e1ae1d6..12c7c8d58 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -87,8 +87,18 @@ def __init__(
         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
+            return None
         assert isinstance(first_der, np.ndarray)
         first_der = first_der.flatten()
         if not first_der.shape == (self.n_vars,):
@@ -100,7 +110,7 @@ def __init__(
         self._order = DerivativeOrder.first
 
         if second_der is None:
-            return
+            return None
         assert isinstance(second_der, np.ndarray)
         if not second_der.shape == (self.n_vars, self.n_vars):
             raise ValueError(
@@ -115,10 +125,20 @@ def n_vars(self) -> int:
         """Number of variables in this hyper-dual"""
         return len(self._symbols)
 
-    def copy(self):
+    def copy(self) -> "VectorHyperDual":
         return deepcopy(self)
 
-    def _check_compatible(self, other: "VectorHyperDual"):
+    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, "
@@ -134,11 +154,12 @@ def _check_compatible(self, other: "VectorHyperDual"):
         return None
 
     @property
-    def value(self):
+    def value(self) -> float:
+        """Return the value of the hyper-dual number"""
         return self._val
 
     @value.setter
-    def value(self, value):
+    def value(self, value: float):
         assert isinstance(value, numeric)
         self._val = float(value)
 
@@ -263,8 +284,11 @@ def __neg__(self) -> "VectorHyperDual":
         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
@@ -423,7 +447,7 @@ def pow(
 
         elif isinstance(num, VectorHyperDual) and isinstance(power, numeric):
             if num.value < 0 and isinstance(power, float):
-                raise ValueError(
+                raise AssertionError(
                     "Math error, can't raise negative number to fractional power"
                 )
             return VectorHyperDual.apply_operation(
@@ -440,7 +464,7 @@ def pow(
             if (isinstance(num, numeric) and num) < 0 or (
                 isinstance(num, VectorHyperDual) and num.value < 0
             ):
-                raise ValueError(
+                raise AssertionError(
                     "Only positive numbers can be used with"
                     " differentiable exponent"
                 )

From 331b6cc8d940964d4e53fafc9514d0049fc397f5 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Fri, 1 Dec 2023 11:04:43 +0000
Subject: [PATCH 28/37] coverage update

---
 autode/opt/coordinates/autodiff.py |  3 ++
 tests/test_opt/test_autodiff.py    | 72 +++++++++++++++++++++++++++---
 2 files changed, 70 insertions(+), 5 deletions(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index 12c7c8d58..4e3e8e712 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -471,6 +471,9 @@ def pow(
             # 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]):
         """Natural logarithm"""
diff --git a/tests/test_opt/test_autodiff.py b/tests/test_opt/test_autodiff.py
index a44c04de8..1af5eaf6b 100644
--- a/tests/test_opt/test_autodiff.py
+++ b/tests/test_opt/test_autodiff.py
@@ -1,5 +1,5 @@
 import math
-
+import numpy as np
 import pytest
 from autode.opt.coordinates.autodiff import (
     get_differentiable_vars,
@@ -40,10 +40,6 @@ def test_hyperdual_sanity_checks():
         )
 
 
-def x_to_the_y(x, y):
-    return x**y
-
-
 def x_to_the_y_derivs(x, y):
     d_dx = y * x ** (y - 1)
     d_dy = x**y * math.log(x)
@@ -88,3 +84,69 @@ def test_math_funcs_work_with_native_types():
     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)),
+    )
+
+
+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

From e78b2b8f05f8e031e833245f5773eb2395f45ec7 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Fri, 1 Dec 2023 11:05:11 +0000
Subject: [PATCH 29/37] coverage update

---
 tests/test_opt/test_autodiff.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/test_opt/test_autodiff.py b/tests/test_opt/test_autodiff.py
index 1af5eaf6b..56006b099 100644
--- a/tests/test_opt/test_autodiff.py
+++ b/tests/test_opt/test_autodiff.py
@@ -127,6 +127,7 @@ def test_hyperdual_order():
         first_der=np.zeros(2),
         second_der=np.zeros((2, 2)),
     )
+    assert x._order == DerivativeOrder.second
 
 
 def test_derivative_not_available():

From 1f831c1defef81fd2f55bd5b58fa9ac9161fdb1a Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 2 Dec 2023 11:29:28 +0000
Subject: [PATCH 30/37] doc update

---
 autode/opt/coordinates/autodiff.py   | 38 ++++++++++++++++++++--------
 autode/opt/coordinates/primitives.py |  4 ++-
 2 files changed, 31 insertions(+), 11 deletions(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index 4e3e8e712..3d3ed6e6b 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -351,7 +351,9 @@ def __pow__(self, power, modulo=None) -> "VectorHyperDual":
         if modulo is not None:
             raise NotImplementedError("Modulo inverse is not implemented")
 
-        return DifferentiableMath.pow(self, power)
+        result = DifferentiableMath.pow(self, power)
+        assert isinstance(result, VectorHyperDual)
+        return result
 
     def __rpow__(self, other):
         return DifferentiableMath.pow(other, self)
@@ -380,6 +382,8 @@ def apply_operation(
         if isinstance(num, numeric):
             return operator(float(num))
 
+        assert isinstance(num, VectorHyperDual)
+
         val = operator(num._val)
 
         if num._order == DerivativeOrder.zeroth:
@@ -409,7 +413,9 @@ class DifferentiableMath:
     """
 
     @staticmethod
-    def sqrt(num: Union[VectorHyperDual, numeric_type]):
+    def sqrt(
+        num: Union[VectorHyperDual, numeric_type]
+    ) -> Union[VectorHyperDual, numeric_type]:
         """Calculate the square root of a hyperdual number"""
 
         if isinstance(num, numeric):
@@ -425,7 +431,9 @@ def sqrt(num: Union[VectorHyperDual, numeric_type]):
         )
 
     @staticmethod
-    def exp(num: Union[VectorHyperDual, numeric_type]):
+    def exp(
+        num: Union[VectorHyperDual, numeric_type]
+    ) -> Union[VectorHyperDual, numeric_type]:
         """Raise e to the power of num"""
 
         return VectorHyperDual.apply_operation(
@@ -439,7 +447,7 @@ def exp(num: Union[VectorHyperDual, numeric_type]):
     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):
@@ -460,7 +468,9 @@ def pow(
                 * math.pow(x0, power - 2),
             )
 
-        elif isinstance(power, VectorHyperDual):
+        elif isinstance(power, VectorHyperDual) and isinstance(
+            num, (numeric, VectorHyperDual)
+        ):
             if (isinstance(num, numeric) and num) < 0 or (
                 isinstance(num, VectorHyperDual) and num.value < 0
             ):
@@ -475,7 +485,9 @@ def pow(
             raise TypeError("Unknown type for exponentiation")
 
     @staticmethod
-    def log(num: Union[VectorHyperDual, numeric_type]):
+    def log(
+        num: Union[VectorHyperDual, numeric_type]
+    ) -> Union[VectorHyperDual, numeric_type]:
         """Natural logarithm"""
 
         if isinstance(num, numeric):
@@ -491,7 +503,9 @@ def log(num: Union[VectorHyperDual, numeric_type]):
         )
 
     @staticmethod
-    def acos(num: Union[VectorHyperDual, numeric_type]):
+    def acos(
+        num: Union[VectorHyperDual, numeric_type]
+    ) -> Union[VectorHyperDual, numeric_type]:
         """Calculate the arccosine of a hyperdual number"""
 
         if isinstance(num, VectorHyperDual):
@@ -508,7 +522,9 @@ def acos(num: Union[VectorHyperDual, numeric_type]):
         )
 
     @staticmethod
-    def atan(num: Union[VectorHyperDual, numeric_type]):
+    def atan(
+        num: Union[VectorHyperDual, numeric_type]
+    ) -> Union[VectorHyperDual, numeric_type]:
         """Calculate the arctangent of a hyperdual number"""
 
         return VectorHyperDual.apply_operation(
@@ -522,7 +538,7 @@ def atan(num: Union[VectorHyperDual, numeric_type]):
     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)
@@ -598,9 +614,11 @@ def norm(self) -> "VectorHyperDual":
         Returns:
             (VectorHyperDual): A scalar number (with derivatives)
         """
-        return DifferentiableMath.sqrt(
+        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"
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 8f14f914d..63e73124d 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -399,7 +399,9 @@ def _evaluate(
         v1 = u_2.cross(u_3)
         v2 = u_1.cross(u_2)
         v3 = u_1 * norm_u2
-        return DifferentiableMath.atan2(v3.dot(v1), v2.dot(v1))
+        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})"

From 10c563cb0fea3228362fb31b9571b444f6849a24 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 2 Dec 2023 11:44:43 +0000
Subject: [PATCH 31/37] test update

---
 autode/opt/coordinates/autodiff.py   |  5 ++--
 autode/opt/coordinates/primitives.py |  2 +-
 tests/test_opt/test_autodiff.py      | 36 ++++++++++++++++++++++++++++
 3 files changed, 40 insertions(+), 3 deletions(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index 3d3ed6e6b..ce46d2493 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -471,7 +471,7 @@ def pow(
         elif isinstance(power, VectorHyperDual) and isinstance(
             num, (numeric, VectorHyperDual)
         ):
-            if (isinstance(num, numeric) and num) < 0 or (
+            if (isinstance(num, numeric) and num < 0) or (
                 isinstance(num, VectorHyperDual) and num.value < 0
             ):
                 raise AssertionError(
@@ -578,10 +578,11 @@ def __init__(self, items: Sequence["VectorHyperDual"]):
         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 = list(items)
+        self._data = items
 
     @staticmethod
     def _check_same_type(other) -> None:
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 63e73124d..91a5e5055 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -7,11 +7,11 @@
     DifferentiableMath,
     DifferentiableVector3D,
     DerivativeOrder,
+    VectorHyperDual,
 )
 
 if TYPE_CHECKING:
     from autode.opt.coordinates import CartesianCoordinates, CartesianComponent
-    from autode.opt.coordinates.autodiff import VectorHyperDual
 
 
 def _get_3d_vecs_from_atom_idxs(
diff --git a/tests/test_opt/test_autodiff.py b/tests/test_opt/test_autodiff.py
index 56006b099..546c8171b 100644
--- a/tests/test_opt/test_autodiff.py
+++ b/tests/test_opt/test_autodiff.py
@@ -6,6 +6,7 @@
     DifferentiableMath,
     VectorHyperDual,
     DerivativeOrder,
+    DifferentiableVector3D,
 )
 
 
@@ -74,6 +75,20 @@ def test_autodiff_exponential_func():
     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)
@@ -151,3 +166,24 @@ def test_derivative_not_available():
     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

From ffaf5b62d69a745f6920beb8f242d87e50415c60 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Sat, 2 Dec 2023 11:56:08 +0000
Subject: [PATCH 32/37] add more tests

---
 autode/opt/coordinates/autodiff.py | 11 ++++++++++-
 tests/test_opt/test_autodiff.py    |  1 +
 2 files changed, 11 insertions(+), 1 deletion(-)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/autodiff.py
index ce46d2493..45f2afdc0 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/autodiff.py
@@ -120,6 +120,15 @@ def _init_deriv_arrays(
         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"""
@@ -195,7 +204,7 @@ def from_variable(
         second_der = None
         n = len(all_symbols)
         idx = list(all_symbols).index(symbol)
-        assert isinstance(order, DerivativeOrder)
+        order = DerivativeOrder(order)
 
         if order == DerivativeOrder.first or order == DerivativeOrder.second:
             first_der = np.zeros(shape=n, dtype=float)
diff --git a/tests/test_opt/test_autodiff.py b/tests/test_opt/test_autodiff.py
index 546c8171b..9120b2056 100644
--- a/tests/test_opt/test_autodiff.py
+++ b/tests/test_opt/test_autodiff.py
@@ -14,6 +14,7 @@ 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

From b2dcc8038006fb69b6c2764be0892b8cd7fee288 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Tue, 5 Dec 2023 10:53:34 +0000
Subject: [PATCH 33/37] update changelog

---
 doc/changelog.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/changelog.rst b/doc/changelog.rst
index 460ed3494..c525d161d 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -7,7 +7,7 @@ Changelog
 
 Functionality improvements
 **************************
-- ...
+- Replaces hand-coded derivatives for primitive internal coordinates with automatic differentiation
 
 Bug Fixes
 *********

From 31d3d299b67ad5bb80bfe202d99fd200d230ae27 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Tue, 5 Dec 2023 10:53:51 +0000
Subject: [PATCH 34/37] update changelog

---
 doc/changelog.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/changelog.rst b/doc/changelog.rst
index c525d161d..a56a51a0a 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -7,7 +7,7 @@ Changelog
 
 Functionality improvements
 **************************
-- Replaces hand-coded derivatives for primitive internal coordinates with automatic differentiation
+- Replaces hard-coded derivatives for primitive internal coordinates with automatic differentiation
 
 Bug Fixes
 *********

From 304d46e27b0db4b4628dbc4ba9cda3972145ec7e Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Tue, 5 Dec 2023 10:55:51 +0000
Subject: [PATCH 35/37] make autodiff private module

---
 autode/opt/coordinates/{autodiff.py => _autodiff.py} | 2 +-
 autode/opt/coordinates/primitives.py                 | 2 +-
 tests/test_opt/test_autodiff.py                      | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)
 rename autode/opt/coordinates/{autodiff.py => _autodiff.py} (99%)

diff --git a/autode/opt/coordinates/autodiff.py b/autode/opt/coordinates/_autodiff.py
similarity index 99%
rename from autode/opt/coordinates/autodiff.py
rename to autode/opt/coordinates/_autodiff.py
index 45f2afdc0..0fa48d437 100644
--- a/autode/opt/coordinates/autodiff.py
+++ b/autode/opt/coordinates/_autodiff.py
@@ -2,7 +2,7 @@
 Automatic differentiation routines in pure Python
 
 References:
-    1. P. Rehner, G. Bauer, Front. Chem. Eng., 2021, 3, 758090
+[1] P. Rehner, G. Bauer, Front. Chem. Eng., 2021, 3, 758090
 """
 from typing import Union, Callable, Sequence, Optional
 from enum import Enum
diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py
index 91a5e5055..764c322b0 100644
--- a/autode/opt/coordinates/primitives.py
+++ b/autode/opt/coordinates/primitives.py
@@ -2,7 +2,7 @@
 
 from abc import ABC, abstractmethod
 from typing import Tuple, TYPE_CHECKING, List
-from autode.opt.coordinates.autodiff import (
+from autode.opt.coordinates._autodiff import (
     get_differentiable_vars,
     DifferentiableMath,
     DifferentiableVector3D,
diff --git a/tests/test_opt/test_autodiff.py b/tests/test_opt/test_autodiff.py
index 9120b2056..eec7065d1 100644
--- a/tests/test_opt/test_autodiff.py
+++ b/tests/test_opt/test_autodiff.py
@@ -1,7 +1,7 @@
 import math
 import numpy as np
 import pytest
-from autode.opt.coordinates.autodiff import (
+from autode.opt.coordinates._autodiff import (
     get_differentiable_vars,
     DifferentiableMath,
     VectorHyperDual,

From aa9238035e1c352da4f5bdddbd015083430d7769 Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Tue, 5 Dec 2023 14:21:16 +0000
Subject: [PATCH 36/37] bug fix dihedral

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

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
 

From 61f93f01a77935a3ddd422930305dbd54529213a Mon Sep 17 00:00:00 2001
From: Shoubhik Maiti <shoubhikrajmaiti@gmail.com>
Date: Tue, 5 Dec 2023 14:21:58 +0000
Subject: [PATCH 37/37] bug fix dihedral

---
 doc/changelog.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/changelog.rst b/doc/changelog.rst
index a56a51a0a..1d2b97732 100644
--- a/doc/changelog.rst
+++ b/doc/changelog.rst
@@ -11,7 +11,7 @@ Functionality improvements
 
 Bug Fixes
 *********
-- ...
+- Fixes triangular rings being incorrectly treated as dihedral angles
 
 
 1.4.1