In [ ]:

    

from __future__ import division
import math
import numpy as np


    

In [ ]:

    

with open('harmine.mol2','r') as ex_file:
    content = ex_file.readlines()
    atom_counter = 0
    for i in range(len(content)):
        if content[i][0:13] == "@<TRIPOS>ATOM":
            a = i+1
            while content[a][0:9] != '@<TRIPOS>':
                atom_counter += 1
                a += 1
    print 'Number of atoms=',atom_counter


    

In [ ]:

    

atoms = []                                     #
bonds = []
x_array = np.zeros((atom_counter))
y_array = np.zeros((atom_counter))             # This part creates arrays with zeros to store coordinates of atoms
z_array = np.zeros((atom_counter))
bond_mat= np.zeros((atom_counter,atom_counter))#

#print x_array
#print y_array
#print z_array


#overwriting coordinates onto zeros array that i created before
with open('harmine.mol2','r') as ex_file:
    content = ex_file.readlines()
    for i in range(len(content)):
        if content[i][0:13] == "@<TRIPOS>ATOM":
            a = i+1
            at_count = 0
            while content[a][0:9] != '@<TRIPOS>':
                split_content = content[a].split()
                atoms.append(split_content[1])
                x_array[at_count] = split_content[2]
                y_array[at_count] = split_content[3]
                z_array[at_count] = split_content[4]
                a += 1
                at_count += 1
                #print split_content
        
        #storing bond information into a list to write into matrix later        
        elif content[i][0:13] == "@<TRIPOS>BOND":
            b = i+1
            #print b,type(b),content[b],type(content[b])
            while b < len(content) and content[b][0:9] != '@<TRIPOS>':
                #print content[b]
                split2_content= content[b].split()
                bonds.append([int(split2_content[1])-1,int(split2_content[2])-1,split2_content[3]])
                b += 1

                
#print x_array
#print y_array
#print z_array

#Changing string bond names to integer bond names in onder to be able write them into matrix
for item in bonds:
    if item[2] == 'am':
        item[2] = int(4)
    
    elif item[2] == 'ar':
        item[2] = int(5)
    
    elif item[2] == 'du':
        item[2] = int(6)
        
    elif item[2] == 'un':
        item[2] = int(7)
        
    elif item[2] == 'nc':
        item[2] = int(8)
        
for item2 in bonds:
    item2[2] = int(item2[2])

    
#print bonds


    

In [ ]:

    

#writing bond information into matrix
for i in range(len(bonds)):
    mi = bonds[i][0]
    mj = bonds[i][1]
    mt = bonds[i][2]
    
    #print mi,mj,mt
    
    bond_mat[mi][mj]=mt
    bond_mat[mj][mi]=mt

print 'bond_matrix: ';print bond_mat


    

In [ ]:

    

ox_a=np.zeros((3)) #empty array for x coordinates of atoms
oy_a=np.zeros((3)) #empty array for y coordinates of atoms
oz_a=np.zeros((3)) #empty array for z coordinates of atoms

with open('h2o.xyz','r') as co2:
    co = 0 #line counter
    for lines in co2: # atom file with xyz coordinates
        lines=lines.split()
        print 'Lines in file:',lines
        ox_a[co] = float(lines[1]) # writing x coordinates into arrays
        oy_a[co] = float(lines[2]) # writing y coordinates into arrays
        oz_a[co] = float(lines[3]) # writing z coordinates into arrays
        co += 1
        
#print 'ox_a',ox_a
#print 'oy_a',oy_a
#print 'oz_a',oz_a

# for h-o-h molecule, there are 2 bond vectors oh1 and oh2
#
#o-h1 and o-h2 bond vectors are creating with taking oxygen as starting point and hydrogens as terminal points
oh1 = np.array([ox_a[1]-ox_a[0],oy_a[1]-oy_a[0],oz_a[1]-oz_a[0]])
oh2 = np.array([ox_a[2]-ox_a[0],oy_a[2]-oy_a[0],oz_a[2]-oz_a[0]])
#print 'oh1',oh1
#print 'oh2',oh2

#calculating magnitude of o-h1 vector
mag_oh1= np.sqrt(oh1[0]**2+oh1[1]**2+oh1[2]**2)
#print 'mag_oh1',mag_oh1

#normalizing o-h1 vector by dividing to magnitude
unit_oh1 = oh1 / mag_oh1
#print 'unit_oh1',unit_oh1

#calculating magnitude of o-h2 vector
mag_oh2= np.sqrt(oh2[0]**2+oh2[1]**2+oh2[2]**2)
#print 'mag_oh2',mag_oh2

#normalizing o-h2 vector by dividing to magnitude
unit_oh2 = oh2 / mag_oh2
#print 'unit_oh2',unit_oh2

#getting dot product of vectors
dot_oh1_oh2 = np.dot(unit_oh1,unit_oh2)
#print 'dot-oh1-oh2:',dot_oh1_oh2

#getting arc cos of dot product to find bond angle between oh1 and oh2
angle_hoh = np.arccos(dot_oh1_oh2)
print 'angle_hoh:',(angle_hoh/np.pi)*180


    

In [ ]:

    

o2x_a=np.zeros((4)) #empty array for x coordinates of atoms
o2y_a=np.zeros((4)) #empty array for y coordinates of atoms
o2z_a=np.zeros((4)) #empty array for z coordinates of atoms

with open('h2o2.xyz','r') as h2o2: # atom file with xyz coordinates
    h2o2c = 0 #line counter
    for lines2 in h2o2:
        lines2=lines2.split()
        print 'Lines in file:',lines2
        o2x_a[h2o2c] = float(lines2[1]) # writing x coordinates into arrays
        o2y_a[h2o2c] = float(lines2[2]) # writing y coordinates into arrays
        o2z_a[h2o2c] = float(lines2[3]) # writing z coordinates into arrays
        h2o2c += 1
        
#print 'o2x_a',o2x_a
#print 'o2y_a',o2y_a
#print 'o2z_a',o2z_a

