1 /*
2 * (c) Schrodinger, Inc.
3 */
4
5 #include "os_std.h"
6
7 #include "CoordSet.h"
8 #include "ObjectMolecule.h"
9 #include "Selector.h"
10 #include "HydrogenAdder.h"
11 #include "Err.h"
12
13 #include <cassert>
14
15 /**
16 * Add coordinates for atom `atm`.
17 *
18 * @param atm Atom index (in ObjectMolecule::AtomInfo)
19 * @param v Atom coordinates (`float[3]`)
20 *
21 * @pre `atm` doesn't have coordinates yet in given coord set
22 */
23 static
AppendAtomVertex(CoordSet * cs,unsigned atm,const float * v)24 void AppendAtomVertex(CoordSet* cs, unsigned atm, const float* v)
25 {
26 assert(cs->atmToIdx(atm) == -1);
27
28 int idx = cs->NIndex++;
29 VLACheck(cs->Coord, float, idx * 3 + 2);
30 VLACheck(cs->IdxToAtm, int, idx);
31
32 cs->IdxToAtm[idx] = atm;
33
34 if (cs->Obj->DiscreteFlag) {
35 cs->Obj->DiscreteAtmToIdx[atm] = idx;
36 cs->Obj->DiscreteCSet[atm] = cs;
37 } else {
38 cs->AtmToIdx[atm] = idx;
39 }
40
41 copy3f(v, cs->coordPtr(idx));
42 }
43
44 /**
45 * If `atm` has a planar (sp2) configuration, then write the plane's normal
46 * vector to the `normal` out pointer and return true.
47 *
48 * @param atm Atom index (in ObjectMolecule::AtomInfo)
49 * @param[out] normal Normal vector (`float[3]` with length 1)
50 * @param h_fix If true, then ignore hydrogen neighbors
51 *
52 * @pre Neighbors up-to-date
53 */
54 static
get_planer_normal_cs(const ObjectMolecule * I,const CoordSet * cs,unsigned atm,float * normal,bool h_fix)55 bool get_planer_normal_cs(
56 const ObjectMolecule* I,
57 const CoordSet* cs, unsigned atm, float *normal,
58 bool h_fix)
59 {
60 int nOcc = 0;
61 float occ[3 * 3];
62
63 if (I->AtomInfo[atm].geom != cAtomInfoPlanar)
64 return false;
65
66 int idx = cs->atmToIdx(atm);
67 if (idx == -1)
68 return false;
69
70 const float* center_coord = cs->coordPtr(idx);
71
72 int neighbor_atm, tmp;
73 ITERNEIGHBORATOMS(I->Neighbor, atm, neighbor_atm, tmp) {
74 if (h_fix && I->AtomInfo[neighbor_atm].isHydrogen())
75 continue;
76
77 // get neighbor coordinate
78 int neighbor_idx = cs->atmToIdx(neighbor_atm);
79 if (neighbor_idx == -1)
80 continue;
81
82 const float* neighbor_coord = cs->coordPtr(neighbor_idx);
83
84 // points away from center
85 float* vvec = occ + 3 * nOcc;
86 subtract3f(neighbor_coord, center_coord, vvec);
87 normalize3f(vvec);
88
89 if (++nOcc == 3)
90 // more doesn't make sence for a planar system
91 break;
92 }
93
94 if (nOcc < 2)
95 return false;
96
97 cross_product3f(occ, occ + 3, normal);
98
99 // avg of all three cross products
100 if (nOcc > 2) {
101 float v2[3];
102 for (int offset = 0; offset < 6; offset += 3) {
103 cross_product3f(occ + offset, occ + 6, v2);
104 if (dot_product3f(normal, v2) < 0) {
105 scale3f(v2, -1.f, v2);
106 }
107 add3f(normal, v2, normal);
108 }
109 }
110
111 normalize3f(normal);
112
113 return true;
114 }
115
116 /**
117 * Calculate plausible coordinates for those neighbors of `atm` which don't
118 * have coordinates yet.
119 *
120 * @param h_fix also reposition hydrogens with existing coordinates.
121 *
122 * @return Number of added/updated coordinates.
123 *
124 * @pre Neighbors up-to-date
125 *
126 * @note Similar to ::ObjectMoleculeFindOpenValenceVector ("Evolutionary
127 * descendant", code duplication event)
128 */
ObjectMoleculeSetMissingNeighborCoords(ObjectMolecule * I,CoordSet * cs,unsigned atm,bool h_fix)129 int ObjectMoleculeSetMissingNeighborCoords(
130 ObjectMolecule* I, CoordSet* cs, unsigned atm, bool h_fix)
131 {
132 auto G = I->G;
133 int n_present = 0;
134 float cbuf[4 * 3];
135 int present_atm = -1;
136 int missing_atm[4];
137 int n_missing = 0;
138
139 const AtomInfoType* ai = I->AtomInfo + atm;
140
141 int idx = cs->atmToIdx(atm);
142 if (idx == -1)
143 return 0;
144
145 const float* center_coord = cs->coordPtr(idx);
146
147 int neighbor_atm, tmp;
148 ITERNEIGHBORATOMS(I->Neighbor, atm, neighbor_atm, tmp) {
149 if (n_present == 4)
150 break;
151
152 // get neighbor coordinate
153 int neighbor_idx = cs->atmToIdx(neighbor_atm);
154 if (neighbor_idx == -1 ||
155 (h_fix && I->AtomInfo[neighbor_atm].isHydrogen())) {
156 missing_atm[n_missing++] = neighbor_atm;
157 continue;
158 }
159
160 const float* neighbor_coord = cs->coordPtr(neighbor_idx);
161
162 // points away from center
163 float* vvec = cbuf + 3 * n_present;
164 subtract3f(neighbor_coord, center_coord, vvec);
165 normalize3f(vvec);
166
167 present_atm = neighbor_atm;
168 ++n_present;
169 }
170
171 if (n_missing == 0)
172 // nothing to do
173 return 0;
174
175 int n_system = n_present;
176 if (n_system == 0) {
177 get_random3f(cbuf);
178 ++n_system;
179 }
180
181 switch (ai->geom) {
182 float t[3], z[3];
183
184 // Tetrahedral system: 109.5 degree angles
185 case cAtomInfoTetrahedral:
186 switch (n_system) {
187 case 1:
188 get_system1f3f(cbuf, t, z);
189 scale3f(cbuf, -0.334F, t); // cos(-109.5)
190 scale3f(z, 0.943F, z); // sin( 109.5)
191 add3f(z, t, cbuf + 3);
192 normalize3f(cbuf + 3);
193 case 2:
194 add3f(cbuf, cbuf + 3, t);
195 normalize3f(t);
196 scale3f(t, -1.0F, t);
197 cross_product3f(cbuf, cbuf + 3, z);
198 normalize3f(z);
199 scale3f(z, 1.41F, z); // tan(109.5 / 2.0)
200 add3f(t, z, cbuf + 6);
201 normalize3f(cbuf + 6);
202 case 3:
203 add3f(cbuf, cbuf + 3, t);
204 add3f(cbuf + 6, t, t);
205 scale3f(t, -1.0F, cbuf + 9);
206 normalize3f(cbuf + 9);
207 }
208 n_system = 4;
209 break;
210
211 // Planar system: 120.0 degree angles
212 case cAtomInfoPlanar:
213 switch (n_system) {
214 case 1:
215 if (present_atm >= 0 &&
216 get_planer_normal_cs(I, cs, present_atm, t, h_fix)) {
217 get_system2f3f(cbuf, t, z);
218 } else {
219 get_system1f3f(cbuf, t, z);
220 }
221 scale3f(cbuf, -0.500F, t);
222 scale3f(z, 0.866F, z); // sin( 120.0)
223 add3f(z, t, cbuf + 3);
224 normalize3f(cbuf + 3);
225 case 2:
226 add3f(cbuf, cbuf + 3, t);
227 scale3f(t, -1.0F, cbuf + 6);
228 normalize3f(cbuf + 6);
229 }
230 n_system = 3;
231 break;
232
233 // Linear system: 180.0 degree angles
234 case cAtomInfoLinear:
235 switch (n_system) {
236 case 1:
237 scale3f(cbuf, -1.0F, cbuf + 3);
238 normalize3f(cbuf + 3);
239 }
240 n_system = 2;
241 break;
242 }
243
244 if (n_missing > n_system - n_present) {
245 n_missing = n_system - n_present;
246 }
247
248 // adding coordinates will invalidate pointer
249 float center_coord_copy[3];
250 copy3f(center_coord, center_coord_copy);
251 center_coord = nullptr;
252
253 for (int i = 0; i < n_missing; ++i) {
254 float bondlength = AtomInfoGetBondLength(G,
255 I->AtomInfo + atm,
256 I->AtomInfo + missing_atm[i]);
257 float* coord = cbuf + (n_present + i) * 3;
258 scale3f(coord, bondlength, coord);
259 add3f(coord, center_coord_copy, coord);
260
261 if (h_fix && (idx = cs->atmToIdx(missing_atm[i])) != -1) {
262 copy3f(coord, cs->coordPtr(idx));
263 } else {
264 AppendAtomVertex(cs, missing_atm[i], coord);
265 }
266 }
267
268 return n_missing;
269 }
270
271 /**
272 * Add hydrogens to selection
273 *
274 * @param sele Valid atom selection
275 * @param state Object state (can be all (-1) or current (-2))
276 *
277 * @return False if `I` has no atoms in the selection or if the chemistry (atom
278 * geometry and valence) can't be determined.
279 */
ObjectMoleculeAddSeleHydrogensRefactored(ObjectMolecule * I,int sele,int state)280 int ObjectMoleculeAddSeleHydrogensRefactored(ObjectMolecule* I, int sele, int state)
281 {
282 auto G = I->G;
283 auto const n_atom_old = I->NAtom;
284
285 bool seleFlag = false;
286 for (unsigned atm = 0; atm < n_atom_old; atm++) {
287 const auto ai = I->AtomInfo + atm;
288 if (SelectorIsMember(G, ai->selEntry, sele)) {
289 seleFlag = true;
290 break;
291 }
292 }
293
294 if (!seleFlag) {
295 return true;
296 }
297
298 if (!ObjectMoleculeVerifyChemistry(I, state)) {
299 ErrMessage(G, " AddHydrogens", "missing chemical geometry information.");
300 return false;
301 }
302
303 ObjectMoleculeUpdateNeighbors(I);
304
305 // add hydrogens (without coordinates)
306 for (unsigned atm = 0; atm < n_atom_old; ++atm) {
307 const auto ai = I->AtomInfo + atm;
308
309 if (ai->isMetal())
310 continue;
311
312 if (!SelectorIsMember(G, ai->selEntry, sele))
313 continue;
314
315 int nneighbors = I->Neighbor[I->Neighbor[atm]];
316 int nimplicit = ai->valence - nneighbors;
317
318 if (nimplicit <= 0)
319 continue;
320
321 I->AtomInfo.reserve(I->NAtom + nimplicit);
322 I->Bond.reserve(I->NBond + nimplicit);
323
324 for (int i = 0; i < nimplicit; ++i) {
325 // bond
326 auto bond = I->Bond + I->NBond++;
327 BondTypeInit2(bond, atm, I->NAtom, 1);
328
329 // atom
330 auto atom = I->AtomInfo + I->NAtom++;
331 atom->protons = cAN_H;
332 atom->geom = cAtomInfoSingle;
333 atom->valence = 1;
334 ObjectMoleculePrepareAtom(I, atm, atom, /* uniquefy */ false);
335 }
336 }
337
338 // grow index arrays
339 for (StateIterator iter(G, nullptr, cSelectorUpdateTableAllStates, I->NCSet);
340 iter.next();) {
341 CoordSet* cs = I->CSet[iter.state];
342 if (cs)
343 cs->extendIndices(I->NAtom);
344 }
345
346 I->invalidate(cRepAll, cRepInvBonds, state);
347 ObjectMoleculeUpdateNeighbors(I);
348
349 AtomInfoUniquefyNames(G,
350 I->AtomInfo, n_atom_old,
351 I->AtomInfo + n_atom_old, nullptr,
352 I->NAtom - n_atom_old);
353
354 // fill coordinates
355 for (StateIterator iter(I, state); iter.next();) {
356 CoordSet* cs = I->CSet[iter.state];
357 if (!cs)
358 continue;
359
360 for (unsigned idx = 0; idx < cs->NIndex; ++idx) {
361 auto atm = cs->IdxToAtm[idx];
362 if (atm >= n_atom_old)
363 continue;
364
365 const auto ai = I->AtomInfo + atm;
366 if (!SelectorIsMember(G, ai->selEntry, sele))
367 continue;
368
369 ObjectMoleculeSetMissingNeighborCoords(I, cs, atm);
370 }
371 }
372
373 I->invalidate(cRepAll, cRepInvAtoms, state);
374 ObjectMoleculeSort(I);
375 ObjectMoleculeUpdateIDNumbers(I);
376
377 return true;
378 }
379