omic/propka
1
0
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases Wiki Activity

De-lint vector_algebra.py

Browse Source
This commit is contained in:
Nathan Baker
2020-05-25 13:53:04 -07:00
parent 2321a4df8a
commit 6702f96730
3 changed files with 287 additions and 193 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
propka/ligand.py
View File

@@ -1,6 +1,6 @@
"""Ligand classes and functions."""
from propka.calculations import squared_distance
from propka.vector_algebra import vector
from propka.vector_algebra import Vector
ALL_SYBYL_TYPES = [
@@ -303,12 +303,12 @@ def are_atoms_planar(atoms):
return False
if len(atoms) < 4:
return False
vec1 = vector(atom1=atoms[0], atom2=atoms[1])
vec2 = vector(atom1=atoms[0], atom2=atoms[2])
vec1 = Vector(atom1=atoms[0], atom2=atoms[1])
vec2 = Vector(atom1=atoms[0], atom2=atoms[2])
norm = (vec1**vec2).rescale(1.0)
margin = PLANARITY_MARGIN
for atom in atoms[3:]:
vec = vector(atom1=atoms[0], atom2=atom).rescale(1.0)
vec = Vector(atom1=atoms[0], atom2=atom).rescale(1.0)
if abs(vec*norm) > margin:
return False
return True

30
propka/protonate.py
View File

@@ -2,7 +2,7 @@
import math
import propka.bonds
import propka.atom
from propka.vector_algebra import rotate_vector_around_an_axis, vector
from propka.vector_algebra import rotate_vector_around_an_axis, Vector
from propka.lib import warning, debug
@@ -202,14 +202,14 @@ class Protonate:
"""
debug('TRIGONAL - %d bonded atoms' % len(atom.bonded_atoms))
rot_angle = math.radians(120.0)
cvec = vector(atom1=atom)
cvec = Vector(atom1=atom)
# 0 bonds
if len(atom.bonded_atoms) == 0:
pass
# 1 bond
if len(atom.bonded_atoms) == 1 and atom.number_of_protons_to_add > 0:
# Add another atom with the right angle to the first one
avec = vector(atom1=atom, atom2=atom.bonded_atoms[0])
avec = Vector(atom1=atom, atom2=atom.bonded_atoms[0])
# use plane of bonded trigonal atom - e.g. arg
self.set_steric_number_and_lone_pairs(atom.bonded_atoms[0])
if (atom.bonded_atoms[0].steric_number == 3
@@ -220,15 +220,15 @@ class Protonate:
for i, bonded_atom in enumerate(atom.bonded_atoms[0].bonded_atoms):
if bonded_atom != atom:
other_atom_indices.append(i)
vec1 = vector(atom1=atom, atom2=atom.bonded_atoms[0])
vec2 = vector(atom1=atom.bonded_atoms[0],
vec1 = Vector(atom1=atom, atom2=atom.bonded_atoms[0])
vec2 = Vector(atom1=atom.bonded_atoms[0],
atom2=atom.bonded_atoms[0]
.bonded_atoms[other_atom_indices[0]])
axis = vec1**vec2
# this is a trick to make sure that the order of atoms doesn't
# influence the final postions of added protons
if len(other_atom_indices) > 1:
vec3 = vector(atom1=atom.bonded_atoms[0],
vec3 = Vector(atom1=atom.bonded_atoms[0],
atom2=atom.bonded_atoms[0]
.bonded_atoms[other_atom_indices[1]])
axis2 = vec1**vec3
@@ -244,8 +244,8 @@ class Protonate:
# 2 bonds
if len(atom.bonded_atoms) == 2 and atom.number_of_protons_to_add > 0:
# Add another atom with the right angle to the first two
avec1 = vector(atom1=atom, atom2=atom.bonded_atoms[0]).rescale(1.0)
avec2 = vector(atom1=atom, atom2=atom.bonded_atoms[1]).rescale(1.0)
avec1 = Vector(atom1=atom, atom2=atom.bonded_atoms[0]).rescale(1.0)
avec2 = Vector(atom1=atom, atom2=atom.bonded_atoms[1]).rescale(1.0)
new_a = -avec1 - avec2
new_a = self.set_bond_distance(new_a, atom.element)
@@ -260,14 +260,14 @@ class Protonate:
debug('TETRAHEDRAL - %d bonded atoms' % len(atom.bonded_atoms))
# TODO - might be good to move tetrahedral angle to constant
rot_angle = math.radians(109.5)
cvec = vector(atom1=atom)
cvec = Vector(atom1=atom)
# 0 bonds
if len(atom.bonded_atoms) == 0:
pass
# 1 bond
if len(atom.bonded_atoms) == 1 and atom.number_of_protons_to_add > 0:
# Add another atom with the right angle to the first one
avec = vector(atom1=atom, atom2=atom.bonded_atoms[0])
avec = Vector(atom1=atom, atom2=atom.bonded_atoms[0])
axis = avec.orthogonal()
avec = rotate_vector_around_an_axis(rot_angle, axis, avec)
avec = self.set_bond_distance(avec, atom.element)
@@ -275,8 +275,8 @@ class Protonate:
# 2 bonds
if len(atom.bonded_atoms) == 2 and atom.number_of_protons_to_add > 0:
# Add another atom with the right angle to the first two
avec1 = vector(atom1=atom, atom2=atom.bonded_atoms[0]).rescale(1.0)
avec2 = vector(atom1=atom, atom2=atom.bonded_atoms[1]).rescale(1.0)
avec1 = Vector(atom1=atom, atom2=atom.bonded_atoms[0]).rescale(1.0)
avec2 = Vector(atom1=atom, atom2=atom.bonded_atoms[1]).rescale(1.0)
axis = avec1 + avec2
new_a = rotate_vector_around_an_axis(math.radians(90), axis,
-avec1)
@@ -284,9 +284,9 @@ class Protonate:
self.add_proton(atom, cvec+new_a)
# 3 bonds
if len(atom.bonded_atoms) == 3 and atom.number_of_protons_to_add > 0:
avec1 = vector(atom1=atom, atom2=atom.bonded_atoms[0]).rescale(1.0)
avec2 = vector(atom1=atom, atom2=atom.bonded_atoms[1]).rescale(1.0)
avec3 = vector(atom1=atom, atom2=atom.bonded_atoms[2]).rescale(1.0)
avec1 = Vector(atom1=atom, atom2=atom.bonded_atoms[0]).rescale(1.0)
avec2 = Vector(atom1=atom, atom2=atom.bonded_atoms[1]).rescale(1.0)
avec3 = Vector(atom1=atom, atom2=atom.bonded_atoms[2]).rescale(1.0)
new_a = -avec1-avec2-avec3
new_a = self.set_bond_distance(new_a, atom.element)
self.add_proton(atom, cvec+new_a)

410
propka/vector_algebra.py
View File

@@ -1,11 +1,21 @@
from __future__ import division
from __future__ import print_function
"""Vector algebra for PROPKA."""
import math
from propka.lib import info
from propka.lib import info, get_sorted_configurations
class vector:
""" Vector """
def __init__(self, xi=0.0, yi=0.0, zi=0.0, atom1 = 0, atom2 = 0):
class Vector:
"""Vector"""
def __init__(self, xi=0.0, yi=0.0, zi=0.0, atom1=None, atom2=None):
"""Initialize vector.
Args:
xi: default x-coordinate
yi: default y-coordinate
zi: default z-coordinate
atom1: atom center used to define default coordinate
atom2: two atom centers used to define vector
"""
self.x = xi
self.y = yi
self.z = zi
@@ -22,56 +32,53 @@ class vector:
self.y = atom2.y - self.y
self.z = atom2.z - self.z
return
def __add__(self, other):
return vector(self.x + other.x,
return Vector(self.x + other.x,
self.y + other.y,
self.z + other.z)
def __sub__(self, other):
return vector(self.x - other.x,
return Vector(self.x - other.x,
self.y - other.y,
self.z - other.z)
def __mul__(self, other):
""" Dot product, scalar and matrix multiplication """
if isinstance(other,vector):
"""Dot product, scalar and matrix multiplication."""
if isinstance(other, Vector):
return self.x * other.x + self.y * other.y + self.z * other.z
elif isinstance(other, matrix4x4):
return vector(
xi = other.a11*self.x + other.a12*self.y + other.a13*self.z + other.a14*1.0,
yi = other.a21*self.x + other.a22*self.y + other.a23*self.z + other.a24*1.0,
zi = other.a31*self.x + other.a32*self.y + other.a33*self.z + other.a34*1.0
elif isinstance(other, Matrix4x4):
return Vector(
xi=other.a11*self.x + other.a12*self.y + other.a13*self.z + other.a14*1.0,
yi=other.a21*self.x + other.a22*self.y + other.a23*self.z + other.a24*1.0,
zi=other.a31*self.x + other.a32*self.y + other.a33*self.z + other.a34*1.0
)
elif type(other) in [int, float]:
return vector(self.x * other, self.y * other, self.z * other)
return Vector(self.x * other, self.y * other, self.z * other)
else:
info('%s not supported' % type(other))
raise TypeError
def __rmul__(self,other):
def __rmul__(self, other):
return self.__mul__(other)
def __pow__(self, other):
""" Cross product """
return vector(self.y * other.z - self.z * other.y,
"""Cross product."""
return Vector(self.y * other.z - self.z * other.y,
self.z * other.x - self.x * other.z,
self.x * other.y - self.y * other.x)
def __neg__(self):
res = vector(xi = -self.x,
yi = -self.y,
zi = -self.z)
res = Vector(xi=-self.x,
yi=-self.y,
zi=-self.z)
return res
def sq_length(self):
"""Return vector squared-length"""
return self.x * self.x + self.y * self.y + self.z * self.z
def length(self):
"""Return vector length."""
return math.sqrt(self.sq_length())
def __str__(self):
@@ -82,164 +89,220 @@ class vector:
def orthogonal(self):
""" Returns a vector orthogonal to self """
res = vector(self.y, -self.x, 0)
res = Vector(self.y, -self.x, 0)
if abs(self.y) < abs(self.z):
res = vector(self.z, 0, -self.x)
res = Vector(self.z, 0, -self.x)
return res
def rescale(self, new_length):
""" Rescale vector to new length while preserving direction """
frac = new_length/(self.length())
res = vector(xi = self.x*frac,
yi = self.y*frac,
zi = self.z*frac)
res = Vector(xi=self.x*frac,
yi=self.y*frac,
zi=self.z*frac)
return res
class matrix4x4:
class Matrix4x4:
"""A 4-by-4 matrix class."""
def __init__(self,
a11i=0.0, a12i=0.0, a13i=0.0, a14i=0.0,
a21i=0.0, a22i=0.0, a23i=0.0, a24i=0.0,
a31i=0.0, a32i=0.0, a33i=0.0, a34i=0.0,
a41i=0.0, a42i=0.0, a43i=0.0, a44i=0.0):
"""Initialize with matrix elements."""
# Row 1
self.a11 = a11i
self.a12 = a12i
self.a13 = a13i
self.a14 = a14i
# Row 2
self.a21 = a21i
self.a22 = a22i
self.a23 = a23i
self.a24 = a24i
# Row 3
self.a31 = a31i
self.a32 = a32i
self.a33 = a33i
self.a34 = a34i
# Row 4
self.a41 = a41i
self.a42 = a42i
self.a43 = a43i
self.a44 = a44i
return
def angle(avec, bvec):
"""Get the angle between two vectors.
Args:
avec: vector 1
bvec: vector 2
Returns:
angle in radians
"""
dot = avec * bvec
return math.acos(dot / (avec.length() * bvec.length()))
def angle_degrees(avec, bvec):
"""Get the angle between two vectors in degrees.
Args:
avec: vector 1
bvec: vector 2
Returns:
angle in degrees
"""
return math.degrees(angle(avec, bvec))
def signed_angle_around_axis(avec, bvec, axis):
"""Get signed angle of two vectors around axis in radians.
# methods working on vectors
Args:
avec: vector 1
bvec: vector 2
axis: axis
Returns:
angle in radians
"""
norma = avec**axis
normb = bvec**axis
ang = angle(norma, normb)
dot_ = bvec*(avec**axis)
if dot_ < 0:
ang = -ang
return ang
def angle(a, b):
dot = a * b
return math.acos(dot / (a.length() * b.length()))
def rotate_vector_around_an_axis(theta, axis, vec):
"""Rotate vector around an axis.
def angle_degrees(a,b):
return math.degrees(angle(a, b))
def signed_angle_around_axis(a,b, axis):
na = a**axis
nb = b**axis
v = angle(na,nb)
d = b*(a**axis)
if d < 0:
v =-v
return v
def signed_angle_degrees(a,b):
return 180/math.pi * signed_angle(a, b)
def rotate_vector_around_an_axis(theta, axis, v):
#print "# 1. rotate space about the z-axis so that the rotation axis lies in the xz-plane"
Args:
theta: rotation angle (in radians)
axis: axis for rotation
vec: vector to rotate
Returns:
rotated vector
"""
gamma = 0.0
if axis.y != 0:
if axis.x != 0:
gamma = -axis.x/abs(axis.x)*math.asin(axis.y/(math.sqrt(axis.x*axis.x + axis.y*axis.y)))
gamma = -axis.x/abs(axis.x)*math.asin(
axis.y/(math.sqrt(axis.x*axis.x + axis.y*axis.y)))
else:
gamma = math.pi/2.0
Rz = rotate_atoms_around_z_axis(gamma)
v = Rz * v
axis = Rz * axis
#print "# 2. rotate space about the y-axis so that the rotation axis lies along the z-axis"
rot_z = rotate_atoms_around_z_axis(gamma)
vec = rot_z * vec
axis = rot_z * axis
beta = 0.0
if axis.x != 0:
beta = -axis.x/abs(axis.x)*math.acos(axis.z/math.sqrt(axis.x*axis.x + axis.z*axis.z))
Ry = rotate_atoms_around_y_axis(beta)
v = Ry * v
axis = Ry *axis
beta = -axis.x/abs(axis.x)*math.acos(
axis.z/math.sqrt(axis.x*axis.x + axis.z*axis.z))
rot_y = rotate_atoms_around_y_axis(beta)
vec = rot_y * vec
axis = rot_y * axis
rot_z = rotate_atoms_around_z_axis(theta)
vec = rot_z * vec
rot_y = rotate_atoms_around_y_axis(-beta)
vec = rot_y * vec
rot_z = rotate_atoms_around_z_axis(-gamma)
vec = rot_z * vec
return vec
#print "# 3. perform the desired rotation by theta about the z-axis"
Rz = rotate_atoms_around_z_axis(theta)
v = Rz * v
#print "# 4. apply the inverse of step 2."
Ry = rotate_atoms_around_y_axis(-beta)
v = Ry * v
def rotate_atoms_around_z_axis(theta):
"""Get rotation matrix for z-axis.
#print "# 5. apply the inverse of step 1."
Rz = rotate_atoms_around_z_axis(-gamma)
v = Rz * v
return v
def rotate_atoms_around_z_axis(angle):
Rz = matrix4x4(
a11i = math.cos(angle), a12i = -math.sin(angle), a13i = 0.0, a14i = 0.0,
a21i = math.sin(angle), a22i = math.cos(angle), a23i = 0.0, a24i = 0.0,
a31i = 0.0 , a32i = 0.0 , a33i = 1.0, a34i = 0.0,
a41i = 0.0 , a42i = 0.0 , a43i = 0.0, a44i = 1.0
Args:
theta: angle of rotation (radians)
Returns:
rotation matrix
"""
return Matrix4x4(
a11i=math.cos(theta),
a12i=-math.sin(theta),
a13i=0.0,
a14i=0.0,
a21i=math.sin(theta),
a22i=math.cos(theta),
a23i=0.0,
a24i=0.0,
a31i=0.0,
a32i=0.0,
a33i=1.0,
a34i=0.0,
a41i=0.0,
a42i=0.0,
a43i=0.0,
a44i=1.0
)
return Rz
def rotate_atoms_around_y_axis(theta):
"""Get rotation matrix for y-axis.
def rotate_atoms_around_y_axis(angle):
Ry = matrix4x4(
a11i = math.cos(angle), a12i = 0.0, a13i = math.sin(angle), a14i = 0.0,
a21i = 0.0 , a22i = 1.0, a23i = 0.0 , a24i = 0.0,
a31i = -math.sin(angle), a32i = 0.0, a33i = math.cos(angle), a34i = 0.0,
a41i = 0.0 , a42i = 0.0, a43i = 0.0 , a44i = 1.0
Args:
theta: angle of rotation (radians)
Returns:
rotation matrix
"""
return Matrix4x4(
a11i=math.cos(theta),
a12i=0.0,
a13i=math.sin(theta),
a14i=0.0,
a21i=0.0,
a22i=1.0,
a23i=0.0,
a24i=0.0,
a31i=-math.sin(theta),
a32i=0.0,
a33i=math.cos(theta),
a34i=0.0,
a41i=0.0,
a42i=0.0,
a43i=0.0,
a44i=1.0
)
return Ry
class MultiVector:
"""Collection of vectors for multiple configurations of atoms.
TODO - this class does not appear to be used or covered by tests
"""
class multi_vector:
def __init__(self, atom1=0, atom2=0):
def __init__(self, atom1=None, atom2=None):
"""Initialize with atom configurations.
Args:
atom1: first atom to define vector
atom2: second atom to define vector
"""
self.vectors = []
self.keys = []
self.result = None
# store vectors for all configurations of atoms
if atom1!=0:
self.keys = lib.get_sorted_configurations(atom1.configurations.keys())
if atom2!=0:
keys2 = lib.get_sorted_configurations(atom2.configurations.keys())
if atom1 is not None:
self.keys = get_sorted_configurations(atom1.configurations.keys())
if atom2 is not None:
keys2 = get_sorted_configurations(atom2.configurations.keys())
if self.keys != keys2:
raise 'Cannot make multi vector: Atomic configurations mismatch for\n %s\n %s\n'%(atom1,atom2)
str_ = ('Cannot make multi vector: Atomic configurations '
'mismatch for\n %s\n %s\n' % (atom1, atom2))
raise KeyError(str_)
for key in self.keys:
atom1.setConfiguration(key)
if atom2!=0:
if atom2 != 0:
atom2.setConfiguration(key)
v = vector(atom1=atom1, atom2=atom2)
self.vectors.append(v)
#info(key,v)
return
vec = Vector(atom1=atom1, atom2=atom2)
self.vectors.append(vec)
def __getattribute__(self,name):
def __getattribute__(self, name):
try:
return object.__getattribute__(self, name)
except AttributeError:
@@ -247,72 +310,103 @@ class multi_vector:
def __str__(self):
res = ''
for i in range(len(self.keys)):
res += '%s %s\n'%(self.keys[i], self.vectors[i])
for i, key in enumerate(self.keys):
res += '%s %s\n' % (key, self.vectors[i])
return res
def do_job(self, job):
#info(job)
self.res = multi_vector()
for i in range(len(self.vectors)):
self.res.vectors.append(eval('self.vectors[%d].%s()'%(i,job)))
self.res.keys.append(self.keys[i])
"""Append vectors to configuration.
Args:
job: name of function to apply to vectors
Returns:
TODO - figure out what this is
"""
self.result = MultiVector()
for i, vector in enumerate(self.vectors):
func = getattr(vector, job)
self.result.vectors.append(func())
self.result.keys.append(self.keys[i])
return self.get_result
@property
def get_result(self):
return self.res
"""Return the latest result."""
return self.result
def generic_operation(self, operation, other):
if self.keys != other.keys:
raise 'Incompatable keys'
"""Perform a generic operation between two MultiVector objects.
self.res = multi_vector()
Args:
operation: operation to perform (string)
other: other MultiVector object
"""
if self.keys != other.keys:
raise 'Incompatible keys'
self.result = MultiVector()
for i in range(len(self.vectors)):
self.res.vectors.append(eval('self.vectors[%d] %s other.vectors[%d]'%(i,operation,i)))
self.res.keys.append(self.keys[i])
return
self.result.vectors.append(
# TODO - eliminate eval() or entire class
eval('self.vectors[%d] %s other.vectors[%d]'
% (i, operation, i)))
self.result.keys.append(self.keys[i])
def __add__(self, other):
self.generic_operation('+',other)
return self.res
self.generic_operation('+', other)
return self.result
def __sub__(self, other):
self.generic_operation('-',other)
return self.res
self.generic_operation('-', other)
return self.result
def __mul__(self, other):
self.generic_operation('*',other)
return self.res
self.generic_operation('*', other)
return self.result
def __pow__(self, other):
self.generic_operation('**',other)
return self.res
self.generic_operation('**', other)
return self.result
def generic_self_operation(self, operation):
@staticmethod
def generic_self_operation(_):
"""TODO - delete this."""
return
def __neg__(self):
self.generic_operation('*',-1.0)
return self.res
self.generic_operation('*', -1.0)
return self.result
def rescale(self, new_length):
self.res = multi_vector()
for i in range(len(self.vectors)):
self.res.vectors.append(self.vectors[i].rescale(new_length))
self.res.keys.append(self.keys[i])
"""Rescale multi-vector to new length.
Args:
new_length: new length for multi-vector
Result:
MultiVector object
"""
self.result = MultiVector()
for i, vector in enumerate(self.vectors):
self.result.vectors.append(vector.rescale(new_length))
self.result.keys.append(self.keys[i])
return self.res
def rotate_multi_vector_around_an_axis(theta, axis, v):
""" both axis ans v must be multi_vectors """
def rotate_multi_vector_around_an_axis(theta, axis, vec):
"""Rotate a multi-vector around an axis.
if axis.keys != v.keys:
raise 'Incompatible keys in rotate multi_vector'
res = multi_vector()
for i in range(len(v.keys)):
res.vectors.append(rotate_vector_around_an_axis(theta, axis.vectors[i], v.vectors[i]))
res.keys.append(v.keys[i])
NOTE - both axis ans v must be MultiVectors.
Args:
theta: angle (in radians)
axis: multi-vector axis
vec: multi-vector vector
"""
if axis.keys != vec.keys:
raise 'Incompatible keys in rotate MultiVector'
res = MultiVector()
for i, key in enumerate(vec.keys):
res.vectors.append(rotate_vector_around_an_axis(theta,
axis.vectors[i],
vec.vectors[i]))
res.keys.append(key)
return res