# for h-o-o-h molecule, there are 4 bonds vectors:
# 
# taking first oxygen as starting point, oh1 which vector ends at h and oo2 which vector end at second oxygen
# taking second oxygen as starting point, oo3 and oh4 with the same manner.
# 
# --oo2 and oo3 reverse of each other-- 
#
#creating oh1 and oo2 bond vectors
oh1 = np.array([o2x_a[2]-o2x_a[0],o2y_a[2]-o2y_a[0],o2z_a[2]-o2z_a[0]])
oo2 = np.array([o2x_a[1]-o2x_a[0],o2y_a[1]-o2y_a[0],o2z_a[1]-o2z_a[0]])
print 'oh1',oh1
print 'oo2',oo2

#calculating magnitude of oh1 vector
mag_oh1= np.sqrt(oh1[0]**2+oh1[1]**2+oh1[2]**2)
#print 'mag_oh1',mag_oh1

#normalizing oh1 by dividing to  its magnitude
unit_oh1 = oh1 / mag_oh1
#print 'unit_oh1',unit_oh1

#calculating magnitude of oo2 vector
mag_oo2= np.sqrt(oo2[0]**2+oo2[1]**2+oo2[2]**2)
#print 'mag_oo2',mag_oo2

#normalizing oo2 by dividing to  its magnitude
unit_oo2 = oo2 / mag_oo2
#print 'unit_oo2',unit_oo2

#getting dot product  of oh1 and oo2 
dot_oh1_oo2 = np.dot(unit_oh1,unit_oo2)
#print 'dot-oh1-oo2:',dot_oh1_oo2

#calculating bond angle between oh1 and oo2
angle_hoo = np.arccos(dot_oh1_oo2)
print 'angle_hoo:',(angle_hoo/np.pi)*180

#creating oo3 oh4 bond vectors
oo3 = np.array([o2x_a[0]-o2x_a[1],o2y_a[0]-o2y_a[1],o2z_a[0]-o2z_a[1]])
oh4 = np.array([o2x_a[3]-o2x_a[1],o2y_a[3]-o2y_a[1],o2z_a[3]-o2z_a[1]])
print 'oo3 :',oo3
print 'oh4 :',oh4

#calculating magnitude of oo3 bond vector
mag_oo3= np.sqrt(oo3[0]**2+oo3[1]**2+oo3[2]**2)
#print 'mag_oo3',mag_oo3

#normalizing oo3 by dividing to its magnitude
unit_oo3 = oo3 / mag_oo3
#print 'unit_oo3',unit_oo3

#calculating magnitude of oh4 bond vector
mag_oh4= np.sqrt(oh4[0]**2+oh4[1]**2+oh4[2]**2)
#print 'mag_oh4',mag_oh4

#normalizing oh4 by dividing to its magnitude
unit_oh4 = oh4 / mag_oh4
#print 'unit_oh4',unit_oh4

# getting dot product of oo3 and  oh4 bond vectos
dot_oo3_oh4 = np.dot(unit_oo3,unit_oh4)
#print 'dot-oo3-oh4:',dot_oo3_oh4

#calculating bond angle between oo3 and oh4
angle_ooh = np.arccos(dot_oo3_oh4)
print 'angle_ooh:',(angle_ooh/np.pi)*180

#calculating torsion angle according to formula given at above in the image
torsion_hooh = np.arccos((np.dot(np.cross(unit_oh1,unit_oo2),np.cross(unit_oo3,unit_oh4)))/(np.sin(angle_hoo)*np.sin(angle_ooh)))

sig = np.dot(np.cross(unit_oh1,unit_oo2),unit_oh4)
sign = 2*float(sig < 0) -1

print 'torsion :' ,(torsion_hooh/np.pi)*180*sign


    

In [ ]:

    

from __future__ import division
import numpy as np

#--------------------------------------------------------------------------
o2x_a=np.zeros((4)) #empty array for x coordinates of atoms
o2y_a=np.zeros((4)) #empty array for y coordinates of atoms
o2z_a=np.zeros((4)) #empty array for z coordinates of atoms                 # This part of code will remove 
atoms = []                                                                  # when functions done

with open('h2o2.xyz','r') as h2o2: # atom file with xyz coordinates          
    h2o2c = 0 #line counter
    for lines2 in h2o2:
        lines2=lines2.split()
        #print 'Lines in file:',lines2
        atoms.extend(lines2[0])
        o2x_a[h2o2c] = float(lines2[1]) # writing x coordinates into arrays
        o2y_a[h2o2c] = float(lines2[2]) # writing y coordinates into arrays
        o2z_a[h2o2c] = float(lines2[3]) # writing z coordinates into arrays
        h2o2c += 1

#--------------------------------------------------------------------------        
        
        
#
# This part of code will add into mol2 read part of script when its working as desired, 
# o2x,o2y,o2z will change to x_array,y_array,z_array relatively
#

def distance(atom1,atom2): #give atom number as input 0 to n-1 for n atom
    x = o2x_a[atom2]-o2x_a[atom1] # x component of vector that goes to atom2 from atom1
    y = o2y_a[atom2]-o2y_a[atom1] # y component of vector that goes to atom2 from atom1
    z = o2z_a[atom2]-o2z_a[atom1] # z component of vector that goes to atom2 from atom1
    r = np.sqrt(x**2+y**2+z**2)
    return [r,[x,y,z]]# [distance between given atoms,[vector 1-2]]



def angle(atom1,atom2,atom3):
    info21 = distance(atom2,atom1)
    info23 = distance(atom2,atom3)
    
    v1 = info21[1] # vector which goes atom2 from atom1
    v1 = np.asarray(v1)
    v2 = info23[1] # vector which goes atom3 from atom1
    v2 = np.asarray(v2)
    
    
    r1 = info21[0] # magnitude of vector that goes to atom2 from atom1
    r2 = info23[0] # magnitude of vector that goes to atom3 from atom1
    
    unit_v1 = v1/r1 # normalized v1
    unit_v2 = v2/r2 # normalized v2
    
    dot12 = np.dot(unit_v1,unit_v2) # dot product of vectors
    
    rad_bond_angle = np.arccos(dot12) # angle as radian
    
    deg_bond_angle = (rad_bond_angle / np.pi)*180 # angle as degree
    
    print 'for atoms %i %i %i\n'%(atom1,atom2,atom3),'#angle_degree,angle_radian,vector,vector\n',[deg_bond_angle,rad_bond_angle,unit_v1,unit_v2]
    return [deg_bond_angle,rad_bond_angle,unit_v1,unit_v2] #[angle_degree,angle_radian,vector,vector]



