get turns, bridges, ladders and helix oop

src/atom.py

@@ -1,4 +1,5 @@
 import math
+from structure import *
 class Atom:
 
     def dist_atoms(self, atom2):
@@ -34,21 +35,225 @@ class Residue:
     def h_bond(self, res2):
         if("H" not in res2.atoms.keys()):
             return(False)
-        # elementary charge in Coulomb
+        # dimensionnal factor, in kcal/mole
-        #e = 1.602176634 * 10*math.exp(-19)
-        # dimensionnal factor f
         f = 332
         # partial charges
         q1 = 0.42
         q2 = 0.20
-        # r, distance between O-N atoms, in angströms
+        # distance between O-N atoms, in angströms
         r_ON = self.atoms["O"].dist_atoms(res2.atoms["N"])
-        # r, distance between C-H atoms, in angströms
+        # distance between C-H atoms, in angströms
         r_CH = self.atoms["C"].dist_atoms(res2.atoms["H"])
-        # r, distance between O-H atoms, in angströms
+        # distance between O-H atoms, in angströms
         r_OH = self.atoms["O"].dist_atoms(res2.atoms["H"])
-        # r, distance between C-N atoms, in angströms
+        # distance between C-N atoms, in angströms
         r_CN = self.atoms["C"].dist_atoms(res2.atoms["N"])
-        # Electrostatic interaction energy, in kcal/mole
+        # electrostatic interaction energy, in kcal/mole
         E = q1*q2*(1/r_ON + 1/r_CH - 1/r_OH - 1/r_CN)*f
         return(E)
