import sys, math
# # CONSTANTS # #
# threshold beyond average of covalent radiii to determine bond cutoff
bond_thresh = 1.2
# covalent (or ionic) radii by atomic element (Angstroms) from
# "Inorganic Chemistry" 3rd ed, Housecroft, Appendix 6, pgs 1013-1014
cov_rads = { 'H' : 0.37, 'C' : 0.77, 'O' : 0.73, 'N' : 0.75, 'F' : 0.71,
'P' : 1.10, 'S' : 1.03, 'Cl': 0.99, 'Br': 1.14, 'I' : 1.33, 'He': 0.30,
'Ne': 0.84, 'Ar': 1.00, 'Li': 1.02, 'Be': 0.27, 'B' : 0.88, 'Na': 1.02,
'Mg': 0.72, 'Al': 1.30, 'Si': 1.18, 'K' : 1.38, 'Ca': 1.00, 'Sc': 0.75,
'Ti': 0.86, 'V' : 0.79, 'Cr': 0.73, 'Mn': 0.67, 'Fe': 0.61, 'Co': 0.64,
'Ni': 0.55, 'Cu': 0.46, 'Zn': 0.60, 'Ga': 1.22, 'Ge': 1.22, 'As': 1.22,
'Se': 1.17, 'Kr': 1.03, 'X' : 0.00}
# # IO FUNCTIONS # #
# read file data into a 2-d array
def get_file_string_array(file_name):
try:
file = open(file_name, "r")
except IOError:
print('Error: file (%s) not found!\n' % (file_name))
sys.exit()
lines = file.readlines()
file.close()
array = []
for line in lines:
array.append(line.split())
return array
# read in geometry from xyz file
def get_geom(xyz_file_name):
xyz_array = get_file_string_array(xyz_file_name)
n_atoms = int(xyz_array[0][0])
at_types = ['' for i in range(n_atoms)]
coords = [[0.0 for j in range(3)] for i in range(n_atoms)]
for i in range(n_atoms):
at_types[ i] = xyz_array[i+2][0]
for j in range(3):
coords[ i][j] = float(xyz_array[i+2][j+1])
geom = [at_types, coords]
return geom
# input syntax and usage warnings
def get_inputs():
if (not len(sys.argv) == 2):
print('Usage: bonds.py XYZ_FILE\n')
print(' XYZ_FILE: coordinates of target molecule\n')
sys.exit()
else:
xyz_file_name = sys.argv[1]
return xyz_file_name
# print geometry to screen
def print_geom(geom, comment):
at_types, coords = geom[0:2]
n_atoms = len(at_types)
print('%i\n%s\n' % (n_atoms, comment), end='')
for i in range(n_atoms):
print('%-2s' % (at_types[ i]), end='')
for j in range(3):
print(' %12.6f' % (coords[ i][j]), end='')
print('\n', end='')
print('\n', end='')
# print bond graph to screen
def print_bond_graph(geom, bond_graph, comment):
at_types = geom[0]
n_atoms = len(at_types)
print('%s\n' % (comment), end='')
for i in range(n_atoms):
print(' %4i %-2s -' % (i+1, at_types[I]), end='')
for j in range(len(bond_graph[ i])):
print(' %i' % (bond_graph[ i ][j] + 1), end='')
print('\n', end='')
print('\n', end='')
# print list of bond lengths to screen
def print_bonds(geom, bonds):
at_types = geom[0]
n_bonds = len(bonds)
print('%i bond(s) found (Angstrom)' % (n_bonds))
for q in range(n_bonds):
n1, n2 = bonds[q][0:2]
r12 = bonds[q][2]
nstr = '%i-%i' % (n1+1, n2+1)
tstr = '(%s-%s) ' % (at_types[n1], at_types[n2])
print(' %-15s %-13s %6.4f\n' % (nstr, tstr, r12), end='')
print('\n', end='')
# # MATH FUNCTIONS # #
# calculate distance between two 3-d cartesian coordinates
def get_r12(coords1, coords2):
r2 = 0.0
for p in range(3):
r2 += (coords2[p] - coords1[p])**2
r = math.sqrt(r2)
return r
# # TOPOLOGY FUNCTIONS # #
# build graph of which atoms are covalently bonded
def get_bond_graph(geom):
at_types, coords = geom[0:2]
n_atoms = len(at_types)
bond_graph = [[] for i in range(n_atoms)]
for i in range(n_atoms):
covrad1 = cov_rads[at_types[I]]
for j in range(i+1, n_atoms):
covrad2 = cov_rads[at_types[j]]
thresh = bond_thresh * (covrad1 + covrad2)
r12 = get_r12(coords[I], coords[j])
if (r12 < thresh):
bond_graph[I].append(j)
bond_graph[j].append(i)
return bond_graph
# determine atoms which are covalently bonded from bond graph
def get_bonds(geom, bond_graph):
at_types, coords = geom[0:2]
n_atoms = len(at_types)
bonds = []
for i in range(n_atoms):
for a in range(len(bond_graph[ i])):
j = bond_graph[I][a]
if (i < j):
r12 = get_r12(coords[I], coords[j])
bonds.append([i, j, r12])
return bonds
# # MAIN BLOCK # #
# read in geometry, determine bonded topology
xyz_file_name = get_inputs()
geom = get_geom(xyz_file_name)
bond_graph = get_bond_graph(geom)
# calculate bond lengths
bonds = get_bonds(geom, bond_graph)
# print resulting values
print_geom(geom, 'initial geometry')
print_bond_graph(geom, bond_graph, 'bond graph')
print_bonds(geom, bonds)
# end of program