1
2 /*
3 A* -------------------------------------------------------------------
4 B* This file contains source code for the PyMOL computer program
5 C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific.
6 D* -------------------------------------------------------------------
7 E* It is unlawful to modify or remove this copyright notice.
8 F* -------------------------------------------------------------------
9 G* Please see the accompanying LICENSE file for further information.
10 H* -------------------------------------------------------------------
11 I* Additional authors of this source file include:
12 -*
13 -*
14 -*
15 Z* -------------------------------------------------------------------
16 */
17 #include"os_python.h"
18 #include"os_predef.h"
19 #include"os_std.h"
20 #include"os_gl.h"
21 #include"OOMac.h"
22 #include"Feedback.h"
23 #include"Util.h"
24 #include"Sculpt.h"
25 #include"SculptCache.h"
26 #include"Scene.h"
27 #include"Vector.h"
28 #include"Word.h"
29 #include"Editor.h"
30 #include "Lex.h"
31 #include "ObjectMolecule.h"
32 #include "CoordSet.h"
33
34 #include"CGO.h"
35
36 #ifndef R_SMALL8
37 #define R_SMALL8 0.00000001
38 #endif
39
40 #define NB_HASH_SIZE 262144
41 #define EX_HASH_SIZE 65536
42
43 #define nb_hash(v) \
44 (((((int)*(v ))>> 2)&0x0003F)|\
45 ((((int)*(v+1))<< 4)&0x00FC0)|\
46 ((((int)*(v+2))<<10)&0x3F000))
47
48 #define nb_hash_off_i0(v0i,d) \
49 ((((d)+v0i)>> 2)&0x0003F)
50
51 #define nb_hash_off_i1(v1i,e) \
52 ((((e)+v1i)<< 4)&0x00FC0)
53
54 #define nb_hash_off_i2(v2i,f) \
55 ((((f)+v2i)<<10)&0x3F000)
56
57 #define nb_hash_off(v,d,e,f) \
58 (((((d)+(int)*(v ))>> 2)&0x0003F)|\
59 ((((e)+(int)*(v+1))<< 4)&0x00FC0)|\
60 ((((f)+(int)*(v+2))<<10)&0x3F000))
61
62
63 /* below are empirically optimized */
64
65 #define ex_hash_i0(a) \
66 (((a)^((a)>>5))&0x00FF)
67
68 #define ex_hash_i1(b) \
69 (( ((b)<<5))&0xFF00)
70
71 #define ex_hash(a,b) \
72 (((((a)^((a)>>5)))&0x00FF)|\
73 ((( ((b)<<5)))&0xFF00))
74
ShakerDoDist(float target,float * v0,float * v1,float * d0to1,float * d1to0,float wt)75 static float ShakerDoDist(float target, float *v0, float *v1, float *d0to1, float *d1to0,
76 float wt)
77 {
78 float d[3], push[3];
79 float len, dev, dev_2, sc, result;
80
81 subtract3f(v0, v1, d);
82 len = (float) length3f(d);
83 dev = target - len;
84 if((result = (float) fabs(dev)) > R_SMALL8) {
85 dev_2 = wt * dev / 2.0F;
86 if(len > R_SMALL8) { /* nonoverlapping */
87 sc = dev_2 / len;
88 scale3f(d, sc, push);
89 add3f(push, d0to1, d0to1);
90 subtract3f(d1to0, push, d1to0);
91 } else { /* overlapping, so just push along X */
92 float rd[3];
93 get_random3f(rd);
94 d0to1[0] -= rd[0] * dev_2;
95 d1to0[0] += rd[0] * dev_2;
96 d0to1[1] -= rd[1] * dev_2;
97 d1to0[1] += rd[1] * dev_2;
98 d0to1[2] -= rd[2] * dev_2;
99 d1to0[2] += rd[2] * dev_2;
100 }
101 } else
102 result = 0.0;
103 return result;
104 }
105
ShakerDoTors(int type,float * v0,float * v1,float * v2,float * v3,float * p0,float * p1,float * p2,float * p3,float tole,float wt)106 static float ShakerDoTors(int type, float *v0, float *v1, float *v2, float *v3,
107 float *p0, float *p1, float *p2, float *p3, float tole,
108 float wt)
109 {
110
111 float push0[3], push3[3];
112 float axis[3], seg0[3], seg1[3], perp0[3], perp1[3];
113 float dir[3];
114 float sc;
115 float sign, dp;
116 float result = 0.0F;
117
118 /* v0 v3
119 \ /
120 v1__v2 */
121
122 subtract3f(v2, v1, axis);
123 subtract3f(v0, v1, seg0);
124 subtract3f(v3, v2, seg1);
125 cross_product3f(seg0, axis, perp0);
126 cross_product3f(axis, seg1, perp1);
127
128 normalize3f(perp0);
129 normalize3f(perp1);
130
131 dp = dot_product3f(perp0, perp1);
132
133 switch (type) {
134 case cShakerTorsSP3SP3:
135 if(dp < -0.5F) {
136 result = ((float) fabs(dp)) - 0.5F;
137 if(result < tole) /* discontinuous low bottom well */
138 result = result / 5.0F;
139 } else if(dp < 0.5) {
140 result = -0.5F - dp;
141 } else {
142 result = 1.0F - dp;
143 }
144 break;
145 case cShakerTorsFlat:
146 if(fabs(dp) < 0.5F) /* don't attempt to resolve when ambiguous */
147 return 0.0F;
148 if(dp > 0.0F) {
149 result = 1.0F - dp;
150 } else {
151 result = -1.0F - dp;
152 }
153 result *= 5.0F; /* emphasize */
154 break;
155 case cShakerTorsAmide:
156 if(dp > -0.7F) { /* highly biased in favor of the input state */
157 result = 1.0F - dp;
158 } else {
159 result = -1.0F - dp;
160 }
161 result *= 50.0F; /* emphasize */
162 break;
163 case cShakerTorsDisulfide:
164 if(fabs(dp) < tole)
165 return 0.0F;
166 result = -dp;
167 if(result < tole)
168 result = result / 25.F;
169 break;
170 }
171
172 cross_product3f(perp0, perp1, dir);
173 sign = dot_product3f(axis, dir);
174
175 if(sign < 0.0F)
176 sc = wt * result;
177 else
178 sc = -wt * result;
179
180 scale3f(perp0, sc, push0);
181 scale3f(perp1, sc, push3);
182
183 add3f(p0, push0, p0);
184 add3f(p3, push3, p3);
185 subtract3f(p1, push0, p1);
186 subtract3f(p2, push3, p2);
187
188 return result;
189
190 }
191
ShakerDoDistLimit(float target,float * v0,float * v1,float * d0to1,float * d1to0,float wt)192 static float ShakerDoDistLimit(float target, float *v0, float *v1, float *d0to1,
193 float *d1to0, float wt)
194 {
195 float d[3], push[3];
196 float len, dev, dev_2, sc;
197
198 subtract3f(v0, v1, d);
199 len = (float) length3f(d);
200 dev = len - target;
201 if(dev > 0.0F) { /* assuming len is non-zero since it is above target */
202 dev_2 = wt * dev * (-0.5F);
203 sc = dev_2 / len;
204 scale3f(d, sc, push);
205 add3f(push, d0to1, d0to1);
206 subtract3f(d1to0, push, d1to0);
207 return dev;
208 } else {
209 return 0.0F;
210 }
211 }
212
ShakerDoDistMinim(float target,float * v0,float * v1,float * d0to1,float * d1to0,float wt)213 static float ShakerDoDistMinim(float target, float *v0, float *v1, float *d0to1,
214 float *d1to0, float wt)
215 {
216 float d[3], push[3];
217 float len, dev, dev_2, sc;
218
219 subtract3f(v0, v1, d);
220 len = (float) length3f(d);
221 dev = len - target;
222
223 if(dev < 0.0F) { /* assuming len is non-zero since it is above target */
224 dev_2 = -wt * dev * 0.5F;
225 sc = dev_2 / len;
226 scale3f(d, sc, push);
227 add3f(push, d0to1, d0to1);
228 subtract3f(d1to0, push, d1to0);
229 return -dev;
230 } else {
231 return 0.0F;
232 }
233 }
234
CSculpt(PyMOLGlobals * G)235 CSculpt::CSculpt (PyMOLGlobals * G)
236 {
237 this->G = G;
238 this->Shaker = pymol::make_unique<CShaker>(G);
239 this->NBList = pymol::vla<int>(150000);
240 this->NBHash = std::vector<int>(NB_HASH_SIZE);
241 this->EXList = pymol::vla<int>(100000);
242 this->EXHash = std::vector<int>(EX_HASH_SIZE);
243 this->Don = pymol::vla<int>(1000);
244 this->Acc = pymol::vla<int>(1000);
245 {
246 int a;
247 for(a = 1; a < 256; a++)
248 this->inverse[a] = 1.0F / a;
249 }
250 }
251
252 typedef struct {
253 const int *neighbor;
254 AtomInfoType *ai;
255 const int *atm2idx1, *atm2idx2;
256 } CountCall;
257
258 typedef struct {
259 PyMOLGlobals *G;
260 CShaker *Shaker;
261 AtomInfoType *ai;
262 const int *atm2idx;
263 const CoordSet *cSet;
264 const CoordSet *const *discCSet;
265 const float *coord;
266 const int *neighbor;
267 int atom0;
268 int min, max, mode;
269 } ATLCall;
270
count_branch(CountCall * CNT,int atom,int limit)271 static int count_branch(CountCall * CNT, int atom, int limit)
272 {
273 AtomInfoType *ai = CNT->ai + atom;
274 int count = 0;
275
276 if(!ai->temp1) {
277 count = (ai->isHydrogen() ? 0 : 1);
278 if(count) {
279 if((CNT->atm2idx1[atom] < 0) || (CNT->atm2idx2[atom] < 0))
280 count = 0;
281 }
282 if(count && (limit > 0)) {
283 int n0 = CNT->neighbor[atom] + 1;
284 int b1;
285 ai->temp1 = true;
286 while((b1 = CNT->neighbor[n0]) >= 0) {
287 count += count_branch(CNT, b1, limit - 1);
288 n0 += 2;
289 }
290 ai->temp1 = false;
291 }
292 }
293 return count;
294 }
295
add_triangle_limits(ATLCall * ATL,int prev,int cur,float dist,int count)296 static void add_triangle_limits(ATLCall * ATL, int prev, int cur, float dist, int count)
297 {
298 ATLCall *I = ATL;
299 int n0;
300 int n1;
301 float dist_limit;
302 int atom1;
303
304 n0 = I->neighbor[cur];
305 if((count >= I->min) && (count > 1)) {
306 int add_flag = false;
307 switch (I->mode) {
308 case 1:
309 add_flag = 1; /* all */
310 break;
311 case 2:
312 add_flag = (count && !(count & 1)); /* evens */
313 break;
314 case 3:
315 add_flag = ((count & (count - 1)) == 0); /* powers of two */
316 break;
317 case 0:
318 default:
319 add_flag = (!I->ai[I->atom0].isHydrogen()); /* all heavies */
320 break;
321 }
322 if(add_flag) {
323 n1 = n0 + 1;
324
325 /* first mark and register */
326 while((atom1 = I->neighbor[n1]) >= 0) {
327 if((!I->ai[atom1].temp1) && (I->atom0 < atom1)) {
328 int ref = prev;
329 if(count & 0x1) { /* odd */
330 ref = cur;
331 }
332 if(((!I->discCSet) ||
333 ((I->cSet == I->discCSet[ref]) && (I->cSet == I->discCSet[atom1]))) &&
334 ((I->mode != 0) || (!I->ai[atom1].isHydrogen()))) {
335 int ia = I->atm2idx[ref];
336 int ib = I->atm2idx[atom1];
337 if((ia >= 0) && (ib >= 0)) {
338 const float *va = I->coord + 3 * ia;
339 const float *vb = I->coord + 3 * ib;
340 dist_limit = dist + diff3f(va, vb);
341 ShakerAddDistCon(I->Shaker, I->atom0, atom1, dist_limit, cShakerDistLimit,
342 1.0F);
343 }
344 }
345 I->ai[atom1].temp1 = 1;
346 }
347 n1 += 2;
348 }
349 }
350 }
351
352 if(count <= I->max) {
353 /* then recurse */
354 n1 = n0 + 1;
355 while((atom1 = I->neighbor[n1]) >= 0) {
356 if(I->ai[atom1].temp1 < 2) {
357 dist_limit = dist;
358 if(!(count & 0x1)) { /* accumulate distances between even atoms only */
359 if((!I->discCSet)
360 || ((I->cSet == I->discCSet[prev]) && (I->cSet == I->discCSet[atom1]))) {
361 int ia = I->atm2idx[prev];
362 int ib = I->atm2idx[atom1];
363 if((ia >= 0) && (ib >= 0)) {
364 const float *va = I->coord + 3 * ia;
365 const float *vb = I->coord + 3 * ib;
366 dist_limit += diff3f(va, vb);
367 }
368 }
369 I->ai[atom1].temp1 = 2;
370 }
371 I->ai[atom1].temp1 = 2;
372 add_triangle_limits(I, cur, atom1, dist_limit, count + 1);
373 }
374 n1 += 2;
375 }
376 }
377 }
378
SculptMeasureObject(CSculpt * I,ObjectMolecule * obj,int state,int match_state,int match_by_segment)379 void SculptMeasureObject(CSculpt * I, ObjectMolecule * obj, int state, int match_state,
380 int match_by_segment)
381 {
382 PyMOLGlobals *G = I->G;
383 int a, a0, a1, a2, a3, b0, b1, b2, b3, b4;
384 const BondType *b;
385 float *v0, *v1, *v2, *v3, d, dummy;
386 CoordSet *cs;
387 int n0, n1, n2, n3;
388 int *planar = NULL;
389 int *linear = NULL;
390 int *single = NULL;
391 int *crdidx = NULL;
392 int nex = 1;
393 int *j, *k, xhash;
394 int ex_type;
395 AtomInfoType *ai, *ai1, *ai2, *obj_atomInfo;
396 int xoffset;
397 int use_cache = 1;
398
399 PRINTFD(G, FB_Sculpt)
400 " SculptMeasureObject-Debug: entered.\n" ENDFD;
401
402 if(match_state < 0)
403 match_state = state;
404 if(state < 0)
405 state = obj->getCurrentState();
406
407 ShakerReset(I->Shaker.get());
408
409 UtilZeroMem(I->NBHash.data(), NB_HASH_SIZE * sizeof(int));
410 UtilZeroMem(I->EXHash.data(), EX_HASH_SIZE * sizeof(int));
411
412 if((state >= 0) && (state < obj->NCSet) && (obj->CSet[state])) {
413 obj_atomInfo = obj->AtomInfo.data();
414
415 VLACheck(I->Don, int, obj->NAtom);
416 VLACheck(I->Acc, int, obj->NAtom);
417 ai = obj_atomInfo;
418 for(a = 0; a < obj->NAtom; a++) {
419 I->Don[a] = false;
420 I->Acc[a] = false;
421 AtomInfoCheckUniqueID(G, ai);
422 ai++;
423 }
424
425 ObjectMoleculeVerifyChemistry(obj, state);
426 ObjectMoleculeUpdateNeighbors(obj);
427
428 cs = obj->CSet[state];
429
430 use_cache = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_sculpt_memory);
431 if(obj->NBond) {
432 const int *neighbor = obj->Neighbor;
433 int n_atom = obj->NAtom;
434
435 planar = pymol::malloc<int>(n_atom);
436 linear = pymol::malloc<int>(n_atom);
437 single = pymol::malloc<int>(n_atom);
438 crdidx = pymol::malloc<int>(n_atom);
439 ai = obj_atomInfo;
440
441 for(a = 0; a < n_atom; a++) {
442 planar[a] = (ai->geom == cAtomInfoPlanar);
443 linear[a] = (ai->geom == cAtomInfoLinear);
444 single[a] = (ai->geom == cAtomInfoSingle);
445
446 a0 = cs->atmToIdx(a);
447 crdidx[a] = a0;
448
449 ai++;
450 }
451
452 /* brain-dead donor/acceptor assignment
453 * REPLACE later on with pattern-based system */
454
455 /* pass 1 */
456
457 b = obj->Bond;
458 for(a = 0; a < obj->NBond; a++) {
459 b1 = b->index[0];
460 b2 = b->index[1];
461
462 ai1 = obj_atomInfo + b1;
463 ai2 = obj_atomInfo + b2;
464
465 /* make blanket assumption that all nitrogens with
466 <3 bonds are donors -- we qualify this below... */
467
468 if(ai1->protons == cAN_N) {
469 n1 = neighbor[b1];
470 if(neighbor[n1] < 3) { /* N with L.P. */
471 I->Don[b1] = true;
472 }
473 }
474
475 if(ai2->protons == cAN_N) {
476 n2 = neighbor[b2];
477 if(neighbor[n2] < 3) { /* N with L.P. */
478 I->Don[b2] = true;
479 }
480 }
481
482 /* assume O is always an acceptor... */
483
484 if(ai1->protons == cAN_O)
485 I->Acc[b1] = true;
486 if(ai2->protons == cAN_O)
487 I->Acc[b2] = true;
488 b++;
489 }
490
491 /* pass 2 */
492 b = obj->Bond;
493 for(a = 0; a < obj->NBond; a++) {
494 b1 = b->index[0];
495 b2 = b->index[1];
496
497 /* nitrogens with lone pairs are acceptors
498 (not donors as assumed above) */
499
500 ai1 = obj_atomInfo + b1;
501 ai2 = obj_atomInfo + b2;
502
503 if(ai1->protons == cAN_N) {
504 if(b->order == 2) {
505 n1 = neighbor[b1];
506 if(neighbor[n1] < 3) { /* N with L.P. */
507 I->Acc[b1] = true;
508 I->Don[b1] = false;
509 }
510 }
511 }
512 if(ai2->protons == cAN_N) {
513 if(b->order == 2) {
514 n2 = neighbor[b2];
515 if(neighbor[n2] < 3) { /* N with L.P. */
516 I->Acc[b2] = true;
517 I->Don[b2] = false;
518 }
519 }
520 }
521 b++;
522 }
523
524 /* pass 3 */
525 b = obj->Bond;
526 for(a = 0; a < obj->NBond; a++) {
527 b1 = b->index[0];
528 b2 = b->index[1];
529
530 ai1 = obj_atomInfo + b1;
531 ai2 = obj_atomInfo + b2;
532
533 /* however, every NH is a donor,
534 even if it's SP2 */
535
536 if(ai1->protons == cAN_H) {
537
538 /* donors: any H attached to O, N */
539 switch (ai2->protons) {
540 case cAN_O:
541 I->Don[b1] = true;
542 I->Don[b2] = true; /* mark heavy atom too... */
543 break;
544 case cAN_N:
545 I->Don[b1] = true;
546 I->Don[b2] = true;
547 break;
548 }
549 } else if(ai2->protons == cAN_H) {
550 switch (ai1->protons) {
551 case cAN_O:
552 I->Don[b1] = true;
553 I->Don[b2] = true; /* mark heavy atom too... */
554 break;
555 case cAN_N:
556 I->Don[b1] = true;
557 I->Don[b2] = true; /* mark heavy atom too... */
558 break;
559 }
560 }
561
562 b++;
563 }
564
565 /* atom pass */
566 ai1 = obj_atomInfo;
567 for(a = 0; a < n_atom; a++) {
568 /* make sure all nonbonded atoms get categorized */
569
570 n0 = neighbor[a];
571 if(neighbor[n0] == 0) { /* nonbonded */
572 if(ai1->protons == cAN_O) {
573 I->Don[a] = true;
574 I->Acc[a] = true;
575 } else if(ai1->protons == cAN_N) {
576 I->Don[a] = true;
577 }
578 }
579 /*
580 if(I->Acc[a]) {
581 printf("ACC %s %s %s\n",ai1->chain,ai1->resi,ai1->name);
582 }
583 if(I->Don[a]) {
584 printf("DON %s %s %s\n",ai1->chain,ai1->resi,ai1->name);
585 } */
586
587 ai1++;
588 }
589
590 /* exclusions */
591 b = obj->Bond;
592 for(a = 0; a < obj->NBond; a++) {
593 b1 = b->index[0];
594 b2 = b->index[1];
595
596 ai1 = obj_atomInfo + b1;
597 ai2 = obj_atomInfo + b2;
598
599 xhash = ((b2 > b1) ? ex_hash(b1, b2) : ex_hash(b2, b1));
600 VLACheck(I->EXList, int, nex + 3);
601 j = I->EXList + nex;
602 *(j++) = I->EXHash[xhash];
603 if(b2 > b1) {
604 *(j++) = b1;
605 *(j++) = b2;
606 } else {
607 *(j++) = b2;
608 *(j++) = b1;
609 }
610 *(j++) = 2; /* 1-2 exclusion */
611 I->EXHash[xhash] = nex;
612 nex += 4;
613
614 a1 = crdidx[b1];
615 a2 = crdidx[b2];
616
617 if((a1 >= 0) && (a2 >= 0)) {
618 v1 = cs->coordPtr(a1);
619 v2 = cs->coordPtr(a2);
620 d = (float) diff3f(v1, v2);
621 if(use_cache) {
622 if(!SculptCacheQuery(G, cSculptBond,
623 obj_atomInfo[b1].unique_id,
624 obj_atomInfo[b2].unique_id, 0, 0, &d))
625 SculptCacheStore(G, cSculptBond,
626 obj_atomInfo[b1].unique_id,
627 obj_atomInfo[b2].unique_id, 0, 0, d);
628 }
629 ShakerAddDistCon(I->Shaker.get(), b1, b2, d, cShakerDistBond, 1.0F);
630 /* NOTE: storing atom indices, not coord. ind.! */
631 }
632 b++;
633 }
634
635 /* triangle relationships */
636 {
637 ATLCall atl;
638 ai1 = obj_atomInfo;
639
640 atl.G = I->G;
641 atl.Shaker = I->Shaker.get();
642 atl.ai = obj_atomInfo;
643 atl.cSet = cs;
644
645 if(obj->DiscreteFlag) {
646 atl.atm2idx = obj->DiscreteAtmToIdx;
647 atl.discCSet = obj->DiscreteCSet;
648 } else {
649 atl.atm2idx = cs->AtmToIdx;
650 atl.discCSet = NULL;
651 }
652 atl.coord = cs->Coord;
653 atl.neighbor = neighbor;
654 atl.min = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_sculpt_tri_min);
655 atl.max = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_sculpt_tri_max);
656 atl.mode =
657 SettingGet_i(G, cs->Setting, obj->Setting, cSetting_sculpt_tri_mode);
658
659 for(a = 0; a < n_atom; a++) {
660
661 atl.atom0 = a;
662
663 /* clear the flag -- TODO replace with array */
664 {
665 int aa;
666 ai = obj_atomInfo;
667 for(aa = 0; aa < n_atom; aa++) {
668 ai->temp1 = false;
669 ai++;
670 }
671 }
672
673 ai1->temp1 = true;
674 add_triangle_limits(&atl, a, a, 0.0F, 1);
675 ai1++;
676 }
677 }
678
679 /* if we have a match state, establish minimum distances */
680 if((match_state >= 0) && (match_state < obj->NCSet) && (!obj->DiscreteFlag)) {
681 CoordSet *cs2 = obj->CSet[match_state];
682 int n_site = 0;
683 if(cs2) {
684 float minim_min =
685 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_min_min);
686 float minim_max =
687 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_min_max);
688 float maxim_min =
689 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_max_min);
690 float maxim_max =
691 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_max_max);
692
693 int *site = pymol::calloc<int>(n_atom);
694 float *weight = pymol::calloc<float>(n_atom);
695 /* first, find candidate atoms with sufficient connectivity */
696 CountCall cnt;
697
698 cnt.ai = obj_atomInfo;
699 cnt.neighbor = neighbor;
700 cnt.atm2idx1 = cs->AtmToIdx;
701 cnt.atm2idx2 = cs2->AtmToIdx;
702
703 {
704 int aa;
705 ai = obj_atomInfo;
706 for(aa = 0; aa < n_atom; aa++) {
707 ai->temp1 = false;
708 ai++;
709 }
710 }
711
712 ai1 = obj_atomInfo;
713 for(b0 = 0; b0 < n_atom; b0++) {
714 int n_qual_branch = 0, cb;
715 int adj_site = false;
716 ai1->temp1 = true;
717 n0 = neighbor[b0] + 1;
718 while((b1 = neighbor[n0]) >= 0) {
719 if(site[b1]) {
720 adj_site = true;
721 break;
722 }
723 cb = count_branch(&cnt, b1, 3);
724 if(cb > 3) {
725 n_qual_branch++;
726 }
727 n0 += 2;
728 }
729 ai1->temp1 = false;
730 if((n_qual_branch > 2) && (!adj_site)) {
731 site[b0] = 10;
732 } else if(!adj_site) {
733 const char * name = LexStr(G, ai1->name);
734 switch (name[0]) {
735 case 'O':
736 if(!name[1])
737 if(AtomInfoKnownPolymerResName(LexStr(G, ai1->resn)))
738 site[b0] = 40; /* main-chain carbonyl */
739 break;
740 case 'C':
741 switch (name[1]) {
742 case 'Z':
743 switch (name[2]) {
744 case 0:
745 if(ai1->resn == G->lex_const.ARG)
746 site[b0] = 20; /* ARG/CZ */
747 else if(ai1->resn == G->lex_const.TYR)
748 site[b0] = 20; /* TYR/CZ */
749 else if(ai1->resn == G->lex_const.PHE)
750 site[b0] = 20; /* PHE/CZ */
751 break;
752 }
753 break;
754 case 'E':
755 switch (name[2]) {
756 case 0:
757 if(ai1->resn == G->lex_const.LYS)
758 site[b0] = 20; /* LYS/CE */
759 break;
760 }
761 break;
762 case 'D':
763 switch (name[2]) {
764 case 0:
765 if(ai1->resn == G->lex_const.GLU)
766 site[b0] = 20; /* GLU/CD */
767 else if(ai1->resn == G->lex_const.GLN)
768 site[b0] = 20; /* GLN/CD */
769 break;
770 }
771 break;
772 case 'G':
773 switch (name[2]) {
774 case 0:
775 if(ai1->resn == G->lex_const.LEU)
776 site[b0] = 20; /* LEU/CG */
777 else if(ai1->resn == G->lex_const.ASP)
778 site[b0] = 20; /* ASP/CG */
779 else if(ai1->resn == G->lex_const.ASN)
780 site[b0] = 20; /* ASN/CG */
781 break;
782 }
783 break;
784 }
785 break;
786 case 'S':
787 switch (name[1]) {
788 case 'D':
789 switch (name[2]) {
790 case 0:
791 if(ai1->resn == G->lex_const.MET)
792 site[b0] = 20; /* MET/SD */
793 break;
794 }
795 break;
796 }
797 break;
798 }
799 }
800 ai1++;
801 }
802
803 for(b0 = 0; b0 < n_atom; b0++) {
804 if(site[b0]) {
805 weight[n_site] = 10.0F / site[b0];
806 site[n_site] = b0;
807 n_site++;
808 }
809 }
810
811 {
812
813 for(a0 = 0; a0 < n_site; a0++) {
814 for(a1 = a0 + 1; a1 < n_site; a1++) {
815 float wt = weight[a0] * weight[a1];
816 b0 = site[a0];
817 b1 = site[a1];
818
819 {
820 int i0a = cs->AtmToIdx[b0];
821 int i1a = cs->AtmToIdx[b1];
822 int i0b = cs2->AtmToIdx[b0];
823 int i1b = cs2->AtmToIdx[b1];
824
825 if((i0a >= 0) && (i1a >= 0) && (i0b >= 0) && (i1b >= 0) &&
826 ((!match_by_segment)
827 || (obj_atomInfo[b0].segi == obj_atomInfo[b1].segi))) {
828 const float *v0a = cs->coordPtr(i0a);
829 const float *v1a = cs->coordPtr(i1a);
830 const float *v0b = cs2->coordPtr(i0b);
831 const float *v1b = cs2->coordPtr(i1b);
832 float dist0, dist1, min_dist, max_dist;
833 dist0 = diff3f(v0a, v1a);
834 dist1 = diff3f(v0b, v1b);
835 min_dist = (dist0 < dist1) ? dist0 : dist1;
836 if((min_dist >= minim_min) && (min_dist <= minim_max)) {
837 ShakerAddDistCon(I->Shaker.get(), b0, b1, min_dist, cShakerDistMinim, wt);
838 }
839 max_dist = (dist0 > dist1) ? dist0 : dist1;
840 if((max_dist >= maxim_min) && (max_dist <= maxim_max)) {
841 ShakerAddDistCon(I->Shaker.get(), b0, b1, max_dist, cShakerDistMaxim, wt);
842 }
843 }
844 }
845 }
846 }
847 }
848 FreeP(weight);
849 FreeP(site);
850 }
851 }
852 /* now pick up those 1-3 interations */
853
854 /* b1-b0-b2 */
855
856 for(b0 = 0; b0 < n_atom; b0++) {
857 n0 = neighbor[b0] + 1;
858 while(neighbor[n0] >= 0) {
859 b1 = neighbor[n0];
860 n1 = n0 + 2;
861 while(neighbor[n1] >= 0) {
862 b2 = neighbor[n1];
863
864 xhash = ((b2 > b1) ? ex_hash(b1, b2) : ex_hash(b2, b1));
865 VLACheck(I->EXList, int, nex + 3);
866 j = I->EXList + nex;
867 *(j++) = I->EXHash[xhash];
868 if(b2 > b1) {
869 *(j++) = b1;
870 *(j++) = b2;
871 } else {
872 *(j++) = b2;
873 *(j++) = b1;
874 }
875 *(j++) = 3; /* 1-3 exclusion */
876 I->EXHash[xhash] = nex;
877 nex += 4;
878
879 a0 = crdidx[b0];
880 a1 = crdidx[b1];
881 a2 = crdidx[b2];
882
883 if((a0 >= 0) && (a1 >= 0) && (a2 >= 0)) {
884 v1 = cs->coordPtr(a1);
885 v2 = cs->coordPtr(a2);
886 d = (float) diff3f(v1, v2);
887 if(use_cache) {
888 if(!SculptCacheQuery(G, cSculptAngl,
889 obj_atomInfo[b0].unique_id,
890 obj_atomInfo[b1].unique_id,
891 obj_atomInfo[b2].unique_id, 0, &d))
892 SculptCacheStore(G, cSculptAngl,
893 obj_atomInfo[b0].unique_id,
894 obj_atomInfo[b1].unique_id,
895 obj_atomInfo[b2].unique_id, 0, d);
896 }
897
898 ShakerAddDistCon(I->Shaker.get(), b1, b2, d, cShakerDistAngle, 1.0F);
899
900 if(linear[b0] && (linear[b1] || linear[b2])) {
901
902 if(use_cache) {
903 if(!SculptCacheQuery(G, cSculptLine,
904 obj_atomInfo[b1].unique_id,
905 obj_atomInfo[b0].unique_id,
906 obj_atomInfo[b2].unique_id, 0, &dummy))
907 SculptCacheStore(G, cSculptLine,
908 obj_atomInfo[b1].unique_id,
909 obj_atomInfo[b0].unique_id,
910 obj_atomInfo[b2].unique_id, 0, 0.0);
911 }
912 ShakerAddLineCon(I->Shaker.get(), b1, b0, b2);
913 }
914 }
915 n1 += 2;
916 }
917 n0 += 2;
918 }
919 }
920
921 /* and record the pyramidal and planar geometries */
922
923 /* b1-b0-b2
924 * |
925 * b3 */
926
927 for(b0 = 0; b0 < n_atom; b0++) {
928 n0 = neighbor[b0] + 1;
929 while(neighbor[n0] >= 0) {
930 b1 = neighbor[n0];
931 n1 = n0 + 2;
932 while(neighbor[n1] >= 0) {
933 b2 = neighbor[n1];
934 n2 = n1 + 2;
935 while(neighbor[n2] >= 0) {
936 b3 = neighbor[n2];
937
938 a0 = crdidx[b0];
939 a1 = crdidx[b1];
940 a2 = crdidx[b2];
941 a3 = crdidx[b3];
942
943 if((a0 >= 0) && (a1 >= 0) && (a2 >= 0) && (a3 >= 0)) {
944 float d2 = 0.0F;
945
946 v0 = cs->coordPtr(a0);
947 v1 = cs->coordPtr(a1);
948 v2 = cs->coordPtr(a2);
949 v3 = cs->coordPtr(a3);
950 d = ShakerGetPyra(&d2, v0, v1, v2, v3);
951
952 if(fabs(d) < 0.05) {
953 planar[b0] = true;
954 }
955 if(planar[b0])
956 d = 0.0;
957 if(use_cache) {
958 if(!SculptCacheQuery(G, cSculptPyra,
959 obj_atomInfo[b1].unique_id,
960 obj_atomInfo[b0].unique_id,
961 obj_atomInfo[b2].unique_id,
962 obj_atomInfo[b3].unique_id, &d))
963 SculptCacheStore(G, cSculptPyra,
964 obj_atomInfo[b1].unique_id,
965 obj_atomInfo[b0].unique_id,
966 obj_atomInfo[b2].unique_id,
967 obj_atomInfo[b3].unique_id, d);
968 if(!SculptCacheQuery(G, cSculptPyra + 1,
969 obj_atomInfo[b1].unique_id,
970 obj_atomInfo[b0].unique_id,
971 obj_atomInfo[b2].unique_id,
972 obj_atomInfo[b3].unique_id, &d2))
973 SculptCacheStore(G, cSculptPyra + 1,
974 obj_atomInfo[b1].unique_id,
975 obj_atomInfo[b0].unique_id,
976 obj_atomInfo[b2].unique_id,
977 obj_atomInfo[b3].unique_id, d2);
978 }
979 ShakerAddPyraCon(I->Shaker.get(), b0, b1, b2, b3, d, d2);
980 }
981 n2 += 2;
982 }
983 n1 += 2;
984 }
985 n0 += 2;
986 }
987 }
988
989 /* b1\b0_b2/b3 */
990
991 for(b0 = 0; b0 < n_atom; b0++) {
992 n0 = neighbor[b0] + 1;
993 while((b1 = neighbor[n0]) >= 0) {
994 n1 = neighbor[b0] + 1;
995 while((b2 = neighbor[n1]) >= 0) {
996 if(b1 != b2) {
997 n2 = neighbor[b2] + 1;
998 while((b3 = neighbor[n2]) >= 0) {
999 if((b3 != b0) && (b3 > b1)) {
1000 if(!(planar[b0] || planar[b2] || linear[b0] || linear[b2])) {
1001 int type;
1002 if((obj_atomInfo[b0].protons == cAN_S) &&
1003 (obj_atomInfo[b2].protons == cAN_S))
1004 type = cShakerTorsDisulfide;
1005 else
1006 type = cShakerTorsSP3SP3;
1007 ShakerAddTorsCon(I->Shaker.get(), b1, b0, b2, b3, type);
1008 }
1009 if(planar[b0] && planar[b2]) {
1010
1011 /* special extra-rigid torsion for hydrogens on
1012 planar acyclic systems (amides, etc.) */
1013
1014 if(((obj_atomInfo[b1].protons == cAN_H) && single[b1] &&
1015 (obj_atomInfo[b3].protons != cAN_H) && planar[b3]) ||
1016 ((obj_atomInfo[b3].protons == cAN_H) && single[b3] &&
1017 (obj_atomInfo[b1].protons != cAN_H) && planar[b1])) {
1018
1019 int cycle = 0;
1020 /* b1\b0_b2/b3-b4-b5-b6-b7... */
1021
1022 int b5, b6, b7, b8, b9, b10;
1023 int n4, n5, n6, n7, n8, n9;
1024 n3 = neighbor[b2] + 1;
1025 while((!cycle) && (b4 = neighbor[n3]) >= 0) {
1026 if(b4 != b0) {
1027 n4 = neighbor[b4] + 1;
1028 while((!cycle) && (b5 = neighbor[n4]) >= 0) {
1029 if(b5 != b2) {
1030 n5 = neighbor[b5] + 1;
1031 while((!cycle) && (b6 = neighbor[n5]) >= 0) {
1032 if(b6 == b0) { /* 4-cycle */
1033 cycle = 4;
1034 } else if((b6 != b4) && (b6 != b2)) {
1035 n6 = neighbor[b6] + 1;
1036 while((!cycle) && (b7 = neighbor[n6]) >= 0) {
1037 if(b7 == b0) { /* 5-cycle */
1038 cycle = 5;
1039 } else if((b7 != b5) && (b7 != b2)) {
1040 n7 = neighbor[b7] + 1;
1041 while((!cycle) && (b8 = neighbor[n7]) >= 0) {
1042 if(b8 == b0) { /* 6-cycle */
1043 cycle = 6;
1044 } else if((b8 != b6) && (b8 != b2)) {
1045 n8 = neighbor[b8] + 1;
1046 while((!cycle) && (b9 = neighbor[n8]) >= 0) {
1047 if(b9 == b0) { /* 7-cycle */
1048 cycle = 7;
1049 } else if((b9 != b7) && (b9 != b2)) {
1050 n9 = neighbor[b9] + 1;
1051 while((!cycle) && (b10 = neighbor[n9]) >= 0) {
1052 if(b10 == b0) { /* 8-cycle */
1053 cycle = 8;
1054 }
1055 n9 += 2;
1056 }
1057 }
1058 n8 += 2;
1059 }
1060 }
1061 n7 += 2;
1062 }
1063 }
1064 n6 += 2;
1065 }
1066 }
1067 n5 += 2;
1068 }
1069 }
1070 n4 += 2;
1071 }
1072 }
1073 n3 += 2;
1074 }
1075 if(!cycle) { /* don't add special amide constraints within small rings */
1076
1077 if(((obj_atomInfo[b1].protons == cAN_H) && single[b1] &&
1078 (obj_atomInfo[b0].protons == cAN_N) &&
1079 (obj_atomInfo[b2].protons == cAN_C) &&
1080 (obj_atomInfo[b3].protons == cAN_O) && planar[b3]) ||
1081 ((obj_atomInfo[b1].protons == cAN_H) && single[b3] &&
1082 (obj_atomInfo[b2].protons == cAN_N) &&
1083 (obj_atomInfo[b0].protons == cAN_C) &&
1084 (obj_atomInfo[b1].protons == cAN_O) && planar[b1])) {
1085 /* biased, asymmetric term for amides */
1086 ShakerAddTorsCon(I->Shaker.get(), b1, b0, b2, b3, cShakerTorsAmide);
1087 } else {
1088 /* biased, symmetric term for all others */
1089 ShakerAddTorsCon(I->Shaker.get(), b1, b0, b2, b3, cShakerTorsFlat);
1090 }
1091 }
1092 }
1093 }
1094 /* check 1-4 exclusion */
1095 xhash = ex_hash(b1, b3);
1096
1097 ex_type = 4;
1098
1099 xoffset = I->EXHash[xhash];
1100 while(xoffset) {
1101 k = I->EXList + xoffset;
1102 if((abs(*(k + 3)) == 4) && (*(k + 1) == b1) && (*(k + 2) == b3)) {
1103 if((b0 != *(k + 4)) && (b2 != *(k + 5))) {
1104 if(planar[b0] && planar[b2] &&
1105 planar[*(k + 4)] && planar[*(k + 5)]) {
1106 /* two planar paths -> likely a planar aromatic system */
1107 *(k + 3) = -4;
1108 }
1109 }
1110 ex_type = 0; /* duplicate, skip */
1111 break;
1112 }
1113 xoffset = *k;
1114 }
1115 if(ex_type) {
1116 VLACheck(I->EXList, int, nex + 5);
1117 j = I->EXList + nex;
1118 *(j++) = I->EXHash[xhash];
1119 *(j++) = b1;
1120 *(j++) = b3;
1121 if(planar[b0] && planar[b2])
1122 *(j++) = -4;
1123 else
1124 *(j++) = ex_type;
1125 *(j++) = b0;
1126 *(j++) = b2;
1127 I->EXHash[xhash]= nex;
1128
1129 nex += 6;
1130 }
1131
1132 /* planarity */
1133
1134 a0 = crdidx[b0];
1135 a1 = crdidx[b1];
1136 a2 = crdidx[b2];
1137 a3 = crdidx[b3];
1138
1139 if((a0 >= 0) && (a1 >= 0) && (a2 >= 0) && (a3 >= 0)) {
1140 v0 = cs->coordPtr(a0);
1141 v1 = cs->coordPtr(a1);
1142 v2 = cs->coordPtr(a2);
1143 v3 = cs->coordPtr(a3);
1144
1145 d = 0.0;
1146 if(planar[b0] && planar[b2]) {
1147 float deg = get_dihedral3f(v1, v0, v2, v3);
1148 if(fabs(deg) < deg_to_rad(10.0))
1149 d = 1.0;
1150 else if(fabs(deg) > deg_to_rad(170))
1151 d = -1.0;
1152
1153 {
1154 int cycle = false;
1155 /* look for 4, 5, 6, 7, or 8 cycle that
1156 connects back to b1 if found, then this
1157 planar system is fixed (either at zero
1158 or 180 -- it can't flip it over) */
1159 /* b1\b0_b2/b3-b4-b5-b6-b7... */
1160
1161 int b5, b6, b7, b8, b9, b10;
1162 int n4, n5, n6, n7, n8, n9;
1163 n3 = neighbor[b2] + 1;
1164 while((!cycle) && (b4 = neighbor[n3]) >= 0) {
1165 if(b4 != b0) {
1166 n4 = neighbor[b4] + 1;
1167 while((!cycle) && (b5 = neighbor[n4]) >= 0) {
1168 if(b5 != b2) {
1169 n5 = neighbor[b5] + 1;
1170 while((!cycle) && (b6 = neighbor[n5]) >= 0) {
1171 if(b6 == b0) { /* 4-cycle */
1172 cycle = 4;
1173 } else if((b6 != b4) && (b6 != b2)) {
1174 n6 = neighbor[b6] + 1;
1175 while((!cycle) && (b7 = neighbor[n6]) >= 0) {
1176 if(b7 == b0) { /* 5-cycle */
1177 cycle = 5;
1178 } else if((b7 != b5) && (b7 != b2)) {
1179 n7 = neighbor[b7] + 1;
1180 while((!cycle) && (b8 = neighbor[n7]) >= 0) {
1181 if(b8 == b0) { /* 6-cycle */
1182 cycle = 6;
1183 } else if((b8 != b6) && (b8 != b2)) {
1184 n8 = neighbor[b8] + 1;
1185 while((!cycle) && (b9 = neighbor[n8]) >= 0) {
1186 if(b9 == b0) { /* 7-cycle */
1187 cycle = 7;
1188 } else if((b9 != b7) && (b9 != b2)) {
1189 n9 = neighbor[b9] + 1;
1190 while((!cycle)
1191 && (b10 = neighbor[n9]) >= 0) {
1192 if(b10 == b0) { /* 8-cycle */
1193 cycle = 8;
1194 }
1195 n9 += 2;
1196 }
1197 }
1198 n8 += 2;
1199 }
1200 }
1201 n7 += 2;
1202 }
1203 }
1204 n6 += 2;
1205 }
1206 }
1207 n5 += 2;
1208 }
1209 }
1210 n4 += 2;
1211 }
1212 }
1213 n3 += 2;
1214 }
1215 /* don't get jacked by pseudo-planar PRO */
1216
1217 if(((obj_atomInfo[b0].protons != cAN_N) ||
1218 (!WordMatchExact(G, obj_atomInfo[b0].resn, G->lex_const.PRO, true))) &&
1219 ((obj_atomInfo[b2].protons != cAN_N) ||
1220 (!WordMatchExact(G, obj_atomInfo[b2].resn, G->lex_const.PRO, true)))) {
1221
1222 if(use_cache) {
1223 if(!SculptCacheQuery(G, cSculptPlan,
1224 obj_atomInfo[b1].unique_id,
1225 obj_atomInfo[b0].unique_id,
1226 obj_atomInfo[b2].unique_id,
1227 obj_atomInfo[b3].unique_id, &d))
1228 SculptCacheStore(G, cSculptPlan,
1229 obj_atomInfo[b1].unique_id,
1230 obj_atomInfo[b0].unique_id,
1231 obj_atomInfo[b2].unique_id,
1232 obj_atomInfo[b3].unique_id, d);
1233 }
1234
1235 ShakerAddPlanCon(I->Shaker.get(), b1, b0, b2, b3, d, cycle);
1236
1237 if(planar[b1] && planar[b3] && ((cycle == 5) || (cycle == 6))) {
1238
1239 /* also add minimum distance constraints to keep small rings from folding */
1240
1241 d = (float) diff3f(v1, v3);
1242
1243 ShakerAddDistCon(I->Shaker.get(), b1, b3, d, cShakerDistBond, 1.0F);
1244
1245 }
1246 }
1247 }
1248 }
1249 }
1250 }
1251 n2 += 2;
1252 }
1253 }
1254 n1 += 2;
1255 }
1256 n0 += 2;
1257 }
1258 }
1259
1260 /* add 1,5 exclusions for hydrogens off arg-like planar systems */
1261
1262 /* b1\b0_b2_b3/b4 */
1263
1264 for(b0 = 0; b0 < n_atom; b0++) {
1265 n0 = neighbor[b0] + 1;
1266 while((b1 = neighbor[n0]) >= 0) {
1267 if(obj_atomInfo[b1].protons == cAN_H) {
1268 n1 = neighbor[b0] + 1;
1269 while((b2 = neighbor[n1]) >= 0) {
1270 if(b1 != b2) {
1271 n2 = neighbor[b2] + 1;
1272 while((b3 = neighbor[n2]) >= 0) {
1273 if(b3 != b0) {
1274 if(planar[b0] && planar[b2] && planar[b3]) {
1275 n3 = neighbor[b3] + 1;
1276 while((b4 = neighbor[n3]) >= 0) {
1277 if((b4 != b2) && (b4 > b1) && (obj_atomInfo[b4].protons == cAN_H)) {
1278
1279 xhash = ex_hash(b1, b4);
1280
1281 ex_type = 5;
1282
1283 xoffset = I->EXHash[xhash];
1284 while(xoffset) {
1285 k = I->EXList + xoffset;
1286 if(((*(k + 3)) == ex_type) &&
1287 (*(k + 1) == b1) && (*(k + 2) == b4)) {
1288 ex_type = 0; /* duplicate, skip */
1289 break;
1290 }
1291 xoffset = *k;
1292 }
1293
1294 if(ex_type) {
1295 VLACheck(I->EXList, int, nex + 6);
1296 j = I->EXList + nex;
1297 *(j++) = I->EXHash[xhash];
1298 *(j++) = b1;
1299 *(j++) = b4;
1300 *(j++) = ex_type;
1301 *(j++) = b0;
1302 *(j++) = b2;
1303 *(j++) = b3;
1304 I->EXHash[xhash] = nex;
1305 nex += 6;
1306 }
1307 }
1308 n3 += 2;
1309 }
1310 }
1311 }
1312 n2 += 2;
1313 }
1314 }
1315 n1 += 2;
1316 }
1317 }
1318 n0 += 2;
1319 }
1320 }
1321
1322 {
1323 /* longer-range exclusions (1-5,1-6,1-7,1-8,1-9) -- only locate & store when needed */
1324
1325 int mask =
1326 SettingGet_i(G, cs->Setting, obj->Setting, cSetting_sculpt_field_mask);
1327 int max_excl =
1328 SettingGet_i(G, cs->Setting, obj->Setting, cSetting_sculpt_avd_excl);
1329 if(max_excl > 9)
1330 max_excl = 9;
1331
1332 if((cSculptAvoid & mask) && (max_excl > 4)) {
1333 int b_stack[10];
1334 int n_stack[10];
1335 int stop_depth = max_excl - 1;
1336 int depth;
1337 int bd, skip;
1338 for(b0 = 0; b0 < n_atom; b0++) {
1339 b_stack[0] = b0;
1340 n_stack[0] = neighbor[b_stack[0]] + 1;
1341 depth = 0;
1342 while(depth >= 0) {
1343 if((bd = neighbor[n_stack[depth]]) < 0) {
1344 depth--;
1345 if(depth >= 0) { /* iterate next atom */
1346 n_stack[depth] += 2;
1347 }
1348 } else {
1349 skip = (depth == stop_depth);
1350 if(!skip) {
1351 for(a = 0; a < depth; a++) {
1352 if(b_stack[a] == bd) {
1353 skip = true;
1354 break;
1355 }
1356 }
1357 }
1358 if(!skip) {
1359 depth++;
1360 b_stack[depth] = bd;
1361 n_stack[depth] = neighbor[bd] + 1;
1362 if((depth > 3) && (b0 < bd)) {
1363
1364 xhash = ex_hash(b0, bd);
1365
1366 VLACheck(I->EXList, int, nex + 3);
1367 j = I->EXList + nex;
1368 *(j++) = I->EXHash[xhash];
1369 *(j++) = b0;
1370 *(j++) = bd;
1371 *(j++) = depth + 1; /* 1-5, 1-6, 1-7 etc. */
1372 I->EXHash[xhash] = nex;
1373 nex += 4;
1374 }
1375 } else {
1376 n_stack[depth] += 2;
1377 }
1378 }
1379 }
1380 }
1381 }
1382 }
1383 FreeP(planar);
1384 FreeP(linear);
1385 FreeP(single);
1386 FreeP(crdidx);
1387 }
1388 }
1389
1390 PRINTFB(G, FB_Sculpt, FB_Blather)
1391 " Sculpt: I->Shaker->NDistCon %d\n", I->Shaker->NDistCon ENDFB(G);
1392 PRINTFB(G, FB_Sculpt, FB_Blather)
1393 " Sculpt: I->Shaker->NPyraCon %d\n", I->Shaker->NPyraCon ENDFB(G);
1394 PRINTFB(G, FB_Sculpt, FB_Blather)
1395 " Sculpt: I->Shaker->NPlanCon %d\n", I->Shaker->NPlanCon ENDFB(G);
1396
1397 PRINTFD(G, FB_Sculpt)
1398 " SculptMeasureObject-Debug: leaving...\n" ENDFD;
1399
1400 }
1401
SculptCheckBump(float * v1,float * v2,float * diff,float * dist,float cutoff)1402 static int SculptCheckBump(float *v1, float *v2, float *diff, float *dist, float cutoff)
1403 {
1404 float d2;
1405 diff[0] = (v1[0] - v2[0]);
1406 diff[1] = (v1[1] - v2[1]);
1407 if(fabs(diff[0]) > cutoff)
1408 return (false);
1409 diff[2] = (v1[2] - v2[2]);
1410 if(fabs(diff[1]) > cutoff)
1411 return (false);
1412 if(fabs(diff[2]) > cutoff)
1413 return (false);
1414 d2 = (diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]);
1415 if(d2 < (cutoff * cutoff)) {
1416 *dist = (float) sqrt(d2);
1417 return (true);
1418 }
1419 return (false);
1420 }
1421
SculptCGOBump(float * v1,float * v2,float vdw1,float vdw2,float cutoff,float min,float mid,float max,float * good_color,float * bad_color,int mode,CGO * cgo)1422 static int SculptCGOBump(float *v1, float *v2,
1423 float vdw1, float vdw2,
1424 float cutoff,
1425 float min, float mid, float max,
1426 float *good_color, float *bad_color, int mode, CGO * cgo)
1427 {
1428 float d2;
1429 float diff[3];
1430 float dist;
1431 float min_cutoff = cutoff - min;
1432 diff[0] = (v1[0] - v2[0]);
1433 diff[1] = (v1[1] - v2[1]);
1434 if(fabs(diff[0]) > min_cutoff)
1435 return (false);
1436 diff[2] = (v1[2] - v2[2]);
1437 if(fabs(diff[1]) > min_cutoff)
1438 return (false);
1439 if(fabs(diff[2]) > min_cutoff)
1440 return (false);
1441 d2 = (diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]);
1442 if(d2 > (min_cutoff * min_cutoff)) {
1443 return false;
1444 } else {
1445 dist = (float) sqrt(d2);
1446 if(dist <= min_cutoff) {
1447 float avg[3], vv1[3], vv2[3], tmp[3], color[3];
1448 float good_bad = cutoff - dist; /* if negative, then good */
1449 float color_factor;
1450 float radius = 0.5 * (good_bad - min);
1451
1452 if(good_bad < mid) {
1453 color_factor = 0.0F;
1454 } else {
1455 color_factor = (good_bad - mid) / max;
1456 if(color_factor > 1.0F)
1457 color_factor = 1.0F;
1458 }
1459
1460 {
1461 float one_minus_color_factor = 1.0F - color_factor;
1462 scale3f(bad_color, color_factor, color);
1463 scale3f(good_color, one_minus_color_factor, tmp);
1464 add3f(tmp, color, color);
1465
1466 switch (mode) {
1467 case 2:
1468 if(good_bad > mid) {
1469 CGOLinewidth(cgo, 1 + color_factor * 3);
1470 CGOColorv(cgo, color);
1471 {
1472 float *vertexVals = CGODrawArrays(cgo, GL_LINES, CGO_VERTEX_ARRAY, 2);
1473 copy3f(v1, vertexVals);
1474 copy3f(v2, &vertexVals[3]);
1475 }
1476 }
1477 break;
1478 case 1:
1479 {
1480 float delta, one_minus_delta;
1481 if(good_bad < 0.0) {
1482 delta = fabs(good_bad);
1483 } else {
1484 delta = 0.5 * (0.01F + fabs(good_bad)) / (cutoff);
1485 }
1486 if(delta < 0.01F)
1487 delta = 0.01F;
1488 if(delta > 0.1F) {
1489 delta = 0.1F;
1490 }
1491 if(radius < 0.01F)
1492 radius = 0.01F;
1493
1494 one_minus_delta = 1.0F - delta;
1495 scale3f(v2, vdw1, avg);
1496 scale3f(v1, vdw2, tmp);
1497 add3f(tmp, avg, avg);
1498 {
1499 float inv = 1.0F / (vdw1 + vdw2);
1500 scale3f(avg, inv, avg);
1501 }
1502 scale3f(v1, delta, vv1);
1503 scale3f(avg, one_minus_delta, tmp);
1504 add3f(tmp, vv1, vv1);
1505 scale3f(v2, delta, vv2);
1506 scale3f(avg, one_minus_delta, tmp);
1507 add3f(tmp, vv2, vv2);
1508 if(good_bad < 0.0F) {
1509 CGOLinewidth(cgo, 1 + color_factor * 3);
1510 CGOResetNormal(cgo, true);
1511 CGOColorv(cgo, color);
1512 {
1513 float *vertexVals = CGODrawArrays(cgo, GL_LINES, CGO_VERTEX_ARRAY, 2);
1514 copy3f(vv1, vertexVals);
1515 copy3f(vv2, &vertexVals[3]);
1516 }
1517 } else {
1518 CGOCustomCylinderv(cgo, vv1, vv2, radius, color, color, 1, 1);
1519 }
1520 }
1521 break;
1522 }
1523 }
1524 }
1525 if(dist > cutoff)
1526 return false;
1527 return (true);
1528 }
1529 }
1530
SculptDoBump(float target,float actual,float * d,float * d0to1,float * d1to0,float wt,float * strain)1531 static int SculptDoBump(float target, float actual, float *d,
1532 float *d0to1, float *d1to0, float wt, float *strain)
1533 {
1534 float push[3];
1535 float dev, dev_2, sc, abs_dev;
1536
1537 dev = target - actual;
1538 if((abs_dev = (float) fabs(dev)) > R_SMALL8) {
1539 dev_2 = wt * dev / 2.0F;
1540 (*strain) += abs_dev;
1541 if(actual > R_SMALL8) { /* nonoverlapping */
1542 sc = dev_2 / actual;
1543 scale3f(d, sc, push);
1544 add3f(push, d0to1, d0to1);
1545 subtract3f(d1to0, push, d1to0);
1546 } else { /* overlapping, so just push along X */
1547 d0to1[0] -= dev_2;
1548 d1to0[0] += dev_2;
1549 }
1550 return 1;
1551 }
1552 return 0;
1553 }
1554
SculptCheckAvoid(float * v1,float * v2,float * diff,float * dist,float avoid,float range)1555 static int SculptCheckAvoid(float *v1, float *v2, float *diff,
1556 float *dist, float avoid, float range)
1557 {
1558 float d2, l2;
1559 float cutoff = avoid + range;
1560 float low_cutoff;
1561 diff[0] = (v1[0] - v2[0]);
1562 diff[1] = (v1[1] - v2[1]);
1563 if(fabs(diff[0]) > cutoff)
1564 return (false);
1565 diff[2] = (v1[2] - v2[2]);
1566 if(fabs(diff[1]) > cutoff)
1567 return (false);
1568 low_cutoff = avoid - range;
1569 if(fabs(diff[2]) > cutoff)
1570 return (false);
1571 l2 = low_cutoff * low_cutoff;
1572 d2 = (diff[0] * diff[0] + diff[1] * diff[1] + diff[2] * diff[2]);
1573 if((d2 < (cutoff * cutoff)) && (d2 > l2)) { /* we are in the avoid range */
1574 *dist = (float) sqrt(d2);
1575 return (true);
1576 }
1577 return (false);
1578 }
1579
SculptDoAvoid(float avoid,float range,float actual,float * d,float * d0to1,float * d1to0,float wt,float * strain)1580 static int SculptDoAvoid(float avoid, float range, float actual, float *d,
1581 float *d0to1, float *d1to0, float wt, float *strain)
1582 {
1583 float push[3];
1584 float dev, dev_2, sc, abs_dev;
1585 float target;
1586 if(actual > avoid) {
1587 target = avoid + range;
1588 } else {
1589 target = avoid - range;
1590 }
1591 dev = target - actual;
1592 if((abs_dev = (float) fabs(dev)) > R_SMALL8) {
1593 dev_2 = wt * dev / 2.0F;
1594 (*strain) += abs_dev;
1595 if(actual > R_SMALL8) { /* nonoverlapping */
1596 sc = dev_2 / actual;
1597 scale3f(d, sc, push);
1598 add3f(push, d0to1, d0to1);
1599 subtract3f(d1to0, push, d1to0);
1600 } else { /* overlapping, so just push along X */
1601 d0to1[0] -= dev_2;
1602 d1to0[0] += dev_2;
1603 }
1604 return 1;
1605 }
1606 return 0;
1607 }
1608
SculptIterateObject(CSculpt * I,ObjectMolecule * obj,int state,int n_cycle,float * center)1609 float SculptIterateObject(CSculpt * I, ObjectMolecule * obj,
1610 int state, int n_cycle, float *center)
1611 {
1612 PyMOLGlobals *G = I->G;
1613 CShaker *shk;
1614 int a0, a1, a2, a3, b0, b3;
1615 int aa;
1616 float *disp = NULL;
1617 float *v, *v0, *v1, *v2, *v3;
1618 float diff[3], len;
1619 int *atm2idx = NULL;
1620 int *cnt = NULL;
1621 int *i;
1622 int hash;
1623 int nb_next;
1624 int h, k, l;
1625 int offset;
1626 float cutoff, vdw_cutoff;
1627 int ex;
1628 int eval_flag;
1629 int mask;
1630 float wt;
1631 float vdw;
1632 float vdw14;
1633 float vdw_wt;
1634 float vdw_wt14;
1635 float bond_wt;
1636 float angl_wt;
1637 float pyra_wt, pyra_inv_wt;
1638 float plan_wt;
1639 float line_wt;
1640 float tors_wt;
1641 float tors_tole;
1642 int active_flag = false;
1643 float hb_overlap, hb_overlap_base;
1644 int *active, n_active;
1645 int *exclude;
1646 const AtomInfoType *ai0, *ai1;
1647 double task_time;
1648 float vdw_magnify, vdw_magnified = 1.0F;
1649 int nb_skip, nb_skip_count;
1650 float total_strain = 0.0F, strain;
1651 int total_count = 1;
1652 CGO *cgo = NULL;
1653 float good_color[3] = { 0.2, 1.0, 0.2 };
1654 float bad_color[3] = { 1.0, 0.2, 0.2 };
1655 int vdw_vis_mode;
1656 float vdw_vis_min = 0.0F, vdw_vis_mid = 0.0F, vdw_vis_max = 0.0F;
1657 float tri_sc, tri_wt;
1658 float min_sc, min_wt;
1659 float max_sc = 1.025F, max_wt = 0.75F;
1660 float *cs_coord;
1661 float solvent_radius;
1662 float avd_wt, avd_gp, avd_rg;
1663 int avd_ex;
1664
1665 PRINTFD(G, FB_Sculpt)
1666 " SculptIterateObject-Debug: entered state=%d n_cycle=%d\n", state, n_cycle ENDFD;
1667 if(!n_cycle)
1668 n_cycle = -1;
1669
1670 auto cs = obj->getCoordSet(state);
1671
1672 if (cs) {
1673
1674 disp = pymol::malloc<float>(3 * obj->NAtom);
1675 atm2idx = pymol::malloc<int>(obj->NAtom);
1676 cnt = pymol::malloc<int>(obj->NAtom);
1677 active = pymol::malloc<int>(obj->NAtom);
1678 exclude = pymol::calloc<int>(obj->NAtom);
1679 shk = I->Shaker.get();
1680
1681 PRINTFD(G, FB_Sculpt)
1682 " SIO-Debug: NDistCon %d\n", shk->NDistCon ENDFD;
1683
1684 cs_coord = cs->Coord.data();
1685
1686 vdw = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_vdw_scale);
1687 vdw14 = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_vdw_scale14);
1688 vdw_wt = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_vdw_weight);
1689 vdw_wt14 =
1690 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_vdw_weight14);
1691 bond_wt = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_bond_weight);
1692 angl_wt = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_angl_weight);
1693 pyra_wt = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_pyra_weight);
1694 pyra_inv_wt =
1695 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_pyra_inv_weight);
1696 plan_wt = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_plan_weight);
1697 line_wt = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_line_weight);
1698 tri_wt = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_tri_weight);
1699 tri_sc = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_tri_scale);
1700
1701 min_wt = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_min_weight);
1702 min_sc = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_min_scale);
1703 max_wt = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_max_weight);
1704 max_sc = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_max_scale);
1705
1706 mask = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_sculpt_field_mask);
1707 hb_overlap =
1708 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_hb_overlap);
1709 hb_overlap_base =
1710 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_hb_overlap_base);
1711 tors_tole =
1712 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_tors_tolerance);
1713 tors_wt = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_tors_weight);
1714 vdw_vis_mode =
1715 SettingGet_i(G, cs->Setting, obj->Setting, cSetting_sculpt_vdw_vis_mode);
1716 solvent_radius =
1717 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_solvent_radius);
1718
1719 avd_wt = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_avd_weight);
1720 avd_gp = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_avd_gap);
1721 avd_rg = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_avd_range);
1722 avd_ex = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_sculpt_avd_excl);
1723 if(avd_gp < 0.0F)
1724 avd_gp = 1.5F * solvent_radius;
1725 if(avd_rg < 0.0F)
1726 avd_rg = solvent_radius;
1727
1728 if(vdw_vis_mode) {
1729 vdw_vis_min =
1730 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_vdw_vis_min);
1731 vdw_vis_mid =
1732 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_vdw_vis_mid);
1733 vdw_vis_max =
1734 SettingGet_f(G, cs->Setting, obj->Setting, cSetting_sculpt_vdw_vis_max);
1735
1736 if(!cs->SculptCGO)
1737 cs->SculptCGO = CGONew(G);
1738 else
1739 CGOReset(cs->SculptCGO);
1740 } else if(cs->SculptCGO) {
1741 CGOReset(cs->SculptCGO);
1742 }
1743 cgo = cs->SculptCGO;
1744
1745 nb_skip = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_sculpt_nb_interval);
1746 if(nb_skip > n_cycle)
1747 nb_skip = n_cycle;
1748 if(nb_skip < 0)
1749 nb_skip = 0;
1750
1751 n_active = 0;
1752 ai0 = obj->AtomInfo;
1753 {
1754 int a;
1755 for(a = 0; a < obj->NAtom; a++) {
1756 if(ai0->flags & cAtomFlag_exclude) {
1757 exclude[a] = true;
1758 a1 = -1;
1759 } else {
1760 a1 = cs->atmToIdx(a);
1761 }
1762 if(a1 >= 0) {
1763 active_flag = true;
1764 active[n_active] = a;
1765 n_active++;
1766 }
1767 atm2idx[a] = a1;
1768 ai0++;
1769 }
1770 }
1771
1772 if(active_flag) {
1773
1774 /* first, create coordinate -> vertex mapping */
1775 /* and count number of constraints */
1776
1777 task_time = UtilGetSeconds(G);
1778 vdw_magnify = 1.0F;
1779 nb_skip_count = 0;
1780
1781 if(center) {
1782 int *a_ptr = active;
1783 int a;
1784 for(aa = 0; aa < n_active; aa++) {
1785 a = *(a_ptr++);
1786 {
1787 AtomInfoType *ai = obj->AtomInfo + a;
1788 if((ai->protekted != 1) && !(ai->flags & cAtomFlag_fix)) {
1789 v2 = cs_coord + 3 * atm2idx[a];
1790 center[4] += *(v2);
1791 center[5] += *(v2 + 1);
1792 center[6] += *(v2 + 2);
1793 center[7] += 1.0F;
1794 }
1795 }
1796 }
1797 }
1798
1799 while(n_cycle--) {
1800
1801 total_strain = 0.0F;
1802 total_count = 0;
1803 /* initialize displacements to zero */
1804
1805 v = disp;
1806 i = cnt;
1807 for(aa = 0; aa < n_active; aa++) {
1808 int a = active[aa];
1809 v = disp + a * 3;
1810 cnt[a] = 0;
1811 *(v) = 0.0F;
1812 *(v + 1) = 0.0F;
1813 *(v + 2) = 0.0F;
1814 }
1815
1816 /* apply distance constraints */
1817
1818 {
1819 const ShakerDistCon *sdc = shk->DistCon.data();
1820 int a,ndc = shk->NDistCon;
1821 for(a = 0; a < ndc; a++) {
1822 int sdc_type = sdc->type;
1823 int b1 = sdc->at0;
1824 int b2 = sdc->at1;
1825
1826 switch (sdc_type) {
1827 case cShakerDistBond:
1828 eval_flag = cSculptBond & mask;
1829 wt = bond_wt;
1830 break;
1831 case cShakerDistAngle:
1832 eval_flag = cSculptAngl & mask;
1833 wt = angl_wt;
1834 break;
1835 case cShakerDistLimit:
1836 eval_flag = cSculptTri & mask;
1837 wt = tri_wt;
1838 break;
1839 case cShakerDistMinim:
1840 eval_flag = cSculptMin & mask;
1841 wt = min_wt * sdc->weight;
1842 break;
1843 case cShakerDistMaxim:
1844 eval_flag = cSculptMax & mask;
1845 wt = max_wt * sdc->weight;
1846 break;
1847 default:
1848 eval_flag = false;
1849 wt = 0.0F;
1850 break;
1851 }
1852
1853 if(eval_flag && !(exclude[b1] || exclude[b2])) {
1854 a1 = atm2idx[b1]; /* coordinate set indices */
1855 a2 = atm2idx[b2];
1856 if((a1 >= 0) && (a2 >= 0)) {
1857 v1 = cs_coord + 3 * a1;
1858 v2 = cs_coord + 3 * a2;
1859 switch (sdc_type) {
1860 case cShakerDistLimit:
1861 strain =
1862 ShakerDoDistLimit(sdc->targ * tri_sc, v1, v2, disp + b1 * 3,
1863 disp + b2 * 3, wt);
1864 if(strain > 0.0F) {
1865 cnt[b1]++;
1866 cnt[b2]++;
1867 total_strain += strain;
1868 total_count++;
1869 }
1870 break;
1871 case cShakerDistMaxim:
1872 strain =
1873 ShakerDoDistLimit(sdc->targ * max_sc, v1, v2, disp + b1 * 3,
1874 disp + b2 * 3, wt);
1875 if(strain > 0.0F) {
1876 cnt[b1]++;
1877 cnt[b2]++;
1878 total_strain += strain;
1879 total_count++;
1880 }
1881 break;
1882 case cShakerDistMinim:
1883 strain =
1884 ShakerDoDistMinim(sdc->targ * min_sc, v1, v2, disp + b1 * 3,
1885 disp + b2 * 3, wt);
1886 if(strain > 0.0F) {
1887 cnt[b1]++;
1888 cnt[b2]++;
1889 total_strain += strain;
1890 total_count++;
1891 }
1892 break;
1893 default:
1894 total_strain +=
1895 ShakerDoDist(sdc->targ, v1, v2, disp + b1 * 3, disp + b2 * 3, wt);
1896 cnt[b1]++;
1897 cnt[b2]++;
1898 total_count++;
1899 }
1900 }
1901 }
1902 sdc++;
1903 }
1904 }
1905 /* apply line constraints */
1906
1907 if(cSculptLine & mask) {
1908 const ShakerLineCon *slc = shk->LineCon.data();
1909 int nlc = shk->NLineCon;
1910 int a,b1,b2;
1911 for(a = 0; a < nlc; a++) {
1912 b0 = slc->at0;
1913 b1 = slc->at1;
1914 b2 = slc->at2;
1915 a0 = atm2idx[b0]; /* coordinate set indices */
1916 a1 = atm2idx[b1];
1917 a2 = atm2idx[b2];
1918
1919 if((a0 >= 0) && (a1 >= 0) && (a2 >= 0)
1920 && !(exclude[b0] || exclude[b1] || exclude[b2])) {
1921 cnt[b0]++;
1922 cnt[b1]++;
1923 cnt[b2]++;
1924 v0 = cs_coord + 3 * a0;
1925 v1 = cs_coord + 3 * a1;
1926 v2 = cs_coord + 3 * a2;
1927 total_strain +=
1928 ShakerDoLine(v0, v1, v2, disp + b0 * 3, disp + b1 * 3, disp + b2 * 3,
1929 line_wt);
1930 total_count++;
1931 }
1932 slc++;
1933 }
1934 }
1935
1936 /* apply pyramid constraints */
1937
1938 if(cSculptPyra & mask) {
1939 const ShakerPyraCon *spc = shk->PyraCon.data();
1940 int npc = shk->NPyraCon;
1941 int a,b1,b2;
1942 for(a = 0; a < npc; a++) {
1943
1944 b0 = spc->at0;
1945 b1 = spc->at1;
1946 b2 = spc->at2;
1947 b3 = spc->at3;
1948 a0 = atm2idx[b0];
1949 a1 = atm2idx[b1];
1950 a2 = atm2idx[b2];
1951 a3 = atm2idx[b3];
1952
1953 if((a0 >= 0) && (a1 >= 0) && (a2 >= 0) && (a3 >= 0)
1954 && !(exclude[b0] || exclude[b1] || exclude[b2] || exclude[b3])) {
1955 v0 = cs_coord + 3 * a0;
1956 v1 = cs_coord + 3 * a1;
1957 v2 = cs_coord + 3 * a2;
1958 v3 = cs_coord + 3 * a3;
1959 total_strain += ShakerDoPyra(spc->targ1,
1960 spc->targ2,
1961 v0, v1, v2, v3,
1962 disp + b0 * 3,
1963 disp + b1 * 3,
1964 disp + b2 * 3,
1965 disp + b3 * 3, pyra_wt, pyra_inv_wt);
1966 total_count++;
1967
1968 cnt[b0]++;
1969 cnt[b1]++;
1970 cnt[b2]++;
1971 cnt[b3]++;
1972 }
1973 spc++;
1974 }
1975 }
1976
1977 if(cSculptPlan & mask) {
1978 const ShakerPlanCon *snc = shk->PlanCon.data();
1979 int npc = shk->NPlanCon;
1980 int a,b1,b2;
1981 /* apply planarity constraints */
1982
1983 for(a = 0; a < npc; a++) {
1984
1985 b0 = snc->at0;
1986 b1 = snc->at1;
1987 b2 = snc->at2;
1988 b3 = snc->at3;
1989 a0 = atm2idx[b0];
1990 a1 = atm2idx[b1];
1991 a2 = atm2idx[b2];
1992 a3 = atm2idx[b3];
1993
1994 if((a0 >= 0) && (a1 >= 0) && (a2 >= 0) && (a3 >= 0)
1995 && !(exclude[b0] || exclude[b1] || exclude[b2] || exclude[b3])) {
1996 v0 = cs_coord + 3 * a0;
1997 v1 = cs_coord + 3 * a1;
1998 v2 = cs_coord + 3 * a2;
1999 v3 = cs_coord + 3 * a3;
2000 total_strain += ShakerDoPlan(v0, v1, v2, v3,
2001 disp + b0 * 3,
2002 disp + b1 * 3,
2003 disp + b2 * 3,
2004 disp + b3 * 3,
2005 snc->target, snc->fixed, plan_wt);
2006 total_count++;
2007 cnt[b0]++;
2008 cnt[b1]++;
2009 cnt[b2]++;
2010 cnt[b3]++;
2011 }
2012
2013 snc++;
2014 }
2015 }
2016
2017 /* apply torsion constraints */
2018
2019 if(cSculptTors & mask) {
2020 const ShakerTorsCon *stc = shk->TorsCon.data();
2021 int ntc = shk->NTorsCon;
2022 int a,b1,b2;
2023 /* apply planarity constraints */
2024
2025 for(a = 0; a < ntc; a++) {
2026
2027 b0 = stc->at0;
2028 b1 = stc->at1;
2029 b2 = stc->at2;
2030 b3 = stc->at3;
2031 a0 = atm2idx[b0];
2032 a1 = atm2idx[b1];
2033 a2 = atm2idx[b2];
2034 a3 = atm2idx[b3];
2035
2036 if((a0 >= 0) && (a1 >= 0) && (a2 >= 0) && (a3 >= 0)
2037 && !(exclude[b0] || exclude[b1] || exclude[b2] || exclude[b3])) {
2038 v0 = cs_coord + 3 * a0;
2039 v1 = cs_coord + 3 * a1;
2040 v2 = cs_coord + 3 * a2;
2041 v3 = cs_coord + 3 * a3;
2042 total_strain += ShakerDoTors(stc->type,
2043 v0, v1, v2, v3,
2044 disp + b0 * 3,
2045 disp + b1 * 3,
2046 disp + b2 * 3,
2047 disp + b3 * 3, tors_tole, tors_wt);
2048 total_count++;
2049 cnt[b0]++;
2050 cnt[b1]++;
2051 cnt[b2]++;
2052 cnt[b3]++;
2053 }
2054 stc++;
2055 }
2056 }
2057 /* apply nonbonded interactions */
2058
2059 if((n_cycle > 0) && (nb_skip_count > 0)) {
2060 /*skip and then weight extra */
2061 nb_skip_count--;
2062 vdw_magnify += 1.0F;
2063 } else {
2064 int nb_off0, nb_off1;
2065 int v0i, v1i, v2i;
2066 int x0i;
2067 int don_b0;
2068 int acc_b0;
2069 int b1;
2070 vdw_magnified = vdw_magnify;
2071 vdw_magnify = 1.0F;
2072
2073 nb_skip_count = nb_skip;
2074 if((cSculptVDW | cSculptVDW14 | cSculptAvoid) & mask) {
2075 /* compute non-bonded interations */
2076
2077 /* construct nonbonded hash */
2078
2079 nb_next = 1;
2080 for(aa = 0; aa < n_active; aa++) {
2081 b0 = active[aa];
2082 a0 = atm2idx[b0];
2083 VLACheck(I->NBList, int, nb_next + 2);
2084 v0 = cs_coord + 3 * a0;
2085 hash = nb_hash(v0);
2086 i = I->NBList + nb_next;
2087 *(i++) = I->NBHash[hash];
2088 *(i++) = hash;
2089 *(i++) = b0;
2090 I->NBHash[hash] = nb_next;
2091 nb_next += 3;
2092 }
2093
2094 /* find neighbors for each atom */
2095 if((cSculptVDW | cSculptVDW14) & mask) {
2096 for(aa = 0; aa < n_active; aa++) {
2097 b0 = active[aa];
2098 a0 = atm2idx[b0];
2099 ai0 = obj->AtomInfo + b0;
2100 v0 = cs_coord + 3 * a0;
2101 don_b0 = I->Don[b0];
2102 acc_b0 = I->Acc[b0];
2103 v0i = (int) (*v0);
2104 v1i = (int) (*(v0 + 1));
2105 v2i = (int) (*(v0 + 2));
2106 x0i = ex_hash_i0(b0);
2107 for(h = -4; h < 5; h += 4) {
2108 nb_off0 = nb_hash_off_i0(v0i, h);
2109 for(k = -4; k < 5; k += 4) {
2110 nb_off1 = nb_off0 | nb_hash_off_i1(v1i, k);
2111 for(l = -4; l < 5; l += 4) {
2112 /* offset = *(I->NBHash+nb_hash_off(v0,h,k,l)); */
2113 offset = I->NBHash[nb_off1 | nb_hash_off_i2(v2i, l)];
2114 while(offset) {
2115 i = I->NBList + offset;
2116 b1 = *(i + 2);
2117 if(b1 > b0) {
2118 /* determine exclusion (if any) */
2119 {
2120 int xoffset;
2121 const int *I_EXList = I->EXList.data();
2122 int ex1;
2123 const int *j;
2124 xoffset = I->EXHash[x0i | ex_hash_i1(b1)];
2125 ex = 10;
2126 while(xoffset) {
2127 xoffset = (*(j = I_EXList + xoffset));
2128 if((*(j + 1) == b0) && (*(j + 2) == b1)) {
2129 ex1 = *(j + 3);
2130 if(ex1 < ex) {
2131 ex = ex1;
2132 }
2133 }
2134 }
2135 }
2136 if(ex > 3) {
2137 ai1 = obj->AtomInfo + b1;
2138 cutoff = ai0->vdw + ai1->vdw;
2139
2140 if(ex == 10) { /* standard interaction -- no exclusion */
2141 if(don_b0 && I->Acc[b1]) { /* h-bond */
2142 if(ai0->protons == cAN_H) {
2143 cutoff -= hb_overlap;
2144 } else {
2145 cutoff -= hb_overlap_base;
2146 }
2147 } else if(acc_b0 && I->Don[b1]) { /* h-bond */
2148 if(ai1->protons == cAN_H) {
2149 cutoff -= hb_overlap;
2150 } else {
2151 cutoff -= hb_overlap_base;
2152 }
2153 }
2154 if(cSculptVDW & mask) {
2155 vdw_cutoff = cutoff * vdw;
2156 wt = vdw_wt * vdw_magnified;
2157 a1 = atm2idx[b1];
2158 v1 = cs_coord + 3 * a1;
2159 if(vdw_vis_mode && cgo && (n_cycle < 1)
2160 && ((!((ai0->protekted && ai1->protekted)
2161 || (ai0->flags & ai1->flags & cAtomFlag_fix))
2162 ) || (ai0->flags & cAtomFlag_study)
2163 || (ai1->flags & cAtomFlag_study))) {
2164 SculptCGOBump(v0, v1, ai0->vdw, ai1->vdw, cutoff,
2165 vdw_vis_min, vdw_vis_mid, vdw_vis_max,
2166 good_color, bad_color, vdw_vis_mode, cgo);
2167 }
2168 if(SculptCheckBump(v0, v1, diff, &len, vdw_cutoff))
2169 if(SculptDoBump(vdw_cutoff, len, diff,
2170 disp + b0 * 3, disp + b1 * 3, wt,
2171 &total_strain)) {
2172 cnt[b0]++;
2173 cnt[b1]++;
2174 total_count++;
2175 }
2176 }
2177 } else if(ex == 4) { /* 1-4 interation */
2178 cutoff *= vdw14;
2179 wt = vdw_wt14 * vdw_magnified;
2180
2181 if(cSculptVDW14 & mask) {
2182 a1 = atm2idx[b1];
2183 v1 = cs_coord + 3 * a1;
2184 if(SculptCheckBump(v0, v1, diff, &len, cutoff)) {
2185 if(SculptDoBump(cutoff, len, diff,
2186 disp + b0 * 3, disp + b1 * 3, wt,
2187 &total_strain)) {
2188 cnt[b0]++;
2189 cnt[b1]++;
2190 total_count++;
2191 }
2192 }
2193 }
2194 } else if(ex == 5) {
2195 /* do nothing */
2196 }
2197 }
2198 }
2199 offset = (*i);
2200 }
2201 }
2202 }
2203 }
2204 }
2205 }
2206
2207 if(cSculptAvoid & mask) {
2208 float target;
2209 float range = solvent_radius * 0.75;
2210 /* tweak nb distances to avoid
2211 sitting in the surface
2212 rendition danger zone for too
2213 long (vdw1+vdw2+0.75*solvent) */
2214 for(aa = 0; aa < n_active; aa++) {
2215 b0 = active[aa];
2216 a0 = atm2idx[b0];
2217 ai0 = obj->AtomInfo + b0;
2218 v0 = cs_coord + 3 * a0;
2219 don_b0 = I->Don[b0];
2220 acc_b0 = I->Acc[b0];
2221 v0i = (int) (*v0);
2222 v1i = (int) (*(v0 + 1));
2223 v2i = (int) (*(v0 + 2));
2224 x0i = ex_hash_i0(b0);
2225 for(h = -8; h < 9; h += 4) {
2226 nb_off0 = nb_hash_off_i0(v0i, h);
2227 for(k = -8; k < 9; k += 4) {
2228 nb_off1 = nb_off0 | nb_hash_off_i1(v1i, k);
2229 for(l = -8; l < 9; l += 4) {
2230 /* offset = *(I->NBHash+nb_hash_off(v0,h,k,l)); */
2231 offset = I->NBHash[nb_off1 | nb_hash_off_i2(v2i, l)];
2232 while(offset) {
2233 i = I->NBList + offset;
2234 b1 = *(i + 2);
2235 if(b1 > b0) {
2236 /* determine exclusion (if any) */
2237 {
2238 int xoffset;
2239 const int *I_EXList = I->EXList.data();
2240 int ex1;
2241 const int *j;
2242 xoffset = I->EXHash[x0i | ex_hash_i1(b1)];
2243 ex = 10;
2244 while(xoffset) {
2245 xoffset = (*(j = I_EXList + xoffset));
2246 if((*(j + 1) == b0) && (*(j + 2) == b1)) {
2247 ex1 = *(j + 3);
2248 if(ex1 < ex) {
2249 ex = ex1;
2250 }
2251 }
2252 }
2253 }
2254 if(ex > avd_ex) { /* either non-covalent or extended chain */
2255 ai1 = obj->AtomInfo + b1;
2256 target = ai0->vdw + ai1->vdw + avd_gp;
2257 a1 = atm2idx[b1];
2258 v1 = cs_coord + 3 * a1;
2259
2260 if(SculptCheckAvoid(v0, v1, diff, &len, target, avd_rg)) {
2261 if(SculptDoAvoid(target, range, len, diff,
2262 disp + b0 * 3, disp + b1 * 3, avd_wt,
2263 &total_strain)) {
2264 cnt[b0]++;
2265 cnt[b1]++;
2266 total_count++;
2267 }
2268 }
2269 }
2270 }
2271 offset = (*i);
2272 }
2273 }
2274 }
2275 }
2276 }
2277 }
2278
2279 /* clean up nonbonded hash */
2280
2281 i = I->NBList + 2;
2282 while(nb_next > 1) {
2283 I->NBHash[*i] = 0;
2284 i += 3;
2285 nb_next -= 3;
2286 }
2287 }
2288 }
2289 /* average the displacements */
2290
2291 if(n_cycle >= 0) {
2292 int cnt_a,a;
2293 float _1 = 1.0F;
2294 float inv_cnt;
2295 int *a_ptr = active;
2296 float *lookup_inverse = I->inverse;
2297 for(aa = 0; aa < n_active; aa++) {
2298 if((cnt_a = cnt[(a = *(a_ptr++))])) {
2299 AtomInfoType *ai = obj->AtomInfo + a;
2300 const RefPosType *cs_refpos = cs->RefPos.data();
2301 int flags;
2302 if(!(ai->protekted || ((flags = ai->flags) & cAtomFlag_fix))) {
2303 v1 = disp + 3 * a;
2304 v2 = cs_coord + 3 * atm2idx[a];
2305
2306 if((flags & cAtomFlag_restrain) && cs_refpos) {
2307 const RefPosType *rp = cs_refpos + atm2idx[a];
2308 if(rp->specified) {
2309 const float *v3 = rp->coord;
2310 cnt_a++;
2311 v1[0] += v3[0] - v2[0];
2312 v1[1] += v3[1] - v2[1];
2313 v1[2] += v3[2] - v2[2];
2314 }
2315 }
2316
2317 if(!(cnt_a & 0xFFFFFF00)) /* don't divide -- too slow! */
2318 inv_cnt = lookup_inverse[cnt_a];
2319 else
2320 inv_cnt = _1 / cnt_a;
2321
2322 *(v2) += (*(v1)) * inv_cnt;
2323 *(v2 + 1) += (*(v1 + 1)) * inv_cnt;
2324 *(v2 + 2) += (*(v1 + 2)) * inv_cnt;
2325 }
2326 }
2327 }
2328 cs->invalidateRep(cRepAll, cRepInvCoord);
2329 } else if(cgo) {
2330 SceneDirty(G);
2331 }
2332 if(n_cycle <= 0) {
2333 int *a_ptr = active;
2334 if(center)
2335 for(aa = 0; aa < n_active; aa++) {
2336 int a = *(a_ptr++);
2337 {
2338 AtomInfoType *ai = obj->AtomInfo + a;
2339 if((ai->protekted != 1) && !(ai->flags & cAtomFlag_fix)) {
2340 v2 = cs_coord + 3 * atm2idx[a];
2341 center[0] += *(v2);
2342 center[1] += *(v2 + 1);
2343 center[2] += *(v2 + 2);
2344 center[3] += 1.0F;
2345 }
2346 }
2347 }
2348 break;
2349 }
2350 }
2351
2352 task_time = UtilGetSeconds(G) - task_time;
2353 PRINTFB(G, FB_Sculpt, FB_Blather)
2354 " Sculpt: %2.5f seconds %8.3f %d %8.3f\n", task_time, total_strain, total_count,
2355 100 * total_strain / total_count ENDFB(G);
2356
2357 if(total_count)
2358 total_strain = (1000 * total_strain) / total_count;
2359 }
2360 FreeP(exclude);
2361 FreeP(active);
2362 FreeP(cnt);
2363 FreeP(disp);
2364 FreeP(atm2idx);
2365 if(cgo) {
2366 CGOStop(cgo);
2367 {
2368 int est = CGOCheckComplex(cgo);
2369 if(est) {
2370 cs->SculptCGO = CGOSimplify(cgo, est);
2371 CGOFree(cgo);
2372 CGOFree(cs->SculptShaderCGO);
2373 }
2374 }
2375 }
2376
2377 EditorDihedralInvalid(G, obj);
2378 }
2379
2380 PRINTFD(G, FB_Sculpt)
2381 " SculptIterateObject-Debug: leaving...\n" ENDFD;
2382
2383 return total_strain;
2384 }
2385