Skip to content

Commit ca7dc55

Browse files
larsbratholmlarsbratholm
committed
small changes
1 parent a082c30 commit ca7dc55

File tree

4 files changed

+53
-69
lines changed

4 files changed

+53
-69
lines changed

‎atomorder/definitions.py

+3-3
Original file line numberDiff line numberDiff line change
@@ -179,12 +179,12 @@ def update(self, args):
179179

180180
self.atomic_sybyl_weight = args.atomic_sybyl_weight
181181
self.bond_weight = args.bond_weight
182-
self.hydrogen_rotation_weight = 0.1
182+
self.hydrogen_rotation_weight = 1.0
183183
self.annealing_method = "multiplication" # multiplication/addition
184184
#self.annealing_method = "addition"
185185
self.initial_inverse_temperature = 1e-1
186-
self.final_inverse_temperature = 1e1
187-
self.max_annealing_iterations = 100
186+
self.final_inverse_temperature = 1e2
187+
self.max_annealing_iterations = 200
188188
self.max_relaxation_iterations = 100
189189
self.max_softassign_iterations = 1000
190190
self.annealing_convergence_threshold = 1e-2

‎atomorder/objectives.py

+9-4
Original file line numberDiff line numberDiff line change
@@ -129,9 +129,10 @@ def analytical_solver(self, m):
129129
squared_distances = np.zeros(self.M.match_matrix.shape, dtype=float)
130130

131131
match[np.ix_(self.reactant_hydrogen_mask, self.product_hydrogen_mask)] *= settings.hydrogen_rotation_weight
132+
all_X = []
132133
for i, reactant_indices in enumerate(self.M.reactant_subset_indices):
134+
X = self.X[reactant_indices]
133135
for j, product_indices in enumerate(self.M.product_subset_indices):
134-
X = self.X[reactant_indices]
135136
Y = self.Y[product_indices]
136137
sub_matrix_indices = np.ix_(reactant_indices, product_indices)
137138
match_sub_matrix = match[sub_matrix_indices]
@@ -152,8 +153,11 @@ def analytical_solver(self, m):
152153
rot, trans = self.transform(r,s)
153154

154155
squared_distances[sub_matrix_indices] = np.sum((Y[None,:,:] - trans[None,None,:] - rot.dot(X.T).T[:,None,:])**2, axis=2)
156+
all_X.append(trans + rot.dot(X.T).T)
155157

156158
squared_distances[np.ix_(self.reactant_hydrogen_mask, self.product_hydrogen_mask)] *= settings.hydrogen_rotation_weight
159+
write_mol2(np.concatenate(all_X), self.M.reaction.reactants.atoms, "reactant%d.mol2" % self.M.it,self.M.reaction.num_reactant_atoms)
160+
write_mol2(self.Y+np.asarray([10,0,0]), self.M.reaction.products.atoms, "product%d.mol2" % self.M.it, self.M.reaction.num_product_atoms, match)
157161

158162
return squared_distances
159163

@@ -390,10 +394,10 @@ def objective(flat_q, m, self):
390394
squared_distances[reactant_indices,:] = np.sum((Y[None,:,:] - xtrans[None,None,:] - xrot.dot(X.T).T[:,None,:])**2, axis=2)
391395
for i, v in enumerate(all_X):
392396
all_X[i] += np.asarray([1*i,0,0])
393-
write_xyz(np.concatenate(all_X), self.M.reactants_elements, "reactant%d.xyz" % self.M.it, str(np.sum(match*squared_distances)))
394-
write_xyz(Y, self.M.products_elements, "product%d.xyz" % self.M.it, str(np.sum(match*squared_distances)))
397+
#write_xyz(np.concatenate(all_X), self.M.reactants_elements, "reactant%d.xyz" % self.M.it, str(np.sum(match*squared_distances)))
398+
#write_xyz(Y, self.M.products_elements, "product%d.xyz" % self.M.it, str(np.sum(match*squared_distances)))
395399
write_mol2(np.concatenate(all_X), self.M.reaction.reactants.atoms, "reactant%d.mol2" % self.M.it,self.M.reaction.num_reactant_atoms)
396-
write_mol2(Y, self.M.reaction.products.atoms, "product%d.mol2" % self.M.it, self.M.reaction.num_product_atoms)
400+
write_mol2(Y+np.asarray([10,0,0]), self.M.reaction.products.atoms, "product%d.mol2" % self.M.it, self.M.reaction.num_product_atoms, match)
397401

