From c4a02113f3f8e86e26804385ab59631787406c4c Mon Sep 17 00:00:00 2001 From: tom Date: Sun, 3 Dec 2023 15:52:35 +0000 Subject: [PATCH 1/5] Template --- doc/changelog.rst | 12 ++++++++++++ setup.py | 2 +- 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/doc/changelog.rst b/doc/changelog.rst index 7d4c49f61..460ed3494 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -1,6 +1,18 @@ Changelog ========= +1.4.2 +------ +------- + +Functionality improvements +************************** +- ... + +Bug Fixes +********* +- ... + 1.4.1 ------ diff --git a/setup.py b/setup.py index 0e243e7b4..300adb952 100644 --- a/setup.py +++ b/setup.py @@ -32,7 +32,7 @@ setup( name="autode", - version="1.4.1", + version="1.4.2", python_requires=">3.7", packages=[ "autode", From 3199ba5dee0429ad1e55af0d1b42bff354e01af0 Mon Sep 17 00:00:00 2001 From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com> Date: Wed, 6 Dec 2023 19:14:03 +0000 Subject: [PATCH 2/5] Automatic differentiation for primitives (#311) * initial commit autodiff * second commit autodiff * third commit autodiff * autodiff update * autodiff update * doc update * autodiff update * test update * type checking and B matrix * fix some tests * fix some tests * add tests * fix test * fix tests 2 * fix tests 2 * doc update * add test * minor refactor autodiff * add more tests * remove legacy derivatives and add test * pr suggestions 1 * Use enum for derivative order * fix tests * doc update and use lambda func * differentiable vector class * doc update * doc update * coverage update * coverage update * doc update * test update * add more tests * update changelog * update changelog * make autodiff private module * bug fix dihedral * bug fix dihedral --- autode/opt/coordinates/_autodiff.py | 714 +++++++++++++++++++++++++++ autode/opt/coordinates/dic.py | 48 +- autode/opt/coordinates/internals.py | 21 +- autode/opt/coordinates/primitives.py | 374 +++++++------- autode/opt/optimisers/crfo.py | 4 + doc/changelog.rst | 4 +- tests/test_opt/test_autodiff.py | 190 +++++++ tests/test_opt/test_coordiantes.py | 178 ++++--- 8 files changed, 1224 insertions(+), 309 deletions(-) create mode 100644 autode/opt/coordinates/_autodiff.py create mode 100644 tests/test_opt/test_autodiff.py diff --git a/autode/opt/coordinates/_autodiff.py b/autode/opt/coordinates/_autodiff.py new file mode 100644 index 000000000..0fa48d437 --- /dev/null +++ b/autode/opt/coordinates/_autodiff.py @@ -0,0 +1,714 @@ +""" +Automatic differentiation routines in pure Python + +References: +[1] P. Rehner, G. Bauer, Front. Chem. Eng., 2021, 3, 758090 +""" +from typing import Union, Callable, Sequence, Optional +from enum import Enum +from copy import deepcopy +import numpy as np +import math + +numeric = (float, int) +numeric_type = Union[float, int] + + +class DerivativeOrder(Enum): + """Order of derivative""" + + zeroth = 0 + first = 1 + second = 2 + + +def get_differentiable_vars( + values: Sequence[numeric_type], + symbols: Sequence[str], + deriv_order: DerivativeOrder = DerivativeOrder.second, +): + """ + Obtain differentiable variables from a series of numbers + + Args: + values: The values of the variables (numbers) + symbols: List of symbols (strings) of the numbers + deriv_order: Order of differentiation + + Returns: + (list[VectorHyperDual]): A list of hyper dual numbers + """ + assert all(isinstance(sym, str) for sym in symbols) + assert len(symbols) == len(values) + symbols = list(symbols) + + hyperduals = [] + for symbol, value in zip(symbols, values): + var = VectorHyperDual.from_variable( + value, symbol, all_symbols=symbols, order=deriv_order + ) + hyperduals.append(var) + + return hyperduals + + +class VectorHyperDual: + """ + Hyper-dual numbers with vector infinitesimals upto the + second order (i.e., upto second partial derivatives) + """ + + def __init__( + self, + value: numeric_type, + symbols: Sequence[str], + first_der: Optional[np.ndarray] = None, + second_der: Optional[np.ndarray] = None, + ): + """ + Create a vector hyper dual number, i.e. a scalar function + with one or more variables + + Args: + value: The scalar value of the hyper-dual + symbols: A list of unique strings representing the variables + first_der: 1D array of first derivatives against variables + second_der: 2D matrix of second derivatives against variables + """ + assert isinstance(value, numeric) + self._val = float(value) + + assert all(isinstance(sym, str) for sym in symbols) + if len(set(symbols)) != len(list(symbols)): + raise RuntimeError("Symbols must be unique!") + self._symbols = tuple(symbols) + + # load the derivatives with sanity checks + self._first_der: Optional[np.ndarray] = None + self._second_der: Optional[np.ndarray] = None + self._order = DerivativeOrder.zeroth + self._init_deriv_arrays(first_der, second_der) + + def _init_deriv_arrays( + self, first_der: Optional[np.ndarray], second_der: Optional[np.ndarray] + ) -> None: + """ + Initialise the derivative matrices, checking they have the + correct shape, and set the derivative order for this + hyper-dual number + """ + if first_der is None: + return None + assert isinstance(first_der, np.ndarray) + first_der = first_der.flatten() + if not first_der.shape == (self.n_vars,): + raise ValueError( + f"Number of symbols ({self.n_vars}) does not match with" + f" shape of derivative array {first_der.shape}" + ) + self._first_der = first_der.astype(float) + self._order = DerivativeOrder.first + + if second_der is None: + return None + assert isinstance(second_der, np.ndarray) + if not second_der.shape == (self.n_vars, self.n_vars): + raise ValueError( + f"Number of symbols ({self.n_vars}) does not match with" + f" shape of second derivative matrix {first_der.shape}" + ) + self._second_der = second_der.astype(float) + self._order = DerivativeOrder.second + + def __repr__(self): + rstring = f"HyperDual({self.value}" + if self._order in [DerivativeOrder.first, DerivativeOrder.second]: + rstring += f", f'[{self.n_vars}]" + if self._order == DerivativeOrder.second: + rstring += f', f"[{self.n_vars}, {self.n_vars}]' + rstring += ")" + return rstring + + @property + def n_vars(self) -> int: + """Number of variables in this hyper-dual""" + return len(self._symbols) + + def copy(self) -> "VectorHyperDual": + return deepcopy(self) + + def _check_compatible(self, other: "VectorHyperDual") -> None: + """ + Check the compatibility of two VectorHyperDual numbers for + any operation that involves the two. + + Args: + other (VectorHyperDual): + + Raises: + (ValueError): If they are incompatible + """ + if self.n_vars != other.n_vars: + raise ValueError( + "Incompatible number of differentiable variables, " + "cannot perform operation" + ) + if self._symbols != other._symbols: + raise ValueError( + "The differentiable variable symbols do not match, " + "cannot perform operation!" + ) + if self._order != other._order: + raise ValueError("The order of derivative do not match!") + return None + + @property + def value(self) -> float: + """Return the value of the hyper-dual number""" + return self._val + + @value.setter + def value(self, value: float): + assert isinstance(value, numeric) + self._val = float(value) + + @classmethod + def from_variable( + cls, + value: float, + symbol: str, + all_symbols: Sequence[str], + order: DerivativeOrder, + ): + """ + Create a hyper-dual number from one variable, requires + list of symbols and the symbol of this specific variable. + Essentially a variable x can be considered as a scalar function + of a list of variables - 1 * x + 0 * y + 0 * z + ... + + Args: + value: The value of the variable (will be converted to float) + symbol: The symbol of the current variable, must be in all_symbols + all_symbols: List of strings indicating all required variables + order: The order of differentiation to consider + + Returns: + (VectorHyperDual): The hyper-dual representing the variable + """ + assert all(isinstance(sym, str) for sym in all_symbols) + assert isinstance(symbol, str) + assert symbol in all_symbols + + val = float(value) + first_der = None + second_der = None + n = len(all_symbols) + idx = list(all_symbols).index(symbol) + order = DerivativeOrder(order) + + if order == DerivativeOrder.first or order == DerivativeOrder.second: + first_der = np.zeros(shape=n, dtype=float) + first_der[idx] = 1.0 + if order == DerivativeOrder.second: + second_der = np.zeros(shape=(n, n), dtype=float) + + return VectorHyperDual(val, all_symbols, first_der, second_der) + + def differentiate_wrt( + self, + symbol1: str, + symbol2: Union[str, None] = None, + ) -> Optional[float]: + """ + Derivative of this hyper-dual number (scalar function) against one + or two variable(s) identified by their string(s). + + Args: + symbol1 (str): + symbol2 (str|None): + + Returns: + (float|None): The derivative value, or None if not available + """ + assert isinstance(symbol1, str) + if symbol1 not in self._symbols: + return None + + if self._order == DerivativeOrder.zeroth: + return None + + idx_1 = self._symbols.index(symbol1) + assert self._first_der is not None + if symbol2 is None: + return self._first_der[idx_1] + + assert isinstance(symbol2, str) + if symbol2 not in self._symbols: + return None + idx_2 = self._symbols.index(symbol2) + # check if second derivs are available + if self._order == DerivativeOrder.first: + return None + assert self._second_der is not None + return self._second_der[idx_1, idx_2] + + def __add__( + self, other: Union["VectorHyperDual", numeric_type] + ) -> "VectorHyperDual": + """Adding a hyper dual number""" + + if isinstance(other, numeric): + new = self.copy() + new._val += float(other) + return new + + # add to another dual number + elif isinstance(other, VectorHyperDual): + self._check_compatible(other) + + val = self._val + other._val + if self._order == DerivativeOrder.zeroth: + return VectorHyperDual(val, self._symbols) + + assert self._first_der is not None + assert other._first_der is not None + first_der = self._first_der + other._first_der + if self._order == DerivativeOrder.first: + return VectorHyperDual(val, self._symbols, first_der) + + assert self._second_der is not None + assert other._second_der is not None + second_der = self._second_der + other._second_der + return VectorHyperDual(val, self._symbols, first_der, second_der) + + else: + raise TypeError("Unknown type for addition") + + def __radd__(self, other): + """Addition is commutative""" + return self.__add__(other) + + def __neg__(self) -> "VectorHyperDual": + """Unary negative operation""" + new = self.copy() + new._val = -new._val + if self._order == DerivativeOrder.first: + assert new._first_der is not None + new._first_der = -new._first_der + elif self._order == DerivativeOrder.second: + assert new._first_der is not None + assert new._second_der is not None + new._first_der = -new._first_der + new._second_der = -new._second_der + return new + + def __sub__(self, other): + """Subtraction of hyper dual numbers""" + return self.__add__(-other) + + def __rsub__(self, other): + """Reverse subtraction""" + return other + (-self) + + def __mul__(self, other) -> "VectorHyperDual": + """Multiply a hyper dual number with float or another hyper dual""" + if isinstance(other, numeric): + new = self.copy() + new._val *= float(other) + return new + + # Product rule for derivatives, Eqn (24) in ref. [1] + elif isinstance(other, VectorHyperDual): + self._check_compatible(other) + + val = self._val * other._val + if self._order == DerivativeOrder.zeroth: + return VectorHyperDual(val, self._symbols) + + assert self._first_der is not None + assert other._first_der is not None + first_der = ( + self._val * other._first_der + other._val * self._first_der + ) + if self._order == DerivativeOrder.first: + return VectorHyperDual(val, self._symbols, first_der) + + assert self._second_der is not None + assert other._second_der is not None + second_der = ( + self._val * other._second_der + + np.outer(self._first_der, other._first_der) + + np.outer(other._first_der, self._first_der) + + other._val * self._second_der + ) + return VectorHyperDual(val, self._symbols, first_der, second_der) + else: + raise TypeError("Unknown type for multiplication") + + def __rmul__(self, other): + return self.__mul__(other) + + def __truediv__(self, other): + """True division, defined by multiplicative inverse""" + return self.__mul__(DifferentiableMath.pow(other, -1)) + + def __rtruediv__(self, other): + """Reverse true division""" + return DifferentiableMath.pow(self, -1).__mul__(other) + + def __pow__(self, power, modulo=None) -> "VectorHyperDual": + if modulo is not None: + raise NotImplementedError("Modulo inverse is not implemented") + + result = DifferentiableMath.pow(self, power) + assert isinstance(result, VectorHyperDual) + return result + + def __rpow__(self, other): + return DifferentiableMath.pow(other, self) + + @staticmethod + def apply_operation( + num: Union["VectorHyperDual", numeric_type], + operator: Callable[[float], float], + operator_first_deriv: Callable[[float], float], + operator_second_deriv: Callable[[float], float], + ) -> Union["VectorHyperDual", numeric_type]: + """ + Perform an operation on the hyperdual (i.e. apply a scalar function), + also compatible with Python numeric (float/int) types. + + Args: + num: Number that is hyper-dual (or float/int) + operator: Function that returns the value (result) of operation + operator_first_deriv: Should return first derivative of operation + operator_second_deriv: Should return second derivative of operation + + Returns: + (VectorHyperDual|float): The result + """ + # pass through numeric types + if isinstance(num, numeric): + return operator(float(num)) + + assert isinstance(num, VectorHyperDual) + + val = operator(num._val) + + if num._order == DerivativeOrder.zeroth: + return VectorHyperDual(val, num._symbols) + + # Eqn (25) in reference [1] + assert num._first_der is not None + f_dash_x0 = operator_first_deriv(num._val) + first_der = num._first_der * f_dash_x0 + + if num._order == DerivativeOrder.first: + return VectorHyperDual(val, num._symbols, first_der) + + assert num._second_der is not None + second_der = np.outer(num._first_der, num._first_der) * ( + operator_second_deriv(num._val) + ) + second_der += num._second_der * f_dash_x0 + return VectorHyperDual(val, num._symbols, first_der, second_der) + + +class DifferentiableMath: + """ + Class defining math functions that can be used on + hyper dual numbers (i.e. differentiable functions), + as well as standard numeric types (float and int) + """ + + @staticmethod + def sqrt( + num: Union[VectorHyperDual, numeric_type] + ) -> Union[VectorHyperDual, numeric_type]: + """Calculate the square root of a hyperdual number""" + + if isinstance(num, numeric): + assert num > 0 + else: + assert num.value > 0 + + return VectorHyperDual.apply_operation( + num, + operator=lambda x0: math.sqrt(x0), + operator_first_deriv=lambda x0: 1 / (2 * math.sqrt(x0)), + operator_second_deriv=lambda x0: -1 / (4 * math.pow(x0, 3 / 2)), + ) + + @staticmethod + def exp( + num: Union[VectorHyperDual, numeric_type] + ) -> Union[VectorHyperDual, numeric_type]: + """Raise e to the power of num""" + + return VectorHyperDual.apply_operation( + num, + operator=lambda x0: math.exp(x0), + operator_first_deriv=lambda x0: math.exp(x0), + operator_second_deriv=lambda x0: math.exp(x0), + ) + + @staticmethod + def pow( + num: Union[VectorHyperDual, numeric_type], + power: Union[VectorHyperDual, numeric_type], + ) -> Union[VectorHyperDual, numeric_type]: + """Exponentiation of one hyperdual to another""" + + if isinstance(num, numeric) and isinstance(power, numeric): + return math.pow(num, power) + + elif isinstance(num, VectorHyperDual) and isinstance(power, numeric): + if num.value < 0 and isinstance(power, float): + raise AssertionError( + "Math error, can't raise negative number to fractional power" + ) + return VectorHyperDual.apply_operation( + num, + operator=lambda x0: math.pow(x0, power), # type: ignore + operator_first_deriv=lambda x0: power # type: ignore + * math.pow(x0, power - 1), + operator_second_deriv=lambda x0: power # type: ignore + * (power - 1) + * math.pow(x0, power - 2), + ) + + elif isinstance(power, VectorHyperDual) and isinstance( + num, (numeric, VectorHyperDual) + ): + if (isinstance(num, numeric) and num < 0) or ( + isinstance(num, VectorHyperDual) and num.value < 0 + ): + raise AssertionError( + "Only positive numbers can be used with" + " differentiable exponent" + ) + # use identity x^y = e^(y log_x) for x > 0 + return DifferentiableMath.exp(power * DifferentiableMath.log(num)) + + else: + raise TypeError("Unknown type for exponentiation") + + @staticmethod + def log( + num: Union[VectorHyperDual, numeric_type] + ) -> Union[VectorHyperDual, numeric_type]: + """Natural logarithm""" + + if isinstance(num, numeric): + assert num > 0 + else: + assert num.value > 0 + + return VectorHyperDual.apply_operation( + num, + operator=lambda x0: math.log(x0), + operator_first_deriv=lambda x0: 1.0 / x0, + operator_second_deriv=lambda x0: -1.0 / (x0**2), + ) + + @staticmethod + def acos( + num: Union[VectorHyperDual, numeric_type] + ) -> Union[VectorHyperDual, numeric_type]: + """Calculate the arccosine of a hyperdual number""" + + if isinstance(num, VectorHyperDual): + assert -1 < num.value < 1 + else: + assert -1 < num < 1 + + return VectorHyperDual.apply_operation( + num, + operator=lambda x0: math.acos(x0), + operator_first_deriv=lambda x0: -1 / math.sqrt(1 - x0**2), + operator_second_deriv=lambda x0: -x0 + / math.pow(1 - x0**2, 3 / 2), + ) + + @staticmethod + def atan( + num: Union[VectorHyperDual, numeric_type] + ) -> Union[VectorHyperDual, numeric_type]: + """Calculate the arctangent of a hyperdual number""" + + return VectorHyperDual.apply_operation( + num, + operator=lambda x0: math.atan(x0), + operator_first_deriv=lambda x0: 1 / (1 + x0**2), + operator_second_deriv=lambda x0: (-2 * x0) / (x0**2 + 1) ** 2, + ) + + @staticmethod + def atan2( + num_y: Union[VectorHyperDual, numeric_type], + num_x: Union[VectorHyperDual, numeric_type], + ) -> Union[VectorHyperDual, numeric_type]: + """Calculate the arctan2 of two hyper dual numbers""" + if isinstance(num_y, numeric) and isinstance(num_x, numeric): + return math.atan2(num_y, num_x) + + # https://en.wikipedia.org/wiki/Atan2 four overlapping half-planes + def atan2_derivs_x_not_0(y, x): + return DifferentiableMath.atan(y / x) + + def atan2_derivs_x_close_0(y, x): + return -DifferentiableMath.atan(x / y) + + x_val = float(num_x) if isinstance(num_x, numeric) else num_x.value + y_val = float(num_y) if isinstance(num_y, numeric) else num_y.value + res_val = math.atan2(y_val, x_val) + + # when atan2(y,x)->pi/2, x->0 or y/x->inf, use other formula for derivs + if math.isclose(abs(res_val), math.pi / 2, abs_tol=0.1): + res = atan2_derivs_x_close_0(num_y, num_x) + res.value = res_val + return res + else: + res = atan2_derivs_x_not_0(num_y, num_x) + res.value = res_val + return res + + +class DifferentiableVector3D: + """ + Convenience class to represent a 3D vector of differentiable + hyper-dual numbers + """ + + def __init__(self, items: Sequence["VectorHyperDual"]): + """ + Initialise the 3D vector from a list of 3 hyperdual numbers + + Args: + items: A list of 3 hyper-dual numbers + """ + items = list(items) + if len(items) != 3: + raise ValueError("A 3D vector must have only 3 components") + assert all(isinstance(item, VectorHyperDual) for item in items) + self._data = items + + @staticmethod + def _check_same_type(other) -> None: + """Check that another object is also a 3D differentiable vector""" + if not isinstance(other, DifferentiableVector3D): + raise ValueError("Operation must be done with another 3D vector!") + return None + + def dot(self, other: "DifferentiableVector3D") -> "VectorHyperDual": + """ + Dot product of two 3D vectors + + Args: + other (DifferentiableVector3D): + + Returns: + (VectorHyperDual): A scalar number (with derivatives) + """ + self._check_same_type(other) + dot = 0 + for k in range(3): + dot = dot + self._data[k] * other._data[k] + assert isinstance(dot, VectorHyperDual) + return dot + + def norm(self) -> "VectorHyperDual": + """ + Euclidean (l2) norm of this 3D vector + + Returns: + (VectorHyperDual): A scalar number (with derivatives) + """ + norm = DifferentiableMath.sqrt( + self._data[0] ** 2 + self._data[1] ** 2 + self._data[2] ** 2 + ) + assert isinstance(norm, VectorHyperDual) + return norm + + def __add__( + self, other: "DifferentiableVector3D" + ) -> "DifferentiableVector3D": + """ + Vector addition in 3D, returns a vector + + Args: + other (DifferentiableVector3D): + + Returns: + (DifferentiableVector3D): + """ + self._check_same_type(other) + return DifferentiableVector3D( + [self._data[k] + other._data[k] for k in range(3)] + ) + + def __neg__(self) -> "DifferentiableVector3D": + """ + Unary negation of a vector, returns another vector + + Returns: + (DifferentiableVector3D): + """ + return DifferentiableVector3D([-self._data[k] for k in range(3)]) + + def __sub__(self, other) -> "DifferentiableVector3D": + """ + Vector subtraction in 3D, defined in terms of addition + and negation + + Args: + other (DifferentiableVector3D): + + Returns: + (DifferentiableVector3D): + """ + return self.__add__(-other) + + def __mul__( + self, other: Union[VectorHyperDual, numeric_type] + ) -> "DifferentiableVector3D": + """ + Multiplication of a 3D vector with a scalar + + Args: + other (VectorHyperDual|float|int): + + Returns: + (DifferentiableVector3D): + """ + assert isinstance(other, numeric) or isinstance(other, VectorHyperDual) + return DifferentiableVector3D( + [self._data[k] * other for k in range(3)] + ) + + def __rmul__(self, other): + """Multiplication of scalar and vector is commutative""" + return self.__mul__(other) + + def cross( + self, other: "DifferentiableVector3D" + ) -> "DifferentiableVector3D": + """ + Cross-product of two 3D vectors, produces another vector + + Args: + other (DifferentiableVector3D): + + Returns: + (DifferentiableVector3D): + """ + return DifferentiableVector3D( + [ + self._data[1] * other._data[2] + - self._data[2] * other._data[1], + self._data[2] * other._data[0] + - self._data[0] * other._data[2], + self._data[0] * other._data[1] + - self._data[1] * other._data[0], + ] + ) diff --git a/autode/opt/coordinates/dic.py b/autode/opt/coordinates/dic.py index b857f2adb..2d5277854 100644 --- a/autode/opt/coordinates/dic.py +++ b/autode/opt/coordinates/dic.py @@ -216,34 +216,42 @@ def iadd(self, value: np.ndarray) -> "OptCoordinates": q_init = self.primitives(x_k) x_1 = self.to("cartesian") + np.matmul(self.B_T_inv, value) + success = False for i in range(1, _max_back_transform_iterations + 1): - x_k = x_k + np.matmul(self.B_T_inv, (s_new - s_k)) + try: + x_k = x_k + np.matmul(self.B_T_inv, (s_new - s_k)) - # Rebuild the primitives & DIC from the back-transformed Cartesians - q_k = self.primitives.close_to(x_k, q_init) - s_k = np.matmul(self.U.T, q_k) - self.B = np.matmul(self.U.T, self.primitives.B) - self.B_T_inv = np.linalg.pinv(self.B) + # Rebuild the primitives & DIC from the back-transformed Cartesians + q_k = self.primitives.close_to(x_k, q_init) + s_k = np.matmul(self.U.T, q_k) + self.B = np.matmul(self.U.T, self.primitives.B) + self.B_T_inv = np.linalg.pinv(self.B) - rms_s = np.sqrt(np.mean(np.square(s_k - s_new))) + rms_s = np.sqrt(np.mean(np.square(s_k - s_new))) + + # for ill-conditioned primitives, there might be math error + except ArithmeticError: + break if rms_s < 1e-10: - logger.info( - f"DIC transformation converged in {i} cycle(s) " - f"in {time() - start_time:.4f} s" - ) + success = True break - if i == _max_back_transform_iterations: - logger.warning( - f"Failed to transform in {i} cycles. " - f"Final RMS(s) = {rms_s:.8f}" + if success: + logger.info( + f"DIC transformation converged in {i} cycle(s) " + f"in {time() - start_time:.4f} s" + ) + else: + logger.warning( + f"Failed to transform in {i} cycles. " + f"Final RMS(s) = {rms_s:.8f}" + ) + x_k = x_1 + if not self.allow_unconverged_back_transform: + raise CoordinateTransformFailed( + "DIC->Cart iterative back-transform did not converge" ) - x_k = x_1 - if not self.allow_unconverged_back_transform: - raise CoordinateTransformFailed( - "DIC->Cart iterative back-transform did not converge" - ) s_k = np.matmul(self.U.T, self.primitives(x_k)) self.B = np.matmul(self.U.T, self.primitives.B) diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py index 675ca1f1c..9e7462e16 100644 --- a/autode/opt/coordinates/internals.py +++ b/autode/opt/coordinates/internals.py @@ -82,6 +82,11 @@ def __init__(self, *args: Any): f"from {args}. Must be primitive internals" ) + def append(self, item: Primitive) -> None: + """Append an item to this set of primitives""" + assert isinstance(item, Primitive), "Must be a Primitive type!" + super().append(item) + @property def B(self) -> np.ndarray: """Wilson B matrix""" @@ -173,22 +178,12 @@ def _calc_B(self, x: np.ndarray) -> None: "primitive internal coordinates" ) - cart_coords = x.reshape((-1, 3)) + cart_coords = x.ravel() - n_atoms, _ = cart_coords.shape - B = np.zeros(shape=(len(self), 3 * n_atoms)) + B = np.zeros(shape=(len(self), len(cart_coords))) for i, primitive in enumerate(self): - for j in range(n_atoms): - B[i, 3 * j + 0] = primitive.derivative( - j, CartesianComponent.x, x=cart_coords - ) - B[i, 3 * j + 1] = primitive.derivative( - j, CartesianComponent.y, x=cart_coords - ) - B[i, 3 * j + 2] = primitive.derivative( - j, CartesianComponent.z, x=cart_coords - ) + B[i] = primitive.derivative(x=cart_coords) self._B = B return None diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py index 22bfabb86..764c322b0 100644 --- a/autode/opt/coordinates/primitives.py +++ b/autode/opt/coordinates/primitives.py @@ -1,12 +1,57 @@ import numpy as np from abc import ABC, abstractmethod -from typing import Tuple, TYPE_CHECKING +from typing import Tuple, TYPE_CHECKING, List +from autode.opt.coordinates._autodiff import ( + get_differentiable_vars, + DifferentiableMath, + DifferentiableVector3D, + DerivativeOrder, + VectorHyperDual, +) if TYPE_CHECKING: from autode.opt.coordinates import CartesianCoordinates, CartesianComponent +def _get_3d_vecs_from_atom_idxs( + *args: int, + x: "CartesianCoordinates", + deriv_order: DerivativeOrder, +) -> List[DifferentiableVector3D]: + """ + Obtain differentiable 3D vectors from the Cartesian components + of each atom, given by atomic indices in order. The symbols are + strings denoting their position in the flat Cartesian coordinate. + + Args: + *args: Integers denoting the atom positions + x: Cartesian coordinate array + deriv_order: Order of derivatives for initialising variables + + Returns: + (list[VectorHyperDual]): A list of differentiable variables + """ + assert all(isinstance(idx, int) and idx >= 0 for idx in args) + # get positions in the flat Cartesian array + _x = x.ravel() + cart_idxs = [] + for atom_idx in args: + for k in range(3): + cart_idxs.append(3 * atom_idx + k) + variables = get_differentiable_vars( + values=[_x[idx] for idx in cart_idxs], + symbols=[str(idx) for idx in cart_idxs], + deriv_order=deriv_order, + ) + atom_vecs = [] + for pos_idx in range(len(args)): + atom_vecs.append( + DifferentiableVector3D(variables[pos_idx * 3 : pos_idx * 3 + 3]) + ) + return atom_vecs + + class Primitive(ABC): """Primitive internal coordinate""" @@ -17,18 +62,35 @@ def __init__(self, *atom_indexes: int): self._atom_indexes = atom_indexes @abstractmethod + def _evaluate( + self, x: "CartesianCoordinates", deriv_order: DerivativeOrder + ): + """ + The function that performs the main evaluation of the PIC, + and optionally returns derivative or second derivatives. + The returned hyper-dual must have the proper cartesian idxs + set. + + Args: + x: Cartesian coordinates + deriv_order: The order of derivatives requested - 0, 1 or 2 + + Returns: + (VectorHyperDual): The result, optionally containing derivatives + """ + def __call__(self, x: "CartesianCoordinates") -> float: """Return the value of this PIC given a set of cartesian coordinates""" + _x = x.ravel() + res = self._evaluate(_x, deriv_order=DerivativeOrder.zeroth) + return res.value - @abstractmethod def derivative( self, - i: int, - component: "CartesianComponent", x: "CartesianCoordinates", - ) -> float: + ) -> np.ndarray: r""" - Calculate the derivative with respect to a cartesian coordinate + Calculate the derivatives with respect to cartesian coordinates .. math:: @@ -40,17 +102,56 @@ def derivative( ----------------------------------------------------------------------- Arguments: - i: Cartesian index to take the derivative with respect to. [0-N), - for N atoms - component: Cartesian component (x, y, z) to take the derivative - with respect to + x: Cartesian coordinate array of shape (N, ) - x: Cartesian coordinates + Returns: + (np.ndarray): Derivative array of shape (N, ) + """ + _x = x.ravel() + res = self._evaluate(_x, deriv_order=DerivativeOrder.first) + derivs = np.zeros_like(_x, dtype=float) + for i in range(_x.shape[0]): + dqdx_i = res.differentiate_wrt(str(i)) + if dqdx_i is not None: + derivs[i] = dqdx_i + + return derivs + + def second_derivative( + self, + x: "CartesianCoordinates", + ) -> np.ndarray: + r""" + Calculate the second derivatives with respect to cartesian coordinates + + .. math:: + + \frac{d^2 q} + {d\boldsymbol{X}_{i, k}^2} {\Bigg\rvert}_{X=X0} + + where :math:`q` is the primitive coordinate and :math:`\boldsymbol{X}` + are the cartesian coordinates. + + ----------------------------------------------------------------------- + Arguments: + + x: Cartesian coordinate array of shape (N, ) Returns: - (float): Derivative + (np.ndarray): Second derivative matrix of shape (N, N) """ + _x = x.ravel() + x_n = _x.shape[0] + res = self._evaluate(_x, deriv_order=DerivativeOrder.second) + derivs = np.zeros(shape=(x_n, x_n), dtype=float) + for i in range(x_n): + for j in range(x_n): + d2q_dx2_ij = res.differentiate_wrt(str(i), str(j)) + if d2q_dx2_ij is not None: + derivs[i, j] = d2q_dx2_ij + + return derivs @abstractmethod def __eq__(self, other): @@ -121,8 +222,8 @@ def __eq__(self, other) -> bool: class PrimitiveInverseDistance(_DistanceFunction): - """ - Inverse distance between to atoms: + r""" + Inverse distance between two atoms: .. math:: @@ -130,40 +231,21 @@ class PrimitiveInverseDistance(_DistanceFunction): {|\boldsymbol{X}_i - \boldsymbol{X}_j|} """ - def derivative( - self, - i: int, - component: "CartesianComponent", - x: "CartesianCoordinates", - ) -> float: - """ - Derivative with respect to Cartesian displacement - - ----------------------------------------------------------------------- - See Also: - :py:meth:`Primitive.derivative ` - """ - - _x = x.reshape((-1, 3)) - k = int(component) - - if i != self.i and i != self.j: - return 0 # Atom does not form part of this distance - - elif i == self.i: - return -(_x[i, k] - _x[self.j, k]) * self(x) ** 3 - - else: # i == self.idx_j: - return (_x[self.i, k] - _x[self.j, k]) * self(x) ** 3 - - def __call__(self, x: "CartesianCoordinates") -> float: + def _evaluate( + self, x: "CartesianCoordinates", deriv_order: DerivativeOrder + ) -> "VectorHyperDual": """1 / |x_i - x_j|""" - _x = x.reshape((-1, 3)) - return 1.0 / np.linalg.norm(_x[self.i] - _x[self.j]) + vec_i, vec_j = _get_3d_vecs_from_atom_idxs( + self.i, self.j, x=x, deriv_order=deriv_order + ) + return 1.0 / (vec_i - vec_j).norm() + + def __repr__(self): + return f"InverseDistance({self.i}-{self.j})" class PrimitiveDistance(_DistanceFunction): - """ + r""" Distance between two atoms: .. math:: @@ -171,33 +253,14 @@ class PrimitiveDistance(_DistanceFunction): q = |\boldsymbol{X}_i - \boldsymbol{X}_j| """ - def derivative( - self, - i: int, - component: "CartesianComponent", - x: "CartesianCoordinates", - ) -> float: - """ - Derivative with respect to Cartesian displacement - - ----------------------------------------------------------------------- - See Also: - :py:meth:`Primitive.derivative ` - """ - _x = x.reshape((-1, 3)) - k = int(component) - - if i != self.i and i != self.j: - return 0 # Atom does not form part of this distance - - val = (_x[self.i, k] - _x[self.j, k]) / self(x) - - return val if i == self.i else -val - - def __call__(self, x: "CartesianCoordinates") -> float: + def _evaluate( + self, x: "CartesianCoordinates", deriv_order: DerivativeOrder + ) -> "VectorHyperDual": """|x_i - x_j|""" - _x = x.reshape((-1, 3)) - return np.linalg.norm(_x[self.i] - _x[self.j]) + vec_i, vec_j = _get_3d_vecs_from_atom_idxs( + self.i, self.j, x=x, deriv_order=deriv_order + ) + return (vec_i - vec_j).norm() def __repr__(self): return f"Distance({self.i}-{self.j})" @@ -230,13 +293,18 @@ def __repr__(self): class PrimitiveBondAngle(Primitive): + """ + Bond angle between three atoms, calculated with the + arccosine of the normalised dot product + """ + def __init__(self, o: int, m: int, n: int): """Bond angle m-o-n""" super().__init__(o, m, n) - self.o = o - self.m = m - self.n = n + self.o = int(o) + self.m = int(m) + self.n = int(n) def __eq__(self, other) -> bool: """Equality of two distance functions""" @@ -247,58 +315,16 @@ def __eq__(self, other) -> bool: and other._ordered_idxs == self._ordered_idxs ) - def __call__(self, x: "CartesianCoordinates") -> float: - _x = x.reshape((-1, 3)) - u = _x[self.m, :] - _x[self.o, :] - v = _x[self.n, :] - _x[self.o, :] - - theta = np.arccos(u.dot(v) / (np.linalg.norm(u) * np.linalg.norm(v))) - return theta - - def derivative( - self, - i: int, - component: "CartesianComponent", - x: "CartesianCoordinates", - ) -> float: - if i not in (self.o, self.m, self.n): - return 0.0 - - k = int(component) - - _x = x.reshape((-1, 3)) - u = _x[self.m, :] - _x[self.o, :] - lambda_u = np.linalg.norm(u) - u /= lambda_u - - v = _x[self.n, :] - _x[self.o, :] - lambda_v = np.linalg.norm(v) - v /= lambda_v - - t0, t1 = np.array([1.0, -1.0, 1.0]), np.array([-1.0, 1.0, 1.0]) - - if not np.isclose(np.abs(np.arccos(u.dot(v))), 1.0): - w = np.cross(u, v) - elif not np.isclose( - np.abs(np.arccos(u.dot(t0))), 1.0 - ) and not np.isclose(np.abs(np.arccos(v.dot(t0))), 1.0): - w = np.cross(u, t0) - else: - w = np.cross(u, t1) - - w /= np.linalg.norm(w) - - dqdx = 0.0 - - if i in (self.m, self.o): - sign = 1 if i == self.m else -1 - dqdx += sign * np.cross(u, w)[k] / lambda_u - - if i in (self.n, self.o): - sign = 1 if i == self.n else -1 - dqdx += sign * np.cross(w, v)[k] / lambda_v - - return dqdx + def _evaluate( + self, x: "CartesianCoordinates", deriv_order: DerivativeOrder + ): + """m - o - n angle""" + vec_m, vec_o, vec_n = _get_3d_vecs_from_atom_idxs( + self.m, self.o, self.n, x=x, deriv_order=deriv_order + ) + u = vec_m - vec_o + v = vec_n - vec_o + return DifferentiableMath.acos(u.dot(v) / (u.norm() * v.norm())) def __repr__(self): return f"Angle({self.m}-{self.o}-{self.n})" @@ -344,22 +370,10 @@ def __init__(self, m: int, o: int, p: int, n: int): """Dihedral angle: m-o-p-n""" super().__init__(m, o, p, n) - self.m = m - self.o = o - self.p = p - self.n = n - - def __call__(self, x: "CartesianCoordinates") -> float: - """Value of the dihedral""" - return self._value(x, return_derivative=False) - - def derivative( - self, - i: int, - component: "CartesianComponent", - x: "CartesianCoordinates", - ) -> float: - return self._value(x, i=i, component=component, return_derivative=True) + self.m = int(m) + self.o = int(o) + self.p = int(p) + self.n = int(n) def __eq__(self, other) -> bool: """Equality of two distance functions""" @@ -368,64 +382,26 @@ def __eq__(self, other) -> bool: or self._atom_indexes == tuple(reversed(other._atom_indexes)) ) - def _value(self, x, i=None, component=None, return_derivative=False): - """Evaluate either the value or the derivative. Shared function - to reuse local variables""" - - _x = x.reshape((-1, 3)) - u = _x[self.m, :] - _x[self.o, :] - lambda_u = np.linalg.norm(u) - u /= lambda_u - - v = _x[self.n, :] - _x[self.p, :] - lambda_v = np.linalg.norm(v) - v /= lambda_v - - w = _x[self.p, :] - _x[self.o, :] - lambda_w = np.linalg.norm(w) - w /= lambda_w - - phi_u = np.arccos(u.dot(w)) - phi_v = np.arccos(w.dot(v)) - - if not return_derivative: - v1 = np.cross(u, w) - v2 = np.cross(-w, v) - return -np.arctan2(np.cross(v1, w).dot(v2), (v1.dot(v2))) - - # are now computing the derivative.. - if i not in self._atom_indexes: - return 0.0 - - k = int(component) - dqdx = 0.0 - - if i in (self.m, self.o): - sign = 1 if i == self.m else -1 - dqdx += sign * ( - np.cross(u, w)[k] / (lambda_u * np.sin(phi_u) ** 2) - ) - - if i in (self.p, self.n): - sign = 1 if i == self.p else -1 - dqdx += sign * ( - np.cross(v, w)[k] / (lambda_v * np.sin(phi_v) ** 2) - ) - - if i in (self.o, self.p): - sign = 1 if i == self.o else -1 - dqdx += sign * ( - ( - (np.cross(u, w)[k] * np.cos(phi_u)) - / (lambda_w * np.sin(phi_u) ** 2) - ) - - ( - (np.cross(v, w)[k] * np.cos(phi_v)) - / (lambda_w * np.sin(phi_v) ** 2) - ) - ) - - return dqdx + def _evaluate( + self, x: "CartesianCoordinates", deriv_order: DerivativeOrder + ) -> "VectorHyperDual": + """Dihedral m-o-p-n""" + # https://en.wikipedia.org/wiki/Dihedral_angle#In_polymer_physics + _x = x.ravel() + vec_m, vec_o, vec_p, vec_n = _get_3d_vecs_from_atom_idxs( + self.m, self.o, self.p, self.n, x=_x, deriv_order=deriv_order + ) + u_1 = vec_o - vec_m + u_2 = vec_p - vec_o + u_3 = vec_n - vec_p + + norm_u2 = u_2.norm() + v1 = u_2.cross(u_3) + v2 = u_1.cross(u_2) + v3 = u_1 * norm_u2 + dihedral = DifferentiableMath.atan2(v3.dot(v1), v2.dot(v1)) + assert isinstance(dihedral, VectorHyperDual) + return dihedral def __repr__(self): return f"Dihedral({self.m}-{self.o}-{self.p}-{self.n})" diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py index f6056aa68..8dceaccf5 100644 --- a/autode/opt/optimisers/crfo.py +++ b/autode/opt/optimisers/crfo.py @@ -238,6 +238,10 @@ def _dihedrals(species): if n == o: continue + # avoid triangular rings like cyclopropane + if m == n: + continue + if np.isclose(species.angle(o, p, n), Angle(np.pi), atol=0.04): continue diff --git a/doc/changelog.rst b/doc/changelog.rst index 460ed3494..1d2b97732 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -7,11 +7,11 @@ Changelog Functionality improvements ************************** -- ... +- Replaces hard-coded derivatives for primitive internal coordinates with automatic differentiation Bug Fixes ********* -- ... +- Fixes triangular rings being incorrectly treated as dihedral angles 1.4.1 diff --git a/tests/test_opt/test_autodiff.py b/tests/test_opt/test_autodiff.py new file mode 100644 index 000000000..eec7065d1 --- /dev/null +++ b/tests/test_opt/test_autodiff.py @@ -0,0 +1,190 @@ +import math +import numpy as np +import pytest +from autode.opt.coordinates._autodiff import ( + get_differentiable_vars, + DifferentiableMath, + VectorHyperDual, + DerivativeOrder, + DifferentiableVector3D, +) + + +def test_hyperdual_sanity_checks(): + a, b = get_differentiable_vars( + values=[1, 2], symbols=["a", "b"], deriv_order=DerivativeOrder.first + ) + assert repr(a) is not None + + c = get_differentiable_vars( + values=[1], symbols=["c"], deriv_order=DerivativeOrder.first + )[0] + with pytest.raises(ValueError, match="Incompatible number"): + _ = a + c + + c, d = get_differentiable_vars( + values=[1, 2], symbols=["c", "d"], deriv_order=DerivativeOrder.first + ) + with pytest.raises(ValueError, match="symbols do not match"): + _ = a * c + + a2, b2 = get_differentiable_vars( + values=[1, 2], symbols=["a", "b"], deriv_order=DerivativeOrder.second + ) + with pytest.raises(ValueError, match="order of derivative do not match!"): + _ = DifferentiableMath.atan2(a2, b) + + with pytest.raises(RuntimeError, match="Symbols must be unique"): + _ = get_differentiable_vars( + values=[1, 2], + symbols=["a", "a"], + deriv_order=DerivativeOrder.first, + ) + + +def x_to_the_y_derivs(x, y): + d_dx = y * x ** (y - 1) + d_dy = x**y * math.log(x) + return d_dx, d_dy + + +def x_to_the_y_second_derivs(x, y): + d2_dx2 = y * (y - 1) * x ** (y - 2) + d2_dxdy = math.log(x) * y * x ** (y - 1) + x**y * (1 / x) + d2_dydx = x ** (y - 1) + y * x ** (y - 1) * math.log(x) + d2_dy2 = math.log(x) * math.log(x) * x**y + return d2_dx2, d2_dxdy, d2_dydx, d2_dy2 + + +def test_autodiff_exponential_func(): + x_val, y_val = 1.0, 2.0 + x, y = get_differentiable_vars( + values=[x_val, y_val], + symbols=["x", "y"], + ) + result = x**y + d_dx, d_dy = x_to_the_y_derivs(x_val, y_val) + assert math.isclose(result.differentiate_wrt("x"), d_dx) + assert math.isclose(result.differentiate_wrt("y"), d_dy) + + d2_dx2, d2_dxdy, d2_dydx, d2_dy2 = x_to_the_y_second_derivs(x_val, y_val) + assert math.isclose(result.differentiate_wrt("x", "x"), d2_dx2) + assert math.isclose(result.differentiate_wrt("x", "y"), d2_dxdy) + assert math.isclose(result.differentiate_wrt("y", "x"), d2_dydx) + assert math.isclose(result.differentiate_wrt("y", "y"), d2_dy2) + # check the second derivatives are symmetric + assert math.isclose(d2_dydx, d2_dxdy) + + +def test_exponential_math_sanity_checks(): + x, y = get_differentiable_vars([-0.1, -0.2], ["x", "y"]) + assert x**2 is not None + assert x**-1 is not None + # negative number raised to fractional power is complex + with pytest.raises(AssertionError): + _ = x**0.2 + # negative number cannot be raised to differentiable power + with pytest.raises(AssertionError): + _ = (-2) ** y + with pytest.raises(AssertionError): + _ = x**y + + +def test_math_funcs_work_with_native_types(): + # python float or int types should be passed through + assert math.isclose(DifferentiableMath.acos(0), math.pi / 2) + (y,) = get_differentiable_vars([1], ["y"]) + res = DifferentiableMath.atan2(y, 0) + assert isinstance(res, VectorHyperDual) + assert math.isclose(res.value, math.pi / 2) + res = DifferentiableMath.atan2(1, 0) + assert isinstance(res, float) + assert math.isclose(DifferentiableMath.pow(0.9, 1.3), math.pow(0.9, 1.3)) + + # however, python's complex type is not supported + with pytest.raises(TypeError, match="Unknown type for addition"): + _ = y + (4 + 1j) + with pytest.raises(TypeError, match="Unknown type for multiplication"): + _ = y * (1 + 2j) + with pytest.raises(TypeError, match="Unknown type for exponentiation"): + _ = y ** (1 + 2j) + + +def test_hyperdual_init_checks(): + with pytest.raises(ValueError): + _ = VectorHyperDual( + 0.1, symbols=["x", "y"], first_der=np.array([0.1, 0.2, 0.3]) + ) + + with pytest.raises(ValueError): + _ = VectorHyperDual( + 0.1, + symbols=["x", "y"], + first_der=np.array([0.1, 0.2]), + second_der=np.zeros(shape=(2, 3)), + ) + + +def test_hyperdual_order(): + symbols = ["x", "y"] + x = VectorHyperDual( + 0.1, + symbols=symbols, + ) + assert x._order == DerivativeOrder.zeroth + x = VectorHyperDual(0.1, symbols=symbols, first_der=np.zeros(2)) + assert x._order == DerivativeOrder.first + # will ignore second derivatives if first is not present + x = VectorHyperDual(0.1, symbols=symbols, second_der=np.zeros((3, 2))) + assert x._order == DerivativeOrder.zeroth + x = VectorHyperDual( + 0.1, + symbols=symbols, + first_der=np.zeros(2), + second_der=np.zeros((2, 2)), + ) + assert x._order == DerivativeOrder.second + + +def test_derivative_not_available(): + x, y = get_differentiable_vars( + [1.0, 2.0], symbols=["x", "y"], deriv_order=DerivativeOrder.first + ) + res = 1 - x**2 + y + assert isinstance(res, VectorHyperDual) + assert math.isclose(res.value, 2) + assert math.isclose(res.differentiate_wrt("x"), -2) + # higher order derivatives are not available + assert res.differentiate_wrt("x", "y") is None + # unknown variables + assert res.differentiate_wrt("z") is None + assert res.differentiate_wrt("y", "z") is None + + # with zero order, only value is present + (x,) = get_differentiable_vars( + [1.0], symbols=["x"], deriv_order=DerivativeOrder.zeroth + ) + res = 1 + x**2 + assert isinstance(res, VectorHyperDual) + assert res.differentiate_wrt("x") is None + + +def test_hyperdual_3d_vector_sanity_checks(): + x, y, z = get_differentiable_vars( + [0.1, 0.2, 0.3], + ["x", "y", "z"], + ) + + with pytest.raises(ValueError): + _ = DifferentiableVector3D([x, y]) + + vec1 = DifferentiableVector3D([x, y, z]) + vec2 = DifferentiableVector3D([x, y, z]) + + with pytest.raises(ValueError): + _ = vec1.dot(2) + + with pytest.raises(AssertionError): + _ = vec1 * vec2 + + assert 2 * vec1 is not None diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py index 2f775e81b..9fe31f72e 100644 --- a/tests/test_opt/test_coordiantes.py +++ b/tests/test_opt/test_coordiantes.py @@ -28,18 +28,12 @@ def test_inv_dist_primitives(): assert np.isclose(inv_dist(x), 0.5) # 1/2.0 = 0.5 Å-1 # Check a couple of derivatives by hand - assert np.isclose( - inv_dist.derivative(0, CartesianComponent.x, x=x), 2 * inv_dist(x) ** 3 - ) - assert np.isclose( - inv_dist.derivative(1, CartesianComponent.x, x=x), - -2 * inv_dist(x) ** 3, - ) + derivs = inv_dist.derivative(x=x) + assert np.isclose(derivs[3 * 0 + 0], 2 * inv_dist(x) ** 3) + assert np.isclose(derivs[3 * 1 + 0], -2 * inv_dist(x) ** 3) # Derivatives with respect to zero components - assert np.isclose(inv_dist.derivative(0, CartesianComponent.y, x=x), 0) - # or those that are not present in the system should be zero - assert np.isclose(inv_dist.derivative(2, CartesianComponent.x, x=x), 0) + assert np.isclose(derivs[3 * 0 + 1], 0) def test_dist_primitives(): @@ -50,18 +44,12 @@ def test_dist_primitives(): inv_dist = PrimitiveDistance(0, 1) assert np.isclose(inv_dist(x), 2.0) - assert np.isclose( - inv_dist.derivative(0, CartesianComponent.x, x=x), -2 / 2 - ) + derivs = inv_dist.derivative(x) + assert np.isclose(derivs[3 * 0 + 0], -2 / 2) + assert np.isclose(derivs[3 * 1 + 0], +2 / 2) - assert np.isclose( - inv_dist.derivative(1, CartesianComponent.x, x=x), +2 / 2 - ) - - for component in (CartesianComponent.y, CartesianComponent.z): - assert np.isclose(inv_dist.derivative(1, component, x=x), 0) - - assert np.isclose(inv_dist.derivative(2, CartesianComponent.x, x=x), 0) + for k in (1, 2): + assert np.isclose(derivs[3 * 1 + k], 0) def test_primitive_equality(): @@ -539,6 +527,14 @@ def test_pic_b_no_primitives(): c._calc_B(np.arange(6, dtype=float).reshape(2, 3)) +def test_pic_append_type_checking(): + c = PIC() + # append should check for primitive type + c.append(PrimitiveDistance(0, 1)) + with pytest.raises(AssertionError): + c.append(3) + + def test_constrained_distance_satisfied(): d = ConstrainedPrimitiveDistance(0, 1, value=1.0) @@ -563,20 +559,15 @@ def numerical_derivative(a, b, h=1e-8): init_coords = m.coordinates.copy() angle = PrimitiveBondAngle(0, 1, 2) - + derivs = angle.derivative(init_coords) for atom_idx in (0, 1, 2): - for component in CartesianComponent: - analytic = angle.derivative(atom_idx, component, init_coords) + for component in (0, 1, 2): + analytic = derivs[3 * atom_idx + component] assert np.isclose( analytic, numerical_derivative(atom_idx, component), atol=1e-6 ) - # Derivative should be zero for an atom not present the bond angle - assert np.isclose( - angle.derivative(3, CartesianComponent.x, init_coords), 0.0 - ) - def test_angle_primitive_equality(): assert PrimitiveBondAngle(0, 1, 2) == PrimitiveBondAngle(0, 2, 1) @@ -605,17 +596,11 @@ def numerical_derivative(a, b, h=1e-8): init_coords = m.coordinates.copy() dihedral = PrimitiveDihedralAngle(2, 0, 1, 3) - + analytic = dihedral.derivative(init_coords) for atom_idx in (0, 1, 2, 3): - for component in CartesianComponent: - analytic = dihedral.derivative(atom_idx, component, init_coords) - numerical = numerical_derivative(atom_idx, component) - - assert np.isclose(analytic, numerical, atol=1e-6) - - assert np.isclose( - dihedral.derivative(5, CartesianComponent.x, init_coords), 0.0 - ) + for k in [0, 1, 2]: + numerical = numerical_derivative(atom_idx, k) + assert np.isclose(analytic[3 * atom_idx + k], numerical, atol=1e-6) def test_dihedral_equality(): @@ -627,40 +612,91 @@ def test_dihedral_equality(): ) -@pytest.mark.parametrize( - "h_coord", - [ - np.array([1.08517, 1.07993, 0.05600]), - np.array([1.28230, -0.63391, -0.54779]), - ], -) -def test_dihedral_primitive_derivative_over_zero(h_coord): - def numerical_derivative(a, b, h=1e-6): - coords = init_coords.copy() - coords[a, int(b)] += h - y_plus = dihedral(coords) - coords[a, int(b)] -= 2 * h - y_minus = dihedral(coords) - return (y_plus - y_minus) / (2 * h) +def test_primitives_consistent_with_mol_values(): + # test that the primitive values are the same as the mol.distance etc. + h2o2 = h2o2_mol() + coords = h2o2.coordinates + dist = PrimitiveDistance(0, 1) + assert np.isclose(dist(coords), h2o2.distance(0, 1), rtol=1e-8) + invdist = PrimitiveInverseDistance(1, 2) + assert np.isclose(invdist(coords), 1 / h2o2.distance(1, 2), rtol=1e-8) + ang = PrimitiveBondAngle(2, 0, 1) # bond is defined in a different way + assert np.isclose(ang(coords), h2o2.angle(0, 2, 1), rtol=1e-8) + dihedral = PrimitiveDihedralAngle(2, 0, 1, 3) + assert np.isclose(dihedral(coords), h2o2.dihedral(2, 0, 1, 3), rtol=1e-8) - m = Molecule( + +# fmt: off +dihedral_mols = [ + Molecule( atoms=[ Atom("C", 0.63365, 0.11934, -0.13163), Atom("C", -0.63367, -0.11938, 0.13153), - Atom("H", 0.0, 0.0, 0.0), + Atom("H", 1.08517, 1.07993, 0.05600), Atom("H", -1.08517, -1.07984, -0.05599), ] - ) - m.atoms[2].coord = h_coord - init_coords = m.coordinates - - dihedral = PrimitiveDihedralAngle(2, 0, 1, 3) - - for atom_idx in (0, 1, 2, 3): - for component in CartesianComponent: - analytic = dihedral.derivative(atom_idx, component, init_coords) - numerical = numerical_derivative(atom_idx, component) - assert np.isclose(analytic, numerical, atol=1e-6) + ), + Molecule( + atoms=[ + Atom("C", 0.63365, 0.11934, -0.13163), + Atom("C", -0.63367, -0.11938, 0.13153), + Atom("H", 1.28230, -0.63391, -0.54779), + Atom("H", -1.08517, -1.07984, -0.05599), + ] + ) # for testing dihedral derivatives over zero +] + +test_mols = [ + h2o2_mol(), h2o2_mol(), water_mol(), + water_mol(), water_mol(), *dihedral_mols +] +test_prims = [ + PrimitiveDihedralAngle(2, 0, 1, 3), PrimitiveBondAngle(2, 0, 1), + PrimitiveBondAngle(0, 1, 2), PrimitiveDistance(0, 1), + PrimitiveInverseDistance(0, 1), PrimitiveDihedralAngle(2, 0, 1, 3), + PrimitiveDihedralAngle(2, 0, 1, 3) +] +# fmt: on + + +@pytest.mark.parametrize("mol,prim", list(zip(test_mols, test_prims))) +def test_primitive_first_derivs(mol, prim): + init_coords = CartesianCoordinates(mol.coordinates) + init_prim = prim(init_coords) + + def numerical_first_deriv(coords, h=1e-8): + coords = coords.flatten() + derivs = np.zeros_like(coords) + for i in range(coords.shape[0]): + coords[i] += h + derivs[i] = (prim(coords) - init_prim) / h + coords[i] -= h + return derivs + + analytic = prim.derivative(init_coords) + numeric = numerical_first_deriv(init_coords) + assert np.allclose(analytic, numeric, atol=1e-6) + + +@pytest.mark.parametrize("mol,prim", list(zip(test_mols, test_prims))) +def test_primitve_second_deriv(mol, prim): + init_coords = CartesianCoordinates(mol.coordinates) + init_first_der = prim.derivative(init_coords) + + def numerical_second_deriv(coords, h=1e-8): + coords = coords.flatten() + derivs = np.zeros((coords.shape[0], coords.shape[0])) + for i in range(coords.shape[0]): + coords[i] += h + derivs[i] = (prim.derivative(coords) - init_first_der) / h + coords[i] -= h + return derivs + + analytic = prim.second_derivative(init_coords) + # second derivative matrix should be symmetric + assert np.allclose(analytic, analytic.T) + numeric = numerical_second_deriv(init_coords) + assert np.allclose(analytic, numeric, atol=1e-6) def test_repr(): @@ -679,14 +715,6 @@ def test_repr(): assert repr(p) is not None -@pytest.mark.parametrize("sign", [1, -1]) -def test_angle_normal(sign): - angle = PrimitiveBondAngle(0, 1, 2) - x = np.array([[0.0, 0.0, 0.0], [1.0, sign * 1.0, 1.0], [1.0, 0.0, 1.0]]) - - assert not np.isinf(angle.derivative(0, 1, x)) - - def test_dic_large_step_allowed_unconverged_back_transform(): x = CartesianCoordinates(water_mol().coordinates) dic = DIC.from_cartesian(x) From 05b2d985b5f5752659184af0644c0cbe97ab6b34 Mon Sep 17 00:00:00 2001 From: Shoubhik Maiti <17470159+shoubhikraj@users.noreply.github.com> Date: Thu, 4 Jan 2024 10:05:09 +0000 Subject: [PATCH 3/5] More comprehensive primitive generation (#314) * initial commit * second commit * second commit * simplify code * linear angle * linear angle * test update, more fixes needed * unifinished update * linear angle update * unfinished update * unfinished update * unifinished update * unfinished update * better linear angles * better linear angles * minor change * fix linear bend code * fix some tests * add tests * fix failing test * add tests * minor dihedral fix * dihedral code fix * add test for linear bend with ref * add test for linear bend with ref * add tests * make bond angle consistent * add tests and fix value bug * add tests * fix test * change backup hessian update scheme for crfo * revert change to hessian update scheme * correction to angle code * assert only for graph * fix tests * pr suggestions 1 * pr suggestions 2 partial * test fix * put methods into anypic class * pr update 3 * pr update 4 * pr update 5 * pr update 6 * pr update 7 * pr update 8 * pr update * pr update, changelog update --- autode/atoms.py | 5 +- autode/mol_graphs.py | 9 + autode/opt/coordinates/__init__.py | 2 +- autode/opt/coordinates/_autodiff.py | 31 ++- autode/opt/coordinates/base.py | 10 - autode/opt/coordinates/internals.py | 402 ++++++++++++++++++++++++++- autode/opt/coordinates/primitives.py | 161 ++++++++++- autode/opt/optimisers/crfo.py | 99 +------ autode/opt/optimisers/trm.py | 21 -- autode/values.py | 4 + doc/changelog.rst | 2 + tests/test_opt/molecules.py | 45 +++ tests/test_opt/test_coordiantes.py | 238 ++++++++++++++-- tests/test_opt/test_crfo.py | 65 ++++- 14 files changed, 910 insertions(+), 184 deletions(-) diff --git a/autode/atoms.py b/autode/atoms.py index 699e7f67f..6e88d8f9f 100644 --- a/autode/atoms.py +++ b/autode/atoms.py @@ -1041,9 +1041,10 @@ def angle(self, i: int, j: int, k: int) -> Angle: f"least one zero vector" ) - value = np.arccos(np.dot(vec1, vec2) / norms) + # Cos(theta) must lie within [-1, 1] + cos_value = np.clip(np.dot(vec1, vec2) / norms, a_min=-1, a_max=1) - return Angle(value) + return Angle(np.arccos(cos_value)) def dihedral(self, w: int, x: int, y: int, z: int) -> Angle: r""" diff --git a/autode/mol_graphs.py b/autode/mol_graphs.py index 9ef148633..521036a90 100644 --- a/autode/mol_graphs.py +++ b/autode/mol_graphs.py @@ -115,6 +115,15 @@ def node_matcher(self): return matcher + @property + def is_connected(self) -> bool: + """Is this graph fully connected (i.e. not separate parts)""" + return nx.is_connected(self) + + def connected_components(self): + """Generate the separate connected components""" + return nx.connected_components(self) + def make_graph( species: "Species", diff --git a/autode/opt/coordinates/__init__.py b/autode/opt/coordinates/__init__.py index 2aadb5de9..36daee459 100644 --- a/autode/opt/coordinates/__init__.py +++ b/autode/opt/coordinates/__init__.py @@ -1,3 +1,3 @@ -from autode.opt.coordinates.base import OptCoordinates, CartesianComponent +from autode.opt.coordinates.base import OptCoordinates from autode.opt.coordinates.cartesian import CartesianCoordinates from autode.opt.coordinates.dic import DIC, DICWithConstraints diff --git a/autode/opt/coordinates/_autodiff.py b/autode/opt/coordinates/_autodiff.py index 0fa48d437..0855a441e 100644 --- a/autode/opt/coordinates/_autodiff.py +++ b/autode/opt/coordinates/_autodiff.py @@ -580,7 +580,9 @@ class DifferentiableVector3D: hyper-dual numbers """ - def __init__(self, items: Sequence["VectorHyperDual"]): + def __init__( + self, items: Sequence[Union["VectorHyperDual", numeric_type]] + ): """ Initialise the 3D vector from a list of 3 hyperdual numbers @@ -590,7 +592,9 @@ def __init__(self, items: Sequence["VectorHyperDual"]): items = list(items) if len(items) != 3: raise ValueError("A 3D vector must have only 3 components") - assert all(isinstance(item, VectorHyperDual) for item in items) + assert all( + isinstance(item, (VectorHyperDual, *numeric)) for item in items + ) self._data = items @staticmethod @@ -600,7 +604,9 @@ def _check_same_type(other) -> None: raise ValueError("Operation must be done with another 3D vector!") return None - def dot(self, other: "DifferentiableVector3D") -> "VectorHyperDual": + def dot( + self, other: "DifferentiableVector3D" + ) -> Union["VectorHyperDual", numeric_type]: """ Dot product of two 3D vectors @@ -611,13 +617,12 @@ def dot(self, other: "DifferentiableVector3D") -> "VectorHyperDual": (VectorHyperDual): A scalar number (with derivatives) """ self._check_same_type(other) - dot = 0 + dot: Union[VectorHyperDual, numeric_type] = 0 for k in range(3): dot = dot + self._data[k] * other._data[k] - assert isinstance(dot, VectorHyperDual) return dot - def norm(self) -> "VectorHyperDual": + def norm(self) -> Union["VectorHyperDual", numeric_type]: """ Euclidean (l2) norm of this 3D vector @@ -627,7 +632,6 @@ def norm(self) -> "VectorHyperDual": norm = DifferentiableMath.sqrt( self._data[0] ** 2 + self._data[1] ** 2 + self._data[2] ** 2 ) - assert isinstance(norm, VectorHyperDual) return norm def __add__( @@ -690,6 +694,18 @@ def __rmul__(self, other): """Multiplication of scalar and vector is commutative""" return self.__mul__(other) + def __truediv__(self, other: Union[VectorHyperDual, numeric_type]): + """ + Division of a 3D vector with a scalar + + Args: + other (VectorHyperDual|float|int): + + Returns: + (DifferentiableVector3D): + """ + return self.__mul__(1 / other) + def cross( self, other: "DifferentiableVector3D" ) -> "DifferentiableVector3D": @@ -702,6 +718,7 @@ def cross( Returns: (DifferentiableVector3D): """ + self._check_same_type(other) return DifferentiableVector3D( [ self._data[1] * other._data[2] diff --git a/autode/opt/coordinates/base.py b/autode/opt/coordinates/base.py index 0f5135ced..747680e99 100644 --- a/autode/opt/coordinates/base.py +++ b/autode/opt/coordinates/base.py @@ -1,7 +1,6 @@ # mypy: disable-error-code="has-type" import numpy as np from copy import deepcopy -from enum import IntEnum, unique from typing import Optional, Union, Sequence, List, TYPE_CHECKING from abc import ABC, abstractmethod @@ -15,15 +14,6 @@ from autode.hessians import Hessian -@unique -class CartesianComponent(IntEnum): - """Cartesian component in 3D space""" - - x = 0 - y = 1 - z = 2 - - class OptCoordinates(ValueArray, ABC): """Coordinates used to perform optimisations""" diff --git a/autode/opt/coordinates/internals.py b/autode/opt/coordinates/internals.py index 9e7462e16..3de919b13 100644 --- a/autode/opt/coordinates/internals.py +++ b/autode/opt/coordinates/internals.py @@ -6,20 +6,29 @@ B : Wilson B matrix q : Primitive internal coordinates G : Spectroscopic G matrix + +Set-up of redundant primitives is based on J. Chem. Phys., 117, 2002, 9160 """ import numpy as np - +import itertools from typing import Any, Optional, Type, List, TYPE_CHECKING from abc import ABC, abstractmethod -from autode.opt.coordinates.base import OptCoordinates, CartesianComponent +from autode.values import Angle, Distance +from autode.opt.coordinates.base import OptCoordinates from autode.opt.coordinates.primitives import ( PrimitiveInverseDistance, Primitive, PrimitiveDistance, + ConstrainedPrimitiveDistance, + PrimitiveBondAngle, + PrimitiveDummyLinearAngle, + PrimitiveLinearAngle, PrimitiveDihedralAngle, + LinearBendType, ) if TYPE_CHECKING: + from autode.species import Species from autode.opt.coordinates.cartesian import CartesianCoordinates from autode.opt.coordinates.primitives import ( ConstrainedPrimitive, @@ -27,12 +36,19 @@ ) +# Angle threshold for linearity (should be in radians) +_lin_thresh = Angle(170, "deg").to("rad") + + class InternalCoordinates(OptCoordinates, ABC): # lgtm [py/missing-equals] def __new__(cls, input_array) -> "InternalCoordinates": """New instance of these internal coordinates""" arr = super().__new__(cls, input_array, units="Å") + arr._x = None + arr.primitives = None + for attr in ("_x", "primitives"): setattr(arr, attr, getattr(input_array, attr, None)) @@ -82,10 +98,18 @@ def __init__(self, *args: Any): f"from {args}. Must be primitive internals" ) + def add(self, item: Primitive) -> None: + """Add a primitive to this set of primitive coordinates""" + assert isinstance(item, Primitive), "Must be a primitive" + # prevent duplication of primitives + if item not in self: + super().append(item) + def append(self, item: Primitive) -> None: """Append an item to this set of primitives""" - assert isinstance(item, Primitive), "Must be a Primitive type!" - super().append(item) + raise NotImplementedError( + "Please use PIC.add() to add new primitives to the set" + ) @property def B(self) -> np.ndarray: @@ -210,7 +234,7 @@ def _populate_all(self, x: np.ndarray): # Add all the unique inverse distances (i < j) for i in range(n_atoms): for j in range(i + 1, n_atoms): - self.append(self._primitive_type(i, j)) + self.add(self._primitive_type(i, j)) return None @@ -234,3 +258,371 @@ def _primitive_type(self): class AnyPIC(PIC): def _populate_all(self, x: np.ndarray) -> None: raise RuntimeError("Cannot populate all on an AnyPIC instance") + + @classmethod + def from_species(cls, mol: "Species") -> "AnyPIC": + """ + Build a set of primitives from the species, using the graph as + a starting point for the connectivity of the species. Also joins + any disjoint parts of the graph, and adds hydrogen bonds to + ensure that the primitives are redundant. + + Args: + mol: The species object + + Returns: + (AnyPIC): The set of primitive internals + """ + pic = cls() + # take a copy of mol as mol.graph might be changed + mol = mol.copy() + _connect_graph_for_species(mol) + pic._add_bonds_from_species(mol) + pic._add_angles_from_species(mol) + pic._add_dihedrals_from_species(mol) + pic._add_chain_dihedrals_from_species(mol) + return pic + + def _add_bonds_from_species( + self, + mol: "Species", + ): + """ + Add bonds to the current set of primitives, from the + connectivity graph of the species + + Args: + mol: The species object + """ + assert mol.graph is not None + + n = 0 + for i, j in sorted(mol.graph.edges): + if ( + mol.constraints.distance is not None + and (i, j) in mol.constraints.distance + ): + r = mol.constraints.distance[(i, j)] + self.add(ConstrainedPrimitiveDistance(i, j, r)) + n += 1 + else: + self.add(PrimitiveDistance(i, j)) + assert n == mol.constraints.n_distance + + return None + + @staticmethod + def _get_ref_for_linear_angle( + mol, + a, + b, + c, + bonded: bool, + dist_thresh=Distance(4, "ang"), + ) -> Optional[int]: + """ + Get a reference atom for describing a linear angle, which + must not itself be linear to the atoms in the angle in + any combination. The linear angle is a--b--c here. + + Args: + mol: + a: + b: + c: + bonded: Whether to look for only atoms bonded to the central + atom (b) for reference + dist_thresh: The distance threshold to check for atoms if + bonded = False + + Returns: + (int|None): The index of the ref. atom if found, else None + """ + + # only check bonded atoms if requested + if bonded: + near_atoms = list(mol.graph.neighbors(b)) + near_atoms.remove(a) + near_atoms.remove(c) + + # otherwise get all atoms in 4 A radius except a, b, c + else: + near_atoms = [ + idx + for idx in range(mol.n_atoms) + if mol.distance(b, idx) < dist_thresh and idx not in (a, b, c) + ] + + # get atoms closest to perpendicular + deviations_from_90 = {} + for atom in near_atoms: + i_b_a = mol.angle(atom, b, a) + if i_b_a > _lin_thresh or i_b_a < (np.pi - _lin_thresh): + continue + i_b_c = mol.angle(atom, b, c) + if i_b_c > _lin_thresh or i_b_c < (np.pi - _lin_thresh): + continue + deviation_a = abs(i_b_a - np.pi / 2) + deviation_b = abs(i_b_c - np.pi / 2) + avg_dev = (deviation_a + deviation_b) / 2 + deviations_from_90[atom] = avg_dev + + if len(deviations_from_90) == 0: + return None + + return min(deviations_from_90, key=deviations_from_90.get) # type: ignore + + def _add_angles_from_species( + self, + mol: "Species", + ) -> None: + """ + Modify the set of primitives in-place by adding angles, from the + connectivity graph supplied + + Args: + mol (Species): The species object + """ + assert mol.graph is not None + + for o in range(mol.n_atoms): + for n, m in itertools.combinations(mol.graph.neighbors(o), r=2): + if mol.angle(m, o, n) < _lin_thresh: + self.add(PrimitiveBondAngle(m=m, o=o, n=n)) + else: + # If central atom is connected to another atom, then the + # linear angle is skipped and instead an out-of-plane + # (improper dihedral) coordinate is used + r = self._get_ref_for_linear_angle( + mol, m, o, n, bonded=True + ) + if r is not None: + self.add(PrimitiveDihedralAngle(m, r, o, n)) + continue + + # Otherwise, we use a nearby (< 4.0 A) reference atom to + # define two orthogonal linear bends + r = self._get_ref_for_linear_angle( + mol, m, o, n, bonded=False + ) + if r is not None: + self.add( + PrimitiveLinearAngle( + m, o, n, r, LinearBendType.BEND + ) + ) + self.add( + PrimitiveLinearAngle( + m, o, n, r, LinearBendType.COMPLEMENT + ) + ) + + # For completely linear molecules (CO2), there will be no such + # reference atoms, so use dummy atoms instead + else: + self.add( + PrimitiveDummyLinearAngle( + m, o, n, LinearBendType.BEND + ) + ) + self.add( + PrimitiveDummyLinearAngle( + m, o, n, LinearBendType.COMPLEMENT + ) + ) + + return None + + def _add_dihedrals_from_species( + self, + mol: "Species", + ) -> None: + """ + Modify the set of primitives in-place by adding dihedrals (torsions), + from the connectivity graph supplied + + Args: + mol (Species): The species + """ + # no dihedrals possible with less than 4 atoms + if mol.n_atoms < 4: + return + + assert mol.graph is not None + + for o, p in list(mol.graph.edges): + for m in mol.graph.neighbors(o): + if m == p: + continue + + for n in mol.graph.neighbors(p): + if n == o: + continue + + # avoid triangle rings like cyclopropane + if n == m: + continue + + if _is_dihedral_well_defined(mol, m, o, p, n): + self.add(PrimitiveDihedralAngle(m, o, p, n)) + + return None + + @staticmethod + def _get_linear_chains(mol: "Species") -> List[List[int]]: + """ + Obtain a list of all the continuous chains of linear atoms + present in the species i.e. A--B--C--D--E... + + Args: + mol: The species + + Returns: + (list[list[int]): A list of lists, each containing indices + of the atoms of linear chains in order + """ + assert mol.graph is not None + + def extend_chain(chain: List[int]): + """Extend a chain in-place""" + for idx in range(mol.n_atoms): + if idx in chain: + continue + + if mol.angle(chain[1], chain[0], idx) > _lin_thresh: + chain.insert(0, idx) + continue + + if mol.angle(chain[-2], chain[-1], idx) > _lin_thresh: + chain.append(idx) + continue + return None + + linear_chains: List[list] = [] + for b in range(mol.n_atoms): + for a, c in itertools.combinations(mol.graph.neighbors(b), r=2): + if any( + a in chain and b in chain and c in chain + for chain in linear_chains + ): + continue + if mol.angle(a, b, c) > _lin_thresh: + chain = [a, b, c] + extend_chain(chain) + linear_chains.append(chain) + + return linear_chains + + def _add_chain_dihedrals_from_species(self, mol: "Species"): + """ + Add extra dihedrals for chain molecules like allene, which + are required to cover all degrees of freedom + + Args: + mol: + + Returns: + + """ + assert mol.graph is not None + linear_chains = self._get_linear_chains(mol) + + for chain in linear_chains: + o, p = chain[0], chain[-1] + for m in mol.graph.neighbors(o): + if m == p: + continue + + if m in chain: + continue + + for n in mol.graph.neighbors(p): + if n == o: + continue + + if n == m: + continue + + if n in chain: + continue + + if _is_dihedral_well_defined(mol, m, o, p, n): + self.add(PrimitiveDihedralAngle(m, o, p, n)) + + return None + + +def _connect_graph_for_species(mol: "Species") -> None: + """ + Creates a fully connected graph from the graph of a species, by + (1) joining disconnected fragments by their shortest distance, + (2) connecting constrained bonds, (3) joining hydrogen bonds, + if present. The molecular graph is modified in-place. + + Args: + mol: A species (must have atoms and graph) + """ + assert mol.graph is not None, "Species must have graph!" + + # join hydrogen bonds + h_bond_x = ["N", "O", "F", "P", "S", "Cl"] + for i, j in itertools.combinations(range(mol.n_atoms), r=2): + if ( + mol.atoms[i].label in h_bond_x + and mol.atoms[j].label == "H" + or mol.atoms[j].label in h_bond_x + and mol.atoms[i].label == "H" + ): + vdw_sum = mol.atoms[i].vdw_radius + mol.atoms[j].vdw_radius + if mol.distance(i, j) < 0.9 * vdw_sum: + if not mol.graph.has_edge(i, j): + mol.graph.add_edge(i, j, pi=False, active=False) + + # join disconnected graph components + if not mol.graph.is_connected: + components = mol.graph.connected_components() + for comp_i, comp_j in itertools.combinations(components, r=2): + min_dist = float("inf") + min_pair = (-1, -1) + for i, j in itertools.product(list(comp_i), list(comp_j)): + if mol.distance(i, j) < min_dist: + min_dist = mol.distance(i, j) + min_pair = (i, j) + mol.graph.add_edge(*min_pair, pi=False, active=False) + + assert mol.graph.is_connected, "Unknown error in connecting graph" + + # The constraints should be counted as bonds + if mol.constraints.distance is not None: + for i, j in mol.constraints.distance: + if not mol.graph.has_edge(i, j): + mol.graph.add_edge(i, j, pi=False, active=False) + + return None + + +def _is_dihedral_well_defined(mol, a, b, c, d) -> bool: + """ + A dihedral a--b--c--d is only well-defined when the + two constituent angles are not linear + + Args: + mol: + a: + b: + c: + d: + + Returns: + (bool): True if well-defined otherwise False + """ + zero_angle_thresh = np.pi - _lin_thresh + is_linear_1 = ( + mol.angle(a, b, c) > _lin_thresh + or mol.angle(a, b, c) < zero_angle_thresh + ) + is_linear_2 = ( + mol.angle(b, c, d) > _lin_thresh + or mol.angle(b, c, d) < zero_angle_thresh + ) + return not (is_linear_1 or is_linear_2) diff --git a/autode/opt/coordinates/primitives.py b/autode/opt/coordinates/primitives.py index 764c322b0..52efca8a1 100644 --- a/autode/opt/coordinates/primitives.py +++ b/autode/opt/coordinates/primitives.py @@ -1,7 +1,8 @@ import numpy as np from abc import ABC, abstractmethod -from typing import Tuple, TYPE_CHECKING, List +from enum import Enum +from typing import Tuple, TYPE_CHECKING, List, Optional from autode.opt.coordinates._autodiff import ( get_differentiable_vars, DifferentiableMath, @@ -11,7 +12,7 @@ ) if TYPE_CHECKING: - from autode.opt.coordinates import CartesianCoordinates, CartesianComponent + from autode.opt.coordinates import CartesianCoordinates def _get_3d_vecs_from_atom_idxs( @@ -238,7 +239,7 @@ def _evaluate( vec_i, vec_j = _get_3d_vecs_from_atom_idxs( self.i, self.j, x=x, deriv_order=deriv_order ) - return 1.0 / (vec_i - vec_j).norm() + return 1.0 / (vec_i - vec_j).norm() # type: ignore def __repr__(self): return f"InverseDistance({self.i}-{self.j})" @@ -260,7 +261,7 @@ def _evaluate( vec_i, vec_j = _get_3d_vecs_from_atom_idxs( self.i, self.j, x=x, deriv_order=deriv_order ) - return (vec_i - vec_j).norm() + return (vec_i - vec_j).norm() # type: ignore def __repr__(self): return f"Distance({self.i}-{self.j})" @@ -298,12 +299,12 @@ class PrimitiveBondAngle(Primitive): arccosine of the normalised dot product """ - def __init__(self, o: int, m: int, n: int): + def __init__(self, m: int, o: int, n: int): """Bond angle m-o-n""" - super().__init__(o, m, n) + super().__init__(m, o, n) - self.o = int(o) self.m = int(m) + self.o = int(o) self.n = int(n) def __eq__(self, other) -> bool: @@ -331,22 +332,22 @@ def __repr__(self): class ConstrainedPrimitiveBondAngle(ConstrainedPrimitive, PrimitiveBondAngle): - def __init__(self, o: int, m: int, n: int, value: float): + def __init__(self, m: int, o: int, n: int, value: float): """ Angle (m-o-n) constrained to a value (in radians) ----------------------------------------------------------------------- Arguments: - o: Atom index - m: Atom index + o: Atom index + n: Atom index value: Required value of the constrained angle """ - super().__init__(o=o, m=m, n=n) + super().__init__(m=m, o=o, n=n) self._theta0 = value @@ -376,7 +377,7 @@ def __init__(self, m: int, o: int, p: int, n: int): self.n = int(n) def __eq__(self, other) -> bool: - """Equality of two distance functions""" + """Equality of two dihedral angles""" return isinstance(other, self.__class__) and ( self._atom_indexes == other._atom_indexes or self._atom_indexes == tuple(reversed(other._atom_indexes)) @@ -405,3 +406,139 @@ def _evaluate( def __repr__(self): return f"Dihedral({self.m}-{self.o}-{self.p}-{self.n})" + + +class LinearBendType(Enum): + """For linear angles, there are two orthogonal directions""" + + BEND = 0 + COMPLEMENT = 1 + + +class LinearAngleBase(Primitive, ABC): + def __init__(self, m: int, o: int, n: int, r: int, axis: LinearBendType): + """Linear Bend: m-o-n""" + super().__init__(m, o, n, r) + self.m = int(m) + self.o = int(o) + self.n = int(n) + self.r = int(r) + + assert isinstance(axis, LinearBendType) + self.axis = axis + + def __eq__(self, other): + """Equality of two linear bend angles""" + return isinstance(other, self.__class__) and ( + self._atom_indexes == other._atom_indexes + and self.axis == other.axis + ) + + def _calc_linear_bend( + self, + m_vec: DifferentiableVector3D, + o_vec: DifferentiableVector3D, + n_vec: DifferentiableVector3D, + r_vec: DifferentiableVector3D, + ): + """ + Evaluate the linear bend from the vector positions of the + atoms involved in the angle m, o, n, and the reference + atom (or dummy atom) r + + Args: + m_vec: + o_vec: + n_vec: + r_vec: + + Returns: + (VectorHyperDual): The value, optionally containing derivatives + """ + # As defined in J. Comput. Chem., 20(10), 1999, 1067 + o_m = m_vec - o_vec + o_n = n_vec - o_vec + o_r = r_vec - o_vec + + # eq.(44) p 1073 + u = o_m.cross(o_r) + u = u / u.norm() + + # eq. (46) and (47) p 1074 + if self.axis == LinearBendType.BEND: + return u.dot(o_n) / o_n.norm() + elif self.axis == LinearBendType.COMPLEMENT: + return u.dot(o_n.cross(o_m)) / (o_n.norm() * o_m.norm()) + else: + raise ValueError("Unknown axis for linear bend") + + +class PrimitiveLinearAngle(LinearAngleBase): + """Linear Angle w.r.t. a reference atom""" + + def _evaluate( + self, x: "CartesianCoordinates", deriv_order: DerivativeOrder + ): + """Linear Bend angle m-o-n against reference atom r""" + + _x = x.ravel() + vec_m, vec_o, vec_n, vec_r = _get_3d_vecs_from_atom_idxs( + self.m, self.o, self.n, self.r, x=_x, deriv_order=deriv_order + ) + + return self._calc_linear_bend(vec_m, vec_o, vec_n, vec_r) + + def __repr__(self): + axis_str = "B" if self.axis == LinearBendType.BEND else "C" + return f"LinearBend{axis_str}({self.m}-{self.o}-{self.n}, {self.r})" + + +class PrimitiveDummyLinearAngle(LinearAngleBase): + """Linear bend with a dummy atom""" + + def __init__(self, m: int, o: int, n: int, axis: LinearBendType): + super().__init__(m, o, n, -1, axis) + + self._vec_r: Optional[DifferentiableVector3D] = None + + def _get_dummy_atom( + self, x: "CartesianCoordinates" + ) -> DifferentiableVector3D: + """Create the dummy atom r""" + _x = x.reshape(-1, 3) + cart_axes = [ + np.array([1.0, 0.0, 0.0]), + np.array([0.0, 1.0, 0.0]), + np.array([0.0, 0.0, 1.0]), + ] + + # choose cartesian axis with the lowest overlap with m-o vector + w = _x[self.m] - _x[self.o] + w /= np.linalg.norm(w) + overlaps = [] + for axis in cart_axes: + overlaps.append(np.dot(w, axis)) + cart_ax = cart_axes[np.argmin(np.abs(overlaps))] + + # place dummy atom perpendicular to m-o bond 1 A away + perp_axis = np.cross(cart_ax, w) + _o = _x[self.o] + dummy_point = _o + perp_axis / np.linalg.norm(perp_axis) + return DifferentiableVector3D(list(dummy_point)) + + def _evaluate( + self, x: "CartesianCoordinates", deriv_order: DerivativeOrder + ): + if self._vec_r is None: + self._vec_r = self._get_dummy_atom(x) + + _x = x.ravel() + vec_m, vec_o, vec_n = _get_3d_vecs_from_atom_idxs( + self.m, self.o, self.n, x=_x, deriv_order=deriv_order + ) + + return self._calc_linear_bend(vec_m, vec_o, vec_n, self._vec_r) + + def __repr__(self): + axis_str = "B" if self.axis == LinearBendType.BEND else "C" + return f"LinearBend{axis_str}({self.m}-{self.o}-{self.n}, D)" diff --git a/autode/opt/optimisers/crfo.py b/autode/opt/optimisers/crfo.py index 8dceaccf5..038f99666 100644 --- a/autode/opt/optimisers/crfo.py +++ b/autode/opt/optimisers/crfo.py @@ -1,4 +1,5 @@ -"""Constrained rational function optimisation +""" +Constrained rational function optimisation Notation follows: [1] J. Baker, J. Comput. Chem., 18, 8 1080 @@ -7,18 +8,14 @@ import numpy as np from typing import Union -from itertools import combinations from autode.log import logger -from autode.values import GradientRMS, Angle, Distance +from autode.values import GradientRMS, Distance from autode.opt.coordinates import CartesianCoordinates, DICWithConstraints -from autode.opt.coordinates.internals import PIC, AnyPIC +from autode.opt.coordinates.internals import AnyPIC from autode.opt.optimisers.rfo import RFOptimiser -from autode.opt.optimisers.hessian_update import BFGSDampedUpdate, NullUpdate -from autode.opt.coordinates.primitives import ( - PrimitiveDistance, - PrimitiveBondAngle, - PrimitiveDihedralAngle, - ConstrainedPrimitiveDistance, +from autode.opt.optimisers.hessian_update import ( + BFGSDampedUpdate, + BFGSSR1Update, ) @@ -41,7 +38,7 @@ def __init__( self.alpha = Distance(init_alpha, units="ang") assert self.alpha > 0 - self._hessian_update_types = [BFGSDampedUpdate, NullUpdate] + self._hessian_update_types = [BFGSDampedUpdate, BFGSSR1Update] def _step(self) -> None: """Partitioned rational function step""" @@ -125,15 +122,7 @@ def _build_internal_coordinates(self): ) cartesian_coords = CartesianCoordinates(self._species.coordinates) - primitives = self._primitives - - if len(primitives) < cartesian_coords.expected_number_of_dof: - logger.info( - "Had an incomplete set of primitives. Adding " - "additional distances" - ) - for i, j in combinations(range(self._species.n_atoms), 2): - primitives.append(PrimitiveDistance(i, j)) + primitives = AnyPIC.from_species(self._species) self._coords = DICWithConstraints.from_cartesian( x=cartesian_coords, primitives=primitives @@ -141,48 +130,6 @@ def _build_internal_coordinates(self): self._coords.zero_lagrangian_multipliers() return None - @property - def _primitives(self) -> PIC: - """Primitive internal coordinates in this molecule""" - assert self._species and self._species.graph - logger.info("Generating primitive internal coordinates") - graph = self._species.graph.copy() - - # Any distance constraints should also count as 'bonds' when forming - # the set of primitive internal coordinates, so that there is a - # single molecule if those distances are approaching dissociation - if self._species.constraints.distance is not None: - logger.info("Adding distance constraints as primitives") - for i, j in self._species.constraints.distance: - graph.add_edge(i, j) - - pic = AnyPIC() - - for i, j in sorted(graph.edges): - if ( - self._species.constraints.distance - and (i, j) in self._species.constraints.distance - ): - r = self._species.constraints.distance[(i, j)] - pic.append(ConstrainedPrimitiveDistance(i, j, value=r)) - - else: - pic.append(PrimitiveDistance(i, j)) - - for o in range(self._species.n_atoms): - for n, m in combinations(graph.neighbors(o), r=2): - pic.append(PrimitiveBondAngle(o=o, m=m, n=n)) - - if self._species.n_atoms > 2 and not self._species.is_planar(): - for dihedral in _dihedrals(self._species): - pic.append(dihedral) - - logger.info( - f"Using {pic.n_constrained} constraints in {len(pic)} " - f"primitive internal coordinates" - ) - return pic - def _lambda_p_from_eigvals_and_gradient( self, b: np.ndarray, f: np.ndarray ) -> float: @@ -218,31 +165,3 @@ def _lambda_n_from_eigvals_and_gradient( eigenvalues = np.linalg.eigvalsh(aug_h) return eigenvalues[0] - - -def _dihedrals(species): - """ - Iterator over the dihedrals in a species. Skipping those that contain - bond angles close to 180 degrees (tolerance <179) - """ - - for o, p in species.graph.edges: - for m in species.graph.neighbors(o): - if m == p: - continue - - if np.isclose(species.angle(m, o, p), Angle(np.pi), atol=0.04): - continue # Don't add potentially ill-defined dihedrals - - for n in species.graph.neighbors(p): - if n == o: - continue - - # avoid triangular rings like cyclopropane - if m == n: - continue - - if np.isclose(species.angle(o, p, n), Angle(np.pi), atol=0.04): - continue - - yield PrimitiveDihedralAngle(m, o, p, n) diff --git a/autode/opt/optimisers/trm.py b/autode/opt/optimisers/trm.py index 32b674143..a44a48f73 100644 --- a/autode/opt/optimisers/trm.py +++ b/autode/opt/optimisers/trm.py @@ -130,27 +130,6 @@ def _initialise_species_and_method( ), "HybridTRMOptimiser cannot work with constraints!" return None - def _build_internal_coordinates(self) -> None: - """ "Build delocalised internal coordinates""" - if self._species is None: - raise RuntimeError( - "Cannot set initial coordinates. No species set" - ) - - cart_coords = CartesianCoordinates(self._species.coordinates) - primitives = self._primitives - - if len(primitives) < cart_coords.expected_number_of_dof: - logger.info( - "Had an incomplete set of primitives. Adding " - "additional distances" - ) - for i, j in combinations(range(self._species.n_atoms), 2): - primitives.append(PrimitiveDistance(i, j)) - - self._coords = DIC.from_cartesian(x=cart_coords, primitives=primitives) - return None - def _step(self) -> None: """ Hybrid RFO/TRM step; if the RFO step is larger than the trust diff --git a/autode/values.py b/autode/values.py index f36196276..ea771c631 100644 --- a/autode/values.py +++ b/autode/values.py @@ -240,6 +240,10 @@ def __radd__(self, other) -> TypeValue: def __sub__(self, other) -> TypeValue: return self.__add__(-other) + def __neg__(self) -> TypeValue: + """Unary negation operation""" + return self._like_self_from_float(-float(self)) + def __floordiv__(self, other) -> Union[float, TypeValue]: x = float(self) // self._other_same_units(other) return x if isinstance(other, Value) else self._like_self_from_float(x) diff --git a/doc/changelog.rst b/doc/changelog.rst index 1d2b97732..3113ebb9d 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -8,6 +8,8 @@ Changelog Functionality improvements ************************** - Replaces hard-coded derivatives for primitive internal coordinates with automatic differentiation +- More comprehensive primitive internal coordinate generation +- Adds the capacity to handle linear molecules in internal coordinate system Bug Fixes ********* diff --git a/tests/test_opt/molecules.py b/tests/test_opt/molecules.py index 92936dfcc..43956fc06 100644 --- a/tests/test_opt/molecules.py +++ b/tests/test_opt/molecules.py @@ -6,6 +6,17 @@ def h_atom(): return Molecule(atoms=[Atom("H")], mult=2) +def acetylene_mol(): + return Molecule( + atoms=[ + Atom("C", 0.3800, -0.2049, -0.4861), + Atom("C", -0.3727, 0.1836, 0.4744), + Atom("H", 1.0156, -0.6055, -1.2353), + Atom("H", -0.9992, 0.5945, 1.1572), + ] + ) + + def methane_mol(): return Molecule( atoms=[ @@ -41,3 +52,37 @@ def h2o2_mol(): Atom("H", 0.58605, 0.91107, 0.59006), ] ) + + +def feco5_mol(): + return Molecule( + atoms=[ + Atom("O", -1.5139, -3.5069, -0.3015), + Atom("C", -0.8455, -2.5982, -0.3019), + Atom("Fe", 0.3545, -0.9664, -0.3020), + Atom("C", 1.4039, 0.5711, -0.2183), + Atom("O", 2.2216, 1.5751, -0.2992), + Atom("C", 1.2470, -1.6232, 1.3934), + Atom("O", 1.7445, -1.9892, 2.3373), + Atom("C", 1.0993, -1.3993, -2.0287), + Atom("O", 1.5042, -1.8102, -3.1145), + Atom("C", -1.2751, 0.2315, -0.1941), + Atom("O", -2.1828, 0.8986, -0.1345), + ] + ) + + +def cumulene_mol(): + return Molecule( + atoms=[ + Atom("C", -1.3968, 3.7829, 0.9582), + Atom("C", -0.8779, 2.7339, 0.6902), + Atom("C", -1.9158, 4.8319, 1.2262), + Atom("C", -2.4770, 5.9661, 1.5160), + Atom("C", -0.3167, 1.5996, 0.4004), + Atom("H", -2.1197, 6.8879, 1.0720), + Atom("H", -3.3113, 6.0085, 2.2063), + Atom("H", 0.1364, 1.4405, -0.5711), + Atom("H", -0.2928, 0.7947, 1.1256), + ] + ) diff --git a/tests/test_opt/test_coordiantes.py b/tests/test_opt/test_coordiantes.py index 9fe31f72e..eb8e21040 100644 --- a/tests/test_opt/test_coordiantes.py +++ b/tests/test_opt/test_coordiantes.py @@ -1,12 +1,26 @@ +import itertools + import pytest import numpy as np -from .molecules import h2, methane_mol, water_mol, h2o2_mol +from .molecules import ( + h2, + methane_mol, + water_mol, + h2o2_mol, + feco5_mol, + cumulene_mol, + acetylene_mol, +) +from autode.utils import work_in_tmp_dir from autode.atoms import Atom from autode.species.molecule import Molecule from autode.values import Angle from autode.exceptions import CoordinateTransformFailed -from autode.opt.coordinates.base import CartesianComponent -from autode.opt.coordinates.internals import PrimitiveInverseDistances, PIC +from autode.opt.coordinates.internals import ( + PrimitiveInverseDistances, + PIC, + AnyPIC, +) from autode.opt.coordinates.cartesian import CartesianCoordinates from autode.opt.coordinates.dic import DIC from autode.opt.coordinates.primitives import ( @@ -16,6 +30,9 @@ PrimitiveBondAngle, ConstrainedPrimitiveBondAngle, PrimitiveDihedralAngle, + PrimitiveLinearAngle, + PrimitiveDummyLinearAngle, + LinearBendType, ) @@ -527,12 +544,21 @@ def test_pic_b_no_primitives(): c._calc_B(np.arange(6, dtype=float).reshape(2, 3)) -def test_pic_append_type_checking(): - c = PIC() - # append should check for primitive type - c.append(PrimitiveDistance(0, 1)) +def test_pic_add_sanity_checking(): + c = AnyPIC() + # pic add should check for primitive type + c.add(PrimitiveDistance(0, 1)) with pytest.raises(AssertionError): - c.append(3) + c.add(3) + + # pic append is disallowed + with pytest.raises(NotImplementedError, match="Please use PIC.add()"): + c.append(PrimitiveDistance(0, 1)) + + # pic should not allow duplicate coordinates to be added + assert len(c) == 1 + c.add(PrimitiveDistance(1, 0)) + assert len(c) == 1 def test_constrained_distance_satisfied(): @@ -558,7 +584,7 @@ def numerical_derivative(a, b, h=1e-8): m = water_mol() init_coords = m.coordinates.copy() - angle = PrimitiveBondAngle(0, 1, 2) + angle = PrimitiveBondAngle(1, 0, 2) derivs = angle.derivative(init_coords) for atom_idx in (0, 1, 2): for component in (0, 1, 2): @@ -570,8 +596,8 @@ def numerical_derivative(a, b, h=1e-8): def test_angle_primitive_equality(): - assert PrimitiveBondAngle(0, 1, 2) == PrimitiveBondAngle(0, 2, 1) - assert PrimitiveBondAngle(0, 1, 2) != PrimitiveBondAngle(2, 1, 0) + assert PrimitiveBondAngle(1, 0, 2) == PrimitiveBondAngle(2, 0, 1) + assert PrimitiveBondAngle(1, 0, 2) != PrimitiveBondAngle(1, 2, 0) def test_dihedral_value(): @@ -612,6 +638,40 @@ def test_dihedral_equality(): ) +def test_linear_angle(): + acetylene = Molecule( + atoms=[ + Atom("C", 0.35540, -0.20370, -0.44810), + Atom("C", -0.37180, 0.21470, 0.40200), + Atom("H", 1.01560, -0.60550, -1.23530), + Atom("H", -0.99920, 0.59450, 1.15720), + ] + ) + x = CartesianCoordinates(acetylene.coordinates) + angle = PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.BEND) + assert angle._vec_r is None + _ = angle(x) + assert angle._vec_r is not None + old_r_vec = angle._vec_r + # the dummy atom should not change after the first call + _ = angle(x) + _ = angle(x) + assert angle._vec_r is old_r_vec + + axis_vec = np.array(np.array(angle._vec_r._data) - x.reshape(-1, 3)[1]) + m_n_vec = acetylene.coordinates[0] - acetylene.coordinates[1] + assert abs(np.dot(axis_vec, m_n_vec)) < 0.001 + + # check that the linear bond complement does not have the same value + angle2 = PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.COMPLEMENT) + assert angle != angle2 + assert not np.isclose(angle(x), angle2(x), rtol=1e-3) + + # for linear angle, swapping the end points changes the definition + angle3 = PrimitiveDummyLinearAngle(3, 1, 0, LinearBendType.BEND) + assert angle3 != angle + + def test_primitives_consistent_with_mol_values(): # test that the primitive values are the same as the mol.distance etc. h2o2 = h2o2_mol() @@ -620,14 +680,14 @@ def test_primitives_consistent_with_mol_values(): assert np.isclose(dist(coords), h2o2.distance(0, 1), rtol=1e-8) invdist = PrimitiveInverseDistance(1, 2) assert np.isclose(invdist(coords), 1 / h2o2.distance(1, 2), rtol=1e-8) - ang = PrimitiveBondAngle(2, 0, 1) # bond is defined in a different way + ang = PrimitiveBondAngle(0, 2, 1) assert np.isclose(ang(coords), h2o2.angle(0, 2, 1), rtol=1e-8) dihedral = PrimitiveDihedralAngle(2, 0, 1, 3) assert np.isclose(dihedral(coords), h2o2.dihedral(2, 0, 1, 3), rtol=1e-8) # fmt: off -dihedral_mols = [ +extra_mols = [ Molecule( atoms=[ Atom("C", 0.63365, 0.11934, -0.13163), @@ -643,18 +703,29 @@ def test_primitives_consistent_with_mol_values(): Atom("H", 1.28230, -0.63391, -0.54779), Atom("H", -1.08517, -1.07984, -0.05599), ] - ) # for testing dihedral derivatives over zero + ), # for testing dihedral derivatives over zero + Molecule( + atoms=[ + Atom("C", 0.35540, -0.20370, -0.44810), + Atom("C", -0.37180, 0.21470, 0.40200), + Atom("H", 1.01560, -0.60550, -1.23530), + Atom("H", -0.99920, 0.59450, 1.15720), + ] + ), + feco5_mol(), # for testing linear angles ] test_mols = [ h2o2_mol(), h2o2_mol(), water_mol(), - water_mol(), water_mol(), *dihedral_mols + water_mol(), water_mol(), *extra_mols ] test_prims = [ - PrimitiveDihedralAngle(2, 0, 1, 3), PrimitiveBondAngle(2, 0, 1), - PrimitiveBondAngle(0, 1, 2), PrimitiveDistance(0, 1), + PrimitiveDihedralAngle(2, 0, 1, 3), PrimitiveBondAngle(0, 2, 1), + PrimitiveBondAngle(1, 0, 2), PrimitiveDistance(0, 1), PrimitiveInverseDistance(0, 1), PrimitiveDihedralAngle(2, 0, 1, 3), - PrimitiveDihedralAngle(2, 0, 1, 3) + PrimitiveDihedralAngle(2, 0, 1, 3), + PrimitiveDummyLinearAngle(0, 1, 3, LinearBendType.BEND), + PrimitiveLinearAngle(2, 3, 4, 8, LinearBendType.BEND), ] # fmt: on @@ -706,9 +777,12 @@ def test_repr(): PrimitiveInverseDistance(0, 1), PrimitiveDistance(0, 1), ConstrainedPrimitiveDistance(0, 1, value=1e-3), - PrimitiveBondAngle(0, 1, 2), - ConstrainedPrimitiveBondAngle(0, 1, 2, value=1.0), + PrimitiveBondAngle(1, 0, 2), + ConstrainedPrimitiveBondAngle(1, 0, 2, value=1.0), PrimitiveDihedralAngle(0, 1, 2, 3), + PrimitiveLinearAngle(0, 1, 2, 3, LinearBendType.BEND), + PrimitiveLinearAngle(0, 1, 2, 3, LinearBendType.COMPLEMENT), + PrimitiveDummyLinearAngle(0, 1, 2, LinearBendType.BEND), ] for p in prims: @@ -733,7 +807,7 @@ def test_dic_large_step_allowed_unconverged_back_transform(): def test_constrained_angle_delta(): - q = ConstrainedPrimitiveBondAngle(0, 1, 2, value=np.pi) + q = ConstrainedPrimitiveBondAngle(1, 0, 2, value=np.pi) mol = water_mol() theta = mol.angle(1, 0, 2) x = CartesianCoordinates(mol.coordinates) @@ -742,8 +816,8 @@ def test_constrained_angle_delta(): def test_constrained_angle_equality(): - a = ConstrainedPrimitiveBondAngle(0, 1, 2, value=np.pi) - b = ConstrainedPrimitiveBondAngle(0, 2, 1, value=np.pi) + a = ConstrainedPrimitiveBondAngle(1, 0, 2, value=np.pi) + b = ConstrainedPrimitiveBondAngle(2, 0, 1, value=np.pi) assert a == b @@ -757,3 +831,121 @@ def test_dics_cannot_be_built_with_incomplete_primitives(): with pytest.raises(RuntimeError): _ = DIC.from_cartesian(x=x, primitives=primitives) + + +def test_pic_generation_linear_angle_ref(): + # Fe(CO)5 with linear Fe-C-O bonds + m = feco5_mol() + pic = AnyPIC.from_species(m) + + # check that there are no duplicates + assert not any(ic1 == ic2 for ic1, ic2 in itertools.combinations(pic, r=2)) + # check that linear bends use reference atoms, not dummy + assert not any(isinstance(ic, PrimitiveDummyLinearAngle) for ic in pic) + assert PrimitiveLinearAngle(4, 3, 2, 8, LinearBendType.BEND) in pic + # for C-Fe-C, only one out-of-plane dihedral should be present + assert PrimitiveDihedralAngle(3, 5, 2, 1) in pic + assert sum(isinstance(ic, PrimitiveDihedralAngle) for ic in pic) == 1 + # check degrees of freedom = 3N - 6 + _ = pic(m.coordinates.flatten()) + assert np.linalg.matrix_rank(pic.B) == 3 * m.n_atoms - 6 + + +def test_pic_generation_linear_angle_dummy(): + # acetylene molecule + mol = acetylene_mol() + pic = AnyPIC.from_species(mol) + + # there should not be any usual bond angles + assert not any(isinstance(ic, PrimitiveBondAngle) for ic in pic) + # there should not be any linear angles with reference atom + assert not any(isinstance(ic, PrimitiveLinearAngle) for ic in pic) + # there should be linear angles with dummy + assert any(isinstance(ic, PrimitiveDummyLinearAngle) for ic in pic) + + # degrees of freedom = 3N - 5 for linear molecules + _ = pic(mol.coordinates.flatten()) + assert np.linalg.matrix_rank(pic.B) == 3 * mol.n_atoms - 5 + + +@work_in_tmp_dir() +def test_pic_generation_disjoint_graph(): + # the algorithm should fully connect the graph + xyz_string = ( + "16\n\n" + "C -0.00247 1.65108 0.05872\n" + "C 1.19010 1.11169 0.27709\n" + "C 1.58519 -0.30014 0.31049\n" + "C 0.05831 -1.54292 -0.45110\n" + "C -1.18798 -1.04262 0.13551\n" + "C -1.28206 0.99883 -0.23631\n" + "H -0.07432 2.73634 0.08639\n" + "H 2.01755 1.78921 0.47735\n" + "H 1.70503 -0.70916 1.30550\n" + "H 2.40398 -0.55376 -0.34855\n" + "H 0.44229 -2.48695 -0.08638\n" + "H 0.15289 -1.41865 -1.51944\n" + "H -1.25410 -1.13318 1.21833\n" + "H -2.09996 -1.35918 -0.36715\n" + "H -2.09462 1.29055 0.41495\n" + "H -1.56001 1.00183 -1.28217\n" + ) + with open("diels_alder_complex.xyz", "w") as fh: + fh.write(xyz_string) + + mol = Molecule("diels_alder_complex.xyz") + assert not mol.graph.is_connected + pic = AnyPIC.from_species(mol) + + # shortest bond is between 4, 5 which should also generate angle, torsion + assert PrimitiveDistance(4, 5) in pic + assert PrimitiveBondAngle(3, 4, 5) in pic + assert PrimitiveBondAngle(4, 5, 0) in pic + assert PrimitiveDihedralAngle(3, 4, 5, 0) in pic + + # the other distance between fragments is 2, 3 which should not be connected + assert PrimitiveDistance(2, 3) not in pic + assert PrimitiveBondAngle(1, 2, 3) not in pic + # check degrees of freedom = 3N - 6 + _ = pic(mol.coordinates.flatten()) + assert np.linalg.matrix_rank(pic.B) == 3 * mol.n_atoms - 6 + + # if the bond between 2, 3 is made into a constraint, it will generate angles + mol.constraints.distance = {(2, 3): mol.distance(2, 3)} + pic = AnyPIC.from_species(mol) + assert ConstrainedPrimitiveDistance(2, 3, mol.distance(2, 3)) in pic + assert PrimitiveBondAngle(1, 2, 3) in pic + + +def test_pic_generation_chain_dihedrals(): + # extra dihedrals are needed for ends of linear chains like allene + cumulene = cumulene_mol() + pic = AnyPIC.from_species(cumulene) + + assert PrimitiveDihedralAngle(5, 3, 4, 8) in pic + assert PrimitiveDihedralAngle(6, 3, 4, 7) in pic + assert PrimitiveDihedralAngle(8, 4, 3, 6) in pic + assert PrimitiveDihedralAngle(7, 4, 3, 6) in pic + + # check that the 3N-6 degrees of freedom are maintained + _ = pic(cumulene.coordinates.flatten()) + assert np.linalg.matrix_rank(pic.B) == 3 * cumulene.n_atoms - 6 + + +def test_pic_generation_square_planar(): + ptcl4 = Molecule( + atoms=[ + Atom("Pt", -0.1467, -0.2594, -0.0294), + Atom("Cl", -0.4597, -2.5963, -0.0523), + Atom("Cl", 2.1804, -0.5689, -0.2496), + Atom("Cl", -2.4738, 0.0501, 0.1908), + Atom("Cl", 0.1663, 2.0776, -0.0066), + ], + charge=-2, + ) + + # for sq planar, out-of-plane dihedrals are needed to have + # all degrees of freedom + pic = AnyPIC.from_species(ptcl4) + _ = pic(ptcl4.coordinates.flatten()) + assert np.linalg.matrix_rank(pic.B) == 3 * ptcl4.n_atoms - 6 diff --git a/tests/test_opt/test_crfo.py b/tests/test_opt/test_crfo.py index 50854b689..f040278a8 100644 --- a/tests/test_opt/test_crfo.py +++ b/tests/test_opt/test_crfo.py @@ -11,7 +11,7 @@ from autode.opt.coordinates import CartesianCoordinates, DICWithConstraints from autode.opt.coordinates.primitives import PrimitiveDihedralAngle from autode.utils import work_in_tmp_dir -from .molecules import h2o2_mol +from .molecules import h2o2_mol, acetylene_mol, feco5_mol, cumulene_mol from ..testutils import requires_working_xtb_install @@ -116,7 +116,8 @@ def test_primitive_projection_discard(): r_initial = optimiser._species.distance(0, 1) x = CartesianCoordinates(optimiser._species.coordinates) - s = DICWithConstraints.from_cartesian(x, optimiser._primitives) + optimiser._build_internal_coordinates() + s = optimiser._coords assert len(s) == 3 # Shift on the first couple of DIC but nothing on the final one @@ -160,7 +161,7 @@ def test_xtb_opt_with_distance_constraint(): assert np.isclose(water.distance(0, 1), 0.99, atol=0.01) - CRFOptimiser.optimise(species=water, method=XTB()) + CRFOptimiser.optimise(species=water, method=XTB(), etol=1e-6) # Optimisation should generate an O-H distance *very* close to 1.1 Å assert np.isclose(water.distance(0, 1).to("Å"), 1.1, atol=1e-4) @@ -199,25 +200,25 @@ def test_baker1997_example(): r1 = prim.ConstrainedPrimitiveDistance(0, 1, value=1.5) r2 = prim.ConstrainedPrimitiveDistance(3, 4, value=2.5) theta = prim.ConstrainedPrimitiveBondAngle( - 1, 0, 5, value=val.Angle(123.0, "º").to("rad") + 0, 1, 5, value=val.Angle(123.0, "º").to("rad") ) pic = PIC(r1, r2, theta) for pair in ((2, 0), (3, 0), (4, 1), (5, 1)): - pic.append(prim.PrimitiveDistance(*pair)) + pic.add(prim.PrimitiveDistance(*pair)) for triple in ( - (1, 0, 3), - (1, 0, 3), - (2, 0, 3), - (4, 1, 0), - (5, 1, 0), - (4, 1, 5), + (0, 1, 3), + (0, 1, 3), + (0, 2, 3), + (1, 4, 0), + (1, 5, 0), + (1, 4, 5), ): - pic.append(prim.PrimitiveBondAngle(*triple)) + pic.add(prim.PrimitiveBondAngle(*triple)) for quadruple in ((4, 1, 0, 2), (4, 1, 0, 3), (5, 1, 0, 2), (5, 1, 0, 3)): - pic.append(prim.PrimitiveDihedralAngle(*quadruple)) + pic.add(prim.PrimitiveDihedralAngle(*quadruple)) dic = DICWithConstraints.from_cartesian( x=CartesianCoordinates(c2h3f.coordinates), primitives=pic @@ -339,3 +340,41 @@ def test_linear_dihedrals_are_removed(): assert not any( isinstance(q, PrimitiveDihedralAngle) for q in dic.primitives ) + + +@requires_working_xtb_install +@work_in_tmp_dir() +def test_optimise_linear_molecule(): + mol = acetylene_mol() + # the two H-C-C angles are almost linear + assert val.Angle(170, "deg") < mol.angle(0, 1, 3) < val.Angle(176, "deg") + assert val.Angle(170, "deg") < mol.angle(2, 0, 1) < val.Angle(176, "deg") + opt = CRFOptimiser(maxiter=10, gtol=1e-4, etol=1e-5) + opt.run(mol, XTB()) + assert opt.converged + assert mol.angle(0, 1, 3) > val.Angle(179, "deg") + assert mol.angle(2, 0, 1) > val.Angle(179, "deg") + + +@requires_working_xtb_install +@work_in_tmp_dir() +def test_optimise_linear_bend_with_ref(): + mol = feco5_mol() + # the Fe-C-O angle are manually deviated + assert val.Angle(170, "deg") < mol.angle(2, 3, 4) < val.Angle(176, "deg") + # large molecule so allow few iters, no need to converge fully + opt = CRFOptimiser(maxiter=15, gtol=1e-5, etol=1e-6) + opt.run(mol, XTB()) + assert mol.angle(2, 3, 4) > val.Angle(178, "deg") + + +@requires_working_xtb_install +@work_in_tmp_dir() +def test_optimise_chain_dihedrals(): + mol = cumulene_mol() + assert abs(mol.dihedral(6, 3, 4, 8)) < val.Angle(40, "deg") + opt = CRFOptimiser(maxiter=15, gtol=1e-4, etol=1e-4) + opt.run(mol, XTB()) + # 5-C chain, should be close to 90 degrees + assert abs(mol.dihedral(6, 3, 4, 8)) > val.Angle(85, "deg") + assert abs(mol.dihedral(6, 3, 4, 8)) < val.Angle(95, "deg") From 7bf8214629d8b21cf94308a5120cfef944de2553 Mon Sep 17 00:00:00 2001 From: Daniel Hollas Date: Tue, 20 Feb 2024 19:44:37 +0000 Subject: [PATCH 4/5] Speedup package import time (#319) * Speedup package import time * Add changelog entry * WIP: Add test --- README.md | 1 + autode/neb/original.py | 5 ++++- autode/path/path.py | 3 ++- autode/pes/pes_nd.py | 5 ++++- doc/changelog.rst | 3 +++ tests/test_import.py | 32 ++++++++++++++++++++++++++++++++ 6 files changed, 46 insertions(+), 3 deletions(-) create mode 100644 tests/test_import.py diff --git a/README.md b/README.md index 5d402b4b3..3b8e65086 100644 --- a/README.md +++ b/README.md @@ -94,3 +94,4 @@ If **autodE** is used in a publication please consider citing the [paper](https: - Domen Pregeljc ([@dpregeljc](https://github.com/dpregeljc)) - Jonathon Vandezande ([@jevandezande](https://github.com/jevandezande)) - Shoubhik Maiti ([@shoubhikraj](https://github.com/shoubhikraj)) +- Daniel Hollas ([@danielhollas](https://github.com/danielhollas)) diff --git a/autode/neb/original.py b/autode/neb/original.py index b0fd5d1e4..58060ddf9 100644 --- a/autode/neb/original.py +++ b/autode/neb/original.py @@ -3,7 +3,6 @@ Henkelman and H. J ́onsson, J. Chem. Phys. 113, 9978 (2000) """ import numpy as np -import matplotlib.pyplot as plt from typing import Optional, Sequence, List, Any, TYPE_CHECKING, Union, Type from copy import deepcopy @@ -361,6 +360,8 @@ def plot_energies( self, save=False, name="None", color=None, xlabel="NEB coordinate" ): """Plot the NEB surface""" + import matplotlib.pyplot as plt + blues = plt.get_cmap("Blues") color = ( @@ -673,6 +674,8 @@ def calculate( etol_per_image: Energy tolerance per image to use in the L-BFGS-B minimisation """ + import matplotlib.pyplot as plt + self.print_geometries(name=f"{name_prefix}neb_init") # Calculate energy on the first and final points as these will not be recalc-ed diff --git a/autode/path/path.py b/autode/path/path.py index 81ab756fe..64939cc6c 100644 --- a/autode/path/path.py +++ b/autode/path/path.py @@ -1,5 +1,4 @@ import numpy as np -import matplotlib.pyplot as plt from autode.species import Species from autode.input_output import atoms_to_xyz_file @@ -147,6 +146,8 @@ def plot_energies( self, save: bool, name: str, color: str, xlabel: str ) -> None: """Plot this path""" + import matplotlib.pyplot as plt + if len(self) == 0 or any(item.energy is None for item in self): logger.error("Could not plot a surface, an energy was None") return diff --git a/autode/pes/pes_nd.py b/autode/pes/pes_nd.py index edc59332c..7741c599e 100644 --- a/autode/pes/pes_nd.py +++ b/autode/pes/pes_nd.py @@ -5,7 +5,6 @@ """ import numpy as np import itertools as it -import matplotlib.pyplot as plt from abc import ABC, abstractmethod from typing import ( @@ -196,6 +195,7 @@ def plot( units: Units of the surface. One of {'Ha', 'eV', 'kcal', 'kJ'} """ + import matplotlib.pyplot as plt if interp_factor < 0: raise ValueError( @@ -559,6 +559,8 @@ def _plot_1d(self, interp_factor: int, units_name: str) -> None: interp_factor: units_name: """ + import matplotlib.pyplot as plt + r_x = self._rs[0] energies, units = self._energies, energy_unit_from_name(units_name) energies = units.times * (energies - np.min(energies)) @@ -605,6 +607,7 @@ def _plot_2d(self, interp_factor: int, units_name: str) -> None: interp_factor: units_name: """ + import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from matplotlib.ticker import FormatStrFormatter diff --git a/doc/changelog.rst b/doc/changelog.rst index 3113ebb9d..af8f7a3c5 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -15,6 +15,9 @@ Bug Fixes ********* - Fixes triangular rings being incorrectly treated as dihedral angles +Usability improvements/Changes +****************************** +- Faster import of autode package by lazily loading matplotlib 1.4.1 ------ diff --git a/tests/test_import.py b/tests/test_import.py new file mode 100644 index 000000000..db90f0fd3 --- /dev/null +++ b/tests/test_import.py @@ -0,0 +1,32 @@ +"""Tests for the import speed of autode.""" +import sys + +import pytest + +SLOW_IMPORTS = ["matplotlib"] + + +@pytest.fixture +def unimport_slow_imports(): + """Remove modules in ``SLOW_IMPORTS`` from ``sys.modules``.""" + for module in SLOW_IMPORTS: + if module in sys.modules: + del sys.modules[module] + + +@pytest.mark.usefixtures("unimport_slow_imports") +def test_slow_imports_during_tab_completion(): + """Check that importing autode does not import certain python modules that would make import slow.""" + + # Let's double check that the undesired imports are not already loaded + for modulename in SLOW_IMPORTS: + assert ( + modulename not in sys.modules + ), f"Module `{modulename}` was not properly unloaded" + + import autode + + for modulename in SLOW_IMPORTS: + assert ( + modulename not in sys.modules + ), f"Detected loaded module {modulename} after autode import" From ea815e78df66fa13784b9ef96cceca5b0b617240 Mon Sep 17 00:00:00 2001 From: Tom Young <39765193+t-young31@users.noreply.github.com> Date: Mon, 26 Feb 2024 08:23:16 +0000 Subject: [PATCH 5/5] Allow configurable ORCA output extensions (#323) * allow configurable * update changelog [skip actions] --- autode/config.py | 3 +++ autode/wrappers/ORCA.py | 2 +- doc/changelog.rst | 1 + 3 files changed, 5 insertions(+), 1 deletion(-) diff --git a/autode/config.py b/autode/config.py index fd534e49b..48f252247 100644 --- a/autode/config.py +++ b/autode/config.py @@ -198,6 +198,9 @@ class ORCA: # # Path can be unset and will be assigned if it can be found in $PATH path = None + # + # File extensions to copy when a calculation completes + copied_output_exts = [".out", ".hess", ".xyz", ".inp", ".pc"] optts_block = ( "\n%geom\n" diff --git a/autode/wrappers/ORCA.py b/autode/wrappers/ORCA.py index 933619229..1732de645 100644 --- a/autode/wrappers/ORCA.py +++ b/autode/wrappers/ORCA.py @@ -220,7 +220,7 @@ def version_in(self, calc: "CalculationExecutor") -> str: def execute(self, calc): @work_in_tmp_dir( filenames_to_copy=calc.input.filenames, - kept_file_exts=(".out", ".hess", ".xyz", ".inp", ".pc"), + kept_file_exts=Config.ORCA.copied_output_exts, ) def execute_orca(): run_external( diff --git a/doc/changelog.rst b/doc/changelog.rst index af8f7a3c5..924a896a0 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -18,6 +18,7 @@ Bug Fixes Usability improvements/Changes ****************************** - Faster import of autode package by lazily loading matplotlib +- ORCA output files copied after a calculation are configurable 1.4.1 ------