Skip to content

Commit

Permalink
Manual sync from SCM git to GitHub
Browse files Browse the repository at this point in the history
David Ormrod Morley: Merge pull request 'Add PLAMS unit tests for atom, bond and kftools SO--' (!14) from DavidOrmrodMorley/SCMSUITE-8674.atom_bond_unit_tests into trunk SO--
  • Loading branch information
dormrod committed Aug 6, 2024
1 parent 119890d commit a4d9060
Show file tree
Hide file tree
Showing 12 changed files with 812 additions and 73 deletions.
2 changes: 1 addition & 1 deletion interfaces/molecule/rdkit.py
Original file line number Diff line number Diff line change
Expand Up @@ -230,7 +230,7 @@ def plams_to_rd_bonds(bo):
raise exc
conf = Chem.Conformer()
for i, atom in enumerate(plams_mol.atoms):
xyz = Geometry.Point3D(atom._getx(), atom._gety(), atom._getz())
xyz = Geometry.Point3D(atom.x, atom.y, atom.z)
conf.SetAtomPosition(i, xyz)
rdmol.AddConformer(conf)
# REB: Assign all stereochemistry, if it wasn't already there
Expand Down
92 changes: 51 additions & 41 deletions mol/atom.py
Original file line number Diff line number Diff line change
@@ -1,13 +1,17 @@
import math

import numpy as np
from typing import Iterable, Union, List, Tuple, TYPE_CHECKING

from scm.plams.core.settings import Settings
from scm.plams.tools.periodic_table import PT
from scm.plams.tools.units import Units

__all__ = ["Atom"]

if TYPE_CHECKING:
from scm.plams.mol.bond import Bond


class Atom:
"""A class representing a single atom in three dimensional space.
Expand Down Expand Up @@ -128,70 +132,67 @@ def __iter__(self):
"""Iteration through atom yields coordinates. Thanks to that instances of |Atom| can be passed to any method requiring as an argument a point or a vector in 3D space."""
return iter(self.coords)

def _setx(self, value):
self.coords = (value, self.coords[1], self.coords[2])

def _sety(self, value):
self.coords = (self.coords[0], value, self.coords[2])

def _setz(self, value):
self.coords = (self.coords[0], self.coords[1], value)

def _getx(self):
@property
def x(self) -> float:
return self.coords[0]

def _gety(self):
@x.setter
def x(self, value: float) -> None:
self.coords = (value, self.coords[1], self.coords[2])

@property
def y(self) -> float:
return self.coords[1]

def _getz(self):
@y.setter
def y(self, value: float) -> None:
self.coords = (self.coords[0], value, self.coords[2])

@property
def z(self) -> float:
return self.coords[2]

x = property(_getx, _setx)
y = property(_gety, _sety)
z = property(_getz, _setz)
@z.setter
def z(self, value: float) -> None:
self.coords = (self.coords[0], self.coords[1], value)

def _getsymbol(self):
@property
def symbol(self) -> str:
if self.atnum == 0:
return self._dummysymbol
else:
return PT.get_symbol(self.atnum)

def _setsymbol(self, symbol):
@symbol.setter
def symbol(self, symbol: str) -> None:
if symbol.lower().capitalize() in PT.dummysymbols:
self.atnum = 0
self._dummysymbol = symbol.lower().capitalize()
else:
self.atnum = PT.get_atomic_number(symbol)
self._dummysymbol = None

symbol = property(_getsymbol, _setsymbol)

def _getmass(self):
@property
def mass(self) -> float:
return PT.get_mass(self.atnum)

mass = property(_getmass)

def _getradius(self):
@property
def radius(self) -> float:
return PT.get_radius(self.atnum)

radius = property(_getradius)

def _getconnectors(self):
@property
def connectors(self) -> int:
return PT.get_connectors(self.atnum)

connectors = property(_getconnectors)

def _ismetallic(self):
@property
def is_metallic(self) -> bool:
return PT.get_metallic(self.atnum)

is_metallic = property(_ismetallic)

def _iselectronegative(self):
@property
def is_electronegative(self) -> bool:
return PT.get_electronegative(self.atnum)

is_electronegative = property(_iselectronegative)

def translate(self, vector, unit="angstrom"):
def translate(self, vector: Iterable[float], unit: str = "angstrom") -> None:
"""Move this atom in space by *vector*, expressed in *unit*.
*vector* should be an iterable container of length 3 (usually tuple, list or numpy array). *unit* describes unit of values stored in *vector*.
Expand All @@ -201,7 +202,7 @@ def translate(self, vector, unit="angstrom"):
ratio = Units.conversion_ratio(unit, "angstrom")
self.coords = tuple(i + j * ratio for i, j in zip(self, vector))