+
+    def get_turns(self, residues, turns):
+        """
+        Get all the turns from a specific residue.
+        """
+        i = residues.index(self)
+        for j in range(3,6):
+                if(i+j<len(residues)):
+                    if(self.h_bond(residues[i+j])<-0.5):
+                        print(j,"TURN", residues[i].resid, residues[i+j].resid)
+                        turns[i] = Turn(j,residues[i].resid)
+
+    def get_bridges(self, residues):
+        bridges = {}
+        bridge = {}
+        strongest_bridge = {}
+        i = residues.index(self)
+        if(i >= 1 and i < len(residues)-4):
+            E_min = 0
+            for j in range(i+2,len(residues)-1):
+                # select triplet with the minimal energy 
+
+                if(residues[i-1].h_bond(residues[j])<-0.5
+                   and residues[j].h_bond(residues[i+1])<-0.5):
+                    bridge = {'res1':residues[i-1].h_bond(residues[j]),
+                              'res2':residues[j].h_bond(residues[i+1]),
+                              'ipos':residues[i].resid,
+                              'jpos':residues[j].resid,
+                              'btype':"para"}
+
+                if(residues[j-1].h_bond(residues[i])<-0.5
+                   and residues[i].h_bond(residues[j+1])<-0.5):
+                    bridge = {'res1':residues[j-1].h_bond(residues[i]),
+                              'res2':residues[i].h_bond(residues[j+1]),
+                              'ipos':residues[i].resid,
+                              'jpos':residues[j].resid,
+                              'btype':"para"}
+
+                if(residues[i].h_bond(residues[j])<-0.5
+                   and residues[j].h_bond(residues[i])<-0.5):
+                    bridge = {'res1':residues[i].h_bond(residues[j]),
+                              'res2':residues[j].h_bond(residues[i]),
+                              'ipos':residues[i].resid,
+                              'jpos':residues[j].resid,
+                              'btype':"anti"}
+
+                if(residues[i-1].h_bond(residues[j+1])<-0.5
+                   and residues[j-1].h_bond(residues[i+1])<-0.5):
+                    bridge = {'res1':residues[i-1].h_bond(residues[j+1]),
+                              'res2':residues[j-1].h_bond(residues[i+1]),
+                              'ipos':residues[i].resid,
+                              'jpos':residues[j].resid,
+                              'btype':"anti"}
+
+                if(bridge):
+                    if(bridge['res1']+bridge['res2']<E_min):
+                        E_min = bridge['res1']+bridge['res2']
+                        strongest_bridge = bridge
+                        bridge = {}
+                        coord_bridge = [i,j]
+            # finally add the strongest bridge at i and j pos
+            if(strongest_bridge):
+                bridges[coord_bridge[0]] = (Bridge(strongest_bridge['btype'],
+                                                            strongest_bridge['ipos'],
+                                                            strongest_bridge['jpos']))
+                bridges[coord_bridge[1]] = (Bridge(strongest_bridge['btype'],
+                                                            strongest_bridge['jpos'],
+                                                            strongest_bridge['ipos']))
+        if(len(bridges)>0):
+            return(bridges[coord_bridge[0]])
+        else:
+            return(False)
+
+    def get_helix(self, residues, turns):
+        """
+        Return if there is an helix at a given residue,
+        as well as its type.
+        """
+        i = residues.index(self)
+        # if there are no turns or it is the first residue, skip
+        if not turns or i == 0:
+            return False
+        if (i <= len(residues)):
+            if(i in turns.keys() and i-1 in turns.keys()):
+                print(turns[i].turn_type,"- HELIX at", residues[i].resid)
+                return(turns[i].turn_type, residues[i].resid)
+        return(False)
+    
+    def get_ladders(self, residues, ladders):
+        i = residues.index(self)
+        if i != 0:
+            if self.get_bridges(residues):
+                if (residues[i-1].get_bridges(residues)):
+                    local_bridge = self.get_bridges(residues)
+                    consec_bridge = residues[i-1].get_bridges(residues)
+                    if local_bridge.bridge_type == consec_bridge.bridge_type:
+                        print("ladder", consec_bridge.res_num, local_bridge.res_num)
+                        ladders[i] = {'start':i-1,
+                                      'end':i,
+                                      'bridges':[consec_bridge, local_bridge]}
+
+    def get_ladders2(self, bridges, residues):
+        ladders = {}
+        i = residues.index(self)
+        if i != 0:
+            k = 1
+            if i in bridges.keys():
+                temp_bridges = [bridges[i]]
+                while ((i+k in bridges.keys()) and
+                       (bridges[i].bridge_type == bridges[i+k].bridge_type)):
+                    temp_bridges.append(bridges[i+k])
+                    k+=1
+            if k>1:
+                print("ladder", bridges[i].res_num, bridges[i+k-1].res_num)
+                ladders[i] = {'start':bridges[i].res_num,
+                              'end':bridges[i+k-1].res_num,
+                              'bridges':temp_bridges}
+                i+=k-1
+            else:
+                i+=1
+        return ladders
+    
+def get_sheets(ladders):
+    """
+    Bridges between ladders.
+    Check if 1 bridge between one ladder and one or more other ladders.
+    Iterate over all residues of one ladder and check if bridge with other residues
+    of the other ladders.
+    """
+    sheets = {}
+    for ladder in ladders:
+        for ladd2 in ladders:
+            for bridge in ladders[ladder]['bridges']:
+                if bridge.res_partner in res_list(ladders[ladd2]):
+                    print("ladder",ladders[ladder]['start'], ladders[ladder]['end'],"bridge",bridge.res_num, bridge.res_partner,
+                    "ladder 2",ladders[ladd2]['start'], ladders[ladd2]['end'])
+            #print("stop ladder 2")
+        print("stop ladder 1")            
+
+def res_list(ladder):
+    # TODO : method in ladder class
+    l=[]
+    for i in range(ladder['start'], ladder['end']):
+        l.append(i)
+    return(l)
+                   
+def get_bends(residues):
+    bends = {}
+    for i in range(2,len(residues)-2):
+        angle = math.degrees(vector_angles(vectors_substr(position_vector(residues[i].atoms["CA"].coords),
+                                                          position_vector(residues[i-2].atoms["CA"].coords)),
+                                           vectors_substr(position_vector(residues[i+2].atoms["CA"].coords),
+                                                          position_vector(residues[i].atoms["CA"].coords))))
+        if(angle>70):
+            print("angle", residues[i].resid, angle)
+            bends[residues[i].resid] = angle
+    return(bends)
+
+def get_chirality(residues):
+
+    for i in range(1,len(residues)-2):
+        chirality = {}
+        angle = calc_dihedral(residues[i-1].atoms["CA"].coords,
+                              residues[i].atoms["CA"].coords,
+                              residues[i+1].atoms["CA"].coords,
+                              residues[i+2].atoms["CA"].coords)
+        
+        if(angle>0 and angle<=180):
+            sign="+"
+            print("chirality", residues[i].resid, "+", angle)
+            
+        if(angle<=0 and angle > -180):
+            sign="-"
+            print("chirality", residues[i].resid, "-", angle)
+            
+        chirality[residues[i].resid] = sign
+        
+    return chirality
+
+
+def get_phi_psi(residues):
+    for i in range(len(residues)):
+        if(i==0):
+            phi = 360.0
+        else:
+            phi = calc_dihedral(residues[i-1].atoms["C"].coords,
+                                residues[i].atoms["N"].coords,
+                                residues[i].atoms["CA"].coords,
+                                residues[i].atoms["C"].coords)
+        if(i==len(residues)-1):
+            psi = 360.0
+        else:
+            psi = calc_dihedral(residues[i].atoms["N"].coords,
+                                residues[i].atoms["CA"].coords,
+                                residues[i].atoms["C"].coords,
+                                residues[i+1].atoms["N"].coords)
+
+        print("ANGLES", i, phi, psi)
+
+def get_TCO(res1, res2):
+    CO_res1 = vector_from_pos(res1.atoms["C"].coords,
+                              res1.atoms["O"].coords)
+    CO_res2 = vector_from_pos(res2.atoms["C"].coords,
+                              res2.atoms["O"].coords)
+    angle = vector_angles(CO_res1, CO_res2)
+    return(math.cos(angle))

