``````

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:
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:
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

1. O
2. O 1 1.45331355946
3. H 1 0.976162545327 2 96.5718345804
4. 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

deg_bond_angle = (rad_bond_angle / np.pi)*180 # angle as degree

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
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
with open('h2o2.mol2','r') as ex_file:
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

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)

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

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
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
with open('h2o2.mol2','r') as ex_file:
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])

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

deg_bond_angle = (rad_bond_angle / np.pi)*180 # angle as degree

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
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
with open('benzene.mol2','r') as ex_file:
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])

# 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

# 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

``````