def move_to(self, point, unit="angstrom"):
def move_to(self, point: Iterable[float], unit: str = "angstrom") -> None:
"""Move this atom to a given *point* in space, expressed in *unit*.
*point* should be an iterable container of length 3 (for example: tuple, |Atom|, list, numpy array). *unit* describes unit of values stored in *point*.
Expand All @@ -211,7 +212,7 @@ def move_to(self, point, unit="angstrom"):
ratio = Units.conversion_ratio(unit, "angstrom")
self.coords = tuple(i * ratio for i in point)

def distance_to(self, point, unit="angstrom", result_unit="angstrom") -> float:
def distance_to(self, point: Iterable[float], unit: str = "angstrom", result_unit: str = "angstrom") -> float:
"""Measure the distance between this atom and *point*.
*point* should be an iterable container of length 3 (for example: tuple, |Atom|, list, numpy array). *unit* describes unit of values stored in *point*. Returned value is expressed in *result_unit*.
Expand All @@ -224,7 +225,9 @@ def distance_to(self, point, unit="angstrom", result_unit="angstrom") -> float:
res += (i - j * ratio) ** 2
return Units.convert(math.sqrt(res), "angstrom", result_unit)

def vector_to(self, point, unit="angstrom", result_unit="angstrom"):
def vector_to(
self, point: Iterable[float], unit: str = "angstrom", result_unit: str = "angstrom"
) -> Tuple[float, float, float]:
"""Calculate a vector from this atom to *point*.
*point* should be an iterable container of length 3 (for example: tuple, |Atom|, list, numpy array). *unit* describes unit of values stored in *point*. Returned value is expressed in *result_unit*.
Expand All @@ -235,7 +238,14 @@ def vector_to(self, point, unit="angstrom", result_unit="angstrom"):
resultratio = Units.conversion_ratio("angstrom", result_unit)
return tuple((i * ratio - j) * resultratio for i, j in zip(point, self))

def angle(self, point1, point2, point1unit="angstrom", point2unit="angstrom", result_unit="radian"):
def angle(
self,
point1: Iterable[float],
point2: Iterable[float],
point1unit: str = "angstrom",
point2unit: str = "angstrom",
result_unit: str = "radian",
) -> float:
"""Calculate an angle between vectors pointing from this atom to *point1* and *point2*.
*point1* and *point2* should be iterable containers of length 3 (for example: tuple, |Atom|, list, numpy array). Values stored in them are expressed in, respectively, *point1unit* and *point2unit*. Returned value is expressed in *result_unit*.
Expand All @@ -246,7 +256,7 @@ def angle(self, point1, point2, point1unit="angstrom", point2unit="angstrom", re
den = self.distance_to(point1, point1unit) * self.distance_to(point2, point2unit)
return Units.convert(math.acos(num / den), "radian", result_unit)

def rotate(self, matrix):
def rotate(self, matrix: Union[Iterable[float], Iterable[Iterable[float]]]) -> None:
"""Rotate this atom according to a rotation *matrix*.
*matrix* should be a container with 9 numerical values. It can be a list (tuple, numpy array etc.) listing matrix elements row-wise, either flat (``[1,2,3,4,5,6,7,8,9]``) or in two-level fashion (``[[1,2,3],[4,5,6],[7,8,9]]``).
Expand All @@ -258,6 +268,6 @@ def rotate(self, matrix):
matrix = np.array(matrix).reshape(3, 3)
self.coords = tuple(np.dot(matrix, np.array(self.coords)))

def neighbors(self):
def neighbors(self) -> List["Bond"]:
"""Return a list of neighbors of this atom within the molecule. The list follows the same order as the ``bonds`` attribute."""
return [b.other_end(self) for b in self.bonds]
18 changes: 11 additions & 7 deletions mol/bond.py
Original file line number Diff line number Diff line change
@@ -1,10 +1,14 @@
import numpy as np
from typing import Optional, Tuple, TYPE_CHECKING

from scm.plams.core.errors import MoleculeError
from scm.plams.core.settings import Settings

__all__ = ["Bond"]

if TYPE_CHECKING:
from scm.plams.mol.atom import Atom


class Bond:
"""A class representing a bond between two atoms.
Expand Down Expand Up @@ -40,15 +44,15 @@ def __iter__(self):
yield self.atom1
yield self.atom2

def is_aromatic(self):
def is_aromatic(self) -> bool:
"""Check if this bond is aromatic."""
return self.order == Bond.AR

def length(self, unit="angstrom"):
def length(self, unit: str = "angstrom") -> float:
"""Return bond length, expressed in *unit*."""
return self.atom1.distance_to(self.atom2, result_unit=unit)

def as_vector(self, start=None, unit="angstrom"):
def as_vector(self, start: Optional["Atom"] = None, unit: str = "angstrom") -> Tuple[float, float, float]:
"""Return a vector between two atoms that form this bond. *start* can be used to indicate which atom should be the beginning of that vector. If not specified, ``self.atom1`` is used. Returned value if a tuple of length 3, expressed in *unit*."""
if start:
if start not in self:
Expand All @@ -58,7 +62,7 @@ def as_vector(self, start=None, unit="angstrom"):
a, b = self.atom1, self.atom2
return a.vector_to(b, result_unit=unit)

def other_end(self, atom):
def other_end(self, atom: "Atom") -> "Atom":
"""Return the atom on the other end of this bond with respect to *atom*. *atom* has to be one of the atoms forming this bond, otherwise an exception is raised."""
if atom is self.atom1:
return self.atom2
Expand All @@ -67,7 +71,7 @@ def other_end(self, atom):
else:
raise MoleculeError("Bond.other_end: invalid atom passed")

def resize(self, moving_atom, length, unit="angstrom"):
def resize(self, moving_atom: "Atom", length: float, unit: str = "angstrom") -> None:
"""Change the length of this bond to *length* expressed in *unit* by moving *moving_atom*.
*moving_atom* should be one of the atoms that form this bond. This atom is moved along the bond axis in such a way that new bond length equals *length*. If this bond is a part of a |Molecule| the whole part connected to *moving_atom* is moved.
Expand All @@ -85,7 +89,7 @@ def resize(self, moving_atom, length, unit="angstrom"):
trans_v = (1 - length / self.length(unit)) * bond_v
moving_atom.translate(trans_v)

def rotate(self, moving_atom, angle, unit="radian"):
def rotate(self, moving_atom: "Atom", angle: float, unit: str = "radian"):
"""Rotate part of the molecule containing *moving_atom* along axis defined by this bond by an *angle* expressed in *unit*.
Calling this method makes sense only if this bond is a part of a |Molecule|. *moving_atom* should be one of the atoms that form this bond and it indicates which part of the molecule is rotated. A positive value of *angle* denotes counterclockwise rotation (when looking along the bond, from the stationary part of the molecule).
Expand All @@ -98,7 +102,7 @@ def rotate(self, moving_atom, angle, unit="radian"):
if self.mol:
self.mol.rotate_bond(self, moving_atom, angle, unit)

def has_cell_shifts(self):
def has_cell_shifts(self) -> bool:
return (
"suffix" in self.properties
and isinstance(self.properties.suffix, str)
Expand Down
10 changes: 9 additions & 1 deletion setup.py
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,15 @@
packages=packages,
package_dir={"scm.plams": "."},
package_data={
"scm.plams": [".flake8", "plams_defaults", "examples/*", "unit_tests/*", "unit_tests/xyz/*", "unit_tests/pdb/*"]
"scm.plams": [
".flake8",
"plams_defaults",
"examples/*",
"unit_tests/*",
"unit_tests/xyz/*",
"unit_tests/pdb/*",
"unit_tests/rkf/*",
]
},
scripts=[opj("scripts", "plams")],
)
17 changes: 11 additions & 6 deletions tools/geometry.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
import numpy as np


try:
from scipy.spatial.distance import cdist

Expand All @@ -19,6 +20,8 @@
"cellvectors_from_shape",
]