def torsion(atom1,atom2,atom3,atom4):

    t1 = angle(atom1,atom2,atom3)
    t2 = angle(atom2,atom3,atom4)
    torsion = np.arccos((np.dot(np.cross(t1[2],t1[3]),np.cross(t2[2],t2[3])))/(np.sin(t1[1])*np.sin(t2[1])))
    
    sign = 2*float(np.dot(np.cross(t1[2],t1[3]),t2[3]) < 0) -1
    res = (torsion/np.pi)*180*sign
    print 'for atoms %i %i %i %i torsion is %f'%(atom1,atom2,atom3,atom4,res)
    return res

a =torsion(2,0,1,3)


    

In [ ]:

    

from __future__ import division
import numpy as np

with open('h2o2.mol2','r') as ex_file:             # count atoms in file
    content = ex_file.readlines()
    atom_counter = 0
    for i in range(len(content)):
        if content[i][0:13] == "@<TRIPOS>ATOM":
            a = i+1
            while content[a][0:9] != '@<TRIPOS>':
                atom_counter += 1
                a += 1
    print 'Number of atoms=',atom_counter
                                                    # create arrays for coordinates
atoms = []
bonds = []
x_array = np.zeros((atom_counter))
y_array = np.zeros((atom_counter))
z_array = np.zeros((atom_counter))
bond_mat= np.zeros((atom_counter,atom_counter))

#print x_array
#print y_array
#print z_array
                                                    # read lines from file
with open('h2o2.mol2','r') as ex_file:
    content = ex_file.readlines()
    for i in range(len(content)):
        #print content[i]
        if content[i][0:13] == "@<TRIPOS>ATOM":
            a = i+1
            at_count = 0                           # write coordinates into arrays
            while content[a][0:9] != '@<TRIPOS>': 
                print content[a]
                split_content = content[a].split()
                atoms.append(split_content[1])
                x_array[at_count] = float(split_content[2])
                y_array[at_count] = float(split_content[3])
                z_array[at_count] = float(split_content[4])
                a += 1
                at_count += 1
                
                                                    # writing bonds
        elif content[i][0:13] == "@<TRIPOS>BOND":
            print 'bonds in file :'
            b = i+1
            #print b,type(b),content[b],type(content[b])
            while b < len(content) and content[b][0:9] != '@<TRIPOS>':
                print content[b]
                split2_content= content[b].split()
                bonds.append([int(split2_content[1])-1,int(split2_content[2])-1,split2_content[3]])
                b += 1

                
#print x_array
#print y_array
#print z_array
                                        #convert bond names into integer counterparts 
for item in bonds:
    if item[2] == 'am':
        item[2] = int(4)
    
    elif item[2] == 'ar':
        item[2] = int(5)
    
    elif item[2] == 'du':
        item[2] = int(6)
        
    elif item[2] == 'un':
        item[2] = int(7)
        
    elif item[2] == 'nc':
        item[2] = int(8)
        
for item2 in bonds:
    item2[2] = int(item2[2])

    
print 'bonds without bond number and taking first atom as atom 0 (every atom number -1)'
for i in bonds:
    print i
print '\n'
                                    # write bonds into matrix
for i in range(len(bonds)):
    mi = bonds[i][0]
    mj = bonds[i][1]
    mt = bonds[i][2]
    
    #print mi,mj,mt
    
    bond_mat[mi][mj]=mt
    bond_mat[mj][mi]=mt

print 'bond_matrix: ';print bond_mat,'\n'


def distance(atom1,atom2): #give atom number as input 0 to n-1 for n atom
    x = x_array[atom2]-x_array[atom1] # x component of vector that goes to atom2 from atom1
    y = y_array[atom2]-y_array[atom1] # y component of vector that goes to atom2 from atom1
    z = z_array[atom2]-z_array[atom1] # z component of vector that goes to atom2 from atom1
    r = np.sqrt(x**2+y**2+z**2)
    return [r,[x,y,z]]# [distance between given atoms,[vector 1-2]]



def angle(atom1,atom2,atom3):
    info21 = distance(atom2,atom1)
    info23 = distance(atom2,atom3)
    
    v1 = info21[1] # vector which goes atom2 from atom1
    v1 = np.asarray(v1)
    v2 = info23[1] # vector which goes atom3 from atom1
    v2 = np.asarray(v2)
    
    
    r1 = info21[0] # magnitude of vector that goes to atom2 from atom1
    r2 = info23[0] # magnitude of vector that goes to atom3 from atom1
    
    unit_v1 = v1/r1 # normalized v1
    unit_v2 = v2/r2 # normalized v2
    
    dot12 = np.dot(unit_v1,unit_v2) # dot product of vectors
    
    rad_bond_angle = np.arccos(dot12) # angle as radian
    
    deg_bond_angle = (rad_bond_angle / np.pi)*180 # angle as degree
    
    print 'angle_degree,angle_radian,vector21,vector23 (for input atom 123)'
    print 'for atoms %i %i %i'%(atom1,atom2,atom3)
    print deg_bond_angle,'\n',rad_bond_angle,'\n',unit_v1,'\n',unit_v2,'\n'
    return [deg_bond_angle,rad_bond_angle,unit_v1,unit_v2] #[angle_degree,angle_radian,vector21,vector23]



def torsion(atom1,atom2,atom3,atom4):

    t1 = angle(atom1,atom2,atom3)
    t2 = angle(atom2,atom3,atom4)
    torsion = np.arccos((np.dot(np.cross(t1[2],t1[3]),np.cross(t2[2],t2[3])))/(np.sin(t1[1])*np.sin(t2[1])))
    
    sign = 2*float(np.dot(np.cross(t1[2],t1[3]),t2[3]) < 0) -1
    res = (torsion/np.pi)*180*sign
    print 'for atoms %i %i %i %i torsion is %f\n'%(atom1,atom2,atom3,atom4,res)
    return res

