-
Notifications
You must be signed in to change notification settings - Fork 0
/
tau-calc.py
291 lines (251 loc) · 9.68 KB
/
tau-calc.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
#!/usr/bin/env python3
#
# Calculation of tau_4, tau_4' and tau5 for 4- or 5-coordinated atoms
# For a deeper explanation have a look at Wikipedia:
#
# https://en.wikipedia.org/wiki/Geometry_index
#
# or the following articles:
#
# For tau_4, please cite:
#
# "Structural variation in copper(i) complexes with pyridylmethylamide ligands:
# structural analysis with a new four-coordinate geometry index, τ4"
#
# Lei Yang, Douglas R. Powell, Robert P. Houser
# Dalton Trans. 2007, 955-964.
# DOI: https://doi.org/10.1039/B617136B
#
# For the improved tau_4', please cite:
#
# "Coordination polymers and molecular structures among complexes of
# mercury(II) halides with selected 1-benzoylthioureas"
#
# Andrzej Okuniewski, Damian Rosiak, Jarosław Chojnacki, Barbara Becker
# Polyhedron 2015, 90, 47–57.
# DOI: https://doi.org/10.1016/j.poly.2015.01.035
#
# For tau_5, please cite:
#
# "Synthesis, structure, and spectroscopic properties of copper(II) compounds containing
# nitrogen–sulphur donor ligands; the crystal and molecular structure of
# aqua[1,7-bis(N-methylbenzimidazol-2′-yl)-2,6-dithiaheptane]copper(II) perchlorate"
#
# Anthony W. Addison, T. Nageswara Rao, Jan Reedijk, Jacobus van Rijn, Gerrit C. Verschoor
# J. Chem. Soc., Dalton Trans. 1984, 1349-1356.
# DOI: https://doi.org/10.1039/DT9840001349
#
# Gemmi:
# https://gemmi.readthedocs.io/en/latest/
# https://github.com/project-gemmi/gemmi
#
import sys #sys
import argparse #argument parser
import re #regular expressions
from gemmi import cif #CIF processing
#regex for bonds and angles --> string to float, no esd
ang_bond_val = re.compile('\d{1,}[\.]?\d{0,}')
#list for the exclusion of atoms in the angle table
list_of_atoms_with_large_bonds=[]
#calculation of tau5
def calc_tau5(beta, alpha):
tau5 = (beta - alpha) / 60
return tau5
#calculation of tau4 classic
def calc_tau4(beta, alpha):
tau4 = (360 - (alpha + beta )) / (360 - 2*109.5)
return tau4
#calculation of tau4 improved
def calc_tau4impr(beta, alpha):
tau4impr = (beta - alpha) / (360 - 109.5) + (180 - beta) / (180 - 109.5)
return tau4impr
#argument parser
parser = argparse.ArgumentParser(prog='tau-calc',
description = "Calculation of tau_4, tau_4' and tau_5 geometry indices.")
#filename is required
parser.add_argument("filename",
help = "filename, CIF; e.g. mystructure.cif")
#atom for tau_x calculation
parser.add_argument("atom_name",
type = str,
help = "atom names; e.g. Co1")
#exclude atoms
parser.add_argument('-e','--exclude',
nargs="+",
type=str,
help='exclude bonds to specified atoms; e.g. -e N1 or -e N1 N2')
#exclude atoms by distance
parser.add_argument('-d','--distance',
type=float,
help='exclude atoms with distances larger than d in Å; e.g. -d2.2')
#parse arguments
args = parser.parse_args()
#load cif
try:
doc = cif.read_file(args.filename)
#file not found
except IOError:
print(f"'{args.filename}'" + " not found")
sys.exit(1)
#not a valid cif
except ValueError:
print(f"'{args.filename}'" + " is not a valid CIF. Exit.\n")
sys.exit(1)
#more than one or no data_block --> exit
if len(doc) != 1:
print("CIF should contain one structure or one 'data_' block. Exit.")
sys.exit(1)
#get the block
block = doc.sole_block()
#check if selected atom is in the CIF
if args.atom_name not in list(block.find_loop('_atom_site_label')):
print("The atom is not part of the CIF. Exit.\n")
sys.exit(1)
#check if excluded atoms are in the CIF
if args.exclude:
if set(args.exclude).issubset(list(block.find_loop('_atom_site_label'))) is False:
print("One or more excluded atoms are not part of the CIF. Exit.\n")
sys.exit(1)
#build a bond table atom1-atom2 bond-length symmetry
bond_table=block.find(['_geom_bond_atom_site_label_1',
'_geom_bond_atom_site_label_2',
'_geom_bond_distance',
'_geom_bond_site_symmetry_2'])
#delete bonds not containing the atom from input
for i in range(len(bond_table)-1, -1, -1):
if args.atom_name not in bond_table[i]:
del bond_table[i]
#exclude symmetry eq. atoms, e.g. for polymeric compounds M1-X-M1'
#comment in case of problems
for i in range(len(bond_table)-1, -1, -1):
if args.atom_name in bond_table[i][1] and "." not in bond_table[i][3]:
del bond_table[i]
#exclude bonds to specified atoms
if args.exclude:
for atom in args.exclude:
for i in range(len(bond_table)-1, -1, -1):
if atom in bond_table[i]:
del bond_table[i]
#exit if to many excluded atoms or print excluded atoms
if len(bond_table) == 0:
print("Warning! To many excluded atoms. Exit.\n")
sys.exit(1)
else:
print(' ')
print("Excluded atoms: " + str(args.exclude))
#exclude bonds outside larger than a specific distance d
if args.distance:
for i in range(len(bond_table)-1, -1, -1):
if float(ang_bond_val.match(bond_table[i][2]).group()) > abs(args.distance):
#add atom names of excluded atoms to list --> for angles
if args.atom_name != bond_table[i][0]: list_of_atoms_with_large_bonds.append(bond_table[i][0])
if args.atom_name != bond_table[i][1]: list_of_atoms_with_large_bonds.append(bond_table[i][1])
del bond_table[i]
#exit if to many excluded atoms or print excluded atoms
if len(bond_table) == 0:
print("Warning! To many excluded atoms. Exit.\n")
sys.exit(1)
#print only if the list was created
elif list_of_atoms_with_large_bonds:
print(' ')
print("Excluded atoms (distance larger than " + str(abs(args.distance)) + " Å): "
+ str(list_of_atoms_with_large_bonds))
#build an angle table atom2-atom1-atom3 angle symmetry symmetry
angle_table=block.find(['_geom_angle_atom_site_label_1',
'_geom_angle_atom_site_label_2',
'_geom_angle_atom_site_label_3',
'_geom_angle',
'_geom_angle_site_symmetry_1',
'_geom_angle_site_symmetry_3'])
#delete angles not containing the atom from input in the center X-M-X
for i in range(len(angle_table)-1, -1, -1):
if args.atom_name not in angle_table[i][1]:
del angle_table[i]
#exclude symmetry eq. atoms, e.g. for polymeric compounds X-M1-M1'
#comment in case of problems
for i in range(len(angle_table)-1, -1, -1):
if args.atom_name in angle_table[i][2] and "." not in angle_table[i][5]:
del angle_table[i]
#exclude angles with specified atoms
if args.exclude:
for atom in args.exclude:
for i in range(len(angle_table)-1, -1, -1):
if atom in angle_table[i]:
del angle_table[i]
#exclude angles with specified atoms from the list of excluded atoms outside a specific distance d
if args.distance:
for atom in list_of_atoms_with_large_bonds:
for i in range(len(angle_table)-1, -1, -1):
if atom in angle_table[i]:
del angle_table[i]
#print bonds
print(' ')
print(args.atom_name + " binds to:")
print('-----------------------------------------------------------------')
for row in bond_table:
print(row[0] + '-' + row[1] + ' ' + row[2] + ' Å ' + row[3])
#calculate coordination number from occurance of atom in list
cn = (list(block.find_loop('_geom_bond_atom_site_label_1')).count(args.atom_name) +
list(block.find_loop('_geom_bond_atom_site_label_2')).count(args.atom_name))
print(' ')
print('The predicted coordination number for ' + args.atom_name + ' is ' + str(cn) + '.')
#exit if cn is < 3
if cn < 3:
print("Warning! the predicted coordination number is < 3. Exit.\n")
sys.exit(1)
#print angles
list_of_angles = []
print(' ')
print(args.atom_name + " angles are:")
print('-----------------------------------------------------------------')
for row in angle_table:
list_of_angles.append(float(ang_bond_val.match(row[3]).group()))
print(row[0] + '-' + row[1] + '-' + row[2] + ' ' +
row[3] + "° " + row[4] + ' ' + row[5])
#if there is only one angle, exit
if len(list_of_angles) < 2:
print(' ')
print("Warning! Number of angles is too low. Exit.\n")
sys.exit(1)
list_of_angles.sort(reverse=True)
beta = list_of_angles[0]
alpha = list_of_angles[1]
print(' ')
print('The two largest angles are beta = ' +
str(beta) + '°' + ' and alpha = ' +
str(alpha) + '°.')
#calculate the cn from number of angles cn=4 number of angles = 6, cn=5 number of angles=10
cn = list(block.find_loop('_geom_angle_atom_site_label_2')).count(args.atom_name)
printmark4 = ('')
printmark5 = ('')
#print '<--' arrow for most probale tau_x assignement
if cn == 6:
printmark4=('<--')
if cn == 10:
printmark5=('<--')
#print warning if cn is not 4 or 5
elif cn != 6 and cn != 10:
print("Warning! The predicted coordination number seems not suitable\n" +
"for the calculation of tau_4 or tau_5.")
print("Calculated values are probably meaningless.")
#calculate and print tau_x values
print(' ')
print(args.atom_name + ' geometry indices:')
print('-----------------------------------------------------------------')
print(f'tau_4 = {calc_tau4(beta, alpha):.2f} {printmark4}')
print(f"tau_4' = {calc_tau4impr(beta, alpha):.2f} {printmark4}")
print(f'tau_5 = {calc_tau5(beta, alpha):.2f} {printmark5}')
#print a table of typical tau_x values
#values different from 0 or 1 and the corresponding geometries have been taken
#from an internet source - don't take it too seriously
print(' ')
print(f"Table of typical geometries and their corresponding tau_x values")
print(f"-----------------------------------------------------------------")
print(f"Coordination number 4:")
print(f"Tetrahedral : tau_4 = 1.00 tau_4' = 1.00")
print(f"Trigonal pyramidal : tau_4 = 0.85 tau_4' = 0.85")
print(f"Seesaw : tau_4 = 0.43 tau_4' = 0.24")
print(f"Square planar : tau_4 = 0.00 tau_4' = 0.00\n")
print(f"Coordination number 5:")
print(f"Trigonal bipyramidal : tau_5 = 1.00 ")
print(f"Square pyramidal : tau_5 = 0.00 \n")