HALF_PI = np.pi / 2


def rotation_matrix(vec1, vec2):
"""
Expand Down Expand Up @@ -106,7 +109,7 @@ def cell_shape(lattice):
Converts lattice vectors to lengths and angles (in radians)
Sets internal cell size data, based on set of cell vectors.
*cellvectors* is list containing three cell vectors (a 3x3 matrix)
*lattice* is list containing three cell vectors (a 3x3 matrix)
"""
lattice = np.asarray(lattice)
a, b, c = np.sqrt((lattice**2).sum(axis=1))
Expand Down Expand Up @@ -165,16 +168,18 @@ def cell_angles(lattice, unit="degree"):

def cellvectors_from_shape(box):
"""
Converts lengths and angles (in radians) of lattice vectors to the lattice vectors
Converts lengths and angles (in radians) of lattice vectors to the lattice vectors.
*box* should be an iterable of length 3, containing the lengths a, b, c, or an iterable of length 6 additionally containing the angles alpha, beta, gamma
"""
a = box[0]
b = box[1]
c = box[2]
alpha, beta, gamma = 90.0, 90.0, 90
alpha, beta, gamma = HALF_PI, HALF_PI, HALF_PI
if len(box) == 6:
alpha = box[3] # *np.pi/180.
beta = box[4] # *np.pi/180.
gamma = box[5] # *np.pi/180.
alpha = box[3]
beta = box[4]
gamma = box[5]

va = [a, 0.0, 0.0]
vb = [b * np.cos(gamma), b * np.sin(gamma), 0.0]
Expand Down
Loading

0 comments on commit a4d9060

Please sign in to comment.