src/dssp.py

@@ -2,6 +2,7 @@ import pdb
 import atom
 import sys
 from structure import *
+from atom import *
 
 def print_dssp():
 
@@ -21,8 +22,8 @@ else:
     #print(get_bridges(pdb_file.residues))
 
 
-    turns = get_turns(pdb_file.residues)
+    #turns = get_turns(pdb_file.residues)
-    get_helix(pdb_file.residues, turns)
+    #get_helix(pdb_file.residues, turns)
     get_bends(pdb_file.residues)
     # get_chirality(pdb_file.residues)
     # bridges = get_bridges(pdb_file.residues)
@@ -47,3 +48,13 @@ else:
     #get_phi_psi(residues)
     #print(residues[2].atoms, residues[0].resid, residues[1].resid)
     print(get_TCO(residues[2],residues[3]))
+
+    turns = {}
+    ladders = {}
+    for i,res in enumerate(residues):
+        # res.get_turns(residues, turns)
+        # res.get_helix(residues, turns)
+        res.get_bridges(residues)
+        res.get_ladders(residues, ladders)
+        
+ 

src/geom.py

@@ -0,0 +1,66 @@
+import math 
+def vector_from_pos(a, b):
+    xAB = b[0]-a[0]
+    yAB = b[1]-a[1]
+    zAB = b[2]-a[2]
+    coord_AB = [xAB,yAB,zAB]
+    return coord_AB
+
+def vector_norm(v):
+    norm = math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
+    return norm
+
+def dot_product(v1,v2):
+    dot_product = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]
+    return(dot_product)
+
+def position_vector(c):
+    vector = vector_from_pos([0,0,0],c)
+    return vector
+
+def vectors_substr(v1, v2):
+    return ([v1[0]-v2[0],
+             v1[1]-v2[1],
+             v1[2]-v2[2]])
+
+def vector_angles(v1,v2):
+    dot_prod = dot_product(v1,v2)
+    norm_v1 = vector_norm(v1)
+    norm_v2 = vector_norm(v2)
+    term = dot_prod/(abs(norm_v1)*abs(norm_v2))
+    rad_angle = math.acos(term)
+    return rad_angle
+
+def calc_dihedral(u1, u2, u3, u4):
+    """ Calculate dihedral angle method. From bioPython.PDB
+    (adapted to np.array)
+    Calculate the dihedral angle between 4 vectors
+    representing 4 connected points. The angle is in
+    [-pi, pi].
+
+    Adapted function of dihedral_angle_from_coordinates.py 
+    by Eric Alcaide.
+    Source : https://gist.github.com/EricAlcaide
+    URL : https://gist.github.com/EricAlcaide/30d6bfdb8358d3a57d010c9a501fda56
+    """
+    #convert coords to numpy arrays
+    u1 = np.array(u1)
+    u2 = np.array(u2)
+    u3 = np.array(u3)
+    u4 = np.array(u4)
+    
+    a1 = u2 - u1
+    a2 = u3 - u2
+    a3 = u4 - u3
+
+    v1 = np.cross(a1, a2)
+    v1 = v1 / (v1 * v1).sum(-1)**0.5
+    v2 = np.cross(a2, a3)
+    v2 = v2 / (v2 * v2).sum(-1)**0.5
+    porm = np.sign((v1 * a3).sum(-1))
+    rad = np.arccos((v1*v2).sum(-1) / ((v1**2).sum(-1) * (v2**2).sum(-1))**0.5)
+    if not porm == 0:
+        rad = rad * porm
+
+    deg_angle = rad*(180/math.pi)
+    return deg_angle