a =torsion(2,0,1,3)


    

In [ ]:

    

# longer versions of same functions with few differances
# torsion can run seperately from angle unlike the above

def distance(atom1,atom2):
    x_diff = x_array[atom2] - x_array[atom1]
    y_diff = y_array[atom2] - y_array[atom1] 
    z_diff = z_array[atom2] - z_array[atom1]
    distance = np.sqrt(x_diff**2+y_diff**2+z_diff**2)
    return distance

def vector(atom1,atom2):
    vector=np.zeros((3))
    vector[0] = (x_array[atom2] - x_array[atom1])
    vector[1] = (y_array[atom2] - y_array[atom1])
    vector[2] = (z_array[atom2] - z_array[atom1])
    return vector

def angle(atom1,atom2,atom3):
    vector_21 = vector(atom2,atom1)
    vector_23 = vector(atom2,atom3)
    dp21_23 = np.dot(vector_21,vector_23)
    distance_21 = distance(atom2,atom1)
    distance_23 = distance(atom2,atom3)
    angle = (np.arccos(dp21_23/(distance_21*distance_23))/np.pi)*180
    return angle

def torsion(atom1,atom2,atom3,atom4):
    vector_21 = vector(atom2,atom1)
    vector_23 = vector(atom2,atom3)
    distance_21 = distance(atom2,atom1)
    distance_23 = distance(atom2,atom3)
    cp_21_23 = np.cross((vector_21/distance_21),(vector_23/distance_23))
    dp21_23 = np.dot(vector_21,vector_23)
    angle_r123 = np.arccos(dp21_23/(distance_21*distance_23))
    sin_123= np.sin(angle_r123)
    
    vector_32 = vector(atom3,atom2)
    vector_34 = vector(atom3,atom4)
    distance_32 = distance(atom3,atom2)
    distance_34 = distance(atom3,atom4)
    cp_32_34 = np.cross((vector_32/distance_32),(vector_34/distance_34))
    dp32_34 = np.dot(vector_32,vector_34)
    angle_r234 = np.arccos(dp32_34/(distance_32*distance_34))
    sin_234= np.sin(angle_r234)
    
    sign = 2 * float(np.dot(cp_21_23, (vector_34/distance_34)) < 0) - 1
    
    torsion = np.arccos((np.dot(cp_21_23,cp_32_34))/(sin_123*sin_234))*(180/np.pi)*sign
    
    return torsion
    
    
    
    
print torsion(2,0,1,3)


    

In [ ]:

    

#searching through bond mat to find i-j couples that m[i][j] is not equal to 0
#than we can say that they are bonded
def bonds():
    print 'bonds in order'
    for i in range(len(bond_mat)): 
        for j in range(len(bond_mat)):
            if bond_mat[i][j] != 0:
                if i<j:
                    print i,'-',j
    print '\n'




# each row and column of matrix represents atom with same index
# if m[i][j] and m[i][j+n] are not 0 that means atom with index i has 2 bonds
# atom i bonded to atom j and atom j+n
# that means there is an angle between  j--i--j+n
def angles():
    print 'angles in order'
    for j in range(len(bond_mat)):
        bond_list=[]
        for a in range(len(bond_mat)):
            if bond_mat[j][a] != 0:
                bond_list.append(a)
        #print angle_list
        if len(bond_list) >= 2:
            for b in range(len(bond_list)):
                i = bond_list[b]
                for c in range(b+1,len(bond_list)):
                    k = bond_list[c]
                    print i,j,k
    print '\n'


#assume that i-j-k-l are different consecutive atoms that are bonded. 
#starting from atom j searches through bond matrix to look for k atoms that bonded to atom j
#if index number of k greater than index number of j, it searches bond matrix to find bonds of atom k
#then picks another atom that bonded to j which is i. if i is not k
# it picks an atom from bonds of k which is l
# if l is an different atom from j or i
# i-j-k-l establish a torsion angle 
def torsions():
    print 'torsion in order'
    for j in range(len(bond_mat)):
        bonds_j = []
        for a in range(len(bond_mat)):
            if (bond_mat[j][a] != 0):
                bonds_j.append(a)
        for b in range(len(bonds_j)):
            k = bonds_j[b]
            if (k<j):
                continue
            bonds_k = []
            for c in range(len(bond_mat)):
                if (bond_mat[k][c] != 0):
                    bonds_k.append(c)
            for d in range(len(bonds_j)):
                i = bonds_j[d]
                if (i==k):
                    continue
                for e in range(len(bonds_k)):
                    l = bonds_k[e]
                    if (l == j or l == i):
                        continue
                    print i,j,k,l
                
            
    
    print '\n'

torsions()


    

In [ ]:

    

from __future__ import division
import numpy as np

with open('h2o2.mol2','r') as ex_file:             # count atoms in file
    content = ex_file.readlines()
    atom_counter = 0
    for i in range(len(content)):
        if content[i][0:13] == "@<TRIPOS>ATOM":
            a = i+1
            while content[a][0:9] != '@<TRIPOS>':
                atom_counter += 1
                a += 1
    print 'Number of atoms=',atom_counter
                                                    # create arrays for coordinates
atoms = []
bond = []
x_array = np.zeros((atom_counter))
y_array = np.zeros((atom_counter))
z_array = np.zeros((atom_counter))
bond_mat= np.zeros((atom_counter,atom_counter))

#print x_array
#print y_array
#print z_array
                                                    # read lines from file
with open('h2o2.mol2','r') as ex_file:
    content = ex_file.readlines()
    for i in range(len(content)):
        #print content[i]
        if content[i][0:13] == "@<TRIPOS>ATOM":
            a = i+1
            at_count = 0                           # write coordinates into arrays
            while content[a][0:9] != '@<TRIPOS>': 
                print content[a]
                split_content = content[a].split()
                atoms.append(split_content[1])
                x_array[at_count] = float(split_content[2])
                y_array[at_count] = float(split_content[3])
                z_array[at_count] = float(split_content[4])
                a += 1
                at_count += 1
                
                                                    # writing bonds
        elif content[i][0:13] == "@<TRIPOS>BOND":
            print 'bonds in file :'
            b = i+1
            #print b,type(b),content[b],type(content[b])
            while b < len(content) and content[b][0:9] != '@<TRIPOS>':
                print content[b]
                split2_content= content[b].split()
                bond.append([int(split2_content[1])-1,int(split2_content[2])-1,split2_content[3]])
                b += 1

                
