

In [ ]:

    
from __future__ import division
import math
import numpy as np

File of interest will open at down. Number of atoms will count for creating arrays in order to stack x,y,z coordinates.



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

After atom number counted, x,y,z coordinates of atoms will write into arrays.Atom names will store in 'atoms' list. Also bonds between atoms write into list for converting string type bond notations to integer type notations.

Single = 1 Double = 2 Triple = 3 Amide = 4
Aromatic = 5 Dummy = 6 Unknown = 7 Not Connected = 8



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

Bonds between atoms will write into a matrix size of (number of atoms,number of atoms)



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

Image above is from a video belong to TMP Chem channel on youtube named 'Computational Chemistry 1.9 - Torison angles' Calculations

From this part, I will try to find bond angle of water molecule which located on h2o.xyz file as an exercise.



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

Bond angle between 3 non-linear molecule calculated above, now ill try to calculate torsion angle h2o2 molecule which is located on h2o2.xyz file.



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

For given h2o2 molecule example z-matrix should look like

O
O 1 1.45331355946
H 1 0.976162545327 2 96.5718345804
H 2 0.976162545327 1 96.5718345804 3 -180

Funtions below are defined to do same calculations above for h2o2 molecule that is written in h2o2.xyz file.



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)

For getting torsions and angles code below is ready to use. However atom index order as input to functions should be figure out.



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)

Atom order was 2-0-1-3 because molecule is H-O-O-H. To calculate angel or torsion that we looking for, molecules need to be in order while giving their numbers to function and in file atoms order was O(0) O(1) H(2) H(3).

To be able to figure out atom index order ill use the script below



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()

So if we put all of our previous work into one cell it will look like at as below



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

As a reminder for angle/torsion image above added Code below improved for more complex molecules. Few mistakes at z-mat function corrected. Now function limited to define new atoms relative to atoms that are defined before. Also z-matrix to cartesian conversion function and necessary new functions added.



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 also other necessary lists
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)









    



@<TRIPOS>MOLECULE

benzene

 12 12 0 0 0

SMALL

GASTEIGER



@<TRIPOS>ATOM

      1 C          -0.2686   -0.4617   -1.4312 C.ar    1  LIG1       -0.0618

      2 C           0.0131    0.0539   -0.1662 C.ar    1  LIG1       -0.0618

      3 C          -0.2400   -0.7125    0.9715 C.ar    1  LIG1       -0.0618

      4 C          -0.7747   -1.9944    0.8440 C.ar    1  LIG1       -0.0618

      5 C          -1.0564   -2.5100   -0.4210 C.ar    1  LIG1       -0.0618

      6 C          -0.8033   -1.7436   -1.5587 C.ar    1  LIG1       -0.0618

      7 H          -1.0227   -2.1453   -2.5443 H       1  LIG1        0.0618

      8 H          -0.0714    0.1354   -2.3176 H       1  LIG1        0.0618

      9 H           0.4297    1.0527   -0.0669 H       1  LIG1        0.0618

     10 H          -0.0206   -0.3108    1.9571 H       1  LIG1        0.0618

     11 H          -0.9719   -2.5915    1.7304 H       1  LIG1        0.0618

     12 H          -1.4730   -3.5087   -0.5203 H       1  LIG1        0.0618

@<TRIPOS>BOND

     1     7     6    1

     2     8     1    1

     3     6     1   ar

     4     6     5   ar

     5     1     2   ar

     6    12     5    1

     7     5     4   ar

     8     2     9    1

     9     2     3   ar

    10     4     3   ar

    11     4    11    1

    12     3    10    1

Number of atoms= 12
bond_matrix: 
[[ 0.  5.  0.  0.  0.  5.  0.  1.  0.  0.  0.  0.]
 [ 5.  0.  5.  0.  0.  0.  0.  0.  1.  0.  0.  0.]
 [ 0.  5.  0.  5.  0.  0.  0.  0.  0.  1.  0.  0.]
 [ 0.  0.  5.  0.  5.  0.  0.  0.  0.  0.  1.  0.]
 [ 0.  0.  0.  5.  0.  5.  0.  0.  0.  0.  0.  1.]
 [ 5.  0.  0.  0.  5.  0.  1.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.]
 [ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.]] 

bonds=
[0, 1]
[0, 5]
[0, 7]
[1, 2]
[1, 8]
[2, 3]
[2, 9]
[3, 4]
[3, 10]
[4, 5]
[4, 11]
[5, 6]


angles=
[1, 0, 5]
[1, 0, 7]
[5, 0, 7]
[0, 1, 2]
[0, 1, 8]
[2, 1, 8]
[1, 2, 3]
[1, 2, 9]
[3, 2, 9]
[2, 3, 4]
[2, 3, 10]
[4, 3, 10]
[3, 4, 5]
[3, 4, 11]
[5, 4, 11]
[0, 5, 4]
[0, 5, 6]
[4, 5, 6]


torsions=
[5, 0, 1, 2]
[5, 0, 1, 8]
[7, 0, 1, 2]
[7, 0, 1, 8]
[1, 0, 5, 4]
[1, 0, 5, 6]
[7, 0, 5, 4]
[7, 0, 5, 6]
[0, 1, 2, 3]
[0, 1, 2, 9]
[8, 1, 2, 3]
[8, 1, 2, 9]
[1, 2, 3, 4]
[1, 2, 3, 10]
[9, 2, 3, 4]
[9, 2, 3, 10]
[2, 3, 4, 5]
[2, 3, 4, 11]
[10, 3, 4, 5]
[10, 3, 4, 11]
[3, 4, 5, 0]
[3, 4, 5, 6]
[11, 4, 5, 0]
[11, 4, 5, 6]


 Z-mat :
C						
C	1	1.39478430232				
C	2	1.39491571788	1	119.997725694		
C	3	1.39478598717	2	119.997263753	1	0.00305771411986
C	4	1.39478430232	3	120.005010482	2	-0.0030579385561
C	1	1.39478598717	2	120.005010482	3	-0.00305795309681
H	6	1.08669527007	1	120.0027387	2	-179.995898042
H	1	1.08679400532	2	119.998183968	3	-179.999342461
H	2	1.08674628594	1	120.002779054	6	-179.999305377
H	3	1.08669527007	2	119.999997539	1	-179.995898295
H	4	1.08679400532	3	119.996805445	2	-179.999342496
H	5	1.08665437928	4	120.003326469	3	-179.998379974


Cartesian coordinates
      1 C           0.0000    0.0000    0.0000
      2 C           0.0000    0.0000    1.3948
      3 C           0.0000    1.2081    2.0922
      4 C          -0.0001    2.4159    1.3947
      5 C          -0.0001    2.4159   -0.0001
      6 C          -0.0001    1.2079   -0.6975
      7 H          -0.0001    1.2078   -1.7842
      8 H          -0.0000   -0.9412   -0.5434
      9 H           0.0001   -0.9411    1.9382
     10 H           0.0001    1.2081    3.1789
     11 H          -0.0001    3.3571    1.9381
     12 H          -0.0001    3.3570   -0.5435