src/geom.py~

@@ -0,0 +1,66 @@
+
+def vector_from_pos(a, b):
+    xAB = b[0]-a[0]
+    yAB = b[1]-a[1]
+    zAB = b[2]-a[2]
+    coord_AB = [xAB,yAB,zAB]
+    return coord_AB
+
+def vector_norm(v):
+    norm = math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
+    return norm
+
+def dot_product(v1,v2):
+    dot_product = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]
+    return(dot_product)
+
+def position_vector(c):
+    vector = vector_from_pos([0,0,0],c)
+    return vector
+
+def vectors_substr(v1, v2):
+    return ([v1[0]-v2[0],
+             v1[1]-v2[1],
+             v1[2]-v2[2]])
+
+def vector_angles(v1,v2):
+    dot_prod = dot_product(v1,v2)
+    norm_v1 = vector_norm(v1)
+    norm_v2 = vector_norm(v2)
+    term = dot_prod/(abs(norm_v1)*abs(norm_v2))
+    rad_angle = math.acos(term)
+    return rad_angle
+
+def calc_dihedral(u1, u2, u3, u4):
+    """ Calculate dihedral angle method. From bioPython.PDB
+    (adapted to np.array)
+    Calculate the dihedral angle between 4 vectors
+    representing 4 connected points. The angle is in
+    [-pi, pi].
+
+    Adapted function of dihedral_angle_from_coordinates.py 
+    by Eric Alcaide.
+    Source : https://gist.github.com/EricAlcaide
+    URL : https://gist.github.com/EricAlcaide/30d6bfdb8358d3a57d010c9a501fda56
+    """
+    #convert coords to numpy arrays
+    u1 = np.array(u1)
+    u2 = np.array(u2)
+    u3 = np.array(u3)
+    u4 = np.array(u4)
+    
+    a1 = u2 - u1
+    a2 = u3 - u2
+    a3 = u4 - u3
+
+    v1 = np.cross(a1, a2)
+    v1 = v1 / (v1 * v1).sum(-1)**0.5
+    v2 = np.cross(a2, a3)
+    v2 = v2 / (v2 * v2).sum(-1)**0.5
+    porm = np.sign((v1 * a3).sum(-1))
+    rad = np.arccos((v1*v2).sum(-1) / ((v1**2).sum(-1) * (v2**2).sum(-1))**0.5)
+    if not porm == 0:
+        rad = rad * porm
+
+    deg_angle = rad*(180/math.pi)
+    return deg_angle

src/structure.py