#print x_array
#print y_array
#print z_array
                                        #convert bond names into integer counterparts 
for item in bond:
    if item[2] == 'am':
        item[2] = int(4)
    
    elif item[2] == 'ar':
        item[2] = int(5)
    
    elif item[2] == 'du':
        item[2] = int(6)
        
    elif item[2] == 'un':
        item[2] = int(7)
        
    elif item[2] == 'nc':
        item[2] = int(8)
        
for item2 in bond:
    item2[2] = int(item2[2])

    
print 'bonds without bond number and taking first atom as atom 0 (every atom number -1)'
for i in bond:
    print i
print '\n'
                                    # write bonds into matrix
for i in range(len(bond)):
    mi = bond[i][0]
    mj = bond[i][1]
    mt = bond[i][2]
    
    #print mi,mj,mt
    
    bond_mat[mi][mj]=mt
    bond_mat[mj][mi]=mt

print 'bond_matrix: ';print bond_mat,'\n'



#searching through bond mat to find i-j couples that m[i][j] is not equal to 0
#than we can say that they are bonded
def bonds():
    bonds_l=[]
    for i in range(len(bond_mat)): 
        for j in range(len(bond_mat)):
            if bond_mat[i][j] != 0:
                if i<j:
                    bonds_l.append([i,j])
    return bonds_l


# each row and column of matrix represents atom with same index
# if m[i][j] and m[i][j+n] are not 0 that means atom with index i has 2 bonds
# atom i bonded to atom j and atom j+n
# that means there is an angle between  j--i--j+n
def angles():
    angles_l=[]
    for j in range(len(bond_mat)):
        bond_list=[]
        for a in range(len(bond_mat)):
            if bond_mat[j][a] != 0:
                bond_list.append(a)
        #print angle_list
        if len(bond_list) >= 2:
            for b in range(len(bond_list)):
                i = bond_list[b]
                for c in range(b+1,len(bond_list)):
                    k = bond_list[c]
                    angles_l.append([i,j,k])
    return angles_l


#assume that i-j-k-l are different consecutive atoms that are bonded. 
#starting from atom j searches through bond matrix to look for k atoms that bonded to atom j
#if index number of k greater than index number of j, it searches bond matrix to find bonds of atom k
#then picks another atom that bonded to j which is i. if i is not k
# it picks an atom from bonds of k which is l
# if l is an different atom from j or i
# i-j-k-l establish a torsion angle 
def torsions():
    torsions_l=[]
    for j in range(len(bond_mat)):
        bonds_j = []
        for a in range(len(bond_mat)):
            if (bond_mat[j][a] != 0):
                bonds_j.append(a)
        for b in range(len(bonds_j)):
            k = bonds_j[b]
            if (k<j):
                continue
            bonds_k = []
            for c in range(len(bond_mat)):
                if (bond_mat[k][c] != 0):
                    bonds_k.append(c)
            for d in range(len(bonds_j)):
                i = bonds_j[d]
                if (i==k):
                    continue
                for e in range(len(bonds_k)):
                    l = bonds_k[e]
                    if (l == j or l == i):
                        continue
                    torsions_l.append([i,j,k,l]) 
    return torsions_l



def distance(atom1,atom2): #give atom number as input 0 to n-1 for n atom
    x = x_array[atom2]-x_array[atom1] # x component of vector that goes to atom2 from atom1
    y = y_array[atom2]-y_array[atom1] # y component of vector that goes to atom2 from atom1
    z = z_array[atom2]-z_array[atom1] # z component of vector that goes to atom2 from atom1
    r = np.sqrt(x**2+y**2+z**2)
    return [r,[x,y,z]]# [distance between given atoms,[vector 1-2]]



def angle(atom1,atom2,atom3):
    info21 = distance(atom2,atom1)
    info23 = distance(atom2,atom3)
    
    v1 = info21[1] # vector which goes atom2 from atom1
    v1 = np.asarray(v1)
    v2 = info23[1] # vector which goes atom3 from atom1
    v2 = np.asarray(v2)
    
    
    r1 = info21[0] # magnitude of vector that goes to atom2 from atom1
    r2 = info23[0] # magnitude of vector that goes to atom3 from atom1
    
    unit_v1 = v1/r1 # normalized v1
    unit_v2 = v2/r2 # normalized v2
    
    dot12 = np.dot(unit_v1,unit_v2) # dot product of vectors
    
    rad_bond_angle = np.arccos(dot12) # angle as radian
    
    deg_bond_angle = (rad_bond_angle / np.pi)*180 # angle as degree
    
    return [deg_bond_angle,rad_bond_angle,unit_v1,unit_v2] #[angle_degree,angle_radian,vector21,vector23]



def torsion(atom1,atom2,atom3,atom4):

    t1 = angle(atom1,atom2,atom3)
    t2 = angle(atom2,atom3,atom4)
    torsion = np.arccos((np.dot(np.cross(t1[2],t1[3]),np.cross(t2[2],t2[3])))/(np.sin(t1[1])*np.sin(t2[1])))
    
    sign = 2*float(np.dot(np.cross(t1[2],t1[3]),t2[3]) < 0) -1
    res = (torsion/np.pi)*180*sign
    return res

def z_mat():
    out = [[[] for row in range(7)] for atom in range(len(bond_mat))]
    for name in range(len(atoms)):
        out[name][0] = atoms[name]
        if (name == 0):
            pass
        elif (name == 1):
            for aa in bonds_l:
                if (aa[-1] == name):
                    out[name][1] = aa[-2]+1
                    out[name][2] = distance(aa[-1],aa[-2])[0]
        elif (name == 2):
            for ab in angles_l:
                if (ab[-1] == name):
                    out[name][1] = ab[-2]+1
                    out[name][2] = distance(ab[-1],ab[-2])[0]
                    out[name][3] = ab[-3]+1
                    out[name][4] = angle(ab[-1],ab[-2],ab[-3])[0]
        elif (name >= 3):
            for ac in torsions_l:
                if (ac[-1] == name):
                    out[name][1] = ac[-2]+1
                    out[name][2] = distance(ac[-1],ac[-2])[0]
                    out[name][3] = ac[-3]+1
                    out[name][4] = angle(ac[-1],ac[-2],ac[-3])[0]
                    out[name][5] = ac[-4]+1
                    out[name][6] = torsion(ac[-1],ac[-2],ac[-3],ac[-4])
                    

    return out