398402
squared_distances[np.ix_(self.reactant_hydrogen_mask, self.product_hydrogen_mask)] *= settings.hydrogen_rotation_weight
399403
E += np.sum(m*squared_distances)
@@ -418,6 +422,7 @@ def objective(flat_q, m, self):
418422
bounds.extend([(None, None)]*3)
419423
bounds.extend([(-1, 1)]*3)
420424
q = self.q.copy().ravel()
425+
#q = np.zeros(q.shape)
421426
opt = scipy.optimize.minimize(objective, q, method="l-bfgs-b", options={"maxiter": 500, "disp": 0, "ftol": 1e-6}, args=(match, self), bounds=bounds)
422427
#opt = scipy.optimize.minimize(objective, q, method="nelder-mead", options={"maxiter": 500, "disp": 0, "ftol": 1e-6}, args=(match, self))
423428
#assert(np.allclose(self.squared_distances, self.analytical_solver(match)))

‎atomorder/ordering.py

+19-3
Original file line numberDiff line numberDiff line change
@@ -309,7 +309,6 @@ def annealing(self):
309309
return True
310310

311311
self.update_inverse_temperature()
312-
print self.match_matrix.diagonal().sum()
313312

314313
# The annealing terminated without a definitive assignment
315314
eprint(3, "WARNING: Annealing algorithm did not converge to %.2g (reached %.2g) in %d iterations and row dominance was %s" % (self.annealing_convergence_threshold, row_l2_deviation, self.max_annealing_iterations, str(self.row_dominance())))
@@ -392,12 +391,29 @@ def finalize(self):
392391
oprint(4, "%d out of %d matched in diagonal" % (sum(self.match_matrix.diagonal() > 0.5), N))
393392
rd = self.row_dominance(True)
394393
oprint(4, "Order: " + str(rd))
395-
M1 = np.random.random((N,N))*1e-9 + np.eye(N, dtype=float)
396-
M1 /= M1.sum(0)
394+
#M1 = np.random.random((N,N))*1e-9 + np.eye(N, dtype=float)
395+
#M1 /= M1.sum(0)
396+
M1 = np.eye(N, dtype=float)
397+
#M1[17,17] = 0
398+
#M1[17,30] = 1
399+
#M1[14,14] = 0
400+
#M1[14,16] = 1
401+
#M1[16,16] = 0
402+
#M1[16,14] = 1
403+
#M1[30,30] = 0
404+
#M1[30,29] = 1
405+
#M1[19,19] = 0
406+
#M1[19,17] = 1
407+
#M1[18,18] = 0
408+
#M1[18,19] = 1
409+
#M1[29,29] = 0
410+
#M1[29,18] = 1
397411
M2 = np.zeros((N,N), dtype=float)
398412
for i in range(N):
399413
M2[i,rd[i]] = 1
400414

415+
print np.sum(abs(M1-M2))
416+
401417
scores1 = [fun.score(M1) for fun in self.obj]
402418
scores2 = [fun.score(M2) for fun in self.obj]
403419
#for i in range(N):

‎atomorder/parsers.py

+22-59
Original file line numberDiff line numberDiff line change
@@ -159,80 +159,43 @@ def write_xyz(coordinates, elements, filename, title = ""):
159159
line = "{0:2s} {1:15.8f} {2:15.8f} {3:15.8f}\n".format(elements[i], *coordinates[i])
160160
f.write(line)
161161

162-
def write_mol2(coordinates, atoms, filename, split):
162+
def write_mol2(coordinates, atoms, filename, split, match = None):
163+
163164
N = coordinates.shape[0]
165+
order = np.arange(N, dtype=int)
166+
rev_order = np.arange(N, dtype=int)
167+
if match != None:
168+
largest_row_indices = match.argmax(1)
169+
unique_indices, counts = np.unique(largest_row_indices, return_counts=True)
170+
if (counts == 1).all() == False:
171+
pass
172+
else:
173+
order = largest_row_indices
174+
rev_order = np.asarray([np.where(order == i)[0][0] for i in range(N)])
175+
#assert(np.allclose(rev_order, match.argmax(0)))
176+
164177
origin = []
165178
for i, n in enumerate(split):
166179
origin += [i+1]*n
167180