@@ -1,5 +1,6 @@
 import math
 import numpy as np
+from geom import *
 
 class Structure:
     def __init__(self, res):
@@ -26,304 +27,3 @@ class Helix(Structure):
         self.res_num = res_num
         Structure.res = res_num
  
-def get_turns(residues):
-    # TODO : prevent redondency of overlapping turns
-    turns = {}
-    for i,res in enumerate(residues):
-        for j in range(3,6):
-            if(i+j<len(residues)):
-                if(res.h_bond(residues[i+j])<-0.5):
-                    print(j,"TURN", residues[i].resid, residues[i+j].resid)
-                    turns[i] = Turn(j,residues[i].resid)
-    return(turns)
-
-def get_bridges(residues):
-    bridges = {}
-    bridge = {}
-    strongest_bridge = {}
-    for i in range(1,len(residues)-4):
-        E_min = 0
-        for j in range(i+2,len(residues)-1):
-            # select triplet with the minimal energy 
-
-            if(residues[i-1].h_bond(residues[j])<-0.5
-               and residues[j].h_bond(residues[i+1])<-0.5):
-                bridge = {'res1':residues[i-1].h_bond(residues[j]),
-                          'res2':residues[j].h_bond(residues[i+1]),
-                          'ipos':residues[i].resid,
-                          'jpos':residues[j].resid,
-                          'btype':"para"}
-                   
-            if(residues[j-1].h_bond(residues[i])<-0.5
-               and residues[i].h_bond(residues[j+1])<-0.5):
-                bridge = {'res1':residues[j-1].h_bond(residues[i]),
-                          'res2':residues[i].h_bond(residues[j+1]),
-                          'ipos':residues[i].resid,
-                          'jpos':residues[j].resid,
-                          'btype':"para"}
-
-            if(residues[i].h_bond(residues[j])<-0.5
-               and residues[j].h_bond(residues[i])<-0.5):
-                bridge = {'res1':residues[i].h_bond(residues[j]),
-                          'res2':residues[j].h_bond(residues[i]),
-                          'ipos':residues[i].resid,
-                          'jpos':residues[j].resid,
-                          'btype':"anti"}
-                   
-            if(residues[i-1].h_bond(residues[j+1])<-0.5
-               and residues[j-1].h_bond(residues[i+1])<-0.5):
-                bridge = {'res1':residues[i-1].h_bond(residues[j+1]),
-                          'res2':residues[j-1].h_bond(residues[i+1]),
-                          'ipos':residues[i].resid,
-                          'jpos':residues[j].resid,
-                          'btype':"anti"}
-
-            if(bridge):
-                if(bridge['res1']+bridge['res2']<E_min):
-                    E_min = bridge['res1']+bridge['res2']
-                    strongest_bridge = bridge
-                    bridge = {}
-                    coord_bridge = [i,j]
-        # finally add the strongest bridge at i and j pos
-        if(strongest_bridge):
-            bridges[coord_bridge[0]] = (Bridge(strongest_bridge['btype'],
-                                                        strongest_bridge['ipos'],
-                                                        strongest_bridge['jpos']))
-            bridges[coord_bridge[1]] = (Bridge(strongest_bridge['btype'],
-                                                        strongest_bridge['jpos'],
-                                                        strongest_bridge['ipos']))
-    if(len(bridges)>0):
-        return(bridges)
-    else:
-        return(False)
-
-# def get_bonds(residues):
-#     for i,res1 in enumerate(residues):
-#         E_min = 0
-#         strongest_bridge = []
-#         for j in range(-5,6):
-#             if(i+j < len(residues)):
-#                 res2 = residues[i+j]
-#                 if(res1.h_bond(res2)<-0.5) and (res1.h_bond(res2)<E_min):
-#                     E_min = res1.h_bond(res2)
-#                     strongest_bridge = [res1, res2]
-#         if(len(strongest_bridge)>0):
-#             diff = strongest_bridge[0].resid - strongest_bridge[1].resid
-#             if( abs (diff) > 1 and abs(diff)<=5):
-#                 print(diff)
-
-def get_bonds(residues):
-    bonds = {}
-    k = 0
-    for i,res in enumerate(residues):
-        E_min = 0
-        for j in range(-5,-2):
-            if(i+j<len(residues)):
-                if(residues[i+j].h_bond(res)<-0.5):
-                    k+=1
-                    #if(res.h_bond(residues[i+j])<E_min):
-                    E_min = res.h_bond(residues[i+j])
-                    res_a = res
-                    res_b = residues[i+j]
-                    if j not in bonds:
-                        bonds[j] = []
-                    bonds[j].append([res_a.resid, res_b.resid])
-                    #turns[residues[i].resid] = Turn(j,residues[i].resid)
-    for key in bonds.keys():
-        print(key, len(bonds[key]))
-
-    print("LH",k)
-    
-def get_helix(residues, turns):
-    i = 1
-    helix = []
-    while i <= len(residues):
-        if(i in turns.keys() and i-1 in turns.keys()):
-            k = 0
-            temp_res = []
-            while(i+k in turns):
-                k+=1
-                temp_res.append(residues[i])
-                last_res_pos = residues[i].resid+k
-            if(k>=1):
-                helix_size=last_res_pos - residues[i].resid
-                print(turns[i].turn_type,"- HELIX at", residues[i].resid)
-                helix.append(Helix(temp_res,residues[i].resid))
-            i = i+k
-        else:
-            i+=1
-    return(helix)
-
-def get_ladders(bridges, residues):
-    ladders = {}
-    i = 1
-    while i < len(residues):
-        k = 1
-        if i in bridges.keys():
-            temp_bridges = [bridges[i]]
-            while ((i+k in bridges.keys()) and
-                   (bridges[i].bridge_type == bridges[i+k].bridge_type)):
-                temp_bridges.append(bridges[i+k])
-                k+=1
-        if k>1:
-            ladders[i] = k-1
-            print("ladder", bridges[i].res_num, bridges[i+k-1].res_num)
-            ladders[i] = {'start':bridges[i].res_num,
-                          'end':bridges[i+k-1].res_num,
-                          'bridges':temp_bridges}
-            i+=k-1
-        else:
-            i+=1
-    return ladders
- 
-def get_sheets(ladders):
-    """
-    Bridges between ladders.
-    Check if 1 bridge between one ladder and one or more other ladders.
-    Iterate over all residues of one ladder and check if bridge with other residues
-    of the other ladders.
-    """
-    sheets = {}
-    for ladder in ladders:
-        for ladd2 in ladders:
-            for bridge in ladders[ladder]['bridges']:
-                if bridge.res_partner in res_list(ladders[ladd2]):
-                    print("ladder",ladders[ladder]['start'], ladders[ladder]['end'],"bridge",bridge.res_num, bridge.res_partner,
-                    "ladder 2",ladders[ladd2]['start'], ladders[ladd2]['end'])
-            #print("stop ladder 2")
-        print("stop ladder 1")            
-
-def res_list(ladder):
-    # TODO : method in ladder class
-    l=[]
-    for i in range(ladder['start'], ladder['end']):
-        l.append(i)
-    return(l)
-                   
-def get_bends(residues):
-    bends = {}
-    for i in range(2,len(residues)-2):
-        angle = math.degrees(vector_angles(vectors_substr(position_vector(residues[i].atoms["CA"].coords),
-                                                          position_vector(residues[i-2].atoms["CA"].coords)),
-                                           vectors_substr(position_vector(residues[i+2].atoms["CA"].coords),
-                                                          position_vector(residues[i].atoms["CA"].coords))))
-        if(angle>70):
-            print("angle", residues[i].resid, angle)
-            bends[residues[i].resid] = angle
-    return(bends)
-
-def vector_from_pos(a, b):
-    xAB = b[0]-a[0]
-    yAB = b[1]-a[1]
-    zAB = b[2]-a[2]
-    coord_AB = [xAB,yAB,zAB]
-    return coord_AB
-
-def vector_norm(v):
-    norm = math.sqrt(v[0]**2 + v[1]**2 + v[2]**2)
-    return norm
-
-def dot_product(v1,v2):
-    dot_product = v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]
-    return(dot_product)
-
-def position_vector(c):
-    vector = vector_from_pos([0,0,0],c)
-    return vector
-
-def vectors_substr(v1, v2):
-    return ([v1[0]-v2[0],
-             v1[1]-v2[1],
-             v1[2]-v2[2]])
-
-def vector_angles(v1,v2):
-    dot_prod = dot_product(v1,v2)
-    norm_v1 = vector_norm(v1)
-    norm_v2 = vector_norm(v2)
-    term = dot_prod/(abs(norm_v1)*abs(norm_v2))
-    rad_angle = math.acos(term)
-    return rad_angle
-
-def calc_dihedral(u1, u2, u3, u4):
-    """ Calculate dihedral angle method. From bioPython.PDB
-    (adapted to np.array)
-    Calculate the dihedral angle between 4 vectors
-    representing 4 connected points. The angle is in
-    [-pi, pi].
-
-    Adapted function of dihedral_angle_from_coordinates.py 
-    by Eric Alcaide.
-    Source : https://gist.github.com/EricAlcaide
-    URL : https://gist.github.com/EricAlcaide/30d6bfdb8358d3a57d010c9a501fda56
-    """
-    #convert coords to numpy arrays
-    u1 = np.array(u1)
-    u2 = np.array(u2)
-    u3 = np.array(u3)
-    u4 = np.array(u4)
-    
-    a1 = u2 - u1
-    a2 = u3 - u2
-    a3 = u4 - u3
-
-    v1 = np.cross(a1, a2)
-    v1 = v1 / (v1 * v1).sum(-1)**0.5
-    v2 = np.cross(a2, a3)
-    v2 = v2 / (v2 * v2).sum(-1)**0.5
-    porm = np.sign((v1 * a3).sum(-1))
-    rad = np.arccos((v1*v2).sum(-1) / ((v1**2).sum(-1) * (v2**2).sum(-1))**0.5)
-    if not porm == 0:
-        rad = rad * porm
-
-    deg_angle = rad*(180/math.pi)
-    return deg_angle
-
-
-def get_chirality(residues):
-
-    for i in range(1,len(residues)-2):
-        chirality = {}
-        angle = calc_dihedral(residues[i-1].atoms["CA"].coords,
-                              residues[i].atoms["CA"].coords,
-                              residues[i+1].atoms["CA"].coords,
-                              residues[i+2].atoms["CA"].coords)
-        
-        if(angle>0 and angle<=180):
-            sign="+"
-            print("chirality", residues[i].resid, "+", angle)
-            
-        if(angle<=0 and angle > -180):
-            sign="-"
-            print("chirality", residues[i].resid, "-", angle)
-            
-        chirality[residues[i].resid] = sign
-        
-    return chirality
-
-
-def get_phi_psi(residues):
-    for i in range(len(residues)):
-        if(i==0):
-            phi = 360.0
-        else:
-            phi = calc_dihedral(residues[i-1].atoms["C"].coords,
-                                residues[i].atoms["N"].coords,
-                                residues[i].atoms["CA"].coords,
-                                residues[i].atoms["C"].coords)
-        if(i==len(residues)-1):
-            psi = 360.0
-        else:
-            psi = calc_dihedral(residues[i].atoms["N"].coords,
-                                residues[i].atoms["CA"].coords,
-                                residues[i].atoms["C"].coords,
-                                residues[i+1].atoms["N"].coords)
-
-        print("ANGLES", i, phi, psi)
-
-def get_TCO(res1, res2):
-    CO_res1 = vector_from_pos(res1.atoms["C"].coords,
-                              res1.atoms["O"].coords)
-    CO_res2 = vector_from_pos(res2.atoms["C"].coords,
-                              res2.atoms["O"].coords)
-    angle = vector_angles(CO_res1, CO_res2)
-    return(math.cos(angle))