"""
A module for manipulating vectors
"""
import math
import numpy as np
from arc.common import logger
from arc.exceptions import VectorsError
from arc.molecule.molecule import Molecule
from arc.species._zmat_kernels import lib as _ck, available as _ck_available
[docs]
def get_normal(v1: list[float],
v2: list[float],
) -> list[float]:
"""
Calculate a normal vector using cross multiplication.
Args:
v1 (list): Vector 1.
v2 (list): Vector 2.
Returns: list
A normal unit vector to v1 and v2.
"""
normal = [v1[1] * v2[2] - v2[1] * v1[2], - v1[0] * v2[2] + v2[0] * v1[2], v1[0] * v2[1] - v2[0] * v1[1]]
return unit_vector(normal)
[docs]
def get_perpendicular_unit_vector(vector: list[float]) -> list[float]:
"""
Get a deterministic unit vector perpendicular to ``vector``.
Crossing ``vector`` with the Cartesian axis it is least aligned with can never be degenerate:
the smallest component of a vector of length L is at most L / sqrt(3), so the cross product has
a length of at least L * sqrt(2/3). Picking the axis by magnitude rather than arbitrarily keeps
the result reproducible for a given input.
Args:
vector (list): The vector to find a perpendicular of. Need not be normalized.
Raises:
VectorsError: If ``vector`` is not a finite, non-zero vector of length three.
ZeroDivisionError: If ``vector`` is non-zero but its squared length underflows to zero,
which happens below a magnitude of about 1e-162.
Returns: list
A unit vector perpendicular to ``vector``.
"""
if len(vector) != 3:
raise VectorsError(f'vector must be of length three, got {len(vector)}.')
if not all(math.isfinite(component) for component in vector):
raise VectorsError(f'vector must be finite, got {vector}.')
if not any(vector):
raise VectorsError(f'vector must be non-zero, got {vector}.')
reference = [0.0, 0.0, 0.0]
reference[min(range(3), key=lambda index: abs(vector[index]))] = 1.0
return get_normal(vector, reference)
[docs]
def get_perpendicular_axes(vector: list[float],
tol: float = 1e-8,
) -> tuple[list[float], list[float]]:
"""
Get the two unit vectors that complete a right-handed orthonormal frame with a given direction.
Writing ``e_x`` for the unit vector along ``vector``, the returned ``e_y`` and ``e_z`` satisfy
``e_x`` x ``e_y`` = ``e_z``. Their common phase about ``e_x`` is arbitrary: it is fixed by
``get_perpendicular_unit_vector()``, so it is reproducible for a given input but carries no
meaning of its own. Angles measured against the returned pair therefore retain their cyclic
order about ``e_x`` while their absolute values do not.
Args:
vector (list): The vector defining the first axis of the frame. Need not be normalized.
tol (float, optional): The minimal accepted length of ``vector``.
Raises:
VectorsError: If ``vector`` is not of length three, is not finite, is shorter than ``tol``,
or is long enough for its squared length to overflow.
Returns: tuple
The ``e_y`` and ``e_z`` unit vectors.
"""
if len(vector) != 3:
raise VectorsError(f'vector must be of length three, got {len(vector)}.')
if not all(math.isfinite(component) for component in vector):
raise VectorsError(f'vector must be finite, got {vector}.')
if get_vector_length(vector) < tol:
raise VectorsError(f'vector must be longer than {tol}, got {vector}.')
e_x = unit_vector(vector)
e_z = get_perpendicular_unit_vector(e_x)
return get_normal(e_z, e_x), e_z
[docs]
def get_angle(v1: list[float],
v2: list[float],
units: str = 'rads',
) -> float:
"""
Calculate the angle between two vectors.
Args:
v1 (list): Vector 1.
v2 (list): Vector 2.
units (str, optional): The desired units, either 'rads' for radians, or 'degs' for degrees.
Raises:
VectorsError: If ``v1`` and ``v2`` are of different lengths.
Returns: float
The angle between ``v1`` and ``v2`` in the desired units.
"""
if len(v1) != len(v2):
raise VectorsError(f'v1 and v2 must be the same length, got {len(v1)} and {len(v2)}.')
v1_u, v2_u = unit_vector(v1), unit_vector(v2)
conversion = 180 / math.pi if 'degs' in units else 1
return float(np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0)) * conversion)
[docs]
def get_dihedral(v1: list[float],
v2: list[float],
v3: list[float],
units: str = 'degs',
tol: float = 1e-8,
) -> float:
"""
Calculate the dihedral angle between three vectors.
``v2`` connects between ``v1`` and ``v3``.
Args:
v1 (list): Vector 1.
v2 (list): Vector 2.
v3 (list): Vector 3.
units (str, optional): The desired units, either 'rads' for radians, or 'degs' for degrees.
tol (float, optional): The tolerance for zero-length vectors.
Raises:
VectorsError: If either ``v1`` or ``v2`` have lengths different from three.
Returns: float
The dihedral angle between ``v1`` and ``v2`` in the desired units.
"""
# 1) length check
if len(v1) != 3 or len(v2) != 3 or len(v3) != 3:
raise VectorsError(f'v1, v2, and v3 must each be length 3, got {len(v1)}, {len(v2)}, {len(v3)}.')
# 2) convert and zero-length check
v1 = np.array(v1, dtype=np.float64)
v2 = np.array(v2, dtype=np.float64)
v3 = np.array(v3, dtype=np.float64)
if np.linalg.norm(v1) == 0.0 or np.linalg.norm(v2) == 0.0 or np.linalg.norm(v3) == 0.0:
raise VectorsError('Input vectors must be non-zero.')
# 3) compute normals
v2_x_v1 = np.cross(v2, v1)
norm1 = np.linalg.norm(v2_x_v1)
if norm1 < tol:
return float('nan') # colinear
v2_x_v1 /= norm1
v3_x_v2 = np.cross(v3, v2)
norm2 = np.linalg.norm(v3_x_v2)
if norm2 < tol:
return float('nan')
v3_x_v2 /= norm2
# 4) angle between normals
cosine = np.clip(np.vdot(v2_x_v1, v3_x_v2), -1.0, 1.0)
dihedral = math.acos(cosine)
# 5) sign via triple product
if np.vdot(v2_x_v1, v3) > 0:
dihedral = 2 * math.pi - dihedral
# 6) unit conversion
conversion = 180 / math.pi if 'degs' in units else 1
return float(dihedral) * conversion
[docs]
def calculate_distance(coords: list | tuple | dict,
atoms: list,
index: int = 0,
) -> float:
"""
Calculate a distance.
Args:
coords (list, tuple, dict): The array-format or tuple-format coordinates, or the xyz dict.
atoms (list): The 2 atoms to calculate defining the vector for which the length will be calculated.
index (int, optional): Whether ``atoms`` is 0-indexed or 1-indexed (values are 0 or 1).
Raises:
VectorsError: If ``index`` is out of range, or ``atoms`` is of wrong length or has repeating indices.
TypeError: If ``coords`` is of wrong type.
Returns: float
The distance in the coords units.
"""
if isinstance(coords, dict) and 'coords' in coords:
coords = coords['coords']
if not isinstance(coords, (list, tuple)):
raise TypeError(f'coords must be a list or a tuple, got\n{coords}\nwhich is a {type(coords)}')
if index not in [0, 1]:
raise VectorsError(f'index must be either 0 or 1, got {index}')
if len(atoms) != 2:
raise VectorsError(f'distance atom list must be of length two, got {len(atoms)}')
if len(set(atoms)) < 2:
raise VectorsError(f'some atoms are repetitive: {atoms}')
new_atoms = list()
for atom in atoms:
if isinstance(atom, str) and 'X' in atom:
new_atoms.append(int(atom[1:]))
else:
new_atoms.append(atom)
if not all([isinstance(a, int) for a in new_atoms]):
raise VectorsError(f'all entries in atoms must be integers, got: {new_atoms} ({[type(a) for a in new_atoms]})')
new_atoms = [a - index for a in new_atoms] # convert 1-index to 0-index
if _ck_available:
c0, c1 = coords[new_atoms[0]], coords[new_atoms[1]]
return _ck.zmat_r_f32(c0[0], c0[1], c0[2], c1[0], c1[1], c1[2])
coords = np.asarray(coords, dtype=np.float32)
vector = coords[new_atoms[1]] - coords[new_atoms[0]]
return get_vector_length(vector)
[docs]
def calculate_angle(coords: list | tuple | dict,
atoms: list,
index: int = 0,
units: str = 'degs',
) -> float:
"""
Calculate an angle.
Args:
coords (list, tuple, dict): The array-format or tuple-format coordinates, or the xyz dict.
atoms (list): The 3 atoms defining the angle.
index (int, optional): Whether ``atoms`` is 0-indexed or 1-indexed (values are 0 or 1).
units (str, optional): The desired units, either 'rads' for radians, or 'degs' for degrees.
Raises:
VectorsError: If ``index`` is out of range, or ``atoms`` is of wrong length or has repeating indices.
TypeError: If ``coords`` is of wrong type.
Returns: float
The angle.
"""
if isinstance(coords, dict) and 'coords' in coords:
coords = coords['coords']
if not isinstance(coords, (list, tuple)):
raise TypeError(f'coords must be a list or a tuple, got\n{coords}\nwhich is a {type(coords)}')
if index not in [0, 1]:
raise VectorsError(f'index must be either 0 or 1, got {index}')
if len(atoms) != 3:
raise VectorsError(f'angle atom list must be of length three, got {len(atoms)}')
if len(set(atoms)) < 3:
raise VectorsError(f'some atoms are repetitive: {atoms}')
new_atoms = list()
for atom in atoms:
if isinstance(atom, str) and 'X' in atom:
new_atoms.append(int(atom[1:]))
else:
new_atoms.append(atom)
if not all([isinstance(a, int) for a in new_atoms]):
raise VectorsError(f'all entries in atoms must be integers, got: {new_atoms} ({[type(a) for a in new_atoms]})')
new_atoms = [a - index for a in new_atoms] # convert 1-index to 0-index
# Computed in float64, not float32. Coordinates used to be truncated to float32 here
# (dd452891, 2020), which cost up to ~8e-4 deg -- the same order as the 1e-3 deg
# DEFAULT_CONSOLIDATION_A_TOL, so the noise alone could decide whether two equal angles
# consolidate into one z-matrix parameter. The kernel below is exact and matches this
# path bit-for-bit; keeping both in float64 is what lets them agree, so that whether
# the .so happens to be compiled cannot change a result.
if _ck_available:
c0, c1 = coords[new_atoms[0]], coords[new_atoms[1]]
c2 = coords[new_atoms[2]]
deg = _ck.zmat_a(c0[0], c0[1], c0[2], c1[0], c1[1], c1[2], c2[0], c2[1], c2[2])
return deg if 'degs' in units else math.radians(deg)
coords = np.asarray(coords, dtype=np.float64)
v1 = coords[new_atoms[1]] - coords[new_atoms[0]]
v2 = coords[new_atoms[1]] - coords[new_atoms[2]]
return get_angle(v1, v2, units=units)
[docs]
def calculate_dihedral_angle(coords: list | tuple | dict,
torsion: list,
index: int = 0,
units: str = 'degs',
) -> float:
"""
Calculate a dihedral angle.
Args:
coords (list, tuple, dict): The array-format or tuple-format coordinates, or the xyz dict.
torsion (list): The 4 atoms defining the dihedral angle.
index (int, optional): Whether ``torsion`` is 0-indexed or 1-indexed (values are 0 or 1).
units (str, optional): The desired units, either 'rads' for radians, or 'degs' for degrees.
Raises:
VectorsError: If ``index`` is out of range, or ``torsion`` is of wrong length or has repeating indices.
TypeError: If ``coords`` is of wrong type.
Returns: float
The dihedral angle in a 0-360 degrees range.
"""
if isinstance(coords, dict) and 'coords' in coords:
coords = coords['coords']
if not isinstance(coords, (list, tuple)):
raise TypeError(f'coords must be a list or a tuple, got\n{coords}\nwhich is a {type(coords)}')
if index not in [0, 1]:
raise VectorsError(f'index must be either 0 or 1, got {index}')
if torsion is None:
raise VectorsError(f'torsion cannot be None')
if len(torsion) != 4:
raise VectorsError(f'torsion atom list must be of length four, got {len(torsion)}')
if len(set(torsion)) < 4:
raise VectorsError(f'some indices in torsion are repetitive: {torsion}')
new_torsion = list()
for atom in torsion:
if isinstance(atom, str) and 'X' in atom:
new_torsion.append(int(atom[1:]))
else:
new_torsion.append(atom)
if not all([isinstance(t, int) for t in new_torsion]):
raise VectorsError(f'all entries in torsion must be integers, got: {new_torsion} '
f'({[type(t) for t in new_torsion]})')
new_torsion = [t - index for t in new_torsion] # convert 1-index to 0-index if needed
if _ck_available:
c0, c1 = coords[new_torsion[0]], coords[new_torsion[1]]
c2, c3 = coords[new_torsion[2]], coords[new_torsion[3]]
deg = _ck.zmat_d_f32(c0[0], c0[1], c0[2], c1[0], c1[1], c1[2],
c2[0], c2[1], c2[2], c3[0], c3[1], c3[2])
if not math.isnan(deg):
return deg if 'degs' in units else math.radians(deg)
# The kernel collapses two distinct degenerate cases onto NaN, while
# get_dihedral() distinguishes them: a zero-length bond vector raises
# VectorsError, and only a collinear (yet non-degenerate) triplet returns
# NaN. Whether the .so happens to be compiled must not change which errors
# ARC reports, so re-classify here -- but do not recompute the angle: the
# two paths draw the collinearity threshold at slightly different places,
# and letting numpy re-answer a borderline triplet with a finite value
# would silently change which reference atoms zmat.py picks.
pts = np.asarray([coords[i] for i in new_torsion], dtype=np.float32)
if any(not np.any(pts[i + 1] - pts[i]) for i in range(3)):
raise VectorsError('Input vectors must be non-zero.')
return deg # NaN, in either unit
coords = np.asarray(coords, dtype=np.float32)
v1 = coords[new_torsion[1]] - coords[new_torsion[0]]
v2 = coords[new_torsion[2]] - coords[new_torsion[1]]
v3 = coords[new_torsion[3]] - coords[new_torsion[2]]
return get_dihedral(v1, v2, v3, units=units)
[docs]
def calculate_param(coords: list | tuple | dict,
atoms: list,
index: int = 0,
) -> float:
"""
Calculate a distance / angle / dihedral angle parameter.
Default units (deg for angles) are used.
Args:
coords (list, tuple, dict): The array-format or tuple-format coordinates, or the xyz dict.
atoms (list): The 2 atoms to calculate defining the vector for which the length will be calculated.
index (int, optional): Whether ``atoms`` is 0-indexed or 1-indexed (values are 0 or 1).
Returns: float
The calculated parameter.
"""
if len(atoms) == 2:
return calculate_distance(coords=coords, atoms=atoms, index=index)
elif len(atoms) == 3:
return calculate_angle(coords=coords, atoms=atoms, index=index)
elif len(atoms) == 4:
return calculate_dihedral_angle(coords=coords, torsion=atoms, index=index)
else:
raise ValueError(f'Expected 2, 3, or 4 indices in coords, got {len(coords)}: {coords}')
[docs]
def unit_vector(vector: list[float]) -> list[float]:
"""
Calculate a unit vector in the same direction as the input vector.
Args:
vector (list): The input vector.
Returns: list
The unit vector.
"""
length = get_vector_length(vector)
return [vi / length for vi in vector]
[docs]
def set_vector_length(vector: list[float],
length: float,
) -> list[float]:
"""
Set the length of a 3D vector.
Args:
vector (list): The vector to process.
length (float): The desired length to set.
Returns: list
A vector with the desired length.
"""
u = unit_vector(vector)
return [u[0] * length, u[1] * length, u[2] * length]
[docs]
def rotate_vector(point_a: list[float],
point_b: list[float],
normal: list[float],
theta: float,
) -> list[float]:
"""
Rotate a vector in 3D space around a given axis by a certain angle.
Inspired by https://stackoverflow.com/questions/6802577/rotation-of-3d-vector
Args:
point_a (list): The 3D coordinates of the starting point (point A) of the vector to be rotated.
point_b (list): The 3D coordinates of the ending point (point B) of the vector to be rotated.
normal (list): The axis to be rotated around.
theta (float): The degree in radians by which to rotate.
Returns: list
The rotated vector (the new coordinates for point B).
"""
normal = np.asarray(normal)
normal = normal / math.sqrt(np.dot(normal, normal)) # should *not* be replaced by an augmented assignment
a = math.cos(theta / 2.0)
b, c, d = -normal * math.sin(theta / 2.0)
aa, bb, cc, dd = a * a, b * b, c * c, d * d
bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
rotation_matrix = np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
[2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
[2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
vector = [b - a for b, a in zip(point_b, point_a)]
new_vector = np.dot(rotation_matrix, vector).tolist()
new_vector = [n_v + a for n_v, a in zip(new_vector, point_a)] # root the vector at the original starting point
return new_vector
[docs]
def apply_rodrigues_rotation(coords: list,
axis_origin: tuple | list,
axis_unit: tuple | list,
angle_rad: float,
indices: list[int],
) -> list:
"""
Rotate atoms at ``indices`` around an axis using Rodrigues' formula.
Modifies ``coords`` in-place and returns it.
Args:
coords: list of (x, y, z) tuples representing all atom positions.
axis_origin: a point on the rotation axis (e.g., coords of the near-side pivot).
axis_unit: unit vector along the rotation axis.
angle_rad: rotation angle in radians (positive = right-hand rule around axis_unit).
indices: 0-indexed atom indices to rotate.
Returns:
list: The same ``coords`` list with rotated atoms.
"""
c, s = math.cos(angle_rad), math.sin(angle_rad)
kx, ky, kz = axis_unit[0], axis_unit[1], axis_unit[2]
ox, oy, oz = axis_origin[0], axis_origin[1], axis_origin[2]
for idx in indices:
vx = coords[idx][0] - ox
vy = coords[idx][1] - oy
vz = coords[idx][2] - oz
kvx = ky * vz - kz * vy
kvy = kz * vx - kx * vz
kvz = kx * vy - ky * vx
kdv = kx * vx + ky * vy + kz * vz
one_minus_c = 1.0 - c
coords[idx] = (vx * c + kvx * s + kx * kdv * one_minus_c + ox,
vy * c + kvy * s + ky * kdv * one_minus_c + oy,
vz * c + kvz * s + kz * kdv * one_minus_c + oz)
return coords
[docs]
def get_vector(pivot: int,
anchor: int,
xyz: dict,
) -> list:
"""
Get a vector between two atoms in the molecule (pointing from pivot to anchor).
Args:
pivot (int): The 0-index of the pivotal atom around which groups are to be translated.
anchor (int): The 0-index of an additional atom in the molecule.
xyz (dict): The 3D coordinates of the molecule with the same atom order as in mol.
Returns: list
A vector pointing from the pivotal atom towards the anchor atom.
"""
px, py, pz = xyz['coords'][pivot]
ax, ay, az = xyz['coords'][anchor]
return [ax - px, ay - py, az - pz]
[docs]
def get_lp_vector(label: str,
mol: Molecule,
xyz: dict,
pivot: int,
) -> list:
"""
Get a vector from the pivotal atom in the molecule towards its lone electron pair (lp).
The approach is to reverse the average of the three unit vectors between the pivotal atom and its neighbors.
Args:
label (str): The species' label.
mol (Molecule): The 2D graph representation of the molecule.
xyz (dict): The 3D coordinates of the molecule with the same atom order as in mol.
pivot (int): The 0-index of the pivotal atom of interest.
Raises:
VectorsError: If the lp vector cannot be attained.
Returns: list
A unit vector pointing from the pivotal (nitrogen) atom towards its lone electron pairs orbital.
"""
neighbors, vectors = list(), list()
if not mol.atoms[pivot].is_nitrogen():
logger.warning(f'The get_lp_vector specializes in nitrogen atoms, got {mol.atoms[pivot].symbol} '
f'(atom number {pivot}) in species {label}.')
for atom in mol.atoms[pivot].edges.keys():
neighbors.append(mol.atoms.index(atom))
if len(neighbors) < 3:
# N will have 3, S may have more.
raise VectorsError(f'Can only get lp vector if the pivotal atom has at least three neighbors. '
f'Atom {mol.atoms[pivot]} in {label} has only {len(neighbors)}.')
for neighbor in neighbors:
# already taken in "reverse", pointing towards the lone pair
vectors.append(unit_vector(get_vector(pivot=neighbor, anchor=pivot, xyz=xyz)))
x = sum(v[0] for v in vectors)
y = sum(v[1] for v in vectors)
z = sum(v[2] for v in vectors)
return unit_vector([x, y, z])
[docs]
def get_vector_length(v: list[float]) -> float:
"""
Get the length of an ND vector
Args:
v (list): The vector.
Returns: float
The vector's length.
"""
return float(np.dot(v, v) ** 0.5)
[docs]
def get_delta_angle(a1: float,
a2: float,
) -> float:
"""
Get the difference between two (dihedral or regular) angles.
Examples::
3 - 1 = 2
1 - 3 = 2
1- 359 = 2
Args:
a1 (float): Angle 1 in degrees.
a2 (float): Angle 2 in degrees.
Returns: float
The difference between the angles in degrees.
"""
a1 %= 360
a2 %= 360
return min(abs(a1 - a2), abs(a1 + 360 - a2), abs(a1 - a2 - 360))