bonds_l = bonds()

angles_l = angles()

torsions_l = torsions()

z_mat_o = z_mat()




print 'torsions=',torsions_l,'\n'
print 'angles=',angles_l,'\n'
print 'bonds=',bonds_l,'\n'

print ' Z-mat :'
for p1 in z_mat_o:
    print p1


    

In [1]:

    

from __future__ import division
import numpy as np

with open('benzene.mol2','r') as ex_file:             # count atoms in file
    content = ex_file.readlines()
    atom_counter = 0
    for i in range(len(content)):
        if content[i][0:13] == "@<TRIPOS>ATOM":
            a = i+1
            while content[a][0:9] != '@<TRIPOS>':
                atom_counter += 1
                a += 1
    
                                                    # create arrays for coordinates
atoms = []
bond = []
info = []
x_array = np.zeros((atom_counter))
y_array = np.zeros((atom_counter))
z_array = np.zeros((atom_counter))
bond_mat= np.zeros((atom_counter,atom_counter))
bon=[]
#print x_array
#print y_array
#print z_array
                                                    # read lines from file
with open('benzene.mol2','r') as ex_file:
    content = ex_file.readlines()
    for i in range(len(content)):
        #print content[i]
        if content[i][0:13] == "@<TRIPOS>ATOM":
            print content[i]
            a = i+1
            at_count = 0                           # write coordinates into arrays
            while content[a][0:9] != '@<TRIPOS>': 
                print content[a]
                split_content = content[a].split()
                atoms.append(split_content[1])
                x_array[at_count] = float(split_content[2])
                y_array[at_count] = float(split_content[3])
                z_array[at_count] = float(split_content[4])
                a += 1
                at_count += 1
                
                                                    # writing bonds into list
        elif content[i][0:13] == "@<TRIPOS>BOND":
            print content[i]
            b = i+1
            #print b,type(b),content[b],type(content[b])
            while b < len(content) and content[b][0:9] != '@<TRIPOS>':
                print content[b]
                bon.append(content[b])
                split2_content= content[b].split()
                bond.append([int(split2_content[1])-1,int(split2_content[2])-1,split2_content[3]])
                b += 1
        
        elif content[i][0:17] == "@<TRIPOS>MOLECULE": #writing atom info to list
            print content[i]
            c = i+1
            while c < len(content) and content[c][0:9] != '@<TRIPOS>':
                print content[c]
                info.append(content[c])
                c += 1

                    
#print x_array
#print y_array
#print z_array
                                        #convert bond names into integer counterparts 
for item in bond:
    if item[2] == 'am':
        item[2] = int(4)
    
    elif item[2] == 'ar':
        item[2] = int(5)
    
    elif item[2] == 'du':
        item[2] = int(6)
        
    elif item[2] == 'un':
        item[2] = int(7)
        
    elif item[2] == 'nc':
        item[2] = int(8)
        
for item2 in bond:
    item2[2] = int(item2[2])

                                    # write bonds into matrix
for i in range(len(bond)):
    mi = bond[i][0]
    mj = bond[i][1]
    mt = bond[i][2]
    
    #print mi,mj,mt
    
    bond_mat[mi][mj]=mt
    bond_mat[mj][mi]=mt





#searching through bond mat to find i-j couples that m[i][j] is not equal to 0
#than we can say that they are bonded
def bonds():
    bonds_l=[]
    for i in range(len(bond_mat)): 
        for j in range(len(bond_mat)):
            if bond_mat[i][j] != 0:
                if i<j:
                    bonds_l.append([i,j])
    return bonds_l


# each row and column of matrix represents atom with same index
# if m[i][j] and m[i][j+n] are not 0 that means atom with index i has 2 bonds
# atom i bonded to atom j and atom j+n
# that means there is an angle between  j--i--j+n
def angles():
    angles_l=[]
    for j in range(len(bond_mat)):
        bond_list=[]
        for a in range(len(bond_mat)):
            if bond_mat[j][a] != 0:
                bond_list.append(a)
        #print angle_list
        if len(bond_list) >= 2:
            for b in range(len(bond_list)):
                i = bond_list[b]
                for c in range(b+1,len(bond_list)):
                    k = bond_list[c]
                    angles_l.append([i,j,k])
    return angles_l


#assume that i-j-k-l are different consecutive atoms that are bonded. 
#starting from atom j searches through bond matrix to look for k atoms that bonded to atom j
#if index number of k greater than index number of j, it searches bond matrix to find bonds of atom k
#then picks another atom that bonded to j which is i. if i is not k
# it picks an atom from bonds of k which is l
# if l is an different atom from j or i
# i-j-k-l establish a torsion angle 
def torsions():
    torsions_l=[]
    for j in range(len(bond_mat)):
        bonds_j = []
        for a in range(len(bond_mat)):
            if (bond_mat[j][a] != 0):
                bonds_j.append(a)
        for b in range(len(bonds_j)):
            k = bonds_j[b]
            if (k<j):
                continue
            bonds_k = []
            for c in range(len(bond_mat)):
                if (bond_mat[k][c] != 0):
                    bonds_k.append(c)
            for d in range(len(bonds_j)):
                i = bonds_j[d]
                if (i==k):
                    continue
                for e in range(len(bonds_k)):
                    l = bonds_k[e]
                    if (l == j or l == i):
                        continue
                    torsions_l.append([i,j,k,l]) 
    return torsions_l



# mathematical equations for the functions below are present in the image that placed around starting point of notebook

