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 <utility>
18
19 #include"Version.h"
20 #include"os_python.h"
21
22 #include"os_predef.h"
23 #include"os_std.h"
24 #include"os_gl.h"
25
26 #include"Base.h"
27 #include"Parse.h"
28 #include"OOMac.h"
29 #include"Vector.h"
30 #include"PConv.h"
31 #include"ObjectMolecule.h"
32 #include"Feedback.h"
33 #include"Util.h"
34 #include"Util2.h"
35 #include"AtomInfo.h"
36 #include"Selector.h"
37 #include"ObjectDist.h"
38 #include"Executive.h"
39 #include"P.h"
40 #include"ObjectCGO.h"
41 #include"Scene.h"
42 #include "Lex.h"
43
44 #include"AtomInfoHistory.h"
45 #include"BondTypeHistory.h"
46
47 #ifdef _PYMOL_IP_PROPERTIES
48 #include"Property.h"
49 #endif
50
51 #include <functional>
52 #include <iostream>
53 #include <map>
54
55 #define ntrim ParseNTrim
56 #define nextline ParseNextLine
57 #define ncopy ParseNCopy
58 #define nskip ParseNSkip
59
60 #define cResvMask 0x7FFF
61
62 #define cMaxOther 6
63
64 #define strstartswith p_strstartswith
65
strstartswithword(const char * s,const char * word)66 static int strstartswithword(const char * s, const char * word) {
67 while(*word)
68 if(*s++ != *word++)
69 return false;
70 switch(*s) {
71 case ' ':
72 case '\t':
73 case '\r':
74 case '\n':
75 case '\0':
76 return true;
77 }
78 return false;
79 }
80
AddOrthoOutputIfMatchesTags(PyMOLGlobals * G,int n_tags,int nAtom,const char * const * tag_start,const char * p,char * cc,int quiet)81 static void AddOrthoOutputIfMatchesTags(PyMOLGlobals * G, int n_tags,
82 int nAtom, const char* const* tag_start, const char *p, char *cc, int quiet)
83 {
84 if(n_tags && !quiet && !(nAtom > 0 && strstartswith(p, "HEADER"))) {
85 // HEADER is the mark for a new object, this is a new object, and
86 // gets processed on the next pass, when nAtom=0
87 int tc = 0;
88 for(; tc < n_tags; tc++) {
89 if(!strstartswithword(p, tag_start[tc]))
90 continue;
91 ParseNTrimRight(cc, p, MAXLINELEN - 1);
92 OrthoAddOutput(G, cc);
93 OrthoNewLine(G, NULL, true);
94 break;
95 }
96 }
97 }
98
99 typedef struct {
100 int n_cyclic_arom, cyclic_arom[cMaxOther];
101 int n_arom, arom[cMaxOther];
102 int n_high_val, high_val[cMaxOther];
103 int n_cyclic, cyclic[cMaxOther];
104 int n_planer, planer[cMaxOther];
105 int n_rest, rest[cMaxOther];
106 int score;
107 } OtherRec;
108
populate_other(OtherRec * other,int at,const AtomInfoType * ai,const BondType * bd,const int * neighbor)109 static int populate_other(OtherRec * other, int at,
110 const AtomInfoType* ai,
111 const BondType* bd,
112 const int* neighbor)
113 {
114 int five_cycle = false;
115 int six_cycle = false;
116
117 {
118 int mem[9], nbr[7];
119 const int ESCAPE_MAX = 500;
120 int escape_count;
121
122 escape_count = ESCAPE_MAX; /* don't get bogged down with structures
123 that have unreasonable connectivity */
124 mem[0] = bd->index[0];
125 mem[1] = bd->index[1];
126 nbr[1] = neighbor[mem[1]] + 1;
127 while(((mem[2] = neighbor[nbr[1]]) >= 0)) {
128 if(mem[2] != mem[0]) {
129 nbr[2] = neighbor[mem[2]] + 1;
130 while(((mem[3] = neighbor[nbr[2]]) >= 0)) {
131 if(mem[3] != mem[1]) {
132 nbr[3] = neighbor[mem[3]] + 1;
133 while(((mem[4] = neighbor[nbr[3]]) >= 0)) {
134 if((mem[4] != mem[2]) && (mem[4] != mem[1]) && (mem[4] != mem[0])) {
135 nbr[4] = neighbor[mem[4]] + 1;
136 while(((mem[5] = neighbor[nbr[4]]) >= 0)) {
137 if(!(escape_count--))
138 goto escape;
139 if((mem[5] != mem[3]) && (mem[5] != mem[2]) && (mem[5] != mem[1])) {
140 if(mem[5] == mem[0]) { /* five-cycle */
141 five_cycle = true;
142 }
143 nbr[5] = neighbor[mem[5]] + 1;
144 while(((mem[6] = neighbor[nbr[5]]) >= 0)) {
145 if((mem[6] != mem[4]) && (mem[6] != mem[3]) && (mem[6] != mem[2])
146 && (mem[6] != mem[1])) {
147 if(mem[6] == mem[0]) { /* six-cycle */
148 six_cycle = true;
149 }
150 }
151 nbr[5] += 2;
152 }
153 }
154 nbr[4] += 2;
155 }
156 }
157 nbr[3] += 2;
158 }
159 }
160 nbr[2] += 2;
161 }
162 }
163 nbr[1] += 2;
164 }
165 }
166 escape:
167
168 if(bd->order == 4) { /* aromatic */
169 if((five_cycle || six_cycle) && (other->n_cyclic_arom < cMaxOther)) {
170 other->cyclic_arom[other->n_cyclic_arom++] = at;
171 if(five_cycle && six_cycle)
172 other->score += 34;
173 else if(five_cycle)
174 other->score += 33;
175 else
176 other->score += 32;
177 return 1;
178 } else if(other->n_arom < cMaxOther) {
179 other->arom[other->n_arom++] = at;
180 other->score += 64;
181 return 1;
182 }
183 }
184 if(bd->order > 1) {
185 if(other->n_high_val < cMaxOther) {
186 other->high_val[other->n_high_val++] = at;
187 other->score += 16;
188 return 1;
189 }
190 }
191 if(five_cycle || six_cycle) {
192 if(other->n_cyclic < cMaxOther) {
193 other->cyclic[other->n_cyclic++] = at;
194 other->score += 8;
195 return 1;
196 }
197 }
198 if(ai->geom == cAtomInfoPlanar) {
199 if(other->n_planer < cMaxOther) {
200 other->planer[other->n_planer++] = at;
201 other->score += 4;
202 return 1;
203 }
204 }
205 if(other->n_rest < cMaxOther) {
206 other->rest[other->n_rest++] = at;
207 other->score += 1;
208 return 1;
209 }
210 return 0;
211 }
212
append_index(int * result,int offset,int a1,int a2,int score,int ar_count)213 static int append_index(int *result, int offset, int a1, int a2, int score, int ar_count)
214 {
215 int c;
216 c = result[a1];
217 while(c < offset) {
218 if(result[c] == a2) { /* already entered */
219 if(result[c + 1] < score) {
220 result[c + 1] = score;
221 result[c + 2] = ar_count;
222 }
223 return offset;
224 }
225 c += 3;
226 }
227 result[offset++] = a2;
228 result[offset++] = score;
229 result[offset++] = ar_count;
230 return offset;
231 }
232
ObjectMoleculeAddPseudoatom(ObjectMolecule * I,int sele_index,const char * name,const char * resn,const char * resi,const char * chain,const char * segi,const char * elem,float vdw,int hetatm,float b,float q,const char * label,const float * pos,int color,int state,int mode,int quiet)233 int ObjectMoleculeAddPseudoatom(ObjectMolecule * I, int sele_index, const char *name,
234 const char *resn, const char *resi, const char *chain,
235 const char *segi, const char *elem, float vdw,
236 int hetatm, float b, float q, const char *label,
237 const float *pos, int color, int state, int mode, int quiet)
238 {
239 PyMOLGlobals *G = I->G;
240 int start_state = 0, stop_state = 0;
241 int extant_only = false;
242 int ai_merged = false;
243 float pos_array[3] = { 0.0F, 0.0F, 0.0F };
244 int ok = true;
245
246 pymol::vla<AtomInfoType> atInfo(1);
247 AtomInfoType* ai = atInfo.data();
248
249 #ifdef _PYMOL_IP_EXTRAS
250 ai->oldid = -1;
251 #endif
252
253 // FIXME this should be -2
254 if (state == -1) {
255 state = I->getCurrentState();
256 }
257
258 if(state >= 0) { /* specific state */
259 start_state = state;
260 stop_state = state + 1;
261 } else { /* all states */
262 if(sele_index >= 0) {
263 start_state = 0;
264 stop_state = SelectorCountStates(G, sele_index);
265 if(state == -3)
266 extant_only = true;
267 } else {
268 start_state = 0;
269 stop_state = ExecutiveCountStates(G, cKeywordAll);
270 if(stop_state < 1)
271 stop_state = 1;
272 }
273 }
274 {
275 /* match existing properties of the old atom */
276 ai->setResi(resi);
277 ai->hetatm = hetatm;
278 ai->geom = cAtomInfoNone;
279 ai->q = q;
280 ai->b = b;
281 ai->chain = LexIdx(G, chain);
282 ai->segi = LexIdx(G, segi);
283 ai->resn = LexIdx(G, resn);
284 ai->name = LexIdx(G, name);
285 strcpy(ai->elem, elem);
286 ai->id = -1;
287 ai->rank = -1;
288 if(vdw >= 0.0F)
289 ai->vdw = vdw;
290 else
291 ai->vdw = 1.0F;
292 if(label[0]) {
293 ai->label = LexIdx(G, label);
294 ai->visRep = cRepLabelBit;
295 } else {
296 ai->visRep = RepGetAutoShowMask(G);
297 }
298
299 ai->flags |= cAtomFlag_inorganic; // suppress auto_show_classified
300
301 if(color < 0) {
302 AtomInfoAssignColors(I->G, ai);
303 if((ai->elem[0] == 'C') && (ai->elem[1] == 0))
304 /* carbons are always colored according to the object color */
305 ai->color = I->Color;
306 } else {
307 ai->color = color;
308 }
309 AtomInfoAssignParameters(I->G, ai);
310 AtomInfoUniquefyNames(I->G, I->AtomInfo, I->NAtom, ai, NULL, 1);
311 if(!quiet) {
312 PRINTFB(G, FB_ObjectMolecule, FB_Actions)
313 " ObjMol: created %s/%s/%s/%s`%d%c/%s\n",
314 I->Name, LexStr(G, ai->segi), LexStr(G, ai->chain),
315 LexStr(G, ai->resn), ai->resv, ai->getInscode(true),
316 LexStr(G, ai->name) ENDFB(G);
317 }
318 }
319
320 CoordSet *cset = CoordSetNew(G);
321 cset->NIndex = 1;
322 cset->Coord = pymol::vla<float>(3);
323 cset->Obj = I;
324 cset->enumIndices();
325
326 for(state = start_state; state < stop_state; state++) {
327
328
329 if((extant_only && (state < I->NCSet) && I->CSet[state]) || !extant_only) {
330
331 if(sele_index >= 0) {
332 ObjectMoleculeOpRec op;
333 ObjectMoleculeOpRecInit(&op);
334 op.code = OMOP_CSetSumVertices;
335 op.cs1 = state;
336
337 ExecutiveObjMolSeleOp(I->G, sele_index, &op);
338
339 if(op.i1) {
340 float factor = 1.0F / op.i1;
341 scale3f(op.v1, factor, pos_array);
342 pos = pos_array;
343
344 if(vdw < 0.0F) {
345 switch (mode) {
346 case 1:
347 ObjectMoleculeOpRecInit(&op);
348 op.code = OMOP_CSetMaxDistToPt;
349 copy3f(pos_array, op.v1);
350 op.cs1 = state;
351 ExecutiveObjMolSeleOp(I->G, sele_index, &op);
352 vdw = op.f1;
353 break;
354 case 2:
355 ObjectMoleculeOpRecInit(&op);
356 op.code = OMOP_CSetSumSqDistToPt;
357 copy3f(pos_array, op.v1);
358 op.cs1 = state;
359 ExecutiveObjMolSeleOp(I->G, sele_index, &op);
360 vdw = sqrt1f(op.d1 / op.i1);
361 break;
362 case 0:
363 default:
364 vdw = 0.5F;
365 break;
366 }
367 if(vdw >= 0.0F)
368 ai->vdw = vdw; /* NOTE: only uses vdw from first state selection... */
369 }
370 } else {
371 pos = NULL; /* skip this state */
372 }
373 } else if(!pos) {
374 SceneGetCenter(I->G, pos_array);
375 pos = pos_array;
376 }
377
378 if(pos) { /* only add coordinate to state if we have position for it */
379
380 float *coord = cset->Coord.data();
381
382 copy3f(pos, coord);
383
384 if(!ai_merged) {
385 if (ok)
386 ok &= ObjectMoleculeMerge(I, std::move(atInfo), cset, false, cAIC_AllMask, true); /* NOTE: will release atInfo */
387 if (ok)
388 ok &= ObjectMoleculeExtendIndices(I, -1);
389 if (ok)
390 ok &= ObjectMoleculeUpdateNeighbors(I);
391 ai_merged = true;
392 }
393 if(state >= I->NCSet) {
394 VLACheck(I->CSet, CoordSet *, state);
395 I->NCSet = state + 1;
396 }
397 if(!I->CSet[state]) {
398 /* new coordinate set */
399 I->CSet[state] = CoordSetCopy(cset);
400 } else {
401 /* merge coordinate set */
402 if (ok)
403 ok &= CoordSetMerge(I, I->CSet[state], cset);
404 }
405 }
406 }
407 }
408
409 cset->fFree();
410
411 if(ai_merged) {
412 if (ok)
413 ok &= ObjectMoleculeSort(I);
414 ObjectMoleculeUpdateIDNumbers(I);
415 ObjectMoleculeUpdateNonbonded(I);
416 I->invalidate(cRepAll, cRepInvAtoms, -1);
417 } else {
418 VLAFreeP(atInfo);
419 }
420 return (ok);
421 }
422
ObjectMoleculeGetPrioritizedOtherIndexList(ObjectMolecule * I,CoordSet * cs)423 int *ObjectMoleculeGetPrioritizedOtherIndexList(ObjectMolecule * I, CoordSet * cs)
424 {
425 int a, b;
426 int b1, b2, a1, a2, a3;
427 OtherRec *o;
428 OtherRec *other = pymol::calloc<OtherRec>(cs->NIndex);
429 int *result = NULL;
430 int offset;
431 int n_alloc = 0;
432 const BondType *bd;
433 int ok = true;
434
435 CHECKOK(ok, other);
436
437 if (ok){
438 ok &= ObjectMoleculeUpdateNeighbors(I);
439 }
440 bd = I->Bond;
441 for(a = 0; ok && a < I->NBond; a++) {
442 b1 = bd->index[0];
443 b2 = bd->index[1];
444 a1 = cs->atmToIdx(b1);
445 a2 = cs->atmToIdx(b2);
446 if((a1 >= 0) && (a2 >= 0)) {
447 n_alloc += populate_other(other + a1, a2, I->AtomInfo + b2, bd, I->Neighbor);
448 n_alloc += populate_other(other + a2, a1, I->AtomInfo + b1, bd, I->Neighbor);
449 }
450 bd++;
451 ok &= !I->G->Interrupt;
452 }
453 if (ok){
454 n_alloc = 3 * (n_alloc + cs->NIndex);
455 o = other;
456 result = pymol::malloc<int>(n_alloc);
457 CHECKOK(ok, result);
458 }
459 if (ok){
460 for(a = 0; a < cs->NIndex; a++) {
461 result[a] = -1;
462 }
463 offset = cs->NIndex;
464 bd = I->Bond;
465 }
466 for(a = 0; ok && a < I->NBond; a++) {
467 b1 = bd->index[0];
468 b2 = bd->index[1];
469 a1 = cs->atmToIdx(b1);
470 a2 = cs->atmToIdx(b2);
471 if((a1 >= 0) && (a2 >= 0)) {
472 if(result[a1] < 0) {
473 o = other + a1;
474 result[a1] = offset;
475 for(b = 0; b < o->n_cyclic_arom; b++) {
476 a3 = o->cyclic_arom[b];
477 offset = append_index(result, offset, a1, a3, 128 + other[a3].score, 1);
478 }
479 for(b = 0; b < o->n_arom; b++) {
480 a3 = o->arom[b];
481 offset = append_index(result, offset, a1, a3, 64 + other[a3].score, 1);
482 }
483 for(b = 0; b < o->n_high_val; b++) {
484 a3 = o->high_val[b];
485 offset = append_index(result, offset, a1, a3, 16 + other[a3].score, 0);
486 }
487 for(b = 0; b < o->n_cyclic; b++) {
488 a3 = o->cyclic[b];
489 offset = append_index(result, offset, a1, a3, 8 + other[a3].score, 0);
490 }
491 for(b = 0; b < o->n_planer; b++) {
492 a3 = o->planer[b];
493 offset = append_index(result, offset, a1, a3, 2 + other[a3].score, 0);
494 }
495 for(b = 0; b < o->n_rest; b++) {
496 a3 = o->rest[b];
497 offset = append_index(result, offset, a1, a3, 1 + other[a3].score, 0);
498 }
499 result[offset++] = -1;
500 }
501
502 if(result[a2] < 0) {
503 o = other + a2;
504 result[a2] = offset;
505 for(b = 0; b < o->n_cyclic_arom; b++) {
506 a3 = o->cyclic_arom[b];
507 offset = append_index(result, offset, a2, a3, 128 + other[a3].score, 1);
508 }
509 for(b = 0; b < o->n_arom; b++) {
510 a3 = o->arom[b];
511 offset = append_index(result, offset, a2, a3, 64 + other[a3].score, 1);
512 }
513 for(b = 0; b < o->n_high_val; b++) {
514 a3 = o->high_val[b];
515 offset = append_index(result, offset, a2, a3, 16 + other[a3].score, 0);
516 }
517 for(b = 0; b < o->n_cyclic; b++) {
518 a3 = o->cyclic[b];
519 offset = append_index(result, offset, a2, a3, 8 + other[a3].score, 0);
520 }
521 for(b = 0; b < o->n_planer; b++) {
522 a3 = o->planer[b];
523 offset = append_index(result, offset, a2, a3, 2 + other[a3].score, 0);
524 }
525 for(b = 0; b < o->n_rest; b++) {
526 a3 = o->rest[b];
527 offset = append_index(result, offset, a2, a3, 1 + other[a3].score, 0);
528 }
529 result[offset++] = -1;
530 }
531
532 }
533 bd++;
534 ok &= !I->G->Interrupt;
535 }
536 FreeP(other);
537 return result;
538 }
539
ObjectMoleculeGetNearestBlendedColor(ObjectMolecule * I,const float * point,float cutoff,int state,float * dist,float * color,int sub_vdw)540 int ObjectMoleculeGetNearestBlendedColor(ObjectMolecule * I, const float *point,
541 float cutoff, int state, float *dist,
542 float *color, int sub_vdw)
543 {
544 int result = -1;
545 float tot_weight = 0.0F;
546 float cutoff2 = cutoff * cutoff;
547 float nearest = -1.0F;
548
549 color[0] = 0.0F;
550 color[1] = 0.0F;
551 color[2] = 0.0F;
552
553 assert(state != -1 /* all states */);
554 auto* cs = I->getCoordSet(state);
555
556 {
557 if(cs) {
558 MapType *map;
559 CoordSetUpdateCoord2IdxMap(cs, cutoff);
560 if(sub_vdw) {
561 cutoff -= MAX_VDW;
562 cutoff2 = cutoff * cutoff;
563 }
564 nearest = cutoff2;
565 if((map = cs->Coord2Idx)) {
566 int a, b, c, d, e, f, j;
567 float test;
568 float *v;
569 MapLocus(map, point, &a, &b, &c);
570 for(d = a - 1; d <= a + 1; d++)
571 for(e = b - 1; e <= b + 1; e++)
572 for(f = c - 1; f <= c + 1; f++) {
573 j = *(MapFirst(map, d, e, f));
574 while(j >= 0) {
575 v = cs->coordPtr(j);
576 test = diffsq3f(v, point);
577 if(sub_vdw) {
578 test = sqrt1f(test);
579 test -= I->AtomInfo[cs->IdxToAtm[j]].vdw;
580 if(test < 0.0F)
581 test = 0.0F;
582 test = test * test;
583 }
584 if(test < cutoff2) {
585 float weight = cutoff - sqrt1f(test);
586 const float *at_col = ColorGet(I->G, I->AtomInfo[cs->IdxToAtm[j]].color);
587 color[0] += at_col[0] * weight;
588 color[1] += at_col[1] * weight;
589 color[2] += at_col[2] * weight;
590 tot_weight += weight;
591 }
592 if(test <= nearest) {
593 result = j;
594 nearest = test;
595 }
596 j = MapNext(map, j);
597 }
598 }
599 } else {
600 int j;
601 float test;
602 const float* v = cs->Coord.data();
603 for(j = 0; j < cs->NIndex; j++) {
604 test = diffsq3f(v, point);
605 if(sub_vdw) {
606 test = sqrt1f(test);
607 test -= I->AtomInfo[cs->IdxToAtm[j]].vdw;
608 if(test < 0.0F)
609 test = 0.0F;
610 test = test * test;
611 }
612 if(test < cutoff2) {
613 float weight = cutoff - sqrt1f(test);
614 const float *at_col = ColorGet(I->G, I->AtomInfo[cs->IdxToAtm[j]].color);
615 color[0] += at_col[0] * weight;
616 color[1] += at_col[1] * weight;
617 color[2] += at_col[2] * weight;
618 tot_weight += weight;
619 }
620 if(test <= nearest) {
621 result = j;
622 nearest = test;
623 }
624 v += 3;
625 }
626 }
627 if(result >= 0)
628 result = cs->IdxToAtm[result];
629 }
630 }
631 if(dist) {
632 if(result >= 0) {
633 *dist = sqrt1f(nearest);
634 if(tot_weight > 0.0F) {
635 color[0] /= tot_weight;
636 color[1] /= tot_weight;
637 color[2] /= tot_weight;
638 }
639 } else {
640 *dist = -1.0F;
641 }
642 }
643 return result;
644 }
645
ObjectMoleculeGetNearestAtomIndex(ObjectMolecule * I,const float * point,float cutoff,int state,float * dist)646 int ObjectMoleculeGetNearestAtomIndex(ObjectMolecule * I, const float *point, float cutoff,
647 int state, float *dist)
648 {
649 int result = -1;
650 float nearest = -1.0F;
651
652 assert(state != -1 /* all states */);
653 auto* cs = I->getCoordSet(state);
654
655 {
656 if(cs) {
657 MapType *map;
658 CoordSetUpdateCoord2IdxMap(cs, cutoff);
659 nearest = cutoff * cutoff;
660 if((map = cs->Coord2Idx)) {
661 int a, b, c, d, e, f, j;
662 float test;
663 float *v;
664 MapLocus(map, point, &a, &b, &c);
665 for(d = a - 1; d <= a + 1; d++)
666 for(e = b - 1; e <= b + 1; e++)
667 for(f = c - 1; f <= c + 1; f++) {
668 j = *(MapFirst(map, d, e, f));
669 while(j >= 0) {
670 v = cs->coordPtr(j);
671 test = diffsq3f(v, point);
672 if(test <= nearest) {
673 result = j;
674 nearest = test;
675 }
676 j = MapNext(map, j);
677 }
678 }
679 } else {
680 int j;
681 float test;
682 const float* v = cs->Coord.data();
683 for(j = 0; j < cs->NIndex; j++) {
684 test = diffsq3f(v, point);
685 if(test <= nearest) {
686 result = j;
687 nearest = test;
688 }
689 v += 3;
690 }
691 }
692 if(result >= 0)
693 result = cs->IdxToAtm[result];
694 }
695 }
696 if(dist) {
697 if(result >= 0) {
698 *dist = sqrt1f(nearest);
699 } else {
700 *dist = -1.0F;
701 }
702 }
703 return result;
704 }
705
ObjectMoleculeGetPrioritizedOther(const int * other,int a1,int a2,int * double_sided)706 int ObjectMoleculeGetPrioritizedOther(const int *other, int a1, int a2, int *double_sided)
707 {
708 int a3 = -1;
709 int lvl = -1, ck, ck_lvl;
710 int offset;
711 int ar_count = 0;
712
713 a3 = -1;
714 lvl = -1;
715 if(a1 >= 0) {
716 offset = other[a1];
717 if(offset >= 0) {
718 while(1) {
719 ck = other[offset];
720 if(ck != a2) {
721 if(ck >= 0) {
722 ck_lvl = other[offset + 1];
723 if(ck_lvl > lvl) {
724 a3 = ck;
725 lvl = ck_lvl;
726 }
727 ar_count += other[offset + 2];
728 } else
729 break;
730 }
731 offset += 3;
732 }
733 }
734 }
735 if(a2 >= 0) {
736 offset = other[a2];
737 if(offset >= 0) {
738 while(1) {
739 ck = other[offset];
740 if(ck != a1) {
741 if(ck >= 0) {
742 ck_lvl = other[offset + 1];
743 if(ck_lvl > lvl) {
744 a3 = ck;
745 lvl = ck_lvl;
746 }
747 ar_count += other[offset + 2];
748 } else
749 break;
750 }
751 offset += 3;
752 }
753 }
754 }
755
756 if(double_sided) {
757 if(ar_count == 4)
758 *double_sided = true;
759 else
760 *double_sided = false;
761 }
762 return a3;
763 }
764
765 /* Check if atom is bonded to an atom with given name
766 *
767 * param same_res:
768 * 0 = must not be in same residue
769 * 1 = must be in same residue
770 * -1 = don't check residue
771 */
ObjectMoleculeIsAtomBondedToName(ObjectMolecule * obj,int a0,const char * name,int same_res)772 int ObjectMoleculeIsAtomBondedToName(ObjectMolecule * obj, int a0, const char *name, int same_res)
773 {
774 int a2, s;
775 PyMOLGlobals * G = obj->G;
776 AtomInfoType *ai2, *ai0 = obj->AtomInfo + a0;
777
778 if(a0 >= 0) {
779 ITERNEIGHBORATOMS(obj->Neighbor, a0, a2, s) {
780 ai2 = obj->AtomInfo + a2;
781 if(WordMatchExact(G, LexStr(G, ai2->name), name, true) &&
782 (same_res < 0 || (same_res == AtomInfoSameResidue(G, ai0, ai2))))
783 return true;
784 }
785 }
786 return false;
787 }
788
ObjectMoleculeAreAtomsBonded2(ObjectMolecule * obj0,int a0,ObjectMolecule * obj1,int a1)789 int ObjectMoleculeAreAtomsBonded2(ObjectMolecule * obj0, int a0, ObjectMolecule * obj1,
790 int a1)
791 {
792 /* assumes neighbor list is current */
793
794 if(obj0 != obj1)
795 return false;
796 else {
797 int a2, s;
798
799 if(a0 >= 0) {
800 s = obj0->Neighbor[a0];
801 s++; /* skip count */
802 while(1) {
803 a2 = obj0->Neighbor[s];
804 if(a2 < 0)
805 break;
806 if(a1 == a2)
807 return true;
808 s += 2;
809 }
810 }
811 }
812 return false;
813 }
814
ObjectMoleculeDoesAtomNeighborSele(ObjectMolecule * I,int index,int sele)815 int ObjectMoleculeDoesAtomNeighborSele(ObjectMolecule * I, int index, int sele)
816 {
817 int result = false;
818 ObjectMoleculeUpdateNeighbors(I);
819 if(index < I->NAtom) {
820 int a1;
821 int n;
822 AtomInfoType *ai;
823
824 n = I->Neighbor[index] + 1;
825 while(1) { /* look for an attached non-hydrogen as a base */
826 a1 = I->Neighbor[n];
827 n += 2;
828 if(a1 < 0)
829 break;
830 ai = I->AtomInfo + a1;
831 if(SelectorIsMember(I->G, ai->selEntry, sele)) {
832 result = true;
833 break;
834 }
835 }
836 }
837 return result;
838 }
839
840 /*
841 * Based on PDB nomenclature (resn, name), do:
842 *
843 * 1) If `ai1` or `ai2` is a known charged PDB atom, assign `formalCharge`
844 * and seth `chemFlag` to false
845 *
846 * 2) If `ai1` and `ai2` are connected by a double bond, set (*bond_order) = 2
847 */
assign_pdb_known_residue(PyMOLGlobals * G,AtomInfoType * ai1,AtomInfoType * ai2,int * bond_order)848 static void assign_pdb_known_residue(PyMOLGlobals * G, AtomInfoType * ai1,
849 AtomInfoType * ai2, int *bond_order)
850 {
851 int order = *(bond_order);
852 const char *name1 = LexStr(G, ai1->name);
853 const char *name2 = LexStr(G, ai2->name);
854 const char *resn1 = LexStr(G, ai1->resn);
855
856 /* nasty high-speed hack to get bond valences and formal charges
857 for standard residues */
858 if(((!name1[1]) && (!name2[1])) &&
859 (((name1[0] == 'C') && (name2[0] == 'O')) ||
860 ((name1[0] == 'O') && (name2[0] == 'C')))) {
861 order = 2;
862 } else if((!name2[1]) && (name2[0] == 'C') && (!strcmp(name1, "OXT"))) {
863 ai1->formalCharge = -1;
864 ai1->chemFlag = false;
865 } else if((!name1[1]) && (name1[0] == 'C') && (!strcmp(name2, "OXT"))) {
866 ai2->formalCharge = -1;
867 ai2->chemFlag = false;
868 } else {
869 switch (resn1[0]) {
870 case 'A':
871 switch (resn1[1]) {
872 case 'R':
873 switch (resn1[2]) {
874 case 'G': /* ARG... */
875 switch (resn1[3]) {
876 case 0:
877 case 'P': /* ARG, ARGP */
878 if(!strcmp(name1, "NH1")) {
879 ai1->formalCharge = 1;
880 ai1->chemFlag = false;
881 } else if(!strcmp(name2, "NH1")) {
882 ai2->formalCharge = 1;
883 ai2->chemFlag = false;
884 }
885 break;
886 }
887 if(((!strcmp(name1, "CZ")) && (!strcmp(name2, "NH1"))) ||
888 ((!strcmp(name2, "CZ")) && (!strcmp(name1, "NH1"))))
889 order = 2;
890 break;
891 }
892 break;
893 case 'S':
894 switch (resn1[2]) {
895 case 'P': /* ASP... */
896 switch (resn1[3]) {
897 case 0:
898 case 'M': /* ASP, ASPM minus assumption */
899 if(!strcmp(name1, "OD2")) {
900 ai1->formalCharge = -1;
901 ai1->chemFlag = false;
902 } else if(!strcmp(name2, "OD2")) {
903 ai2->formalCharge = -1;
904 ai2->chemFlag = false;
905 }
906 break;
907 }
908 if(((!strcmp(name1, "CG")) && (!strcmp(name2, "OD1"))) ||
909 ((!strcmp(name2, "CG")) && (!strcmp(name1, "OD1"))))
910 order = 2;
911 break;
912 case 'N': /* ASN */
913 if(((!strcmp(name1, "CG")) && (!strcmp(name2, "OD1"))) ||
914 ((!strcmp(name2, "CG")) && (!strcmp(name1, "OD1"))))
915 order = 2;
916 break;
917 }
918 break;
919 case 0:
920 if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
921 ai1->formalCharge = -1;
922 ai1->chemFlag = false;
923 } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
924 ai2->formalCharge = -1;
925 ai2->chemFlag = false;
926 }
927 if(((!strcmp(name1, "C8")) && (!strcmp(name2, "N7"))) ||
928 ((!strcmp(name2, "C8")) && (!strcmp(name1, "N7"))))
929 order = 2;
930 else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "C5"))) ||
931 ((!strcmp(name2, "C4")) && (!strcmp(name1, "C5"))))
932 order = 2;
933
934 else if(((!strcmp(name1, "C6")) && (!strcmp(name2, "N1"))) ||
935 ((!strcmp(name2, "C6")) && (!strcmp(name1, "N1"))))
936 order = 2;
937 else if(((!strcmp(name1, "C2")) && (!strcmp(name2, "N3"))) ||
938 ((!strcmp(name2, "C2")) && (!strcmp(name1, "N3"))))
939 order = 2;
940 else
941 if(((!strcmp(name1, "P"))
942 && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
943 || ((!strcmp(name2, "P"))
944 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
945 order = 2;
946 break;
947 }
948 break;
949 case 'C':
950 if(resn1[1] == 0) {
951 if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
952 ai1->formalCharge = -1;
953 ai1->chemFlag = false;
954 } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
955 ai2->formalCharge = -1;
956 ai2->chemFlag = false;
957 }
958 if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
959 ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
960 order = 2;
961 else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "N3"))) ||
962 ((!strcmp(name2, "C4")) && (!strcmp(name1, "N3"))))
963 order = 2;
964
965 else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
966 ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
967 order = 2;
968 else
969 if(((!strcmp(name1, "P"))
970 && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
971 || ((!strcmp(name2, "P"))
972 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
973 order = 2;
974 }
975 break;
976 case 'D': /* deoxy nucleic acids */
977 switch (resn1[1]) {
978 case 'A':
979 if(resn1[2] == 0) {
980 if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
981 ai1->formalCharge = -1;
982 ai1->chemFlag = false;
983 } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
984 ai2->formalCharge = -1;
985 ai2->chemFlag = false;
986 }
987 if(((!strcmp(name1, "C8")) && (!strcmp(name2, "N7"))) ||
988 ((!strcmp(name2, "C8")) && (!strcmp(name1, "N7"))))
989 order = 2;
990 else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "C5"))) ||
991 ((!strcmp(name2, "C4")) && (!strcmp(name1, "C5"))))
992 order = 2;
993
994 else if(((!strcmp(name1, "C6")) && (!strcmp(name2, "N1"))) ||
995 ((!strcmp(name2, "C6")) && (!strcmp(name1, "N1"))))
996 order = 2;
997 else if(((!strcmp(name1, "C2")) && (!strcmp(name2, "N3"))) ||
998 ((!strcmp(name2, "C2")) && (!strcmp(name1, "N3"))))
999 order = 2;
1000 else
1001 if(((!strcmp(name1, "P"))
1002 && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
1003 || ((!strcmp(name2, "P"))
1004 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
1005 order = 2;
1006 }
1007 break;
1008 case 'C':
1009 if(resn1[2] == 0) {
1010 if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
1011 ai1->formalCharge = -1;
1012 ai1->chemFlag = false;
1013 } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
1014 ai2->formalCharge = -1;
1015 ai2->chemFlag = false;
1016 }
1017 if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
1018 ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
1019 order = 2;
1020 else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "N3"))) ||
1021 ((!strcmp(name2, "C4")) && (!strcmp(name1, "N3"))))
1022 order = 2;
1023
1024 else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
1025 ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
1026 order = 2;
1027 else
1028 if(((!strcmp(name1, "P"))
1029 && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
1030 || ((!strcmp(name2, "P"))
1031 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
1032 order = 2;
1033 }
1034 break;
1035 case 'T':
1036 if(resn1[2] == 0) {
1037 if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2"))))
1038 ai1->formalCharge = -1;
1039 else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2"))))
1040 ai2->formalCharge = -1;
1041
1042 if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
1043 ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
1044 order = 2;
1045 else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "O4"))) ||
1046 ((!strcmp(name2, "C4")) && (!strcmp(name1, "O4"))))
1047 order = 2;
1048
1049 else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
1050 ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
1051 order = 2;
1052 else
1053 if(((!strcmp(name1, "P"))
1054 && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
1055 || ((!strcmp(name2, "P"))
1056 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
1057 order = 2;
1058 }
1059 break;
1060 case 'G':
1061 if(resn1[2] == 0) {
1062 if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
1063 ai1->formalCharge = -1;
1064 ai1->chemFlag = false;
1065 } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
1066 ai2->formalCharge = -1;
1067 ai2->chemFlag = false;
1068 }
1069 if(((!strcmp(name1, "C6")) && (!strcmp(name2, "O6"))) ||
1070 ((!strcmp(name2, "C6")) && (!strcmp(name1, "O6"))))
1071 order = 2;
1072 else if(((!strcmp(name1, "C2")) && (!strcmp(name2, "N3"))) ||
1073 ((!strcmp(name2, "C2")) && (!strcmp(name1, "N3"))))
1074 order = 2;
1075 else if(((!strcmp(name1, "C8")) && (!strcmp(name2, "N7"))) ||
1076 ((!strcmp(name2, "C8")) && (!strcmp(name1, "N7"))))
1077 order = 2;
1078 else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "C5"))) ||
1079 ((!strcmp(name2, "C4")) && (!strcmp(name1, "C5"))))
1080 order = 2;
1081 else
1082 if(((!strcmp(name1, "P"))
1083 && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
1084 || ((!strcmp(name2, "P"))
1085 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
1086 order = 2;
1087 }
1088 break;
1089 case 'U':
1090 if(resn1[2] == 0) {
1091 if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
1092 ai1->formalCharge = -1;
1093 ai1->chemFlag = false;
1094 } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
1095 ai2->formalCharge = -1;
1096 ai2->chemFlag = false;
1097 }
1098
1099 if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
1100 ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
1101 order = 2;
1102 else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "O4"))) ||
1103 ((!strcmp(name2, "C4")) && (!strcmp(name1, "O4"))))
1104 order = 2;
1105
1106 else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
1107 ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
1108 order = 2;
1109 else
1110 if(((!strcmp(name1, "P"))
1111 && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
1112 || ((!strcmp(name2, "P"))
1113 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
1114 order = 2;
1115 }
1116 break;
1117 }
1118 break;
1119 case 'G':
1120 switch (resn1[1]) {
1121 case 'L':
1122 switch (resn1[2]) {
1123 case 'U': /* GLU missing GLUN, GLUH, GLH handling */
1124 switch (resn1[3]) {
1125 case 0:
1126 case 'M': /* minus */
1127 if(!strcmp(name1, "OE2")) {
1128 ai1->formalCharge = -1;
1129 ai1->chemFlag = false;
1130 } else if(!strcmp(name2, "OE2")) {
1131 ai2->formalCharge = -1;
1132 ai2->chemFlag = false;
1133 }
1134 break;
1135 }
1136 if(((!strcmp(name1, "CD")) && (!strcmp(name2, "OE1"))) ||
1137 ((!strcmp(name2, "CD")) && (!strcmp(name1, "OE1"))))
1138 order = 2;
1139 break;
1140 case 'N': /* GLN or GLU */
1141 if(((!strcmp(name1, "CD")) && (!strcmp(name2, "OE1"))) ||
1142 ((!strcmp(name2, "CD")) && (!strcmp(name1, "OE1"))))
1143 order = 2;
1144 break;
1145 }
1146 break;
1147 case 0:
1148 if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
1149 ai1->formalCharge = -1;
1150 ai1->chemFlag = false;
1151 } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
1152 ai2->formalCharge = -1;
1153 ai2->chemFlag = false;
1154 }
1155
1156 if(((!strcmp(name1, "C6")) && (!strcmp(name2, "O6"))) ||
1157 ((!strcmp(name2, "C6")) && (!strcmp(name1, "O6"))))
1158 order = 2;
1159 else if(((!strcmp(name1, "C2")) && (!strcmp(name2, "N3"))) ||
1160 ((!strcmp(name2, "C2")) && (!strcmp(name1, "N3"))))
1161 order = 2;
1162 else if(((!strcmp(name1, "C8")) && (!strcmp(name2, "N7"))) ||
1163 ((!strcmp(name2, "C8")) && (!strcmp(name1, "N7"))))
1164 order = 2;
1165 else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "C5"))) ||
1166 ((!strcmp(name2, "C4")) && (!strcmp(name1, "C5"))))
1167 order = 2;
1168 else
1169 if(((!strcmp(name1, "P"))
1170 && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
1171 || ((!strcmp(name2, "P"))
1172 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
1173 order = 2;
1174 break;
1175 }
1176 break;
1177 case 'H':
1178 switch (resn1[1]) {
1179 case 'I':
1180 switch (resn1[2]) {
1181 case 'P':
1182 if(!strcmp(name1, "ND1")) {
1183 ai1->formalCharge = 1;
1184 ai1->chemFlag = false;
1185 } else if(!strcmp(name2, "ND1")) {
1186 ai2->formalCharge = 1;
1187 ai2->chemFlag = false;
1188 }
1189 if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
1190 ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
1191 order = 2;
1192 else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "ND1"))) ||
1193 ((!strcmp(name2, "CE1")) && (!strcmp(name1, "ND1"))))
1194 order = 2;
1195 break;
1196 case 'S':
1197 switch (resn1[3]) {
1198 case 'A': /* HISA Gromacs */
1199 case 'D': /* HISD Quanta */
1200 if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
1201 ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
1202 order = 2;
1203 else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "NE2"))) ||
1204 ((!strcmp(name2, "CE1")) && (!strcmp(name1, "NE2"))))
1205 order = 2;
1206 break;
1207 case 0: /* plain HIS */
1208 case 'B': /* HISB Gromacs */
1209 case 'E': /* HISE Quanta */
1210 if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
1211 ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
1212 order = 2;
1213 else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "ND1"))) ||
1214 ((!strcmp(name2, "CE1")) && (!strcmp(name1, "ND1"))))
1215 order = 2;
1216 break;
1217 case 'H': /* HISH Gromacs */
1218 case 'P': /* HISP Quanta */
1219 if(!strcmp(name1, "ND1")) {
1220 ai1->formalCharge = 1;
1221 ai1->chemFlag = false;
1222 } else if(!strcmp(name2, "ND1")) {
1223 ai2->formalCharge = 1;
1224 ai2->chemFlag = false;
1225 }
1226 if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
1227 ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
1228 order = 2;
1229 else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "ND1"))) ||
1230 ((!strcmp(name2, "CE1")) && (!strcmp(name1, "ND1"))))
1231 order = 2;
1232 break;
1233 }
1234 break;
1235 case 'E': /* HIE */
1236 if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
1237 ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
1238 order = 2;
1239 else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "ND1"))) ||
1240 ((!strcmp(name2, "CE1")) && (!strcmp(name1, "ND1"))))
1241 order = 2;
1242 break;
1243 case 'D': /* HID */
1244 if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD2"))) ||
1245 ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD2"))))
1246 order = 2;
1247 else if(((!strcmp(name1, "CE1")) && (!strcmp(name2, "NE2"))) ||
1248 ((!strcmp(name2, "CE1")) && (!strcmp(name1, "NE2"))))
1249 order = 2;
1250 break;
1251 }
1252 break;
1253 }
1254 break;
1255 case 'I':
1256 if(resn1[1] == 0) {
1257 if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
1258 ai1->formalCharge = -1;
1259 ai1->chemFlag = false;
1260 } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
1261 ai2->formalCharge = -1;
1262 ai2->chemFlag = false;
1263 }
1264 if(((!strcmp(name1, "C8")) && (!strcmp(name2, "N7"))) ||
1265 ((!strcmp(name2, "C8")) && (!strcmp(name1, "N7"))))
1266 order = 2;
1267 else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "C5"))) ||
1268 ((!strcmp(name2, "C4")) && (!strcmp(name1, "C5"))))
1269 order = 2;
1270
1271 else if(((!strcmp(name1, "C6")) && (!strcmp(name2, "N1"))) ||
1272 ((!strcmp(name2, "C6")) && (!strcmp(name1, "N1"))))
1273 order = 2;
1274 else if(((!strcmp(name1, "C2")) && (!strcmp(name2, "N3"))) ||
1275 ((!strcmp(name2, "C2")) && (!strcmp(name1, "N3"))))
1276 order = 2;
1277 else
1278 if(((!strcmp(name1, "P"))
1279 && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
1280 || ((!strcmp(name2, "P"))
1281 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
1282 order = 2;
1283 }
1284 break;
1285 case 'P':
1286 switch (resn1[1]) {
1287 case 'H': /* PHE */
1288 if(resn1[2] == 'E') {
1289 if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD1"))) ||
1290 ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD1"))))
1291 order = 2;
1292 else if(((!strcmp(name1, "CZ")) && (!strcmp(name2, "CE1"))) ||
1293 ((!strcmp(name2, "CZ")) && (!strcmp(name1, "CE1"))))
1294 order = 2;
1295
1296 else if(((!strcmp(name1, "CE2")) && (!strcmp(name2, "CD2"))) ||
1297 ((!strcmp(name2, "CE2")) && (!strcmp(name1, "CD2"))))
1298 order = 2;
1299 break;
1300 }
1301 }
1302 break;
1303 case 'L':
1304 switch (resn1[1]) {
1305 case 'Y':
1306 switch (resn1[2]) {
1307 case 'S': /* LYS. */
1308 switch (resn1[3]) {
1309 case 0:
1310 case 'P': /* LYS, LYSP */
1311 if(!strcmp(name1, "NZ")) {
1312 ai1->formalCharge = 1;
1313 ai1->chemFlag = false;
1314 } else if(!strcmp(name2, "NZ")) {
1315 ai2->formalCharge = 1;
1316 ai2->chemFlag = false;
1317 }
1318 break;
1319 }
1320 break;
1321 }
1322 break;
1323 }
1324 break;
1325 case 'T':
1326 switch (resn1[1]) {
1327 case 'Y': /* TYR */
1328 if(resn1[2] == 'R') {
1329 if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD1"))) ||
1330 ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD1"))))
1331 order = 2;
1332 else if(((!strcmp(name1, "CZ")) && (!strcmp(name2, "CE1"))) ||
1333 ((!strcmp(name2, "CZ")) && (!strcmp(name1, "CE1"))))
1334 order = 2;
1335
1336 else if(((!strcmp(name1, "CE2")) && (!strcmp(name2, "CD2"))) ||
1337 ((!strcmp(name2, "CE2")) && (!strcmp(name1, "CD2"))))
1338 order = 2;
1339 break;
1340 }
1341 break;
1342 case 'R':
1343 if(resn1[2] == 'P') {
1344 if(((!strcmp(name1, "CG")) && (!strcmp(name2, "CD1"))) ||
1345 ((!strcmp(name2, "CG")) && (!strcmp(name1, "CD1"))))
1346 order = 2;
1347 else if(((!strcmp(name1, "CZ3")) && (!strcmp(name2, "CE3"))) ||
1348 ((!strcmp(name2, "CZ3")) && (!strcmp(name1, "CE3"))))
1349 order = 2;
1350 else if(((!strcmp(name1, "CZ2")) && (!strcmp(name2, "CH2"))) ||
1351 ((!strcmp(name2, "CZ2")) && (!strcmp(name1, "CH2"))))
1352 order = 2;
1353 else if(((!strcmp(name1, "CE2")) && (!strcmp(name2, "CD2"))) ||
1354 ((!strcmp(name2, "CE2")) && (!strcmp(name1, "CD2"))))
1355 order = 2;
1356 break;
1357 }
1358 break;
1359 case 0:
1360 if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
1361 ai1->formalCharge = -1;
1362 ai1->chemFlag = false;
1363 } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
1364 ai2->formalCharge = -1;
1365 ai2->chemFlag = false;
1366 }
1367
1368 if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
1369 ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
1370 order = 2;
1371 else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "O4"))) ||
1372 ((!strcmp(name2, "C4")) && (!strcmp(name1, "O4"))))
1373 order = 2;
1374
1375 else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
1376 ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
1377 order = 2;
1378 else
1379 if(((!strcmp(name1, "P"))
1380 && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
1381 || ((!strcmp(name2, "P"))
1382 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
1383 order = 2;
1384 break;
1385 }
1386 break;
1387 case 'U':
1388 if(resn1[1] == 0) {
1389 if(((!strcmp(name1, "O2P")) || (!strcmp(name1, "OP2")))) {
1390 ai1->formalCharge = -1;
1391 ai1->chemFlag = false;
1392 } else if(((!strcmp(name2, "O2P")) || (!strcmp(name2, "OP2")))) {
1393 ai2->formalCharge = -1;
1394 ai2->chemFlag = false;
1395 }
1396 if(((!strcmp(name1, "C2")) && (!strcmp(name2, "O2"))) ||
1397 ((!strcmp(name2, "C2")) && (!strcmp(name1, "O2"))))
1398 order = 2;
1399 else if(((!strcmp(name1, "C4")) && (!strcmp(name2, "O4"))) ||
1400 ((!strcmp(name2, "C4")) && (!strcmp(name1, "O4"))))
1401 order = 2;
1402
1403 else if(((!strcmp(name1, "C5")) && (!strcmp(name2, "C6"))) ||
1404 ((!strcmp(name2, "C5")) && (!strcmp(name1, "C6"))))
1405 order = 2;
1406 else
1407 if(((!strcmp(name1, "P"))
1408 && (((!strcmp(name2, "O1P")) || (!strcmp(name2, "OP1")))))
1409 || ((!strcmp(name2, "P"))
1410 && (((!strcmp(name1, "O1P")) || (!strcmp(name1, "OP1"))))))
1411 order = 2;
1412 }
1413 break;
1414 }
1415 }
1416 *(bond_order) = order;
1417 }
1418
ObjectMoleculeFixChemistry(ObjectMolecule * I,int sele1,int sele2,int invalidate)1419 void ObjectMoleculeFixChemistry(ObjectMolecule * I, int sele1, int sele2, int invalidate)
1420 {
1421 int b;
1422 int flag = false;
1423 int s1, s2;
1424 AtomInfoType *ai1, *ai2;
1425 int order;
1426 BondType* bond = I->Bond.data();
1427 for(b = 0; b < I->NBond; b++) {
1428 flag = false;
1429 ai1 = I->AtomInfo + bond->index[0];
1430 ai2 = I->AtomInfo + bond->index[1];
1431 s1 = ai1->selEntry;
1432 s2 = ai2->selEntry;
1433
1434 if((SelectorIsMember(I->G, s1, sele1) &&
1435 SelectorIsMember(I->G, s2, sele2)) ||
1436 (SelectorIsMember(I->G, s2, sele1) && SelectorIsMember(I->G, s1, sele2))) {
1437 order = -1;
1438 if(strlen(LexStr(I->G, ai1->resn)) < 4) { /* Standard disconnected PDB residue */
1439 if(AtomInfoSameResidue(I->G, ai1, ai2)) {
1440 assign_pdb_known_residue(I->G, ai1, ai2, &order);
1441 }
1442 }
1443 if(order > 0) {
1444 bond->order = order;
1445 ai1->chemFlag = false;
1446 ai2->chemFlag = false;
1447 flag = true;
1448 } else if(invalidate) {
1449 ai1->chemFlag = false;
1450 ai2->chemFlag = false;
1451 flag = true;
1452 }
1453 }
1454 bond++;
1455 }
1456 if(flag) {
1457 I->invalidate(cRepAll, cRepInvAll, -1);
1458 SceneChanged(I->G);
1459 }
1460 }
1461
ObjMolPairwiseInit(ObjMolPairwise * pairwise)1462 void ObjMolPairwiseInit(ObjMolPairwise * pairwise)
1463 {
1464 UtilZeroMem((char *) pairwise, sizeof(ObjMolPairwise));
1465 pairwise->trg_vla = VLAlloc(int, 10);
1466 pairwise->mbl_vla = VLAlloc(int, 10);
1467 }
1468
ObjMolPairwisePurge(ObjMolPairwise * pairwise)1469 void ObjMolPairwisePurge(ObjMolPairwise * pairwise)
1470 {
1471 VLAFreeP(pairwise->trg_vla);
1472 VLAFreeP(pairwise->mbl_vla);
1473 }
1474
ObjectMoleculeConvertIDsToIndices(ObjectMolecule * I,int * id,int n_id)1475 int ObjectMoleculeConvertIDsToIndices(ObjectMolecule * I, int *id, int n_id)
1476 {
1477 /* return true if all IDs are unique, false if otherwise */
1478
1479 int min_id, max_id, range, *lookup = NULL;
1480 int unique = true;
1481
1482 /* this routine only works if IDs cover a reasonable range --
1483 should rewrite using a hash table */
1484
1485 if(I->NAtom) {
1486
1487 /* determine range */
1488
1489 {
1490 int a, cur_id;
1491 cur_id = I->AtomInfo[0].id;
1492 min_id = cur_id;
1493 max_id = cur_id;
1494 for(a = 1; a < I->NAtom; a++) {
1495 cur_id = I->AtomInfo[a].id;
1496 if(min_id > cur_id)
1497 min_id = cur_id;
1498 if(max_id < cur_id)
1499 max_id = cur_id;
1500 }
1501 }
1502
1503 /* create cross-reference table */
1504
1505 {
1506 int a, offset;
1507
1508 range = max_id - min_id + 1;
1509 lookup = pymol::calloc<int>(range);
1510 for(a = 0; a < I->NAtom; a++) {
1511 offset = I->AtomInfo[a].id - min_id;
1512 if(!lookup[offset])
1513 lookup[offset] = a + 1;
1514 else
1515 unique = false;
1516 }
1517 }
1518
1519 /* iterate through IDs and replace with indices or -1 */
1520
1521 {
1522 int i, offset, lkup;
1523
1524 for(i = 0; i < n_id; i++) {
1525 offset = id[i] - min_id;
1526 if((offset >= 0) && (offset < range)) {
1527 lkup = lookup[offset];
1528 if(lkup > 0) {
1529 id[i] = lkup - 1;
1530 } else {
1531 id[i] = -1; /* negative means no match */
1532 }
1533 } else
1534 id[i] = -1;
1535 }
1536 }
1537 }
1538
1539 FreeP(lookup);
1540 return unique;
1541
1542 }
1543
check_next_pdb_object(const char * p,int skip_to_next)1544 static const char *check_next_pdb_object(const char *p, int skip_to_next)
1545 {
1546 const char *start = p;
1547 while(*p) {
1548 if(strstartswith(p, "HEADER")) {
1549 if(skip_to_next)
1550 return p;
1551 return start;
1552 } else if(strstartswith(p, "ATOM ") || strstartswith(p, "HETATM")) {
1553 return start;
1554 } else if(skip_to_next && strcmp("END", p) == 0) {
1555 /* if we pass over the END of the current PDB file, then reset start */
1556 start = p;
1557 }
1558 p = nextline(p);
1559 }
1560 return NULL;
1561 }
1562
get_multi_object_status(const char * p)1563 static int get_multi_object_status(const char *p)
1564 { /* expensive -- only call
1565 this if there is no
1566 other way to determine
1567 this information */
1568 int seen_header = 0;
1569 while(*p) {
1570 if(strstartswith(p, "HEADER")) {
1571 if(seen_header) {
1572 return 1;
1573 } else {
1574 seen_header = true;
1575 }
1576 }
1577 p = nextline(p);
1578 }
1579 return -1;
1580 }
1581
1582 /*
1583 * If any atom in I->AtomInfo contains the wildcard character (from
1584 * "atom_name_wildcard" or "wildcard" setting), then set the object-level
1585 * "atom_name_wildcard" setting to " " (disables wildcard matching).
1586 */
ObjectMoleculeAutoDisableAtomNameWildcard(ObjectMolecule * I)1587 int ObjectMoleculeAutoDisableAtomNameWildcard(ObjectMolecule * I)
1588 {
1589 PyMOLGlobals *G = I->G;
1590 char wildcard = 0;
1591 int found_wildcard = false;
1592
1593 {
1594 const char *tmp = SettingGet_s(G, NULL, I->Setting, cSetting_atom_name_wildcard);
1595 if(tmp && tmp[0]) {
1596 wildcard = *tmp;
1597 } else {
1598 tmp = SettingGet_s(G, NULL, I->Setting, cSetting_wildcard);
1599 if(tmp) {
1600 wildcard = *tmp;
1601 }
1602 }
1603 if(wildcard == 32)
1604 wildcard = 0;
1605
1606 }
1607
1608 if(wildcard) {
1609 int a;
1610 const char *p;
1611 char ch;
1612 const AtomInfoType *ai = I->AtomInfo;
1613
1614 for(a = 0; a < I->NAtom; a++) {
1615 p = LexStr(G, ai->name);
1616 while((ch = *(p++))) {
1617 if(ch == wildcard) {
1618 found_wildcard = true;
1619 break;
1620 }
1621 }
1622 ai++;
1623 }
1624 if(found_wildcard) {
1625 ExecutiveSetObjSettingFromString(G, cSetting_atom_name_wildcard, " ",
1626 I, -1, true, true);
1627 }
1628 }
1629 return found_wildcard;
1630 }
1631
1632
1633 /*========================================================================*/
1634 #define PDB_MAX_TAGS 64
1635
ObjectMoleculePDBStr2CoordSetPASS1(PyMOLGlobals * G,int * ok,const char ** restart_model,const char * p,int n_tags,const char * const * tag_start,int * nAtom,char cc[],int quiet,int * bogus_name_alignment,int * ssFlag,const char ** next_pdb,PDBInfoRec * info,int only_read_one_model,int * ignore_conect,int * bondFlag,int * have_bond_order)1636 static void ObjectMoleculePDBStr2CoordSetPASS1(PyMOLGlobals * G, int *ok,
1637 const char **restart_model, const char *p, int n_tags, const char* const* tag_start,
1638 int *nAtom, char cc[], int quiet, int *bogus_name_alignment,
1639 int *ssFlag, const char **next_pdb, PDBInfoRec *info, int only_read_one_model,
1640 int *ignore_conect, int *bondFlag, int *have_bond_order) {
1641 int seen_end_of_atoms = false;
1642 *restart_model = NULL;
1643 while(*ok && *p) {
1644 AddOrthoOutputIfMatchesTags(G, n_tags, *nAtom, tag_start, p, cc, quiet);
1645 if((strstartswith(p, "ATOM ") ||
1646 strstartswith(p, "HETATM")) && (!*restart_model)) {
1647 if(!seen_end_of_atoms)
1648 (*nAtom)++;
1649 if(*bogus_name_alignment) {
1650 ncopy(cc, nskip(p, 12), 4); /* copy the atom field */
1651 if((cc[0] == 32) && (cc[1] != 32)) { /* check to see if indentation was followed correctly, 32 - space */
1652 *bogus_name_alignment = false;
1653 }
1654 }
1655 } else if(strstartswith(p, "HELIX ")){
1656 *ssFlag = true;
1657 }else if(strstartswith(p, "SHEET ")){
1658 *ssFlag = true;
1659 }else if(strstartswith(p, "HEADER")) {
1660 if(*nAtom > 0) { /* if we've already found atom records, then this must be a new pdb */
1661 (*next_pdb) = p;
1662 break;
1663 }
1664 } else if(strstartswith(p, "REMARK")) {
1665 // char * start = p; // USED TO SAVE REMARKS (TODO)
1666 do {
1667 if(info && strncmp(" GENERATED BY TRJCONV", p + 6, 24) == 0)
1668 info->ignore_header_names = true;
1669 p = nextline(p);
1670 AddOrthoOutputIfMatchesTags(G, n_tags, *nAtom, tag_start, p, cc, quiet);
1671 } while(strstartswith(p, "REMARK"));
1672 // REMARKS are string from start to p, but not saved currently (TODO)
1673 continue;
1674 } else if(strstartswith(p, "ENDMDL") && (!*restart_model)) {
1675 *restart_model = nextline(p);
1676 seen_end_of_atoms = true;
1677 if(only_read_one_model)
1678 break;
1679 } else if(strstartswith(p, "END")) { /* stop parsing after END */
1680 ntrim(cc, p, 6);
1681 if(strcmp("END", cc) == 0) { /* END */
1682 seen_end_of_atoms = true;
1683 if(next_pdb) {
1684 p = nextline(p);
1685 ncopy(cc, p, 6);
1686 if(strcmp("HEADER", cc) == 0) {
1687 (*next_pdb) = p; /* found another PDB file after this one... */
1688 } else if(strcmp("ENDMDL", cc) == 0) {
1689 seen_end_of_atoms = false;
1690 }
1691 }
1692 break;
1693 }
1694 } else if(strstartswith(p, "CONECT")) {
1695 if(!*ignore_conect)
1696 *bondFlag = true;
1697 } else if(strstartswith(p, "USER") && (!*restart_model)) {
1698 }
1699 p = nextline(p);
1700 }
1701 }
1702
1703 /*
1704 * Datastructure for efficient array-based secondary structure lookup.
1705 */
1706 typedef struct {
1707 int n_ss; // number of ss_list items
1708 int* ss[256]; // one array for each chain identifier
1709 SSEntry *ss_list; // VLA
1710 } SSHash;
1711
sshash_free(SSHash * hash)1712 static void sshash_free(SSHash *hash) {
1713 int a;
1714 if(!hash)
1715 return;
1716 for(a = 0; a <= 255; a++)
1717 FreeP(hash->ss[a]);
1718 VLAFreeP(hash->ss_list);
1719 FreeP(hash);
1720 }
1721
sshash_new()1722 static SSHash * sshash_new() {
1723 SSHash *hash = pymol::calloc<SSHash>(1);
1724 ok_assert(1, hash);
1725 hash->n_ss = 1;
1726 hash->ss_list = VLAlloc(SSEntry, 50);
1727 ok_assert(1, hash->ss_list);
1728 return hash;
1729 ok_except1:
1730 sshash_free(hash);
1731 return NULL;
1732 }
1733
1734 /*
1735 * Insert a secondary structure record into the hash table.
1736 */
sshash_register_rec(SSHash * hash,unsigned char ss_chain1,int ss_resv1,char ss_inscode1,unsigned char ss_chain2,int ss_resv2,char ss_inscode2,char SSCode)1737 static int sshash_register_rec(SSHash * hash,
1738 unsigned char ss_chain1, int ss_resv1, char ss_inscode1,
1739 unsigned char ss_chain2, int ss_resv2, char ss_inscode2,
1740 char SSCode) {
1741 /* pretty confusing how this works... the following efficient (i.e. array-based)
1742 secondary structure lookup even when there are multiple insertion codes
1743 and where there may be multiple SS records for the residue using different
1744 insertion codes */
1745
1746 unsigned char chain;
1747 int ss_found = false, ssi = 0, a, b, index;
1748 SSEntry *sst;
1749
1750 // up to two iterations:
1751 // 1) assume chain1==chain2
1752 // 2) chains are not the same (undefined in PDB spec?)
1753 for (a = 0, chain = ss_chain1; a < 2; a++, chain = ss_chain2) {
1754 // allocate new array for chain if necc.
1755 if(!hash->ss[chain]) {
1756 ok_assert(1, hash->ss[chain] = pymol::calloc<int>(cResvMask + 1));
1757 }
1758
1759 sst = NULL;
1760 // iterate over all residues indicated
1761 for(b = ss_resv1; b <= ss_resv2; b++) {
1762 index = b & cResvMask;
1763
1764 if(hash->ss[chain][index]) {
1765 // make a unique copy in the event of multiple entries for one resv
1766 sst = NULL;
1767 }
1768
1769 if(!sst) {
1770 VLACheck(hash->ss_list, SSEntry, hash->n_ss);
1771 ok_assert(1, hash->ss_list);
1772 ssi = (hash->n_ss)++;
1773 sst = hash->ss_list + ssi;
1774 sst->resv1 = ss_resv1;
1775 sst->resv2 = ss_resv2;
1776 sst->chain1 = ss_chain1;
1777 sst->chain2 = ss_chain2;
1778 sst->type = SSCode;
1779 sst->inscode1 = ss_inscode1;
1780 sst->inscode2 = ss_inscode2;
1781 ss_found = true;
1782 }
1783 sst->next = hash->ss[chain][index];
1784 hash->ss[chain][index] = ssi;
1785 if(sst->next)
1786 sst = NULL; /* force another unique copy */
1787 }
1788 }
1789 return ss_found;
1790 ok_except1:
1791 return false;
1792 }
1793
1794 /*
1795 * Assign ai->ssType
1796 */
sshash_lookup(SSHash * hash,AtomInfoType * ai,unsigned char ss_chain1)1797 static void sshash_lookup(SSHash *hash, AtomInfoType *ai, unsigned char ss_chain1) {
1798 int index, ssi;
1799 SSEntry *sst = NULL;
1800
1801 index = ai->resv & cResvMask;
1802 if(hash->ss[ss_chain1]) {
1803 ssi = hash->ss[ss_chain1][index];
1804 while(ssi) {
1805 sst = hash->ss_list + ssi;
1806 /* contains shared entry, or unique linked list for each residue */
1807 if( ai->resv >= sst->resv1
1808 && ai->resv <= sst->resv2
1809 && (ai->resv != sst->resv1 || ai->inscode >= sst->inscode1)
1810 && (ai->resv != sst->resv2 || ai->inscode <= sst->inscode2))
1811 {
1812 ai->ssType[0] = sst->type;
1813 return;
1814 }
1815 ssi = sst->next;
1816 }
1817 }
1818 }
1819
1820 /*
1821 * PQR atom line parsing
1822 *
1823 * Try to parse columns white space delimited (10 columns with optional
1824 * chain, 9 without).
1825 *
1826 * Where PQR files come from:
1827 *
1828 * pdb2pqr -> writes PDB-like fixed colums. APBS will fail to read those files
1829 * if columns are too wide.
1830 *
1831 * pdb2pqr --whitespace -> puts extra whitespace, starting at column 2, but
1832 * not between chain and resi! PyMOL <= 2.1 fails to read those files.
1833 *
1834 * APBS Tools Plugin by Michael Lerner adds extra whitespace before negative
1835 * coordinates (assuming -100.0 is the most likely source of error).
1836 *
1837 * @param[in,out] p Pointer to parse from, points after the ATOM field.
1838 * Will move the pointer to the end of the line if
1839 * parsing was successful.
1840 * @param[out] ai Atom to populate with data
1841 * @param[out] coord Coordinates to populate
1842 * @return true on success, false otherwise.
1843 */
parse_pqr_atom_line(PyMOLGlobals * G,const char * & p,AtomInfoType * ai,float * coord)1844 static bool parse_pqr_atom_line(PyMOLGlobals * G,
1845 const char * &p,
1846 AtomInfoType * ai,
1847 float * coord)
1848 {
1849 auto p_eol = nskip(p, 999); // end of line pointer
1850 std::string cc(p, p_eol); // line (starting after ATOM field)
1851
1852 // whitespace splitting
1853 auto columns = strsplit(cc);
1854
1855 // insert chain column if missing
1856 if (columns.size() == 9) {
1857 columns.insert(columns.begin() + 3, "");
1858
1859 // check for concatenated chain + resi
1860 if (columns[4].size() > 4 && !isdigit(columns[4][0])) {
1861 columns[3] = columns[4].substr(0, 1);
1862 columns[4] = columns[4].substr(1);
1863 }
1864 }
1865
1866 // for validation: if number parsing consumes the entire string, then dummy
1867 // will never be populated (sscanf(...) == 1, not 2)
1868 char dummy[2];
1869
1870 // validate numeric fields and populate atom info and coordinates
1871 if (columns.size() == 10 &&
1872 sscanf(columns[0].c_str(), "%d%1s", &ai->id, dummy) == 1 &&
1873 sscanf(columns[5].c_str(), "%f%1s", coord + 0, dummy) == 1 &&
1874 sscanf(columns[6].c_str(), "%f%1s", coord + 1, dummy) == 1 &&
1875 sscanf(columns[7].c_str(), "%f%1s", coord + 2, dummy) == 1 &&
1876 sscanf(columns[8].c_str(), "%f%1s", &ai->partialCharge, dummy) == 1 &&
1877 sscanf(columns[9].c_str(), "%f%1s", &ai->elec_radius, dummy) == 1) {
1878 LexAssign(G, ai->name, columns[1].c_str());
1879 LexAssign(G, ai->resn, columns[2].c_str());
1880 LexAssign(G, ai->chain, columns[3].c_str());
1881 ai->setResi(columns[4].c_str());
1882
1883 // move parser to next line
1884 p = p_eol;
1885
1886 return true;
1887 }
1888
1889 return false;
1890 }
1891
ObjectMoleculePDBStr2CoordSet(PyMOLGlobals * G,const char * buffer,AtomInfoType ** atInfoPtr,const char ** restart_model,char * segi_override,char * pdb_name,const char ** next_pdb,PDBInfoRec * info,int quiet,int * model_number)1892 CoordSet *ObjectMoleculePDBStr2CoordSet(PyMOLGlobals * G,
1893 const char *buffer,
1894 AtomInfoType ** atInfoPtr,
1895 const char **restart_model,
1896 char *segi_override,
1897 char *pdb_name,
1898 const char **next_pdb,
1899 PDBInfoRec * info, int quiet, int *model_number)
1900 {
1901
1902 const char *p;
1903 int nAtom;
1904 int a;
1905 float *coord = NULL;
1906 CoordSet *cset = NULL;
1907 AtomInfoType *atInfo = NULL, *ai;
1908 int AFlag;
1909 char SSCode;
1910 int atomCount;
1911 int bondFlag = false;
1912 BondType *bond = NULL, *ii1, *ii2;
1913 int *idx;
1914 int nBond = 0;
1915 int b1, b2, nReal, maxAt;
1916 CSymmetry *symmetry = NULL;
1917 int auto_show = RepGetAutoShowMask(G);
1918 int reformat_names = SettingGetGlobal_i(G, cSetting_pdb_reformat_names_mode);
1919 int truncate_resn = SettingGetGlobal_b(G, cSetting_pdb_truncate_residue_name);
1920 const char *tags_in = SettingGetGlobal_s(G, cSetting_pdb_echo_tags), *tag_start[PDB_MAX_TAGS];
1921 int n_tags = 0;
1922 int foundNextModelFlag = false;
1923 int ssFlag = false;
1924 int ss_resv1 = 0, ss_resv2 = 0;
1925 char ss_inscode1 = '\0', ss_inscode2 = '\0';
1926 unsigned char ss_chain1 = 0, ss_chain2 = 0;
1927 SSHash *ss_hash = NULL;
1928 char cc[MAXLINELEN], tags[MAXLINELEN];
1929 int ignore_pdb_segi = 0;
1930 int ss_valid, ss_found = false;
1931 int only_read_one_model = false;
1932 int ignore_conect = SettingGetGlobal_b(G, cSetting_pdb_ignore_conect);
1933 int have_bond_order = false;
1934 int seen_model, in_model = false;
1935 int is_end_of_object = false;
1936 int literal_names = SettingGetGlobal_b(G, cSetting_pdb_literal_names);
1937 int bogus_name_alignment = true;
1938 AtomName literal_name = "";
1939 int ok = true;
1940 lexidx_t segi_override_idx = LexIdx(G, segi_override);
1941
1942 if(tags_in && (!quiet) && (!*restart_model)) {
1943 char *p = tags;
1944 strcpy(tags, tags_in);
1945
1946 while(*p) {
1947 while(*p == ' ') /* skip spaces */
1948 p++;
1949 if(*p) {
1950 tag_start[n_tags] = p;
1951 n_tags++;
1952 while(*p) {
1953 if(*p != ',') {
1954 if(*p == ' ')
1955 *p = 0;
1956 p++;
1957 } else
1958 break;
1959 }
1960 if(*p) { /* terminate tag */
1961 *p = 0;
1962 p++;
1963 }
1964 }
1965 }
1966 }
1967
1968 if(literal_names)
1969 reformat_names = 0;
1970
1971 ignore_pdb_segi = SettingGetGlobal_b(G, cSetting_ignore_pdb_segi);
1972
1973 p = buffer;
1974 nAtom = 0;
1975 if(atInfoPtr)
1976 atInfo = *atInfoPtr;
1977
1978 if(!atInfo)
1979 ErrFatal(G, "PDBStr2CoordSet", "need atom information record!"); /* failsafe for old version.. */
1980
1981 if(buffer == *restart_model)
1982 only_read_one_model = true;
1983 else if(info && (info->multiplex > 0)) {
1984 only_read_one_model = true;
1985 if(!info->multi_object_status) { /* figure out if this is a multi-object (multi-HEADER) pdb file */
1986 info->multi_object_status = get_multi_object_status(p);
1987 }
1988 }
1989 /* PASS 1 */
1990 ObjectMoleculePDBStr2CoordSetPASS1(G, &ok, restart_model, p, n_tags,
1991 tag_start, &nAtom, cc, quiet, &bogus_name_alignment, &ssFlag, next_pdb,
1992 info, only_read_one_model, &ignore_conect, &bondFlag, &have_bond_order);
1993 /* INPUT USED:
1994 only_read_one_model - only reads one pdb model if set
1995 n_tags, tag_start - tags defined to report in log
1996 cc - just temporary char * buffer that is used, no input/output
1997 OUTPUT USED:
1998 bogus_name_alignment - whether all ATOM/HETATM has correct names in PDB (12-16, starting with space
1999 ssFlag - if PDB has HELIX and/or SHEET records
2000 both INPUT/OUTPUT:
2001 nAtom - returns the number of atoms read in, also uses it to detect multiple PDBs
2002 ignore_conect - do not set bondFlag if set
2003
2004 END PASS 1 */
2005
2006 *restart_model = NULL;
2007 if (ok){
2008 coord = VLAlloc(float, 3 * nAtom);
2009 CHECKOK(ok, coord);
2010 }
2011
2012 if(ok && atInfo){
2013 VLACheck(atInfo, AtomInfoType, nAtom);
2014 CHECKOK(ok, atInfo);
2015 if (ok)
2016 *atInfoPtr = atInfo;
2017 }
2018 if(ok && bondFlag) {
2019 nBond = 0;
2020 bond = VLACalloc(BondType, 6 * nAtom);
2021 CHECKOK(ok, bond);
2022 }
2023 p = buffer;
2024 PRINTFB(G, FB_ObjectMolecule, FB_Blather)
2025 " ObjectMoleculeReadPDB: Found %i atoms...\n", nAtom ENDFB(G);
2026
2027 if(ok && ssFlag) {
2028 ss_hash = sshash_new();
2029 }
2030
2031 a = 0; /* WATCHOUT */
2032 atomCount = 0;
2033
2034 /* PASS 2 */
2035 seen_model = false;
2036
2037 while(ok && *p) {
2038 /* printf("%c%c%c%c%c%c\n",p[0],p[1],p[2],p[3],p[4],p[5]); */
2039 AFlag = false;
2040 SSCode = 0;
2041 if(strstartswith(p, "ATOM "))
2042 AFlag = 1;
2043 else if(strstartswith(p, "HETATM"))
2044 AFlag = 2;
2045 else if(strstartswith(p, "HEADER")) {
2046
2047 if(pdb_name) {
2048 if(atomCount > 0) {
2049 /* have we already read atoms? then this is probably a new PDB file */
2050 (*next_pdb) = p; /* found another PDB file after this one... */
2051 break;
2052 } else if((!info) || (!info->ignore_header_names)) {
2053 /* if this isn't an MD trajectory... */
2054 const char *pp;
2055 pp = nskip(p, 62); /* is there a four-letter PDB code? */
2056 pp = ntrim(pdb_name, pp, 4);
2057 if(pdb_name[0] == 0) { /* if not, is there a plain name (for MERCK!) */
2058 p = nskip(p, 6);
2059 p = ntrim(cc, p, 44);
2060 UtilNCopy(pdb_name, cc, WordLength);
2061 } else {
2062 p = pp;
2063 }
2064 }
2065 }
2066 } else if(strstartswith(p, "MODEL ")) {
2067 if(model_number) {
2068 int tmp;
2069 p = nskip(p, 10);
2070 p = ncopy(cc, p, 5);
2071 if(sscanf(cc, "%d", &tmp) == 1)
2072 *model_number = tmp;
2073 }
2074 seen_model = true;
2075 in_model = true;
2076 } else if(strstartswith(p, "HELIX ")) {
2077 ss_valid = true;
2078
2079 p = nskip(p, 19);
2080 ss_chain1 = (*p);
2081 p = nskip(p, 2);
2082 p = ncopy(cc, p, 4);
2083 if(!sscanf(cc, "%d", &ss_resv1))
2084 ss_valid = false;
2085 ss_inscode1 = makeInscode(*p);
2086
2087 p = nskip(p, 6);
2088 ss_chain2 = (*p);
2089 p = nskip(p, 2);
2090 p = ncopy(cc, p, 4);
2091
2092 if(!sscanf(cc, "%d", &ss_resv2))
2093 ss_valid = false;
2094 ss_inscode2 = makeInscode(*p);
2095
2096 if(ss_valid) {
2097 PRINTFB(G, FB_ObjectMolecule, FB_Blather)
2098 " ObjectMolecule: read HELIX %c %d%.1s %c %d%.1s\n",
2099 ss_chain1, ss_resv1, &ss_inscode1, ss_chain2, ss_resv2, &ss_inscode2 ENDFB(G);
2100 SSCode = 'H';
2101 }
2102
2103 if(ss_chain1 == ' ')
2104 ss_chain1 = 0;
2105 if(ss_chain2 == ' ')
2106 ss_chain2 = 0;
2107
2108 } else if(strstartswith(p, "SHEET ")) {
2109 ss_valid = true;
2110 p = nskip(p, 21);
2111 ss_chain1 = (*p);
2112 p = nskip(p, 1);
2113 p = ncopy(cc, p, 4);
2114 if(!sscanf(cc, "%d", &ss_resv1))
2115 ss_valid = false;
2116 ss_inscode1 = makeInscode(*p);
2117 p = nskip(p, 6);
2118 ss_chain2 = (*p);
2119 p = nskip(p, 1);
2120 p = ncopy(cc, p, 4);
2121 if(!sscanf(cc, "%d", &ss_resv2))
2122 ss_valid = false;
2123 ss_inscode2 = makeInscode(*p);
2124
2125 if(ss_valid) {
2126 PRINTFB(G, FB_ObjectMolecule, FB_Blather)
2127 " ObjectMolecule: read SHEET %c %d%.1s %c %d%.1s\n",
2128 ss_chain1, ss_resv1, &ss_inscode1, ss_chain2, ss_resv2, &ss_inscode2 ENDFB(G);
2129 SSCode = 'S';
2130 }
2131
2132 if(ss_chain1 == ' ')
2133 ss_chain1 = 0;
2134 if(ss_chain2 == ' ')
2135 ss_chain2 = 0;
2136
2137 } else if(strstartswith(p, "ENDMDL")) {
2138 if(*restart_model)
2139 in_model = false;
2140 else {
2141 *restart_model = nextline(p);
2142 in_model = false;
2143 if(only_read_one_model) {
2144 const char *pp;
2145 pp = nextline(p);
2146 if(strstartswith(pp, "END")) { /* END we're going to be starting a new object... */
2147 (*next_pdb) = check_next_pdb_object(nextline(pp), true);
2148 if(info && (info->multiplex == 0)) {
2149 /* multiplex == 0: FORCED multimodel behavior with concatenated PDB files */
2150 (*restart_model) = (*next_pdb);
2151 (*next_pdb) = NULL;
2152 foundNextModelFlag = true;
2153 info->multi_object_status = -1;
2154 } else {
2155 is_end_of_object = true;
2156 }
2157 } else if(strstartswith(pp, "MODEL")) { /* not a new object...just a new state (model) */
2158 if(info && (info->multiplex > 0)) { /* end object if we're multiplexing */
2159 (*next_pdb) = check_next_pdb_object(pp, true);
2160 (*restart_model) = NULL;
2161 } else
2162 is_end_of_object = false;
2163 } else {
2164 if(pp[0] > 32) /* more content follows... */
2165 (*next_pdb) = check_next_pdb_object(pp, true);
2166 else
2167 (*next_pdb) = NULL; /* at end of file */
2168 }
2169 break;
2170 }
2171 }
2172 } else if(strstartswith(p, "END")) {
2173 ntrim(cc, p, 6);
2174 if((strcmp("END", cc) == 0) && (!in_model)) { /* stop parsing here... */
2175 const char *pp;
2176 pp = nextline(p); /* ...unless we're in MODEL or next line is itself ENDMDL */
2177 p = ncopy(cc, p, 6);
2178 if(!((cc[0] == 'E') && (cc[1] == 'N') && (cc[2] == 'D') && (cc[3] == 'M') && (cc[4] == 'D') && (cc[5] == 'L'))) { /* NOTE: this test seems unnecessary given strcmp above... */
2179 if(!*next_pdb) {
2180 (*next_pdb) = check_next_pdb_object(pp, false);
2181 }
2182 if((*next_pdb) && info && (!info->multiplex) && !(*restart_model)) {
2183 /* multiplex == 0: FORCED multimodel behavior with concatenated PDB files */
2184 (*restart_model) = (*next_pdb);
2185 (*next_pdb) = NULL;
2186 foundNextModelFlag = true;
2187 info->multi_object_status = -1;
2188 is_end_of_object = false;
2189 break;
2190 }
2191 if(*next_pdb) { /* we've found another object... */
2192 if(*restart_model)
2193 is_end_of_object = false; /* however, if we're parsing multi-models, then we're not yet at the end */
2194 else
2195 is_end_of_object = true;
2196 break;
2197 } else if(!seen_model)
2198 break;
2199 }
2200 }
2201 } else if(strstartswith(p, "CRYST1") && (!*restart_model)) {
2202 if(!symmetry){
2203 symmetry = new CSymmetry(G);
2204 CHECKOK(ok, symmetry);
2205 }
2206 if(symmetry) {
2207 int symFlag = true;
2208 PRINTFB(G, FB_ObjectMolecule, FB_Blather)
2209 " PDBStrToCoordSet: Attempting to read symmetry information\n" ENDFB(G);
2210 p = nskip(p, 6);
2211 symFlag = true;
2212 p = ncopy(cc, p, 9);
2213 if(sscanf(cc, "%f", &symmetry->Crystal.Dim[0]) != 1)
2214 symFlag = false;
2215 p = ncopy(cc, p, 9);
2216 if(sscanf(cc, "%f", &symmetry->Crystal.Dim[1]) != 1)
2217 symFlag = false;
2218 p = ncopy(cc, p, 9);
2219 if(sscanf(cc, "%f", &symmetry->Crystal.Dim[2]) != 1)
2220 symFlag = false;
2221 p = ncopy(cc, p, 7);
2222 if(sscanf(cc, "%f", &symmetry->Crystal.Angle[0]) != 1)
2223 symFlag = false;
2224 p = ncopy(cc, p, 7);
2225 if(sscanf(cc, "%f", &symmetry->Crystal.Angle[1]) != 1)
2226 symFlag = false;
2227 p = ncopy(cc, p, 7);
2228 if(sscanf(cc, "%f", &symmetry->Crystal.Angle[2]) != 1)
2229 symFlag = false;
2230 p = nskip(p, 1);
2231 p = ncopy(symmetry->SpaceGroup, p, 11);
2232 UtilCleanStr(symmetry->SpaceGroup);
2233 p = ncopy(cc, p, 3);
2234 if(sscanf(cc, "%d", &symmetry->PDBZValue) != 1)
2235 symmetry->PDBZValue = 1;
2236 if(!symFlag) {
2237 ErrMessage(G, "PDBStrToCoordSet", "Error reading CRYST1 record\n");
2238 SymmetryFree(symmetry);
2239 symmetry = NULL;
2240 }
2241 }
2242 } else if(strstartswith(p, "SCALE") && (!*restart_model) && info) { /* SCALEn */
2243 switch (p[5]) {
2244 case '1':
2245 case '2':
2246 case '3':
2247 {
2248 int flag = (p[5] - '1');
2249 int offset = flag * 4;
2250 int scale_flag = true;
2251 p = nskip(p, 10);
2252 p = ncopy(cc, p, 10);
2253 if(sscanf(cc, "%f", &info->scale.matrix[offset]) != 1)
2254 scale_flag = false;
2255 p = ncopy(cc, p, 10);
2256 if(sscanf(cc, "%f", &info->scale.matrix[offset + 1]) != 1)
2257 scale_flag = false;
2258 p = ncopy(cc, p, 10);
2259 if(sscanf(cc, "%f", &info->scale.matrix[offset + 2]) != 1)
2260 scale_flag = false;
2261 p = nskip(p, 5);
2262 p = ncopy(cc, p, 10);
2263 if(sscanf(cc, "%f", &info->scale.matrix[offset + 3]) != 1)
2264 scale_flag = false;
2265 if(scale_flag)
2266 info->scale.flag[flag] = true;
2267 PRINTFB(G, FB_ObjectMolecule, FB_Blather)
2268 " PDBStrToCoordSet: SCALE%d %8.4f %8.4f %8.4f %8.4f\n", flag + 1,
2269 info->scale.matrix[offset],
2270 info->scale.matrix[offset + 1],
2271 info->scale.matrix[offset + 2], info->scale.matrix[offset + 3]
2272 ENDFB(G);
2273 }
2274 break;
2275 }
2276 } else if(strstartswith(p, "CONECT") &&
2277 bondFlag && (!ignore_conect) && ((!*restart_model) || (!in_model))) {
2278 p = nskip(p, 6);
2279 p = ncopy(cc, p, 5);
2280 if(sscanf(cc, "%d", &b1) == 1)
2281 while(1) {
2282 p = ncopy(cc, p, 5);
2283 if(sscanf(cc, "%d", &b2) != 1)
2284 break;
2285 else {
2286 if((b1 >= 0) && (b2 >= 0) && (b1 != b2)) { /* IDs must be positive and distinct */
2287 VLACheck(bond, BondType, nBond);
2288 CHECKOK(ok, bond);
2289 if (ok){
2290 if(b1 <= b2) {
2291 bond[nBond].index[0] = b1; /* temporarily store the atom indexes */
2292 bond[nBond].index[1] = b2;
2293 bond[nBond].order = 1;
2294 bond[nBond].stereo = 0;
2295 } else {
2296 bond[nBond].index[0] = b2;
2297 bond[nBond].index[1] = b1;
2298 bond[nBond].order = 1;
2299 bond[nBond].stereo = 0;
2300 }
2301 nBond++;
2302 }
2303 }
2304 }
2305 }
2306 } else if(strstartswith(p, "USER") && (!*restart_model)) {
2307 } else if(strstartswith(p, "ANISOU") && (!*restart_model) && (atomCount)) {
2308 ai = atInfo + atomCount - 1;
2309
2310 /* TODO: check atom identifier match */
2311
2312 {
2313 int dummy;
2314 p = nskip(p, 6);
2315 p = ncopy(cc, p, 5);
2316 if(!sscanf(cc, "%d", &dummy))
2317 dummy = 0;
2318 if(dummy == ai->id) { /* ATOM ID must match */
2319 int dummy;
2320 float * anisou = ai->get_anisou();
2321 p = nskip(p, 17);
2322 for (int i = 0; i < 6; ++i) {
2323 p = ncopy(cc, p, 7);
2324 if(sscanf(cc, "%d", &dummy))
2325 anisou[i] = dummy / 10000.0F;
2326 }
2327 }
2328 }
2329 }
2330
2331 /* END KEYWORDS */
2332
2333 /* Secondary structure records */
2334
2335 if(ok && SSCode) {
2336 ss_found = sshash_register_rec(ss_hash,
2337 ss_chain1, ss_resv1, ss_inscode1,
2338 ss_chain2, ss_resv2, ss_inscode2, SSCode);
2339 }
2340 /* Atom records */
2341
2342 if(ok && AFlag && (!*restart_model)) {
2343 ai = atInfo + atomCount;
2344 p = nskip(p, 6);
2345
2346 ai->rank = atomCount;
2347
2348 if(info && info->is_pqr_file()) {
2349 if (parse_pqr_atom_line(G, p, ai, coord + a)) {
2350 goto pqr_done;
2351 }
2352 }
2353
2354 p = ncopy(cc, p, 5);
2355 if(!sscanf(cc, "%d", &ai->id))
2356 ai->id = 0;
2357
2358 p = nskip(p, 1); /* to 12 */
2359 p = ncopy(literal_name, p, 4);
2360 if(literal_names) {
2361 LexAssign(G, ai->name, literal_name);
2362 } else {
2363 ParseNTrim(cc, literal_name, 4);
2364 LexAssign(G, ai->name, cc);
2365 }
2366
2367 p = ncopy(cc, p, 1);
2368 if(*cc == 32)
2369 ai->alt[0] = 0;
2370 else {
2371 ai->alt[0] = *cc;
2372 ai->alt[1] = 0;
2373 }
2374
2375 p = ntrim(cc, p, 4); /* now allowing for 4-letter residues */
2376 if (truncate_resn) /* unless specifically disabled */
2377 cc[3] = 0;
2378
2379 LexAssign(G, ai->resn, cc);
2380
2381 if(ai->name) {
2382 const char * ai_name = LexStr(G, ai->name);
2383 int name_len = strlen(ai_name);
2384 char name[5];
2385 switch (reformat_names) {
2386 case 1: /* pdb compliant: HH12 becomes 2HH1, etc. */
2387 if(name_len > 3) {
2388 if((ai_name[0] >= 'A') && ((ai_name[0] <= 'Z')) &&
2389 isdigit(ai_name[3])) {
2390 if(!(((ai_name[1] >= 'a') && (ai_name[1] <= 'z')) ||
2391 ((ai_name[0] == 'C') && (ai_name[1] == 'L')) || /* try to be smart about */
2392 ((ai_name[0] == 'B') && (ai_name[1] == 'R')) || /* distinguishing common atoms */
2393 ((ai_name[0] == 'C') && (ai_name[1] == 'A')) || /* in all-caps from typical */
2394 ((ai_name[0] == 'F') && (ai_name[1] == 'E')) || /* nonatomic abbreviations */
2395 ((ai_name[0] == 'C') && (ai_name[1] == 'U')) ||
2396 ((ai_name[0] == 'N') && (ai_name[1] == 'A')) ||
2397 ((ai_name[0] == 'N') && (ai_name[1] == 'I')) ||
2398 ((ai_name[0] == 'M') && (ai_name[1] == 'G')) ||
2399 ((ai_name[0] == 'M') && (ai_name[1] == 'N')) ||
2400 ((ai_name[0] == 'H') && (ai_name[1] == 'G')) ||
2401 ((ai_name[0] == 'S') && (ai_name[1] == 'E')) ||
2402 ((ai_name[0] == 'S') && (ai_name[1] == 'I')) ||
2403 ((ai_name[0] == 'Z') && (ai_name[1] == 'N'))
2404 )) {
2405 strncpy(name + 1, ai_name, 3);
2406 name[0] = ai_name[3];
2407 name[4] = 0;
2408 LexAssign(G, ai->name, name);
2409 }
2410 }
2411 } else if(name_len == 3) {
2412 if((ai_name[0] == 'H') &&
2413 (ai_name[1] >= 'A') && ((ai_name[1] <= 'Z')) &&
2414 isdigit(ai_name[2])) {
2415 AtomInfoGetPDB3LetHydroName(G, LexStr(G, ai->resn), ai_name, name);
2416 LexAssign(G, ai->name, (name[0] == ' ') ? (name + 1) : name);
2417 }
2418 }
2419 break;
2420 case 2: /* amber compliant: 2HH1 becomes HH12 */
2421 case 3: /* pdb compliant, but use IUPAC within PyMOL */
2422 if(ai_name[0]) {
2423 if(isdigit(ai_name[0]) && ai_name[1] && (!isdigit(ai_name[1]))) {
2424 if (1 < name_len && name_len < 5) {
2425 strcpy(name, ai_name + 1);
2426 name[name_len - 1] = ai_name[0];
2427 name[name_len] = 0;
2428 LexAssign(G, ai->name, name);
2429 }
2430 break;
2431 default: /* AS IS */
2432 break;
2433 }
2434 }
2435 break;
2436 case 4: /* simply read trim and write back out with 3-letter names starting from the
2437 second column, and four-letter names starting in the first */
2438 ntrim(cc, ai_name, 4);
2439 LexAssign(G, ai->name, cc);
2440 break;
2441 }
2442 }
2443
2444 p = ncopy(cc, p, 1);
2445 if (ai->chain){
2446 LexDec(G, ai->chain);
2447 }
2448 if(*cc == ' ') {
2449 ss_chain1 = 0;
2450 ai->chain = 0;
2451 } else {
2452 ss_chain1 = *cc;
2453 ai->chain = LexIdx(G, cc);
2454 }
2455
2456 p = ncopy(cc, p, 4);
2457 if(!sscanf(cc, "%d", &ai->resv))
2458 ai->resv = 0;
2459 ai->setInscode(*p);
2460 p = nskip(p, 1);
2461
2462 if(ssFlag) { /* get secondary structure information (if avail) */
2463 sshash_lookup(ss_hash, ai, ss_chain1);
2464 } else {
2465 ai->cartoon = cCartoon_tube;
2466 }
2467
2468 {
2469 p = nskip(p, 3);
2470 p = ncopy(cc, p, 8);
2471 sscanf(cc, "%f", coord + a);
2472 p = ncopy(cc, p, 8);
2473 sscanf(cc, "%f", coord + (a + 1));
2474 p = ncopy(cc, p, 8);
2475 sscanf(cc, "%f", coord + (a + 2));
2476 }
2477
2478 if((!info) || (!info->is_pqr_file())) { /* standard PDB file */
2479 p = ncopy(cc, p, 6);
2480 if(!sscanf(cc, "%f", &ai->q))
2481 ai->q = 1.0;
2482
2483 p = ncopy(cc, p, 6);
2484 if(!sscanf(cc, "%f", &ai->b))
2485 ai->b = 0.0;
2486
2487 if (info->variant == PDB_VARIANT_PDBQT) {
2488 ignore_pdb_segi = true;
2489 p = nskip(p, 4);
2490 p = ncopy(cc, p, 6);
2491 if(!sscanf(cc, "%f", &ai->partialCharge))
2492 ai->partialCharge = 0.0;
2493
2494 // type is 78-79 in pdbqt, 77-78 in pdb
2495 p = nskip(p, 1);
2496 } else {
2497 p = nskip(p, 6);
2498 p = ncopy(cc, p, 4);
2499 }
2500
2501 if(!ignore_pdb_segi) {
2502 if(!segi_override_idx) {
2503 if(cc[3] == '1' && atomCount && strncmp(p, "0000", 4) == 0) {
2504 /* atom ID overflow? (nonstandard use...)... */
2505 LexAssign(G, segi_override_idx, (ai - 1)->segi);
2506 LexAssign(G, ai->segi, (ai - 1)->segi);
2507 } else {
2508 UtilCleanStr(cc);
2509 LexAssign(G, ai->segi, cc);
2510 }
2511 } else {
2512 LexAssign(G, ai->segi, segi_override_idx);
2513 }
2514 } else {
2515 LexAssign(G, ai->segi, 0);
2516 }
2517
2518 p = ncopy(cc, p, 2);
2519 if(!sscanf(cc, "%s", ai->elem))
2520 ai->elem[0] = 0;
2521 else if(!((((ai->elem[0] >= 'a') && (ai->elem[0] <= 'z')) || /* don't get confused by PDB misuse */
2522 ((ai->elem[0] >= 'A') && (ai->elem[0] <= 'Z'))) &&
2523 (((ai->elem[1] == 0) ||
2524 ((ai->elem[1] >= 'a') && (ai->elem[1] <= 'z')) ||
2525 ((ai->elem[1] >= 'A') && (ai->elem[1] <= 'Z'))))))
2526 ai->elem[0] = 0;
2527 else if (info->variant == PDB_VARIANT_PDBQT) {
2528 if (strcmp(ai->elem, "A") == 0) {
2529 // aromatic carbon
2530 ai->elem[0] = 'C';
2531 } else if (isupper(ai->elem[1])) {
2532 // h-bond donor or acceptor
2533 ai->elem[1] = 0;
2534 }
2535 }
2536
2537 if(!ai->elem[0]) {
2538 if(((literal_name[0] == ' ') || ((literal_name[0] >= '0') && (literal_name[0] <= '9'))) && (literal_name[1] >= 'A') && (literal_name[1] <= 'Z')) { /* infer element from name column */
2539 ai->elem[0] = literal_name[1];
2540 ai->elem[1] = 0;
2541 } else if(((literal_name[0] >= 'A') && (literal_name[0] <= 'Z')) && (((literal_name[1] >= 'A') && (literal_name[1] <= 'Z')) || ((literal_name[1] >= 'a') && (literal_name[1] <= 'z')))) { /* infer element from name column */
2542 ai->elem[0] = literal_name[0];
2543 ai->elem[2] = 0;
2544 if((literal_name[1] >= 'A') && (literal_name[1] <= 'Z')) { /* second letter is capitalized */
2545 if(bogus_name_alignment) {
2546 /* if other atom names aren't properly aligned */
2547 ai->elem[1] = 0; /* kill 2nd letter */
2548 } else if(literal_name[0] == 'H') {
2549 /* or if this is an ultra-bogus PDB with inconsistent
2550 indendentation, and this is likely a hydrogen */
2551 ai->elem[1] = 0; /* kill 2nd letter */
2552 } else {
2553 ai->elem[1] = tolower(literal_name[1]);
2554 }
2555 } else
2556 ai->elem[1] = literal_name[1];
2557 }
2558 }
2559
2560 p = ncopy(cc, p, 2);
2561 if((cc[1] == '-') || (cc[1] == '+')) {
2562 /* only read formal charge when sign is present */
2563 char ctmp = cc[0];
2564 cc[0] = cc[1];
2565 cc[1] = ctmp;
2566 if(!sscanf(cc, "%hhi", &ai->formalCharge))
2567 ai->formalCharge = 0;
2568 }
2569
2570 /* end normal PDB */
2571 } else if(info && info->is_pqr_file()) {
2572 p = ParseWordNumberCopy(cc, p, MAXLINELEN - 1);
2573 if(!sscanf(cc, "%f", &ai->partialCharge))
2574 ai->partialCharge = 0.0F;
2575
2576 p = ParseWordNumberCopy(cc, p, MAXLINELEN - 1);
2577 if(sscanf(cc, "%f", &ai->elec_radius) != 1)
2578 ai->elec_radius = 0.0F;
2579 }
2580
2581 pqr_done:
2582
2583 ai->visRep = auto_show;
2584
2585 if(AFlag == 1)
2586 ai->hetatm = 0;
2587 else {
2588 ai->hetatm = 1;
2589 ai->flags = cAtomFlag_ignore;
2590 }
2591
2592 AtomInfoAssignParameters(G, ai);
2593 AtomInfoAssignColors(G, ai);
2594
2595 PRINTFD(G, FB_ObjectMolecule)
2596 "%s %s %d%c %s %8.3f %8.3f %8.3f %6.2f %6.2f %s\n",
2597 LexStr(G, ai->name), LexStr(G, ai->resn), ai->resv, ai->getInscode(true), LexStr(G, ai->chain),
2598 *(coord + a), *(coord + a + 1), *(coord + a + 2), ai->b, ai->q, LexStr(G, ai->segi) ENDFD;
2599
2600 if(atomCount < nAtom) { /* safety */
2601 a += 3;
2602 atomCount++;
2603 }
2604 }
2605 p = nextline(p);
2606 }
2607
2608 /* END PASS 2 */
2609
2610 if(ok && bondFlag) {
2611 UtilSortInPlace(G, bond, nBond, sizeof(BondType), (UtilOrderFn *) BondInOrder);
2612 if(nBond) {
2613 if(!have_bond_order) { /* handle PDB bond-order kludge */
2614 ii1 = bond;
2615 ii2 = bond + 1;
2616 nReal = 1;
2617 ii1->order = 1;
2618 a = nBond - 1;
2619 while(a) {
2620 if((ii1->index[0] == ii2->index[0]) && (ii1->index[1] == ii2->index[1])) {
2621 ii1->order++; /* count dup */
2622 } else {
2623 ii1++; /* non-dup, make copy */
2624 ii1->index[0] = ii2->index[0];
2625 ii1->index[1] = ii2->index[1];
2626 ii1->order = ii2->order;
2627 ii1->stereo = ii2->stereo;
2628 nReal++;
2629 }
2630 ii2++;
2631 a--;
2632 }
2633 nBond = nReal;
2634 }
2635 /* now, find atoms we're looking for */
2636
2637 /* determine ranges */
2638 maxAt = nAtom;
2639 ii1 = bond;
2640 for(a = 0; a < nBond; a++) {
2641 if(ii1->index[0] > maxAt)
2642 maxAt = ii1->index[0];
2643 if(ii1->index[1] > maxAt)
2644 maxAt = ii1->index[1];
2645 ii1++;
2646 }
2647 for(a = 0; a < nAtom; a++)
2648 if(maxAt < atInfo[a].id)
2649 maxAt = atInfo[a].id;
2650 /* build index */
2651 maxAt++;
2652 idx = pymol::malloc<int>(maxAt + 1);
2653 CHECKOK(ok, idx);
2654 if (ok){
2655 for(a = 0; a < maxAt; a++) {
2656 idx[a] = -1;
2657 }
2658 for(a = 0; a < nAtom; a++)
2659 idx[atInfo[a].id] = a;
2660
2661 /* convert indices to bonds */
2662 ii1 = bond;
2663 ii2 = bond;
2664 nReal = 0;
2665 }
2666 if (ok) {
2667 int unbond_cations = SettingGetGlobal_i(G, cSetting_pdb_unbond_cations);
2668 int flag;
2669
2670 for(a = 0; a < nBond; a++) {
2671
2672 if((ii1->index[0] >= 0) && ((ii1->index[1]) >= 0)) {
2673 if((idx[ii1->index[0]] >= 0) && (idx[ii1->index[1]] >= 0)) { /* in case PDB file has bad bonds */
2674 ii2->index[0] = idx[ii1->index[0]];
2675 ii2->index[1] = idx[ii1->index[1]];
2676 ii2->order = ii1->order;
2677 if((ii2->index[0] >= 0) && (ii2->index[1] >= 0)) {
2678
2679 if(!have_bond_order) { /* handle PDB bond order kludge */
2680 if(ii1->order <= 2)
2681 ii2->order = 1;
2682 else if(ii1->order <= 4)
2683 ii2->order = 2;
2684 else if(ii1->order <= 6)
2685 ii2->order = 3;
2686 else
2687 ii2->order = 4;
2688 }
2689 flag = true;
2690 if(unbond_cations) {
2691 if(AtomInfoIsFreeCation(G, atInfo + ii2->index[0]))
2692 flag = false;
2693 else if(AtomInfoIsFreeCation(G, atInfo + ii2->index[1]))
2694 flag = false;
2695 }
2696 if(flag) {
2697 atInfo[ii2->index[0]].bonded = true;
2698 atInfo[ii2->index[1]].bonded = true;
2699 nReal++;
2700 ii2++;
2701 }
2702 }
2703 }
2704 }
2705 ii1++;
2706 }
2707 }
2708 nBond = nReal;
2709 FreeP(idx);
2710 }
2711 }
2712 if(ss_found && !quiet) {
2713 PRINTFB(G, FB_ObjectMolecule, FB_Details)
2714 " ObjectMolecule: Read secondary structure assignments.\n" ENDFB(G);
2715 }
2716 if(symmetry && !quiet && (!only_read_one_model)) {
2717 PRINTFB(G, FB_ObjectMolecule, FB_Details)
2718 " ObjectMolecule: Read crystal symmetry information.\n" ENDFB(G);
2719 }
2720 PRINTFB(G, FB_ObjectMolecule, FB_Blather)
2721 " PDBStr2CoordSet: Read %d bonds from CONECT records (%p).\n", nBond,
2722 (void *) bond ENDFB(G);
2723
2724 if (ok){
2725 cset = CoordSetNew(G);
2726 CHECKOK(ok, cset);
2727
2728 cset->NIndex = nAtom;
2729 cset->Coord = pymol::vla_take_ownership(coord);
2730 cset->TmpBond = pymol::vla_take_ownership(bond);
2731 cset->NTmpBond = nBond;
2732 if(symmetry)
2733 cset->Symmetry = std::unique_ptr<CSymmetry>(symmetry);
2734 if(atInfoPtr)
2735 *atInfoPtr = atInfo;
2736
2737 if((*restart_model) && (!foundNextModelFlag)) {
2738 /* if plan on need to reading another model into this object,
2739 make sure there is another model to read... */
2740 p = *restart_model;
2741 while(*p) {
2742 if(strstartswith(p, "HEADER")) {
2743 /* seeing HEADER means we're off the end of the existing file */
2744 break;
2745 } else if(strstartswith(p, "MODEL ") ||
2746 strstartswith(p, "ENDMDL")) {
2747 foundNextModelFlag = true;
2748 break;
2749 }
2750 p = nextline(p);
2751 }
2752 if(!foundNextModelFlag) {
2753 *restart_model = NULL;
2754 }
2755 }
2756 }
2757 if (!ok){
2758 if (cset){
2759 cset->fFree();
2760 cset = NULL;
2761 } else {
2762 VLAFreeP(coord);
2763 VLAFreeP(bond);
2764 }
2765 }
2766 sshash_free(ss_hash);
2767
2768 if (ok){
2769 if(!seen_model)
2770 *model_number = 1;
2771
2772 if((*restart_model) && (*next_pdb)) { /* if we're mixing multistate objects and
2773 trajectories, then enforce sanity by
2774 reading the models first... */
2775 if(is_end_of_object)
2776 (*restart_model) = NULL;
2777 else if((*next_pdb) < (*restart_model))
2778 (*next_pdb) = NULL;
2779 }
2780 }
2781
2782 LexDec(G, segi_override_idx);
2783 return (cset);
2784 }
2785
2786
2787 /*========================================================================*/
2788
ObjectMoleculeInitHBondCriteria(PyMOLGlobals * G,HBondCriteria * hbc)2789 void ObjectMoleculeInitHBondCriteria(PyMOLGlobals * G, HBondCriteria * hbc)
2790 {
2791 hbc->maxAngle = SettingGet_f(G, NULL, NULL, cSetting_h_bond_max_angle);
2792 hbc->maxDistAtMaxAngle = SettingGet_f(G, NULL, NULL, cSetting_h_bond_cutoff_edge);
2793 hbc->maxDistAtZero = SettingGet_f(G, NULL, NULL, cSetting_h_bond_cutoff_center);
2794 hbc->power_a = SettingGet_f(G, NULL, NULL, cSetting_h_bond_power_a);
2795 hbc->power_b = SettingGet_f(G, NULL, NULL, cSetting_h_bond_power_b);
2796 hbc->cone_dangle =
2797 (float) cos(PI * 0.5 * SettingGet_f(G, NULL, NULL, cSetting_h_bond_cone) / 180.0F);
2798 if(hbc->maxDistAtMaxAngle != 0.0F) {
2799 hbc->factor_a = (float) (0.5 / pow(hbc->maxAngle, hbc->power_a));
2800 hbc->factor_b = (float) (0.5 / pow(hbc->maxAngle, hbc->power_b));
2801 }
2802 }
2803
ObjectMoleculeTestHBond(float * donToAcc,float * donToH,float * hToAcc,float * accPlane,HBondCriteria * hbc)2804 static int ObjectMoleculeTestHBond(float *donToAcc, float *donToH, float *hToAcc,
2805 float *accPlane, HBondCriteria * hbc)
2806 {
2807 float nDonToAcc[3], nDonToH[3], nAccPlane[3], nHToAcc[3];
2808 double angle;
2809 double cutoff;
2810 double curve;
2811 double dist;
2812 float dangle;
2813
2814 /* A ~~ H - D */
2815
2816 normalize23f(donToAcc, nDonToAcc);
2817 normalize23f(hToAcc, nHToAcc);
2818 if(accPlane) { /* remember, plane need not exist if it's water... */
2819 normalize23f(accPlane, nAccPlane);
2820 if(dot_product3f(nHToAcc, nAccPlane) > (-hbc->cone_dangle)) /* don't allow H behind Acceptor plane */
2821 return 0;
2822 }
2823
2824 normalize23f(donToH, nDonToH);
2825 normalize23f(donToAcc, nDonToAcc);
2826
2827 dangle = dot_product3f(nDonToH, nDonToAcc);
2828 if((dangle < 1.0F) && (dangle > 0.0F))
2829 angle = 180.0 * acos((double) dangle) / PI;
2830 else if(dangle > 0.0F)
2831 angle = 0.0;
2832 else
2833 angle = 90.0;
2834
2835 if(angle > hbc->maxAngle)
2836 return 0;
2837
2838 /* interpolate cutoff based on ADH angle */
2839
2840 if(hbc->maxDistAtMaxAngle != 0.0F) {
2841 curve = (pow(angle, (double) hbc->power_a) * hbc->factor_a +
2842 pow(angle, (double) hbc->power_b) * hbc->factor_b);
2843
2844 cutoff = (hbc->maxDistAtMaxAngle * curve) + (hbc->maxDistAtZero * (1.0 - curve));
2845 } else {
2846 cutoff = hbc->maxDistAtZero;
2847 }
2848
2849 /*
2850 printf("angle %8.3f curve %8.3f %8.3f %8.3f %8.3f\n",angle,
2851 curve,cutoff,hbc->maxDistAtMaxAngle,hbc->maxDistAtZero);
2852 */
2853
2854 dist = length3f(donToAcc);
2855
2856 if(dist > cutoff)
2857 return 0;
2858 else
2859 return 1;
2860
2861 }
2862
2863
2864 /*========================================================================*/
2865
ObjectMoleculeFindBestDonorH(ObjectMolecule * I,int atom,int state,float * dir,float * best,AtomInfoType ** h_real)2866 static int ObjectMoleculeFindBestDonorH(ObjectMolecule * I,
2867 int atom,
2868 int state, float *dir, float *best,
2869 AtomInfoType **h_real)
2870 {
2871 int result = 0;
2872 CoordSet *cs;
2873 int n, nn;
2874 int a1;
2875 float cand[3], cand_dir[3];
2876 float best_dot = 0.0F, cand_dot;
2877
2878 ObjectMoleculeUpdateNeighbors(I);
2879
2880 if((state >= 0) && (state < I->NCSet) && (cs = I->CSet[state]) && (atom < I->NAtom)) {
2881
2882 auto idx = cs->atmToIdx(atom);
2883
2884 if(idx >= 0) {
2885
2886 const float* orig = cs->coordPtr(idx);
2887
2888 /* do we need to add any new hydrogens? */
2889
2890 n = I->Neighbor[atom];
2891 nn = I->Neighbor[n++];
2892
2893 /* printf("nn %d valence %d %s\n",nn,
2894 I->AtomInfo[atom].valence,I->AtomInfo[atom].name);
2895 {
2896 int i;
2897 for(i=0;i<nn;i++) {
2898 printf("%d \n",I->Neighbor[n+2*i]);
2899 }
2900 }
2901 */
2902
2903 if((nn < I->AtomInfo[atom].valence) || I->AtomInfo[atom].hb_donor) { /* is there an implicit hydrogen? */
2904 if(ObjectMoleculeFindOpenValenceVector(I, state, atom, best, dir, -1)) {
2905 result = true;
2906 best_dot = dot_product3f(best, dir);
2907 add3f(orig, best, best);
2908 if(h_real)
2909 *h_real = NULL;
2910 }
2911 }
2912 /* iterate through real hydrogens looking for best match
2913 with desired direction */
2914
2915 while(1) { /* look for an attached non-hydrogen as a base */
2916 a1 = I->Neighbor[n];
2917 n += 2;
2918 if(a1 < 0)
2919 break;
2920 if(I->AtomInfo[a1].protons == 1) { /* hydrogen */
2921 if(ObjectMoleculeGetAtomVertex(I, state, a1, cand)) { /* present */
2922
2923 subtract3f(cand, orig, cand_dir);
2924 normalize3f(cand_dir);
2925 cand_dot = dot_product3f(cand_dir, dir);
2926 if(result) { /* improved */
2927 if((best_dot < cand_dot) || (h_real && !*h_real)) {
2928 best_dot = cand_dot;
2929 copy3f(cand, best);
2930 if(h_real)
2931 *h_real = I->AtomInfo + a1;
2932 }
2933 } else { /* first */
2934 result = true;
2935 copy3f(cand, best);
2936 best_dot = cand_dot;
2937 if(h_real)
2938 *h_real = I->AtomInfo + a1;
2939 }
2940 }
2941 }
2942 }
2943 }
2944 }
2945 return result;
2946 }
2947
2948
2949 /*========================================================================*/
2950
ObjectMoleculeGetCheckHBond(AtomInfoType ** h_real,float * h_crd_ret,ObjectMolecule * don_obj,int don_atom,int don_state,ObjectMolecule * acc_obj,int acc_atom,int acc_state,HBondCriteria * hbc)2951 int ObjectMoleculeGetCheckHBond(AtomInfoType **h_real,
2952 float *h_crd_ret,
2953 ObjectMolecule * don_obj,
2954 int don_atom,
2955 int don_state,
2956 ObjectMolecule * acc_obj,
2957 int acc_atom,
2958 int acc_state,
2959 HBondCriteria * hbc)
2960 {
2961 int result = 0;
2962 const CoordSet *csD, *csA;
2963 float donToAcc[3];
2964 float donToH[3];
2965 float bestH[3];
2966 float hToAcc[3];
2967 float accPlane[3], *vAccPlane = NULL;
2968
2969 /* first, check for existence of coordinate sets */
2970
2971 if((don_state >= 0) &&
2972 (don_state < don_obj->NCSet) &&
2973 (csD = don_obj->CSet[don_state]) &&
2974 (acc_state >= 0) &&
2975 (acc_state < acc_obj->NCSet) &&
2976 (csA = acc_obj->CSet[acc_state]) &&
2977 (don_atom < don_obj->NAtom) && (acc_atom < acc_obj->NAtom)) {
2978
2979 /* now check for coordinates of these actual atoms */
2980
2981 auto idxD = csD->atmToIdx(don_atom);
2982 auto idxA = csA->atmToIdx(acc_atom);
2983
2984 if((idxA >= 0) && (idxD >= 0)) {
2985
2986 /* now get local geometries, including
2987 real or virtual hydrogen atom positions */
2988
2989 const float* vDon = csD->coordPtr(idxD);
2990 const float* vAcc = csA->coordPtr(idxA);
2991
2992 subtract3f(vAcc, vDon, donToAcc);
2993
2994 if(ObjectMoleculeFindBestDonorH(don_obj,
2995 don_atom, don_state, donToAcc, bestH, h_real)) {
2996
2997 subtract3f(bestH, vDon, donToH);
2998 subtract3f(vAcc, bestH, hToAcc);
2999
3000 if(ObjectMoleculeGetAvgHBondVector(acc_obj, acc_atom,
3001 acc_state, accPlane, hToAcc) > 0.1) {
3002 vAccPlane = &accPlane[0];
3003 }
3004 result = ObjectMoleculeTestHBond(donToAcc, donToH, hToAcc, vAccPlane, hbc);
3005 if(result && h_crd_ret && h_real && *h_real)
3006 copy3f(bestH, h_crd_ret);
3007 }
3008 }
3009 }
3010
3011 return (result);
3012 }
3013
3014
3015 /*========================================================================*/
3016
ObjectMoleculeGetMaxVDW(ObjectMolecule * I)3017 float ObjectMoleculeGetMaxVDW(ObjectMolecule * I)
3018 {
3019 float max_vdw = 0.0F;
3020 int a;
3021 const AtomInfoType *ai;
3022 if(I->NAtom) {
3023 ai = I->AtomInfo;
3024 for(a = 0; a < I->NAtom; a++) {
3025 if(max_vdw < ai->vdw)
3026 max_vdw = ai->vdw;
3027 ai++;
3028 }
3029 }
3030 return (max_vdw);
3031 }
3032
3033
3034 /*========================================================================*/
ObjectMoleculeCSetAsPyList(ObjectMolecule * I)3035 static PyObject *ObjectMoleculeCSetAsPyList(ObjectMolecule * I)
3036 {
3037 PyObject *result = NULL;
3038 int a;
3039 result = PyList_New(I->NCSet);
3040 for(a = 0; a < I->NCSet; a++) {
3041 if(I->CSet[a]) {
3042 PyList_SetItem(result, a, CoordSetAsPyList(I->CSet[a]));
3043 } else {
3044 PyList_SetItem(result, a, PConvAutoNone(Py_None));
3045 }
3046 }
3047 return (PConvAutoNone(result));
3048 }
3049
3050
3051 /*static PyObject *ObjectMoleculeDiscreteCSetAsPyList(ObjectMolecule *I)
3052 {
3053 PyObject *result = NULL;
3054 return(PConvAutoNone(result));
3055 }*/
ObjectMoleculeCSetFromPyList(ObjectMolecule * I,PyObject * list)3056 static int ObjectMoleculeCSetFromPyList(ObjectMolecule * I, PyObject * list)
3057 {
3058 int ok = true;
3059 int a;
3060 if(ok)
3061 ok = PyList_Check(list);
3062 if(ok) {
3063 VLACheck(I->CSet, CoordSet *, I->NCSet);
3064 for(a = 0; a < I->NCSet; a++) {
3065 if(ok)
3066 ok = CoordSetFromPyList(I->G, PyList_GetItem(list, a), &I->CSet[a]);
3067 PRINTFB(I->G, FB_ObjectMolecule, FB_Debugging)
3068 " %s: ok %d after CoordSet %d\n", __func__, ok, a ENDFB(I->G);
3069
3070 if(ok)
3071 if(I->CSet[a]) /* WLD 030205 */
3072 I->CSet[a]->Obj = I;
3073 }
3074 }
3075 return (ok);
3076 }
3077
ObjectMoleculeBondAsPyList(ObjectMolecule * I)3078 static PyObject *ObjectMoleculeBondAsPyList(ObjectMolecule * I)
3079 {
3080 PyObject *result = NULL;
3081 PyObject *bond_list;
3082 const BondType *bond;
3083 int a;
3084
3085 #ifndef PICKLETOOLS
3086 PyMOLGlobals *G = I->G;
3087 int pse_export_version = SettingGetGlobal_f(I->G, cSetting_pse_export_version) * 1000;
3088
3089 if (SettingGetGlobal_b(G, cSetting_pse_binary_dump) && (!pse_export_version || pse_export_version >= 1765)){
3090 /* For the pse_binary_dump, save entire Bond array to a binary string array
3091 */
3092
3093 // supported versions
3094 auto version = (!pse_export_version || pse_export_version >= 1810) ?
3095 181 : (pse_export_version >= 1770) ? 177 : 176;
3096
3097 auto blob = Copy_To_BondType_Version(version, I->Bond.data(), I->NBond);
3098 auto blobsize = VLAGetByteSize(blob);
3099
3100 result = PyList_New(2);
3101 PyList_SetItem(result, 0, PyInt_FromLong(version));
3102 PyList_SetItem(result, 1, PyBytes_FromStringAndSize(reinterpret_cast<const char*>(blob), blobsize));
3103
3104 VLAFreeP(blob);
3105
3106 return result;
3107 }
3108 #endif
3109 result = PyList_New(I->NBond);
3110 bond = I->Bond;
3111 for(a = 0; a < I->NBond; a++) {
3112 bond_list = PyList_New(7);
3113 PyList_SetItem(bond_list, 0, PyInt_FromLong(bond->index[0]));
3114 PyList_SetItem(bond_list, 1, PyInt_FromLong(bond->index[1]));
3115 PyList_SetItem(bond_list, 2, PyInt_FromLong(bond->order));
3116 PyList_SetItem(bond_list, 3, PyInt_FromLong(bond->id));
3117 PyList_SetItem(bond_list, 4, PyInt_FromLong(bond->stereo));
3118 PyList_SetItem(bond_list, 5, PyInt_FromLong(bond->unique_id));
3119 PyList_SetItem(bond_list, 6, PyInt_FromLong(bond->has_setting));
3120 PyList_SetItem(result, a, bond_list);
3121 bond++;
3122 }
3123
3124 return (PConvAutoNone(result));
3125 }
3126
ObjectMoleculeBondFromPyList(ObjectMolecule * I,PyObject * list)3127 static int ObjectMoleculeBondFromPyList(ObjectMolecule * I, PyObject * list)
3128 {
3129 PyMOLGlobals *G = I->G;
3130 int ok = true;
3131 int a;
3132 int stereo, ll = 0;
3133 PyObject *bond_list = NULL;
3134 BondType *bond;
3135
3136 if(ok)
3137 ok = PyList_Check(list);
3138 if(ok)
3139 ll = PyList_Size(list);
3140
3141 bool pse_binary_dump = false;
3142
3143 if (ll >= 2) {
3144 // checking if from pse_binary_dump
3145 // pse_binary_dump saves 2 values: bondInfo_version, BondType binary
3146 CPythonVal *val1 = CPythonVal_PyList_GetItem(G, list, 1);
3147 pse_binary_dump = PyBytes_Check(val1);
3148 CPythonVal_Free(val1);
3149 }
3150 if (pse_binary_dump){
3151 CPythonVal *verobj = CPythonVal_PyList_GetItem(G, list, 0);
3152 int bondInfo_version;
3153 ok = PConvPyIntToInt(verobj, &bondInfo_version);
3154
3155 CPythonVal *strobj = CPythonVal_PyList_GetItem(G, list, 1);
3156 auto strval = PyBytes_AsSomeString(strobj);
3157
3158 if(ok)
3159 ok = bool((I->Bond = pymol::vla<BondType>(I->NBond)));
3160
3161 Copy_Into_BondType_From_Version(strval.data(), bondInfo_version, I->Bond.data(), I->NBond);
3162
3163 CPythonVal_Free(verobj);
3164 CPythonVal_Free(strobj);
3165 } else {
3166 if(ok)
3167 ok = bool((I->Bond = pymol::vla<BondType>(I->NBond)));
3168 bond = I->Bond.data();
3169 for(a = 0; a < I->NBond; a++) {
3170 bond_list = NULL;
3171 if(ok)
3172 bond_list = PyList_GetItem(list, a);
3173 if(ok)
3174 ok = PyList_Check(bond_list);
3175 if(ok)
3176 ll = PyList_Size(bond_list);
3177 if(ok)
3178 ok = PConvPyIntToInt(PyList_GetItem(bond_list, 0), &bond->index[0]);
3179 if(ok)
3180 ok = PConvPyIntToInt(PyList_GetItem(bond_list, 1), &bond->index[1]);
3181 if(ok)
3182 if((ok = CPythonVal_PConvPyIntToInt_From_List(I->G, bond_list, 2, &stereo)))
3183 bond->order = stereo;
3184 if(ok)
3185 ok = PConvPyIntToInt(PyList_GetItem(bond_list, 3), &bond->id);
3186 if(ok)
3187 ok = PConvPyIntToInt(PyList_GetItem(bond_list, 4), &stereo);
3188 if(ok)
3189 bond->stereo = (short int) stereo;
3190 if(ok && (ll > 5)) {
3191 int has_setting;
3192 if(ok)
3193 ok = PConvPyIntToInt(PyList_GetItem(bond_list, 5), &bond->unique_id);
3194 if(ok)
3195 ok = PConvPyIntToInt(PyList_GetItem(bond_list, 6), &has_setting);
3196 if(ok)
3197 bond->has_setting = (short int) has_setting;
3198 if(ok && bond->unique_id) { /* reserve existing IDs */
3199 bond->unique_id = SettingUniqueConvertOldSessionID(G, bond->unique_id);
3200 }
3201 }
3202 bond++;
3203 }
3204 }
3205 PRINTFB(G, FB_ObjectMolecule, FB_Debugging)
3206 " %s: ok %d after restore\n", __func__, ok ENDFB(G);
3207
3208 return (ok);
3209 }
3210
3211 /**
3212 * Extract an atom property "column" as a Python list.
3213 *
3214 * As an optimization, trailing None values are removed, so the returned list
3215 * may be shorter than the atom array or even be empty.
3216 *
3217 * @param mol Object molecule which provides the atom array
3218 * @param func Function which extracts a property from an atom
3219 * @return List of properties
3220 */
AtomColumnAsPyList(const ObjectMolecule & mol,std::function<PyObject * (const AtomInfoType &)> func)3221 static PyObject* AtomColumnAsPyList(const ObjectMolecule& mol,
3222 std::function<PyObject*(const AtomInfoType&)> func)
3223 {
3224 auto list = PyList_New(mol.NAtom);
3225 PyObject* prev = Py_None;
3226 int pruned_size = 0;
3227
3228 for (int i = 0; i < mol.NAtom; ++i) {
3229 PyObject* value = func(mol.AtomInfo[i]);
3230
3231 #ifndef PICKLETOOLS
3232 // Simple optimization: Repeated property lists will reference the same
3233 // Python object
3234 if (prev != value && PyObject_RichCompareBool(prev, value, Py_EQ)) {
3235 Py_INCREF(prev);
3236 Py_DECREF(value);
3237 value = prev;
3238 } else {
3239 prev = value;
3240 }
3241
3242 if (value != Py_None) {
3243 pruned_size = i + 1;
3244 }
3245 #endif
3246
3247 PyList_SetItem(list, i, value);
3248 }
3249
3250 #ifndef PICKLETOOLS
3251 assert(pruned_size <= mol.NAtom);
3252
3253 // Simple optimization: Prune "None" tail
3254 PyList_SetSlice(list, pruned_size, mol.NAtom, nullptr);
3255
3256 assert(PyList_Size(list) == pruned_size);
3257 #endif
3258
3259 return list;
3260 }
3261
3262 /**
3263 * Restore an atom property "column" from a Python list.
3264 * @param mol Object molecule to update atoms in
3265 * @param list Property list (may be shorter than atom array)
3266 * @param func Function which sets a property for an atom
3267 */
AtomColumnFromPyList(ObjectMolecule & mol,PyObject * list,std::function<void (AtomInfoType &,PyObject *)> func)3268 static void AtomColumnFromPyList(ObjectMolecule& mol, PyObject* list,
3269 std::function<void(AtomInfoType&, PyObject*)> func)
3270 {
3271 if (!list || !PyList_Check(list)) {
3272 return;
3273 }
3274
3275 int size = PyList_Size(list);
3276 assert(size <= mol.NAtom);
3277
3278 for (int i = 0; i < size; ++i) {
3279 func(mol.AtomInfo[i], PyList_GetItem(list, i));
3280 }
3281 }
3282
ObjectMoleculeAtomAsPyList(ObjectMolecule * I)3283 static PyObject *ObjectMoleculeAtomAsPyList(ObjectMolecule * I)
3284 {
3285 PyMOLGlobals *G = I->G;
3286 PyObject *result = NULL;
3287 const AtomInfoType *ai;
3288 int a;
3289 #ifndef PICKLETOOLS
3290 int pse_export_version = SettingGetGlobal_f(I->G, cSetting_pse_export_version) * 1000;
3291
3292 if (SettingGetGlobal_b(G, cSetting_pse_binary_dump) && (!pse_export_version || pse_export_version >= 1765)){
3293 /* For the pse_binary_dump, record all strings in lex and
3294 write them into separate binary string
3295 */
3296 AtomInfoTypeConverter converter(G, I->NAtom);
3297 std::set<lexidx_t> lexIDs;
3298 int totalstlen = 0;
3299 ai = I->AtomInfo.data();
3300 for(a = 0; a < I->NAtom; a++) {
3301 if (ai->textType) lexIDs.insert(ai->textType);
3302 if (ai->chain) lexIDs.insert(ai->chain);
3303 if (ai->label) lexIDs.insert(ai->label);
3304 if (ai->custom) lexIDs.insert(ai->custom);
3305 if (ai->segi) lexIDs.insert(ai->segi);
3306 if (ai->resn) lexIDs.insert(ai->resn);
3307 if (ai->name) lexIDs.insert(ai->name);
3308 ++ai;
3309 }
3310 for (const auto& lexID : lexIDs) { // need to calculate totalstlen so we can allocate
3311 const char *lexstr = LexStr(G, lexID);
3312 int lexlen = strlen(lexstr);
3313 totalstlen += lexlen + 1;
3314 }
3315 int strinfolen = totalstlen + sizeof(int) * (lexIDs.size() + 1);
3316 void *strinfo = pymol::malloc<unsigned char>(strinfolen);
3317 int *strval = (int*)strinfo;
3318 *(strval++) = lexIDs.size(); // first write number of strings into binary data string
3319 char *strpl = (char*)((char*)strinfo + (1 + lexIDs.size()) * sizeof(int));
3320 /* write map of lex ids and strings into binary data string as an array of ids
3321 and null-terminated strings */
3322 for (const auto& lexID : lexIDs) {
3323 *(strval++) = converter.to_lexidx_int(lexID);
3324 const char *strptr = LexStr(G, lexID);
3325 strcpy(strpl, strptr);
3326 strpl += strlen(strptr) + 1;
3327 }
3328
3329 auto version = AtomInfoVERSION;
3330 if (pse_export_version && pse_export_version < 1810) {
3331 if (pse_export_version < 1770) {
3332 version = 176;
3333 } else {
3334 version = 177;
3335 }
3336 }
3337
3338 auto blob = converter.allocCopy(version, I->AtomInfo);
3339 auto blobsize = VLAGetByteSize(blob);
3340
3341 // PyMOL versions up to 2.3.5 can only restore list size 3
3342 size_t result_size = 3;
3343
3344 // Atom properties (not binary)
3345 PyObject* prop_list = nullptr;
3346 #ifdef _PYMOL_IP_PROPERTIES
3347 if (pse_export_version > 2399) {
3348 prop_list = AtomColumnAsPyList(*I, [G](const AtomInfoType& atom) {
3349 return PConvAutoNone(
3350 atom.prop_id ? PropertyAsPyList(G, atom.prop_id, true) : nullptr);
3351 });
3352
3353 result_size = 4;
3354 }
3355 #endif
3356
3357 result = PyList_New(result_size);
3358 PyList_SetItem(result, 0, PyInt_FromLong(version));
3359 PyList_SetItem(result, 1, PyBytes_FromStringAndSize(reinterpret_cast<const char*>(blob), blobsize));
3360 PyList_SetItem(result, 2, PyBytes_FromStringAndSize(reinterpret_cast<const char*>(strinfo), strinfolen));
3361
3362 if (result_size > 3) {
3363 PyList_SetItem(result, 3, PConvAutoNone(prop_list));
3364 }
3365
3366 VLAFreeP(blob);
3367 FreeP(strinfo);
3368 return result;
3369 }
3370 #endif
3371 result = PyList_New(I->NAtom);
3372 ai = I->AtomInfo;
3373 for(a = 0; a < I->NAtom; a++) {
3374 PyList_SetItem(result, a, AtomInfoAsPyList(I->G, ai));
3375 ai++;
3376 }
3377 return (PConvAutoNone(result));
3378 }
3379
ObjectMoleculeAtomFromPyList(ObjectMolecule * I,PyObject * list)3380 static int ObjectMoleculeAtomFromPyList(ObjectMolecule * I, PyObject * list)
3381 {
3382 PyMOLGlobals *G = I->G;
3383 int ok = true;
3384 int a, ll = 0;
3385 AtomInfoType *ai;
3386
3387 if(ok)
3388 ok = PyList_Check(list);
3389 if (ok)
3390 ll = PyList_Size(list);
3391
3392 bool pse_binary_dump = false;
3393
3394 if (ll >= 3) {
3395 // checking if from pse_binary_dump
3396 // pse_binary_dump saves 3 values: atomInfo_version, AtomInfo binary, and strings array
3397 CPythonVal *val1 = CPythonVal_PyList_GetItem(G, list, 1);
3398 CPythonVal *val2 = CPythonVal_PyList_GetItem(G, list, 2);
3399 pse_binary_dump = PyBytes_Check(val1) && PyBytes_Check(val2);
3400 CPythonVal_Free(val1);
3401 CPythonVal_Free(val2);
3402 }
3403 if (pse_binary_dump){
3404 CPythonVal *verobj = CPythonVal_PyList_GetItem(G, list, 0);
3405 int atomInfo_version;
3406 ok = PConvPyIntToInt(verobj, &atomInfo_version);
3407
3408 CPythonVal *strlookupobj = CPythonVal_PyList_GetItem(G, list, 2);
3409 auto strval_1 = PyBytes_AsSomeString(strlookupobj);
3410 int *strval = (int*)strval_1.data();
3411
3412 AtomInfoTypeConverter converter(G, I->NAtom);
3413
3414 auto& oldIDtoLexID = converter.lexidxmap;
3415 int nstrings = *(strval++);
3416 char *strpl = (char*)(strval + nstrings);
3417 int strcnt = nstrings;
3418 int stlen;
3419 // populate oldIDtoLexID with nstrings from binary string data (3rd entry in list)
3420 while (strcnt){
3421 lexidx_t idx = LexIdx(G, strpl); // increments ref count, need to take into account
3422 int oldidx = *(strval++);
3423 oldIDtoLexID[oldidx] = idx;
3424 stlen = strlen(strpl);
3425 strpl += stlen + 1;
3426 strcnt--;
3427 }
3428
3429 CPythonVal *strobj = CPythonVal_PyList_GetItem(G, list, 1);
3430 auto strval_2 = PyBytes_AsSomeString(strobj);
3431
3432 VLACheck(I->AtomInfo, AtomInfoType, I->NAtom + 1);
3433 converter.copy(I->AtomInfo.data(), strval_2.data(), atomInfo_version);
3434
3435 // go through AtomInfo array, swap new strings, convert colors, convert settings
3436 // (everything that AtomInfoFromPyList does except set properties, which are currently
3437 // not saved for pse_binary_dump)
3438 AtomInfoType *ai = I->AtomInfo.data();
3439 for(a = 0; a < I->NAtom; ++a, ++ai) {
3440 ai->color = ColorConvertOldSessionIndex(G, ai->color);
3441 if (ai->unique_id){
3442 ai->unique_id = SettingUniqueConvertOldSessionID(G, ai->unique_id);
3443 }
3444 }
3445 // need to decrement since we call LexIdx() above on each
3446 for (auto it = oldIDtoLexID.begin(); it != oldIDtoLexID.end(); ++it){
3447 LexDec(G, it->second);
3448 }
3449 CPythonVal_Free(verobj);
3450 CPythonVal_Free(strobj);
3451 CPythonVal_Free(strlookupobj);
3452
3453 #ifdef _PYMOL_IP_PROPERTIES
3454 if (ll > 3) {
3455 // Restore atom properties
3456 AtomColumnFromPyList(*I, PyList_GetItem(list, 3), //
3457 [G](AtomInfoType& atom, PyObject* value) {
3458 assert(atom.prop_id == 0);
3459 atom.prop_id = PropertyFromPyList(G, value);
3460 });
3461 }
3462 #endif
3463
3464 } else {
3465 // The old slow way of loading in AtomInfo, using python lists
3466 if (ok)
3467 VLACheck(I->AtomInfo, AtomInfoType, I->NAtom + 1);
3468 CHECKOK(ok, I->AtomInfo);
3469 ai = I->AtomInfo.data();
3470 for(a = 0; ok && a < I->NAtom; a++) {
3471 PyObject *val = PyList_GetItem(list, a);
3472 ok &= AtomInfoFromPyList(I->G, ai, val);
3473 ai++;
3474 }
3475 }
3476 PRINTFB(I->G, FB_ObjectMolecule, FB_Debugging)
3477 " %s: ok %d \n", __func__, ok ENDFB(I->G);
3478 return (ok);
3479 }
3480
ObjectMoleculeNewFromPyList(PyMOLGlobals * G,PyObject * list,ObjectMolecule ** result)3481 int ObjectMoleculeNewFromPyList(PyMOLGlobals * G, PyObject * list,
3482 ObjectMolecule ** result)
3483 {
3484 int ok = true;
3485 ObjectMolecule *I = NULL;
3486 int discrete_flag = 0;
3487 (*result) = NULL;
3488
3489 if(ok)
3490 ok = PyList_Check(list);
3491 /* TO SUPPORT BACKWARDS COMPATIBILITY...
3492 Always check ll when adding new PyList_GetItem's */
3493 if(ok)
3494 ok = PConvPyIntToInt(PyList_GetItem(list, 8), &discrete_flag);
3495
3496 if (ok)
3497 I = new ObjectMolecule(G, discrete_flag);
3498 CHECKOK(ok, I);
3499
3500 if(ok){
3501 PyObject *val = PyList_GetItem(list, 0);
3502 ok = ObjectFromPyList(G, val, I);
3503 }
3504 if(ok)
3505 ok = PConvPyIntToInt(PyList_GetItem(list, 1), &I->NCSet);
3506 if(ok)
3507 ok = PConvPyIntToInt(PyList_GetItem(list, 2), &I->NBond);
3508 if(ok)
3509 ok = PConvPyIntToInt(PyList_GetItem(list, 3), &I->NAtom);
3510 if(ok)
3511 ok = ObjectMoleculeCSetFromPyList(I, PyList_GetItem(list, 4));
3512 if(ok){
3513 ok = CoordSetFromPyList(G, PyList_GetItem(list, 5), &I->CSTmpl);
3514
3515 if(I->CSTmpl)
3516 I->CSTmpl->Obj = I;
3517 }
3518 if(ok){
3519 CPythonVal *val = CPythonVal_PyList_GetItem(G, list, 6);
3520 ok = ObjectMoleculeBondFromPyList(I, val);
3521 CPythonVal_Free(val);
3522 }
3523 if (!ok && I)
3524 I->NBond = 0;
3525 if(ok){
3526 CPythonVal *val = CPythonVal_PyList_GetItem(G, list, 7);
3527 ok = ObjectMoleculeAtomFromPyList(I, val);
3528 CPythonVal_Free(val);
3529 }
3530 if (!ok && I)
3531 I->NAtom = 0;
3532 if(ok){
3533 CPythonVal *val = CPythonVal_PyList_GetItem(G, list, 10);
3534 I->Symmetry = SymmetryNewFromPyList(G, val);
3535 CPythonVal_Free(val);
3536 }
3537 /* 11 was CurCSet */
3538 if(ok)
3539 ok = PConvPyIntToInt(PyList_GetItem(list, 12), &I->BondCounter);
3540 if(ok)
3541 ok = PConvPyIntToInt(PyList_GetItem(list, 13), &I->AtomCounter);
3542
3543 I->updateAtmToIdx();
3544
3545 if (ok)
3546 I->invalidate(cRepAll, cRepInvAll, -1);
3547 if(ok)
3548 (*result) = I;
3549 else {
3550 /* cleanup */
3551 if (I)
3552 DeleteP(I);
3553 (*result) = NULL;
3554 }
3555 return (ok);
3556 }
3557
3558
3559 /*========================================================================*/
ObjectMoleculeAsPyList(ObjectMolecule * I)3560 PyObject *ObjectMoleculeAsPyList(ObjectMolecule * I)
3561 {
3562 PyObject *result = NULL;
3563
3564 /* first, dump the atoms */
3565
3566 result = PyList_New(16);
3567 PyList_SetItem(result, 0, ObjectAsPyList(I));
3568 PyList_SetItem(result, 1, PyInt_FromLong(I->NCSet));
3569 PyList_SetItem(result, 2, PyInt_FromLong(I->NBond));
3570 PyList_SetItem(result, 3, PyInt_FromLong(I->NAtom));
3571 PyList_SetItem(result, 4, ObjectMoleculeCSetAsPyList(I));
3572 PyList_SetItem(result, 5, CoordSetAsPyList(I->CSTmpl));
3573 PyList_SetItem(result, 6, ObjectMoleculeBondAsPyList(I));
3574 PyList_SetItem(result, 7, ObjectMoleculeAtomAsPyList(I));
3575 PyList_SetItem(result, 8, PyInt_FromLong(I->DiscreteFlag));
3576 PyList_SetItem(result, 9, PyInt_FromLong(I->DiscreteFlag ? I->NAtom : 0 /* NDiscrete */));
3577 PyList_SetItem(result, 10, SymmetryAsPyList(I->Symmetry));
3578 PyList_SetItem(result, 11, PyInt_FromLong(0 /* CurCSet */));
3579 PyList_SetItem(result, 12, PyInt_FromLong(I->BondCounter));
3580 PyList_SetItem(result, 13, PyInt_FromLong(I->AtomCounter));
3581
3582 float pse_export_version = SettingGetGlobal_f(I->G, cSetting_pse_export_version);
3583
3584 if(I->DiscreteFlag
3585 && (pse_export_version || !SettingGetGlobal_b(I->G, cSetting_pse_binary_dump))
3586 && pse_export_version < 1.7699) {
3587 int *dcs;
3588 int a;
3589 CoordSet *cs;
3590
3591 /* get coordinate set indices */
3592
3593 for(a = 0; a < I->NCSet; a++) {
3594 cs = I->CSet[a];
3595 if(cs) {
3596 cs->tmp_index = a;
3597 }
3598 }
3599
3600 dcs = pymol::malloc<int>(I->NAtom);
3601
3602 for(a = 0; a < I->NAtom; a++) {
3603 cs = I->DiscreteCSet[a];
3604 if(cs)
3605 dcs[a] = cs->tmp_index;
3606 else
3607 dcs[a] = -1;
3608 }
3609
3610 PyList_SetItem(result, 14, PConvIntArrayToPyList(I->DiscreteAtmToIdx, I->NAtom));
3611 PyList_SetItem(result, 15, PConvIntArrayToPyList(dcs, I->NAtom));
3612 FreeP(dcs);
3613 } else {
3614 PyList_SetItem(result, 14, PConvAutoNone(NULL));
3615 PyList_SetItem(result, 15, PConvAutoNone(NULL));
3616 }
3617
3618 return (PConvAutoNone(result));
3619 }
3620
3621 /*========================================================================*/
3622
3623 static
connect_cutoff_adjustment(const AtomInfoType * ai1,const AtomInfoType * ai2)3624 float connect_cutoff_adjustment(
3625 const AtomInfoType * ai1,
3626 const AtomInfoType * ai2)
3627 {
3628 if (ai1->isHydrogen() || ai2->isHydrogen())
3629 return -0.2f;
3630
3631 if (ai1->protons == cAN_S || ai2->protons == cAN_S)
3632 return 0.2f;
3633
3634 return 0.f;
3635 }
3636
3637 /*
3638 * True if two atoms should be bonded
3639 */
3640 static
is_distance_bonded(PyMOLGlobals * G,const CoordSet * cs,const AtomInfoType * ai1,const AtomInfoType * ai2,const float * v1,const float * v2,float cutoff,int connect_mode,int discrete_chains,bool connect_bonded,bool unbond_cations)3641 bool is_distance_bonded(
3642 PyMOLGlobals * G,
3643 const CoordSet * cs,
3644 const AtomInfoType * ai1,
3645 const AtomInfoType * ai2,
3646 const float * v1,
3647 const float * v2,
3648 float cutoff,
3649 int connect_mode,
3650 int discrete_chains,
3651 bool connect_bonded,
3652 bool unbond_cations)
3653 {
3654 auto dst = (float) diff3f(v1, v2);
3655 dst -= (ai1->vdw + ai2->vdw) / 2;
3656
3657 cutoff += connect_cutoff_adjustment(ai1, ai2);
3658
3659 if (dst > cutoff)
3660 return false;
3661
3662 if (ai1->isHydrogen() && ai2->isHydrogen())
3663 return false;
3664
3665 if (discrete_chains > 0 && ai1->chain != ai2->chain)
3666 return false;
3667
3668 if (!connect_bonded && ai1->bonded && ai2->bonded)
3669 return false;
3670
3671 bool water_flag = (
3672 AtomInfoKnownWaterResName(G, LexStr(G, ai1->resn)) ||
3673 AtomInfoKnownWaterResName(G, LexStr(G, ai2->resn)));
3674
3675 if (connect_mode != 3 &&
3676 cs->TmpBond && /* connectivity information present in file */
3677 ai1->hetatm &&
3678 ai2->hetatm &&
3679 !water_flag &&
3680 !(AtomInfoKnownPolymerResName(LexStr(G, ai1->resn)) &&
3681 AtomInfoKnownPolymerResName(LexStr(G, ai2->resn))))
3682 return false;
3683
3684 // don't connect water atoms in different residues
3685 if (water_flag && !AtomInfoSameResidue(G, ai1, ai2))
3686 return false;
3687
3688 // don't connect atoms with different, non-NULL alternate conformations
3689 if (ai1->alt[0] != ai2->alt[0] && ai1->alt[0] && ai2->alt[0])
3690 return false;
3691
3692 // if either is a cation, unbond is user wants
3693 if (unbond_cations &&
3694 (AtomInfoIsFreeCation(G, ai1) ||
3695 AtomInfoIsFreeCation(G, ai2)))
3696 return false;
3697
3698 return true;
3699 }
3700
3701 /**
3702 * Do bonding of atoms in `I`, using distances and/or temporary bonds in `cs`.
3703 *
3704 * Incorporates `cs->TmpBond` unless `connect_mode` is 2.
3705 * Incorporates `cs->TmpLinkBond`.
3706 *
3707 * @param I Molecule to modify
3708 * @param cs Coordinates and temporary bonds to consider
3709 * @param bondSearchMode If false and `connect_mode` != 2, do not search for new
3710 * bonds (only use TmpBond/TmpLinkBond).
3711 * @param connectModeOverride Overrides `connect_mode` setting if not -1
3712 *
3713 * `connect_mode` options:
3714 * 0 = distance-based (excluding HETATM to HETATM) and CONECT records (default)
3715 * 1 = CONECT records
3716 * 2 = distance-based, ignores CONECT records
3717 * 3 = distance-based (including HETATM to HETATM) and CONECT records
3718 * 4 = same as `connect_mode` = 0 (special meaning during mmCIF loading)
3719 */
ObjectMoleculeConnect(ObjectMolecule * I,CoordSet * cs,bool bondSearchMode,int connectModeOverride)3720 bool ObjectMoleculeConnect(ObjectMolecule* I, CoordSet* cs, bool bondSearchMode,
3721 int connectModeOverride)
3722 {
3723 return ObjectMoleculeConnect(
3724 I, I->NBond, I->Bond, cs, bondSearchMode, connectModeOverride);
3725 }
3726
3727 /*========================================================================*/
ObjectMoleculeConnect(ObjectMolecule * I,int & nBond,pymol::vla<BondType> & bondvla,struct CoordSet * cs,int bondSearchMode,int connectModeOverride)3728 bool ObjectMoleculeConnect(ObjectMolecule* I, int& nBond, pymol::vla<BondType>& bondvla,
3729 struct CoordSet *cs, int bondSearchMode,
3730 int connectModeOverride)
3731 {
3732 #define cMULT 1
3733 PyMOLGlobals *G = I->G;
3734 int a, b, c, d, e, f, i, j;
3735 int a1, a2;
3736 float *v1, *v2;
3737 int maxBond;
3738 MapType *map;
3739 BondType *ii1;
3740 const BondType* ii2;
3741 int flag;
3742 int order;
3743 AtomInfoType* const ai = I->AtomInfo.data();
3744 AtomInfoType *ai1, *ai2;
3745 /* Sulfur cutoff */
3746 float cutoff_s;
3747 float cutoff_v;
3748 float max_cutoff;
3749 int repeat = true;
3750 int discrete_chains = SettingGetGlobal_i(G, cSetting_pdb_discrete_chains);
3751 int connect_bonded = SettingGetGlobal_b(G, cSetting_connect_bonded);
3752 int connect_mode = SettingGetGlobal_i(G, cSetting_connect_mode);
3753 int unbond_cations = SettingGetGlobal_i(G, cSetting_pdb_unbond_cations);
3754 int ok = true;
3755 cutoff_v = SettingGetGlobal_f(G, cSetting_connect_cutoff);
3756 cutoff_s = cutoff_v + 0.2F;
3757 max_cutoff = cutoff_s;
3758
3759 if(connectModeOverride >= 0)
3760 connect_mode = connectModeOverride;
3761
3762 if(connect_mode == 2) { /* force use of distance-based connectivity,
3763 ignoring that provided with file */
3764 bondSearchMode = true;
3765 cs->NTmpBond = 0;
3766 VLAFreeP(cs->TmpBond);
3767 } else if (connect_mode == 4) {
3768 // mmCIF specific, fall back to default to get any bonds for PDB, XYZ, etc.
3769 connect_mode = 0;
3770 }
3771
3772 /* FeedbackMask[FB_ObjectMolecule]=0xFF; */
3773 nBond = 0;
3774 maxBond = cs->NIndex * 8;
3775 bondvla.resize(maxBond);
3776 CHECKOK(ok, bool(bondvla));
3777 while(ok && repeat) {
3778 repeat = false;
3779
3780 /*
3781 * BOND SEARCH MODE
3782 */
3783 if(cs->NIndex && bondSearchMode) { /* &&(!I->DiscreteFlag) WLD 010527 */
3784 /* if there are atoms, and we need to search for bonds, instead of using
3785 * (possibly) supplied CONECT records... */
3786
3787 PRINTFB(G, FB_ObjectMolecule, FB_Blather)
3788 " %s: Searching for bonds amongst %d coordinates.\n", __func__,
3789 cs->NIndex ENDFB(G);
3790 if(Feedback(G, FB_ObjectMolecule, FB_Debugging)) {
3791 for(a = 0; a < cs->NIndex; a++)
3792 printf(" %s: coord %d %8.3f %8.3f %8.3f\n", __func__,
3793 a, cs->Coord[a * 3], cs->Coord[a * 3 + 1], cs->Coord[a * 3 + 2]);
3794 }
3795
3796 switch (connect_mode) {
3797 /* DISTANCE BASED CONNECTIONS */
3798 case 0: /* distance-based and explicit (not HETATM to HETATM) */
3799 case 3: /* distance-based and explicit (even HETATM to HETATM) */
3800 case 2: /* distance-based only */ {
3801 /* distance-based bond location */
3802 int violations = 0;
3803 int *cnt = pymol::malloc<int>(cs->NIndex);
3804 int valcnt;
3805
3806 CHECKOK(ok, cnt);
3807
3808 if (ok){
3809 /* initialize each atom's (max) expected valence */
3810 for(i = 0; i < cs->NIndex; i++) {
3811 valcnt = AtomInfoGetExpectedValence(G, ai + cs->IdxToAtm[i]);
3812 if(valcnt < 0)
3813 valcnt = 6;
3814 cnt[i] = valcnt;
3815 }
3816 }
3817
3818 /* make a map of the local neighborhood in space */
3819 if (ok)
3820 map = MapNew(G, max_cutoff + MAX_VDW, cs->Coord, cs->NIndex, NULL);
3821 CHECKOK(ok, map);
3822 if(ok) {
3823 int dim12 = map->D1D2;
3824 int dim2 = map->Dim[2];
3825
3826 for(i = 0; ok && i < cs->NIndex; i++) {
3827 if(nBond > maxBond)
3828 break;
3829 /* atom i's position in space */
3830 v1 = cs->coordPtr(i);
3831
3832 a1 = cs->IdxToAtm[i];
3833 ai1 = ai + a1;
3834
3835 MapLocus(map, v1, &a, &b, &c);
3836 /* d = [a-1, a, a+1] */
3837 for(d = a - 1; ok && d <= a + 1; d++) {
3838 int *j_ptr1 = map->Head + d * dim12 + (b - 1) * dim2;
3839 /* e = [b-1, b, b+1] */
3840 for(e = b - 1; ok && e <= b + 1; e++) {
3841 int *j_ptr2 = j_ptr1 + c - 1;
3842 j_ptr1 += dim2;
3843 /* f = [c-1, c, c+1] */
3844 for(f = c - 1; ok && f <= c + 1; f++) {
3845 j = *(j_ptr2++); /* *MapFirst(map,d,e,f)); */
3846 while(ok && j >= 0) {
3847 if(i < j) {
3848 /* position in space for atom 2 */
3849 v2 = cs->coordPtr(j);
3850 a2 = cs->IdxToAtm[j];
3851 ai2 = ai + a2;
3852
3853 flag = is_distance_bonded(G, cs, ai1, ai2, v1, v2,
3854 cutoff_v, connect_mode, discrete_chains,
3855 connect_bonded, unbond_cations);
3856
3857 {
3858
3859 /* we have a bond, now process it */
3860 if(flag) {
3861 auto bnd = bondvla.check(nBond);
3862 CHECKOK(ok, bool(bondvla));
3863 if (ok){
3864 BondTypeInit2(bnd, a1, a2);
3865 order = 1;
3866 }
3867 /* if we allow bonds between chains and it screws up the
3868 * bonding, disallow inter-chain bonds */
3869 if(ok && discrete_chains < 0) { /* if we're allowing bonds between chains,
3870 then make sure things don't get out of hand */
3871 if(cnt[i] == -1)
3872 violations++;
3873 if(cnt[j] == -1)
3874 violations++;
3875 /* decrement free valences, since we have a bond */
3876 cnt[i]--;
3877 cnt[j]--;
3878 if(violations > (cs->NIndex >> 3)) {
3879 /* if more than 12% of the structure has excessive #'s of bonds... */
3880 PRINTFB(G, FB_ObjectMolecule, FB_Blather)
3881 " %s: Assuming chains are discrete...\n", __func__
3882 ENDFB(G);
3883 discrete_chains = 1;
3884 repeat = true;
3885 goto do_it_again;
3886 }
3887 }
3888
3889 if(!ai1->hetatm || ai1->resn == G->lex_const.MSE) {
3890 if(AtomInfoSameResidue(I->G, ai1, ai2)) {
3891 /* hookup standard disconnected PDB residue */
3892 assign_pdb_known_residue(G, ai1, ai2, &order);
3893 }
3894 }
3895 bnd->order = -order; /* store tentative valence as negative */
3896 nBond++;
3897 }
3898 }
3899 }
3900 j = MapNext(map, j);
3901 }
3902 }
3903 }
3904 }
3905 }
3906 do_it_again:
3907 MapFree(map);
3908 FreeP(cnt);
3909 }
3910 }
3911 case 1: /* only use explicit connectivity from file (don't do anything) */
3912 break;
3913 case 4: /* dictionary-based connectivity */
3914 /* TODO */
3915 break;
3916 }
3917 PRINTFB(G, FB_ObjectMolecule, FB_Blather)
3918 " %s: Found %d bonds.\n", __func__, nBond ENDFB(G);
3919 }
3920 if(repeat) {
3921 nBond = 0;
3922 }
3923 }
3924 /* if we have explicit connectivity, determine if we need to set check_conect_all */
3925 if(ok && cs->NTmpBond && cs->TmpBond) {
3926 int check_conect_all = false;
3927 int pdb_conect_all = false;
3928 PRINTFB(G, FB_ObjectMolecule, FB_Blather)
3929 " %s: incorporating explicit bonds. %d %d\n", __func__,
3930 nBond, cs->NTmpBond ENDFB(G);
3931 if((nBond == 0) && (cs->NTmpBond > 0) &&
3932 bondSearchMode && (connect_mode == 0) && cs->NIndex) {
3933 /* if no bonds were found, and we have explicit connectivity,
3934 * try to determine if we need to set pdb_conect_mode */
3935 for(i = 0; i < cs->NIndex; i++) {
3936 a1 = cs->IdxToAtm[i];
3937 ai1 = ai + a1;
3938 if(ai1->bonded && (!ai1->hetatm)) {
3939 /* apparent PDB ATOM record with explicit bonding... */
3940 check_conect_all = true;
3941 break;
3942 }
3943 }
3944 }
3945
3946 bondvla.check(nBond + cs->NTmpBond);
3947 CHECKOK(ok, bool(bondvla));
3948 if (ok){
3949 ii1 = bondvla.data() + nBond;
3950 ii2 = cs->TmpBond;
3951 }
3952 if (ok) {
3953 int n_atom = I->NAtom;
3954 for(a = 0; a < cs->NTmpBond; a++) {
3955 a1 = cs->IdxToAtm[ii2->index[0]]; /* convert bonds from index space */
3956 a2 = cs->IdxToAtm[ii2->index[1]]; /* to atom space */
3957 if((a1 >= 0) && (a2 >= 0) && (a1 < n_atom) && (a2 < n_atom)) {
3958 if(check_conect_all) {
3959 if((!ai[a1].hetatm) && (!ai[a2].hetatm)) {
3960 /* found bond between non-HETATMs -- so tell PyMOL to CONECT all ATOMs
3961 * when writing out a PDB file */
3962 pdb_conect_all = true;
3963 }
3964 }
3965 ai[a1].bonded = true;
3966 ai[a2].bonded = true;
3967 (*ii1) = (*ii2); /* note this copies owned ids and thus settings etc. */
3968 ii1->index[0] = a1;
3969 ii1->index[1] = a2;
3970 ii1++;
3971 ii2++;
3972
3973 }
3974 }
3975 }
3976
3977 if (ok){
3978 nBond = nBond + cs->NTmpBond;
3979 VLAFreeP(cs->TmpBond);
3980 cs->NTmpBond = 0;
3981
3982 if(pdb_conect_all) {
3983 int dummy;
3984 if(!SettingGetIfDefined_b(G, I->Setting, cSetting_pdb_conect_all, &dummy)) {
3985 {
3986 CSetting** handle = I->getSettingHandle(-1);
3987 if(handle) {
3988 SettingCheckHandle(G, handle);
3989 SettingSet_b(*handle, cSetting_pdb_conect_all, true);
3990 }
3991 }
3992 }
3993 }
3994 }
3995 }
3996
3997 /* Link b/t ligand and residue? */
3998 if(ok && cs->NTmpLinkBond && cs->TmpLinkBond) {
3999 PRINTFB(G, FB_ObjectMolecule, FB_Blather)
4000 "%s: incorporating linkage bonds. %d %d\n", __func__,
4001 nBond, cs->NTmpLinkBond ENDFB(G);
4002 bondvla.check(nBond + cs->NTmpLinkBond);
4003 CHECKOK(ok, bool(bondvla));
4004 if (ok){
4005 ii1 = bondvla.data() + nBond;
4006 ii2 = cs->TmpLinkBond;
4007 for(a = 0; a < cs->NTmpLinkBond; a++) {
4008 a1 = ii2->index[0]; /* first atom is in object */
4009 a2 = cs->IdxToAtm[ii2->index[1]]; /* second is in the cset */
4010 ai[a1].bonded = true;
4011 ai[a2].bonded = true;
4012 (*ii1) = (*ii2); /* note this copies owned ids and thus settings etc. */
4013 ii1->index[0] = a1;
4014 ii1->index[1] = a2;
4015 ii1++;
4016 ii2++;
4017 }
4018 nBond = nBond + cs->NTmpLinkBond;
4019 }
4020 VLAFreeP(cs->TmpLinkBond);
4021 cs->NTmpLinkBond = 0;
4022 }
4023
4024 PRINTFD(G, FB_ObjectMolecule)
4025 " %s: elminating duplicates with %d bonds...\n", __func__, nBond ENDFD;
4026
4027 if(ok && !I->DiscreteFlag) {
4028 UtilSortInPlace(G, bondvla.data(), nBond, sizeof(BondType), (UtilOrderFn *) BondInOrder);
4029 if(nBond) { /* eliminate duplicates */
4030 ii1 = bondvla.data() + 1;
4031 ii2 = bondvla.data() + 1;
4032 a = nBond - 1;
4033 nBond = 1;
4034 if(a > 0)
4035 while(a--) {
4036 if((ii2->index[0] != (ii1 - 1)->index[0]) ||
4037 (ii2->index[1] != (ii1 - 1)->index[1])) {
4038 *(ii1++) = *(ii2++); /* copy bond */
4039 nBond++;
4040 } else {
4041 if((ii2->order > 0) && ((ii1 - 1)->order < 0))
4042 (ii1 - 1)->order = ii2->order; /* use most certain valence */
4043 ii2++; /* skip bond */
4044 }
4045 }
4046 bondvla.resize(nBond);
4047 CHECKOK(ok, bool(bondvla));
4048 }
4049 }
4050 /* restore bond oder positivity */
4051
4052 if (ok){
4053 ii1 = bondvla.data();
4054 for(a = 0; a < nBond; a++) {
4055 if(ii1->order < 0)
4056 ii1->order = -ii1->order;
4057 ii1++;
4058 }
4059 }
4060
4061 PRINTFD(G, FB_ObjectMolecule)
4062 " %s: leaving with %d bonds...\n", __func__, nBond ENDFD;
4063 return ok;
4064 }
4065
4066
4067 /*========================================================================*/
ObjectMoleculeSort(ObjectMolecule * I)4068 int ObjectMoleculeSort(ObjectMolecule * I)
4069 { /* sorts atoms and bonds */
4070 int *index;
4071 int *outdex = NULL;
4072 int a, b;
4073 CoordSet *cs;
4074 int ok = true;
4075 if(!I->DiscreteFlag) { /* currently, discrete objects are never sorted */
4076 int n_bytes = sizeof(int) * I->NAtom;
4077 int already_in_order = true;
4078 int i_NAtom = I->NAtom;
4079 index = AtomInfoGetSortedIndex(I->G, I, I->AtomInfo, i_NAtom, &outdex);
4080 CHECKOK(ok, index);
4081 if (ok){
4082 for(a = 0; a < i_NAtom; a++) {
4083 if(index[a] != a) {
4084 already_in_order = false;
4085 break;
4086 }
4087 }
4088 }
4089 if(ok && !already_in_order) { /* if we aren't already in perfect order */
4090
4091 for(a = 0; a < I->NBond; a++) { /* bonds */
4092 I->Bond[a].index[0] = outdex[I->Bond[a].index[0]];
4093 I->Bond[a].index[1] = outdex[I->Bond[a].index[1]];
4094 }
4095
4096 for(a = -1; a < I->NCSet; a++) { /* coordinate set mapping */
4097 if(a < 0) {
4098 cs = I->CSTmpl;
4099 } else {
4100 cs = I->CSet[a];
4101 }
4102
4103 if(cs) {
4104 int cs_NIndex = cs->NIndex;
4105 int *cs_IdxToAtm = cs->IdxToAtm.data();
4106 int *cs_AtmToIdx = cs->AtmToIdx.data();
4107 for(b = 0; b < cs_NIndex; b++)
4108 cs_IdxToAtm[b] = outdex[cs_IdxToAtm[b]];
4109 if(cs_AtmToIdx) {
4110 memset(cs_AtmToIdx, -1, n_bytes);
4111 /* for(b=0;b<i_NAtom;b++)
4112 cs_AtmToIdx[b]=-1; */
4113 for(b = 0; b < cs_NIndex; b++) {
4114 cs_AtmToIdx[cs_IdxToAtm[b]] = b;
4115 }
4116 }
4117 }
4118 }
4119
4120 ExecutiveUniqueIDAtomDictInvalidate(I->G);
4121
4122 pymol::vla<AtomInfoType> atInfo(i_NAtom);
4123 CHECKOK(ok, atInfo);
4124 if (ok){
4125 /* autozero here is important */
4126 for(a = 0; a < i_NAtom; a++)
4127 atInfo[a] = std::move(I->AtomInfo[index[a]]);
4128 }
4129 VLAFreeP(I->AtomInfo);
4130 std::swap(I->AtomInfo, atInfo);
4131 }
4132 AtomInfoFreeSortedIndexes(I->G, &index, &outdex);
4133 if (ok){
4134 UtilSortInPlace(I->G, I->Bond.data(), I->NBond, sizeof(BondType),
4135 (UtilOrderFn *) BondInOrder);
4136 /* sort...important! */
4137 I->invalidate(cRepAll, cRepInvAtoms, -1); /* important */
4138 }
4139 }
4140 return ok;
4141 }
4142
ObjectMoleculeSetGeometry(PyMOLGlobals * G,ObjectMolecule * I,int sele,int geom,int valence)4143 int ObjectMoleculeSetGeometry(
4144 PyMOLGlobals* G, ObjectMolecule* I, int sele, int geom, int valence)
4145 {
4146 int count = 0;
4147 for (int a = 0; a < I->NAtom; ++a) {
4148 auto s = I->AtomInfo[a].selEntry;
4149
4150 if(SelectorIsMember(G, s, sele)) {
4151 auto& ai = I->AtomInfo[a];
4152 ai.geom = geom;
4153 ai.valence = valence;
4154 ai.chemFlag = true;
4155 count++;
4156 break;
4157 }
4158 }
4159 return count;
4160 }
4161
4162