168181
lines = []
169-
for i,atom in enumerate(atoms):
170-
line = "\t{:d} {:s}\t{:.4f}\t{:.4f}\t{:4f} {:s}\t {:d} LIG{:d}\t0.0\n".format(i+1, atom.element_symbol, coordinates[i,0], coordinates[i,1], coordinates[i,2], atom.sybyl, origin[i], origin[i])
182+
for i in range(N):
183+
j = order[i]
184+
line = "\t{:d} {:s}\t{:.4f}\t{:.4f}\t{:4f} {:s}\t {:d} LIG{:d}\t0.0\n".format(i+1, atoms[j].element_symbol, coordinates[j,0], coordinates[j,1], coordinates[j,2], atoms[j].sybyl, origin[j], origin[j])
171185
lines.append(line)
172186
lines.append("@<TRIPOS>BOND\n")
173187
count = 0
174-
for atom1 in atoms:
188+
for i in range(N):
189+
j = order[i]
190+
atom1 = atoms[j]
175191
for atom2 in atom1.bonded_atoms:
176-
if atom1.mixture_index >= atom2.mixture_index:
192+
if i >= rev_order[atom2.mixture_index]:
177193
continue
178194
count += 1
179195
# TODO print sybyl bond
180-
line = "\t{:d}\t{:d}\t{:d}\t1\n".format(count, atom1.mixture_index+1, atom2.mixture_index+1)
196+
line = "\t{:d}\t{:d}\t{:d}\t1\n".format(count, i+1, rev_order[atom2.mixture_index]+1)
181197
lines.append(line)
182198
with open(filename, "w") as f:
183199
f.write("@<TRIPOS>MOLECULE\n\n%d %d 0 0 0\nSMALL\nGASTEIGER\n\n@<TRIPOS>ATOM\n" % (N,count))
184200
for line in lines:
185201
f.write(line)
186-
# 23 23 0 0 0
187-
#SMALL
188-
#GASTEIGER
189-
#
190-
#@<TRIPOS>ATOM
191-
# 1 C -1.7175 1.4267 -0.0723 C.ar 1 LIG1 0.0594
192-
# 2 C -2.9787 0.8162 -0.1439 C.ar 1 LIG1 -0.0477
193-
# 3 C -3.1199 -0.5684 -0.1047 C.ar 1 LIG1 -0.0609
194-
# 4 C -1.9971 -1.3773 0.0127 C.ar 1 LIG1 -0.0582
195-
# 5 C -0.7362 -0.7967 0.0861 C.ar 1 LIG1 -0.0185
196-
# 6 C -0.5838 0.5925 0.0336 C.ar 1 LIG1 0.1420
197-
# 7 C -1.6765 2.9418 -0.1207 C.2 1 LIG1 0.1636
198-
# 8 C -0.3674 3.7031 -0.0636 C.3 1 LIG1 -0.0013
199-
# 9 O -2.7031 3.5888 -0.2123 O.2 1 LIG1 -0.2922
200-
# 10 H -3.8773 1.4220 -0.2316 H 1 LIG1 0.0626
201-
# 11 H -4.1094 -1.0157 -0.1620 H 1 LIG1 0.0618
202-
# 12 H -2.1033 -2.4587 0.0498 H 1 LIG1 0.0619
203-
# 13 H 0.1370 -1.4366 0.1875 H 1 LIG1 0.0654
204-
# 14 O 0.6987 1.0989 0.1355 O.3 1 LIG1 -0.4253
205-
# 15 H 0.2699 3.4283 -0.9065 H 1 LIG1 0.0309
206-
# 16 H 0.1580 3.4909 0.8698 H 1 LIG1 0.0309
207-
# 17 H -0.5477 4.7796 -0.1129 H 1 LIG1 0.0309
208-
# 18 C 1.7296 0.7539 -0.6667 C.2 1 LIG1 0.3092
209-
# 19 O 1.6365 0.0338 -1.6412 O.2 1 LIG1 -0.2507
210-
# 20 C 3.0060 1.4639 -0.2836 C.3 1 LIG1 0.0336
211-
# 21 H 2.8360 2.5412 -0.2266 H 1 LIG1 0.0342
212-
# 22 H 3.7897 1.2702 -1.0194 H 1 LIG1 0.0342
213-
# 23 H 3.3467 1.1170 0.6943 H 1 LIG1 0.0342
214-
#@<TRIPOS>BOND
215-
# 1 19 18 2
216-
# 2 22 20 1
217-
# 3 15 8 1
218-
# 4 18 20 1
219-
# 5 18 14 1
220-
# 6 20 21 1
221-
# 7 20 23 1
222-
# 8 10 2 1
223-
# 9 9 7 2
224-
# 10 11 3 1
225-
# 11 2 3 ar
226-
# 12 2 1 ar
227-
# 13 7 1 1
228-
# 14 7 8 1
229-
# 15 17 8 1
230-
# 16 3 4 ar
231-
# 17 1 6 ar
232-
# 18 8 16 1
233-
# 19 4 12 1
234-
# 20 4 5 ar
235-
# 21 6 5 ar
236-
# 22 6 14 1
237-
# 23 5 13 1
238-

0 commit comments

Comments
 (0)