# returns distance betweem atoms
def distance(atom1,atom2):
    x_diff = x_array[atom2] - x_array[atom1]
    y_diff = y_array[atom2] - y_array[atom1] 
    z_diff = z_array[atom2] - z_array[atom1]
    distance = np.sqrt(x_diff**2+y_diff**2+z_diff**2)
    return distance

# returns a vector that defined by extracting coordinates of atom 2 from atom 1
def vector(atom1,atom2):
    vector=np.zeros((3))
    vector[0] = (x_array[atom2] - x_array[atom1])
    vector[1] = (y_array[atom2] - y_array[atom1])
    vector[2] = (z_array[atom2] - z_array[atom1])
    return vector

#returns angle between 3 given consecutive atom
def angle(atom1,atom2,atom3):
    vector_21 = vector(atom2,atom1)
    vector_23 = vector(atom2,atom3)
    dp21_23 = np.dot(vector_21,vector_23)
    distance_21 = distance(atom2,atom1)
    distance_23 = distance(atom2,atom3)
    angle = (np.arccos(dp21_23/(distance_21*distance_23))/np.pi)*180
    return angle

#returns torsion angle between 4 given consecutive atom
def torsion(atom1,atom2,atom3,atom4):
    vector_21 = vector(atom2,atom1)
    vector_23 = vector(atom2,atom3)
    
    distance_21 = distance(atom2,atom1)
    distance_23 = distance(atom2,atom3)
    
    cp_21_23 = np.cross((vector_21/distance_21),(vector_23/distance_23))
    dp21_23 = np.dot(vector_21,vector_23)
    angle_r123 = np.arccos(dp21_23/(distance_21*distance_23))
    sin_123= np.sin(angle_r123)
    
    vector_32 = vector(atom3,atom2)
    vector_34 = vector(atom3,atom4)
    
    distance_32 = distance(atom3,atom2)
    distance_34 = distance(atom3,atom4)
    
    cp_32_34 = np.cross((vector_32/distance_32),(vector_34/distance_34))
    dp32_34 = np.dot(vector_32,vector_34)
    angle_r234 = np.arccos(dp32_34/(distance_32*distance_34))
    sin_234= np.sin(angle_r234)
    
    sign = 2 * float(np.dot(cp_21_23, (vector_34/distance_34)) < 0) - 1
    
    torsion = np.arccos((np.dot(cp_21_23,cp_32_34))/(sin_123*sin_234))*(180/np.pi)*sign
    
    return torsion
    
    
    
# builds and returns z matrix

def z_mat():
    out = [['' for row in range(7)] for atom in range(len(bond_mat))] #creating empty list (n_atoms)x7 for writing zmat
    for name in range(len(atoms)):
        out[name][0] = atoms[name] #writing bond names into first coloumn
        if (name == 0): #for first atom do nothing
            pass
        elif (name == 1): # for second atom, find a bond that ends with 1 from bond list and gets distance 
            for aa in bonds_l:
                if (aa[-1] == name):
                    out[name][1] = aa[-2]+1
                    out[name][2] = distance(aa[-1],aa[-2])
                    break
            
        elif (name == 2): #for third atom, find a triplet that ends or starts with 2 from angle list  
            for ab in angles_l: # and gets distance , angle respective to triplet order in list
                if (ab[-1] == name): # check if all other indices less than current
                    if (int(ab[-1]) > int(ab[-2])) and (int(ab[-1]) > int(ab[-3])):
                        out[name][1] = ab[-2]+1
                        out[name][2] = distance(ab[-1],ab[-2])
                        out[name][3] = ab[-3]+1
                        out[name][4] = angle(ab[-1],ab[-2],ab[-3])
                        break
                elif (ab[0] == name):
                    if (int(ab[0]) > int(ab[1])) and (int(ab[0]) > int(ab[2])):
                        out[name][1] = ab[1]+1
                        out[name][2] = distance(ab[0],ab[1])
                        out[name][3] = ab[2]+1
                        out[name][4] = angle(ab[0],ab[1],ab[2])
                        break

                    
        elif (name >= 3): #for any atom index >=3 atom, find a quartet ends or starts with that index from torsion list
            for ac in torsions_l: # and gets distance , angle , torsion respective to triplet order in list
                if (ac[-1] == name): # check if all other indices less than current
                    if (int(ac[-1]) > int(ac[-2])) and (int(ac[-1]) > int(ac[-3])) and (int(ac[-1]) > int(ac[-4])):
                        out[name][1] = ac[-2]+1
                        out[name][2] = distance(ac[-1],ac[-2])
                        out[name][3] = ac[-3]+1
                        out[name][4] = angle(ac[-1],ac[-2],ac[-3])
                        out[name][5] = ac[-4]+1
                        out[name][6] = torsion(ac[-1],ac[-2],ac[-3],ac[-4])
                        break
                elif (ac[0] == name):
                    if (int(ac[0]) > int(ac[1])) and (int(ac[0]) > int(ac[2])) and (int(ac[0]) > int(ac[3])):
                        out[name][1] = ac[1]+1
                        out[name][2] = distance(ac[0],ac[1])
                        out[name][3] = ac[2]+1
                        out[name][4] = angle(ac[0],ac[1],ac[2])
                        out[name][5] = ac[3]+1
                        out[name][6] = torsion(ac[0],ac[1],ac[2],ac[3])
                        break
                    

    return out


#priting all  of the output
def printout():
    print 'Number of atoms=',atom_counter
    print 'bond_matrix: ';print bond_mat,'\n'
    print 'bonds='
    for i in bonds_l:
        print i
    print '\n'
    print 'angles='
    for i in angles_l:
        print i
    print '\n'
    print 'torsions='
    for i in torsions_l:
        print i
    print '\n'
    print ' Z-mat :'
    z_mat_string = z_mat_o
    for i in range(len(atoms)):
        for j in range(7):
            z_mat_string[i][j] = str(z_mat_string[i][j])
    for i in z_mat_string:
        i = '\t'.join(i)
        print i
    print '\n'

# running functions
bonds_l = bonds()

angles_l = angles()

torsions_l = torsions()

z_mat_o = z_mat() 

printout()

##########################################################################################################
#Converting back to cartesian from z mat
##########################################################################################################

