From 6702f967300bb54e2183627f8341d5a5755c0cb5 Mon Sep 17 00:00:00 2001 From: Nathan Baker Date: Mon, 25 May 2020 13:53:04 -0700 Subject: [PATCH] De-lint vector_algebra.py --- propka/ligand.py | 8 +- propka/protonate.py | 30 +-- propka/vector_algebra.py | 442 ++++++++++++++++++++++++--------------- 3 files changed, 287 insertions(+), 193 deletions(-) diff --git a/propka/ligand.py b/propka/ligand.py index e97efc6..2be158d 100644 --- a/propka/ligand.py +++ b/propka/ligand.py @@ -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 diff --git a/propka/protonate.py b/propka/protonate.py index d9c0c09..e3d96ff 100644 --- a/propka/protonate.py +++ b/propka/protonate.py @@ -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) diff --git a/propka/vector_algebra.py b/propka/vector_algebra.py index 202d9c5..8d874b5 100644 --- a/propka/vector_algebra.py +++ b/propka/vector_algebra.py @@ -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) + 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): + def __mul__(self, other): + """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): - return self.__mul__(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" - 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))) - 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" - beta = 0.0 + 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: - 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 + 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 + 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)) + 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 - - #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(theta): + """Get rotation matrix for z-axis. -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(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 +def rotate_atoms_around_y_axis(theta): + """Get rotation matrix for y-axis. + + 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. -class multi_vector: - def __init__(self, atom1=0, atom2=0): + TODO - this class does not appear to be used or covered by tests + """ + + 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 - - def __getattribute__(self,name): + vec = Vector(atom1=atom1, atom2=atom2) + self.vectors.append(vec) + + 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 - - def generic_operation(self, operation, other): - if self.keys != other.keys: - raise 'Incompatable keys' + """Return the latest result.""" + return self.result - self.res = multi_vector() + def generic_operation(self, operation, other): + """Perform a generic operation between two MultiVector objects. + + 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 """ - - 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]) - +def rotate_multi_vector_around_an_axis(theta, axis, vec): + """Rotate a multi-vector around an axis. + + 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