coord_mat=[[] for p in range(len(z_mat_o))] # creating a epmty list for appending with atom number and coordinates

# conveting spherical coordinates to cartesian coordinates for each atom>3
def bond_vector(distance, angle, torsion):
    x = distance*np.sin(angle)*np.sin(torsion)
    y = distance*np.sin(angle)*np.cos(torsion)
    z = distance*np.cos(angle) 
    return np.array([x,y,z])



#the first two lines in get_local_axes are to build the vectors r⃗ 12 and r⃗ 23. The next line builds the cross product
#r⃗ 23×r⃗ 12|r⃗ 23×r⃗ 12|,
#The next vector is
#r⃗ 12×(r⃗ 23×r⃗ 12)|r⃗ 12×(r⃗ 23×r⃗ 12)|,
#This is a vector perpendicular to r⃗ 12, and lies in the plane where r⃗ 12 and r⃗ 23 lie.
#The next lines are defining the three axes of the system. e⃗ z
#is just r⃗ 12 normalized, e⃗ y is chosen so that yz is the plane 
#where the three points (atoms) locate. and e⃗ x is chosen so that (e⃗ x,e⃗ y,e⃗ z)
#forms a orthonormal set and are the versors of a right-handed coordinate system. 
#This system is the local axis system.
#Explanation belongs to "Weijun Zhou" '
#https://chemistry.stackexchange.com/questions/89620/conversion-of-z-matrix-to-cartesian-coordinates'    
def local_axes(coords_1,coords_2,coords_3):
    v12 = np.array([coords_2[0]-coords_1[0],coords_2[1]-coords_1[1],coords_2[2]-coords_1[2]])
    norm12= v12/np.sqrt(v12[0]**2+v12[1]**2+v12[2]**2)
    
    v23 = np.array([coords_3[0]-coords_2[0],coords_3[1]-coords_2[1],coords_3[2]-coords_2[2]])
    norm23= v23/np.sqrt(v23[0]**2+v23[1]**2+v23[2]**2)

    cross1 = np.cross(norm23,norm12)/ np.sin(np.arccos(np.dot(norm23,norm12)))
    
    cross2 = np.cross(norm12,cross1)/ np.sin(np.arccos(np.dot(norm12,cross1)))

    z = norm12
    y = cross2
    x = np.cross(y, z)/ np.sin(np.arccos(np.dot(y,z)))
    return np.array([x,y,z])

    
#[atom name, bonded atom, bond length,angled atom,angle value,torsined atom, torison value]
#reading z-mat as input and converting back to cartesian row by row
def coor_mat(z_mat_o):
    for row in range(len(z_mat_o)):
        if row == 0:#atom 0
            coord_mat[row].append(row+1)#atom num
            coord_mat[row].append(0)#x component
            coord_mat[row].append(0)#y component
            coord_mat[row].append(0)#z component
        elif row == 1:#atom 1 setting to axe. Value equal to distace (0,0,z)
            z = float(z_mat_o[row][2])
            coord_mat[row].append(row+1)
            coord_mat[row].append(0)
            coord_mat[row].append(0)
            coord_mat[row].append(z)
        elif row == 2:# atom 2
            y = np.sin(float(z_mat_o[row][4])*np.pi/180)*float(z_mat_o[row][2])
            # y value = r * sin (angle 2 - 1 - 0)
            z1 = float(coord_mat[int(z_mat_o[row][1])-1][3]) #z value of atom 1
            z2 = (1-2*float(int(z_mat_o[row][1])-1==1)) # sign depending to if atom 2 bonded to atom 1 or to atom 0
                                                        # to add or subsract new z value from  atom 1 z value
            z3 = float(z_mat_o[row][2])*np.cos(float(z_mat_o[row][4])*np.pi/180) # r * cos (angle 2 1 0) 
            z = z1+z2*z3 # new z value
            
            
            coord_mat[row].append(row+1)
            coord_mat[row].append(0)
            coord_mat[row].append(y)
            coord_mat[row].append(z)
        elif row >= 3: # atom 3 or with index greater than 3
            bn, an, tn = int(z_mat_o[row][1])-1 , int(z_mat_o[row][3])-1 , int(z_mat_o[row][5])-1
            # bond atom number , angle atom number, torsion atom number
            
            coords1= coord_mat[bn][1:4] # coordinates of bond atom number
            coords2= coord_mat[an][1:4] # coordinates of angle atom number
            coords3= coord_mat[tn][1:4] # coordinates of torsion atom number
            localax = local_axes(coords1,coords2,coords3) # calculating local axes from 3 coordinates
            
            bondvec = bond_vector(float(z_mat_o[row][2]),float(z_mat_o[row][4])*np.pi/180,float(z_mat_o[row][6])*np.pi/180)
            #calculating bond vector from (r,agle,torsion) covert sherical coordinates vector with x,y,z values
            # finding direction
            
            displacement = np.dot(bondvec,localax) # (1x1).(3x3) calculating displacement vector
            new_coords = displacement + np.asarray(coord_mat[bn][1:4]) # adding diplacement to coordinates of bonded atom
                                                                        # to find current coordinates
            coord_mat[row].append(row+1)
            coord_mat[row].append(new_coords[0])
            coord_mat[row].append(new_coords[1])
            coord_mat[row].append(new_coords[2])

#running z mat to cartesian converter function            
coor_mat(z_mat_o)


#printing cartesian coordinates
print 'Cartesian coordinates'
for i in range(len(coord_mat)):
    print '%7i %-4s %13.4f %9.4f %9.4f'%(coord_mat[i][0],atoms[i],coord_mat[i][1],coord_mat[i][2],coord_mat[i][3])

# writing new coordinates into examle.mol2 file
with open('example.mol2','w') as f:
    f.write("@<TRIPOS>MOLECULE\n")
    for item in info:
        f.write(item)
    f.write("@<TRIPOS>ATOM\n")
    for i in range(len(coord_mat)):
        f.write('%7i %-4s %13.4f %9.4f %9.4f\n'%(coord_mat[i][0],atoms[i],coord_mat[i][1],coord_mat[i][2],coord_mat[i][3]))
    f.write("@<TRIPOS>BOND\n")
    for item in bon:
        f.write(item)