1 
2 /*
3 A* -------------------------------------------------------------------
4 B* This file contains source code for the PyMOL computer program
5 C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific.
6 D* -------------------------------------------------------------------
7 E* It is unlawful to modify or remove this copyright notice.
8 F* -------------------------------------------------------------------
9 G* Please see the accompanying LICENSE file for further information.
10 H* -------------------------------------------------------------------
11 I* Additional authors of this source file include:
12 -*
13 -*
14 -*
15 Z* -------------------------------------------------------------------
16 */
17 #include"os_python.h"
18 #include"os_numpy.h"
19 
20 #include"os_std.h"
21 
22 #include <algorithm>
23 
24 #include"Base.h"
25 #include"OOMac.h"
26 #include"MemoryDebug.h"
27 #include"Err.h"
28 #include"Scene.h"
29 #include"CoordSet.h"
30 #include"Color.h"
31 #include"PConv.h"
32 #include"P.h"
33 #include"Matrix.h"
34 #include"Sphere.h"
35 #include"Util.h"
36 #include"Feedback.h"
37 #include"RepWireBond.h"
38 #include"RepCylBond.h"
39 #include"RepDot.h"
40 #include"RepMesh.h"
41 #include"RepSphere.h"
42 #include"RepSphereImmediate.h"
43 #include"RepRibbon.h"
44 #include"RepCartoon.h"
45 #include"RepSurface.h"
46 #include"RepLabel.h"
47 #include"RepNonbonded.h"
48 #include"RepNonbondedSphere.h"
49 #include"RepEllipsoid.h"
50 
51 #include"PyMOLGlobals.h"
52 #include"PyMOLObject.h"
53 #include "Executive.h"
54 #include "Lex.h"
55 
56 #ifdef _PYMOL_IP_PROPERTIES
57 #include "Property.h"
58 #endif
59 
60 
61 /*========================================================================*/
62 /*
63  * Get coordinate index for given atom index
64  */
atmToIdx(int atm) const65 int CoordSet::atmToIdx(int atm) const {
66   if (Obj->DiscreteFlag) {
67     if (this == Obj->DiscreteCSet[atm])
68       return Obj->DiscreteAtmToIdx[atm];
69     return -1;
70   }
71   return AtmToIdx[atm];
72 }
73 
74 /*========================================================================*/
75 static char sATOM[] = "ATOM  ";
76 static char sHETATM[] = "HETATM";
77 
CoordSetValidateRefPos(CoordSet * I)78 int CoordSetValidateRefPos(CoordSet * I)
79 {
80   if(I->RefPos) {
81     VLACheck(I->RefPos, RefPosType, I->NIndex);
82     return true;
83   } else {
84     int ok = true && (I->RefPos = pymol::vla<RefPosType>(I->NIndex));
85     if(ok) {
86       int a;
87       for(a = 0; a < I->NIndex; a++) {
88         const float* src = I->coordPtr(a);
89         copy3f(src, I->RefPos[a].coord);
90         I->RefPos[a].specified = true;
91       }
92     }
93     return ok;
94   }
95 }
96 
97 
98 /*========================================================================*/
BondCompare(BondType * a,BondType * b)99 int BondCompare(BondType * a, BondType * b)
100 {
101   int ai0 = a->index[0];
102   int bi0 = b->index[0];
103   if(ai0 == bi0) {
104     int ai1 = a->index[1];
105     int bi1 = b->index[1];
106     if(ai1 == bi1) {
107       return 0;
108     } else if(ai1 > bi1) {
109       return 1;
110     } else {
111       return -1;
112     }
113   } else if(ai0 > bi0) {
114     return 1;
115   } else {
116     return -1;
117   }
118 }
119 
120 
121 /*========================================================================*/
BondInOrder(BondType * a,int b1,int b2)122 int BondInOrder(BondType * a, int b1, int b2)
123 {
124   return (BondCompare(a + b1, a + b2) <= 0);
125 }
126 
CoordSetFromPyList(PyMOLGlobals * G,PyObject * list,CoordSet ** cs)127 int CoordSetFromPyList(PyMOLGlobals * G, PyObject * list, CoordSet ** cs)
128 {
129   CoordSet *I = NULL;
130   int ok = true;
131   int ll = 0;
132 
133   if(*cs) {
134     (*cs)->fFree();
135     *cs = NULL;
136   }
137 
138   if(list == Py_None) {         /* allow None for CSet */
139     *cs = NULL;
140   } else {
141 
142     if(ok)
143       I = CoordSetNew(G);
144     if(ok)
145       ok = (I != NULL);
146     if(ok)
147       ok = (list != NULL);
148     if(ok)
149       ok = PyList_Check(list);
150     if(ok)
151       ll = PyList_Size(list);
152     /* TO SUPPORT BACKWARDS COMPATIBILITY...
153        Always check ll when adding new PyList_GetItem's */
154     if(ok)
155       ok = PConvPyIntToInt(PyList_GetItem(list, 0), &I->NIndex);
156     if(ok)
157       ok = PConvPyIntToInt(PyList_GetItem(list, 1), &I->NAtIndex);
158     if(ok)
159       ok = PConvPyListToFloatVLA(PyList_GetItem(list, 2), &I->Coord);
160     if(ok)
161       ok = PConvPyListToIntVLA(PyList_GetItem(list, 3), &I->IdxToAtm);
162     if(ok && (ll > 5))
163       ok = CPythonVal_PConvPyStrToStr_From_List(G, list, 5, I->Name, sizeof(WordType));
164     if(ok && (ll > 6)){
165       CPythonVal *val = CPythonVal_PyList_GetItem(G, list, 6);
166       ok = ObjectStateFromPyList(G, val, I);
167       CPythonVal_Free(val);
168     }
169     if(ok && (ll > 7)){
170       CPythonVal *val = CPythonVal_PyList_GetItem(G, list, 7);
171       I->Setting = SettingNewFromPyList(G, val);
172       CPythonVal_Free(val);
173     }
174     if(ok && (ll > 8)){
175       CPythonVal *val = CPythonVal_PyList_GetItem(G, list, 8);
176       ok = CPythonVal_PConvPyListToLabPosVLA(G, val, &I->LabPos);
177       CPythonVal_Free(val);
178     }
179 
180 #ifdef _PYMOL_IP_PROPERTIES
181 #endif
182 
183     if(ok && (ll > 10)){
184       CPythonVal *val = CPythonVal_PyList_GetItem(G, list, 10);
185       if (!CPythonVal_IsNone(val))
186 	I->SculptCGO = CGONewFromPyList(G, val, 0, 1);
187       else {
188 	I->SculptShaderCGO = I->SculptCGO = NULL;
189       }
190       CPythonVal_Free(val);
191     }
192     if(ok){
193       if(ll > 11){
194 	// load in atom-state settings
195 	CPythonVal *val = CPythonVal_PyList_GetItem(G, list, 11);
196 	if (!CPythonVal_IsNone(val)){
197 	  int a;
198 	  I->atom_state_setting_id = pymol::vla<int>(I->NIndex);
199 	  I->has_atom_state_settings = pymol::vla<char>(I->NIndex);
200 	  for (a=0; a<I->NIndex; a++){
201 	    CPythonVal *val2 = CPythonVal_PyList_GetItem(G, val, a);
202 	    if (!CPythonVal_IsNone(val2)){
203 	      CPythonVal_PConvPyIntToInt(val2, &I->atom_state_setting_id[a]);
204 	      I->has_atom_state_settings[a] = (I->atom_state_setting_id[a]) ? 1 : 0;
205 	      if (I->atom_state_setting_id[a]) {
206 	        I->atom_state_setting_id[a] = SettingUniqueConvertOldSessionID(G, I->atom_state_setting_id[a]);
207 	      }
208 	    }
209 	    CPythonVal_Free(val2);
210 	  }
211 	} else {
212 	  I->atom_state_setting_id = NULL;
213 	  I->has_atom_state_settings = NULL;
214 	}
215 	CPythonVal_Free(val);
216       } else {
217 	// check I->LabPos.offset to see if its set
218 	if (I->LabPos){
219 	  int a;
220 	  for (a=0; a<I->NIndex; a++){
221 	    if(length3f(I->LabPos[a].offset) > R_SMALL4) {
222 	      SettingSet(cSetting_label_placement_offset, I->LabPos[a].offset, I, a);
223 	    }
224 	  }
225 	}
226       }
227     }
228     if(!ok) {
229       if(I)
230         I->fFree();
231       *cs = NULL;
232     } else {
233       *cs = I;
234     }
235   }
236   return (ok);
237 }
238 
239 /*
240  * Coord set as numpy array
241  */
CoordSetAsNumPyArray(CoordSet * cs,short copy)242 PyObject *CoordSetAsNumPyArray(CoordSet * cs, short copy)
243 {
244 #ifndef _PYMOL_NUMPY
245   PRINTFB(cs->G, FB_CoordSet, FB_Errors)
246     "No numpy support\n" ENDFB(cs->G);
247   return NULL;
248 #else
249 
250   PyObject *result = NULL;
251   const int base_size = sizeof(float);
252   int typenum = -1;
253   npy_intp dims[2] = {0, 3};
254 
255   import_array1(NULL);
256 
257   switch(base_size) {
258     case 4: typenum = NPY_FLOAT32; break;
259     case 8: typenum = NPY_FLOAT64; break;
260   }
261 
262   if(typenum == -1) {
263     printf("error: no typenum for float size %d\n", base_size);
264     return NULL;
265   }
266 
267   dims[0] = cs->NIndex;
268 
269   if(copy) {
270     if((result = PyArray_SimpleNew(2, dims, typenum)))
271       memcpy(PyArray_DATA((PyArrayObject *)result), cs->Coord, cs->NIndex * 3 * base_size);
272   } else {
273     result = PyArray_SimpleNewFromData(2, dims, typenum, cs->Coord.data());
274   }
275 
276   return result;
277 #endif
278 }
279 
CoordSetAsPyList(CoordSet * I)280 PyObject *CoordSetAsPyList(CoordSet * I)
281 {
282   PyObject *result = NULL;
283 
284   if(I) {
285     auto G = I->G;
286     int pse_export_version = SettingGet<float>(G, cSetting_pse_export_version) * 1000;
287     bool dump_binary = SettingGet<bool>(G, cSetting_pse_binary_dump) && (!pse_export_version || pse_export_version >= 1765);
288     result = PyList_New(12);
289     PyList_SetItem(result, 0, PyInt_FromLong(I->NIndex));
290     PyList_SetItem(result, 1, PyInt_FromLong(I->NAtIndex));
291     PyList_SetItem(result, 2, PConvFloatArrayToPyList(I->Coord, I->NIndex * 3, dump_binary));
292     PyList_SetItem(result, 3, PConvIntArrayToPyList(I->IdxToAtm, I->NIndex, dump_binary));
293     if(I->AtmToIdx
294         && pse_export_version < 1770)
295       PyList_SetItem(result, 4, PConvIntArrayToPyList(I->AtmToIdx, I->NAtIndex, dump_binary));
296     else
297       PyList_SetItem(result, 4, PConvAutoNone(NULL));
298     PyList_SetItem(result, 5, PyString_FromString(I->Name));
299     PyList_SetItem(result, 6, ObjectStateAsPyList(I));
300     PyList_SetItem(result, 7, SettingAsPyList(I->Setting));
301     PyList_SetItem(result, 8, PConvLabPosVLAToPyList(I->LabPos, I->NIndex));
302 
303     PyList_SetItem(result, 9,
304 #ifdef _PYMOL_IP_PROPERTIES
305 #endif
306         PConvAutoNone(Py_None));
307 
308     if(I->SculptCGO) {
309       PyList_SetItem(result, 10, CGOAsPyList(I->SculptCGO));
310     } else {
311       PyList_SetItem(result, 10, PConvAutoNone(NULL));
312     }
313     if (I->has_atom_state_settings){
314       int a;
315       PyObject *settings_list = NULL;
316       settings_list = PyList_New(I->NIndex);
317       for (a=0; a<I->NIndex; a++){
318 	if (I->has_atom_state_settings[a]){
319 	  PyList_SetItem(settings_list, a, PyInt_FromLong(I->atom_state_setting_id[a]));
320 	} else {
321 	  PyList_SetItem(settings_list, a, PConvAutoNone(NULL));
322 	}
323       }
324       PyList_SetItem(result, 11, settings_list);
325     } else {
326       PyList_SetItem(result, 11, PConvAutoNone(NULL));
327     }
328     /* TODO symmetry, spheroid, periodic box ... */
329   }
330   return (PConvAutoNone(result));
331 }
332 
CoordSetAdjustAtmIdx(CoordSet * I,int * lookup,int nAtom)333 void CoordSetAdjustAtmIdx(CoordSet * I, int *lookup, int nAtom)
334 
335 /* performs second half of removal */
336 {
337   /* NOTE: only works in a compressive mode, where lookup[a]<=a  */
338   int a, a0, atm;
339   char *new_has_atom_state_settings_by_atom = NULL;
340   int *new_atom_state_setting_id_by_atom = NULL;
341   int nIndex = I->NIndex;
342 
343   if (I->has_atom_state_settings){
344     new_has_atom_state_settings_by_atom = VLACalloc(char, nIndex);
345     new_atom_state_setting_id_by_atom = VLACalloc(int, nIndex);
346   }
347 
348   for(a = 0; a < nIndex; a++) {
349     atm = I->IdxToAtm[a];
350     a0 = lookup[atm];
351     if (a0 < 0){
352       if (I->has_atom_state_settings){
353 	if (I->has_atom_state_settings[a]){
354 	  SettingUniqueDetachChain(I->G, I->atom_state_setting_id[a]);
355 	  I->has_atom_state_settings[a] = 0;
356 	  I->atom_state_setting_id[a] = 0;
357 	}
358       }
359     } else if (new_has_atom_state_settings_by_atom){
360       new_has_atom_state_settings_by_atom[a0] = I->has_atom_state_settings[a];
361       new_atom_state_setting_id_by_atom[a0] = I->atom_state_setting_id[a];
362     }
363   }
364 
365   if (I->AtmToIdx){
366     for(a = 0; a < I->NAtIndex; a++) {
367       a0 = lookup[a];
368       if(a0 >= 0) {
369 	I->AtmToIdx[a0] = I->AtmToIdx[a];
370       }
371     }
372   }
373   I->NAtIndex = nAtom;
374   if (I->AtmToIdx){
375     VLASize(I->AtmToIdx, int, nAtom);
376   }
377   for(a = 0; a < I->NIndex; a++) {
378     atm = I->IdxToAtm[a] = lookup[I->IdxToAtm[a]];
379     if (new_has_atom_state_settings_by_atom){
380       I->has_atom_state_settings[a] = new_has_atom_state_settings_by_atom[atm];
381       I->atom_state_setting_id[a] = new_atom_state_setting_id_by_atom[atm];
382     }
383   }
384   if (new_has_atom_state_settings_by_atom){
385     VLAFreeP(new_has_atom_state_settings_by_atom);
386     VLAFreeP(new_atom_state_setting_id_by_atom);
387   }
388   PRINTFD(I->G, FB_CoordSet)
389     " CoordSetAdjustAtmIdx-Debug: leaving... NAtIndex: %d NIndex %d\n",
390     I->NAtIndex, I->NIndex ENDFD;
391 
392 }
393 
CoordSetCopy(const CoordSet * src)394 CoordSet* CoordSetCopy(const CoordSet* src)
395 {
396   if(!src) {
397     return nullptr;
398   }
399   return new CoordSet(*src);
400 }
401 
402 
403 /*========================================================================*/
CoordSetMerge(ObjectMolecule * OM,CoordSet * I,CoordSet * cs)404 int CoordSetMerge(ObjectMolecule *OM, CoordSet * I, CoordSet * cs)
405 {                               /* must be non-overlapping */
406   int nIndex;
407   int a, i0;
408   int ok = true;
409   /* calculate new size and make room for new data */
410   nIndex = I->NIndex + cs->NIndex;
411   VLASize(I->IdxToAtm, int, nIndex);
412   CHECKOK(ok, I->IdxToAtm);
413   if (ok)
414     VLACheck(I->Coord, float, nIndex * 3);
415   CHECKOK(ok, I->Coord);
416   if (ok){
417     for(a = 0; a < cs->NIndex; a++) {
418       i0 = a + I->NIndex;
419       I->IdxToAtm[i0] = cs->IdxToAtm[a];
420       if (OM->DiscreteFlag){
421 	int idx = cs->IdxToAtm[a];
422 	OM->DiscreteAtmToIdx[idx] = i0;
423 	OM->DiscreteCSet[idx] = I;
424       } else {
425 	I->AtmToIdx[cs->IdxToAtm[a]] = i0;
426       }
427       copy3f(cs->coordPtr(a), I->coordPtr(i0));
428     }
429   }
430   if (ok){
431     if(cs->LabPos) {
432       if(!I->LabPos)
433 	I->LabPos = pymol::vla<LabPosType>(nIndex);
434       else
435 	VLACheck(I->LabPos, LabPosType, nIndex);
436       if(I->LabPos) {
437 	UtilCopyMem(I->LabPos + I->NIndex, cs->LabPos, sizeof(LabPosType) * cs->NIndex);
438       }
439     } else if(I->LabPos) {
440       VLACheck(I->LabPos, LabPosType, nIndex);
441     }
442   }
443   if (ok){
444     if(cs->RefPos) {
445       if(!I->RefPos)
446 	I->RefPos = pymol::vla<RefPosType>(nIndex);
447       else
448 	VLACheck(I->RefPos, RefPosType, nIndex);
449       if(I->RefPos) {
450 	UtilCopyMem(I->RefPos + I->NIndex, cs->RefPos, sizeof(RefPosType) * cs->NIndex);
451       }
452     } else if(I->RefPos) {
453       VLACheck(I->RefPos, RefPosType, nIndex);
454     }
455     I->invalidateRep(cRepAll, cRepInvAll);
456   }
457   I->NIndex = nIndex;
458 
459   return ok;
460 }
461 
462 
463 /*========================================================================*/
CoordSetPurge(CoordSet * I)464 void CoordSetPurge(CoordSet * I)
465 
466 /* performs first half of removal  */
467 {
468   auto G = I->G;
469   int offset = 0;
470   int a, a1, ao;
471   AtomInfoType *ai;
472   ObjectMolecule *obj;
473   float *c0, *c1;
474   LabPosType *l0, *l1;
475   RefPosType *r0, *r1;
476   int *atom_state0, *atom_state1;
477   char *has_atom_state0, *has_atom_state1;
478   obj = I->Obj;
479 
480   PRINTFD(G, FB_CoordSet)
481     " CoordSetPurge-Debug: entering..." ENDFD;
482 
483   c0 = c1 = I->Coord.data();
484   r0 = r1 = I->RefPos.data();
485   l0 = l1 = I->LabPos.data();
486   atom_state0 = atom_state1 = I->atom_state_setting_id.data();
487   has_atom_state0 = has_atom_state1 = I->has_atom_state_settings.data();
488 
489   /* This loop slides down the atoms that are not deleted (deleteFlag)
490      it moves the Coord, RefPos, and LabPos */
491   for(a = 0; a < I->NIndex; a++) {
492     a1 = I->IdxToAtm[a];
493     ai = obj->AtomInfo + a1;
494     if(ai->deleteFlag) {
495       offset--;
496       c0 += 3;
497       if(l0)
498         l0++;
499       if(r0)
500         r0++;
501       if (has_atom_state0){
502 	atom_state0++;
503 	has_atom_state0++;
504       }
505     } else if(offset) {
506       ao = a + offset;
507       *(c1++) = *(c0++);
508       *(c1++) = *(c0++);
509       *(c1++) = *(c0++);
510       if(r1) {
511         *(r1++) = *(r0++);
512       }
513       if(l0) {
514         *(l1++) = *(l0++);
515       }
516       if (has_atom_state0){
517 	*(atom_state1++) = *(atom_state0++);
518 	*(has_atom_state1++) = *(has_atom_state0++);
519       }
520       if (I->AtmToIdx)
521 	I->AtmToIdx[a1] = ao;
522       I->IdxToAtm[ao] = a1;     /* no adjustment of these indexes yet... */
523       if (I->Obj->DiscreteFlag){
524 	I->Obj->DiscreteAtmToIdx[a1] = ao;
525 	I->Obj->DiscreteCSet[a1] = I;
526       }
527     } else {
528       c0 += 3;
529       c1 += 3;
530       if(r1) {
531         r0++;
532         r1++;
533       }
534       if(l0) {
535         l0++;
536         l1++;
537       }
538       if (has_atom_state0){
539 	atom_state0++; atom_state1++;
540 	has_atom_state0++; has_atom_state1++;
541       }
542     }
543   }
544   if(offset) {
545     /* If there were deleted atoms, (offset < 0), then
546        re-adjust the array sizes */
547     I->NIndex += offset;
548     VLASize(I->Coord, float, I->NIndex * 3);
549     if(I->LabPos) {
550       VLASize(I->LabPos, LabPosType, I->NIndex);
551     }
552     if(I->RefPos) {
553       VLASize(I->RefPos, RefPosType, I->NIndex);
554     }
555     if(I->has_atom_state_settings) {
556       VLASize(I->has_atom_state_settings, char, I->NIndex);
557       VLASize(I->atom_state_setting_id, int, I->NIndex);
558     }
559     VLASize(I->IdxToAtm, int, I->NIndex);
560     PRINTFD(G, FB_CoordSet)
561       " CoordSetPurge-Debug: I->IdxToAtm shrunk to %d\n", I->NIndex ENDFD;
562     I->invalidateRep(cRepAll, cRepInvAtoms);      /* this will free Color */
563   }
564   PRINTFD(G, FB_CoordSet)
565     " CoordSetPurge-Debug: leaving NAtIndex %d NIndex %d...\n",
566     I->NAtIndex, I->NIndex ENDFD;
567 
568 #ifdef _PYMOL_IP_PROPERTIES
569 #endif
570 }
571 
572 
573 /*========================================================================*/
CoordSetTransformAtomTTTf(CoordSet * I,int at,const float * TTT)574 int CoordSetTransformAtomTTTf(CoordSet * I, int at, const float *TTT)
575 {
576   int a1 = I->atmToIdx(at);
577   if(a1 < 0)
578     return false;
579 
580   auto* v1 = I->coordPtr(a1);
581   MatrixTransformTTTfN3f(1, v1, TTT, v1);
582   return true;
583 }
584 
585 
586 /*========================================================================*/
CoordSetTransformAtomR44f(CoordSet * I,int at,const float * matrix)587 int CoordSetTransformAtomR44f(CoordSet * I, int at, const float *matrix)
588 {
589   int a1 = I->atmToIdx(at);
590   if(a1 < 0)
591     return false;
592 
593   auto* v1 = I->coordPtr(a1);
594   MatrixTransformR44fN3f(1, v1, matrix, v1);
595   return true;
596 }
597 
598 
599 /*========================================================================*/
CoordSetRecordTxfApplied(CoordSet * I,const float * matrix,int homogenous)600 void CoordSetRecordTxfApplied(CoordSet * I, const float *matrix, int homogenous)
601 {
602   double temp[16];
603 
604   if(!homogenous) {
605     convertTTTfR44d(matrix, temp);
606   } else {
607     convert44f44d(matrix, temp);
608   }
609 
610   ObjectStateLeftCombineMatrixR44d(I, temp);
611 }
612 
613 
614 /*========================================================================*/
CoordSetMoveAtom(CoordSet * I,int at,const float * v,int mode)615 int CoordSetMoveAtom(CoordSet * I, int at, const float *v, int mode)
616 {
617   int a1 = I->atmToIdx(at);
618   if(a1 < 0)
619     return false;
620 
621   auto* v1 = I->coordPtr(a1);
622   if(mode) {
623     add3f(v, v1, v1);
624   } else {
625     copy3f(v, v1);
626   }
627   return true;
628 }
629 
630 
631 /*========================================================================*/
CoordSetMoveAtomLabel(CoordSet * I,int at,const float * v,const float * diff)632 int CoordSetMoveAtomLabel(CoordSet * I, int at, const float *v, const float *diff)
633 {
634   auto G = I->G;
635   ObjectMolecule *obj = I->Obj;
636   int a1 = I->atmToIdx(at);
637   int result = 0;
638 
639   /* if label is valid, get the label offset
640    * and set the new position relative to that */
641   if(a1 >= 0) {
642     float at_offset[3];
643     const float * at_offset_ptr;
644     int at_label_relative_mode = 0;
645     AtomInfoType *ai = obj->AtomInfo + at;
646 
647     AtomStateGetSetting_i(G, obj, I, a1, ai, cSetting_label_relative_mode, &at_label_relative_mode);
648     switch (at_label_relative_mode){
649     case 0:
650       AtomStateGetSetting(G, obj, I, a1, ai, cSetting_label_placement_offset, &at_offset_ptr);
651       add3f(v, at_offset_ptr, at_offset);
652       SettingSet(cSetting_label_placement_offset, at_offset, I, a1);
653       break;
654     case 1: // screen relative
655     case 2: // screen pixel space
656       {
657 	float voff[3];
658 	int width, height;
659 	SceneGetWidthHeight(G, &width, &height);
660 	if (at_label_relative_mode==1){
661 	  voff[0] = 2.f * diff[0] / width;
662 	  voff[1] = 2.f * diff[1] / height;
663 	} else {
664 	  voff[0] = diff[0];
665 	  voff[1] = diff[1];
666 	}
667 	voff[2] = 0.f;
668 	AtomStateGetSetting(G, obj, I, a1, ai, cSetting_label_screen_point, &at_offset_ptr);
669 	add3f(voff, at_offset_ptr, at_offset);
670 	SettingSet(cSetting_label_screen_point, at_offset, I, a1);
671       }
672       break;
673     }
674   }
675 
676   return (result);
677 }
678 
679 
680 /*========================================================================*/
CoordSetGetAtomVertex(const CoordSet * I,int at,float * v)681 int CoordSetGetAtomVertex(const CoordSet * I, int at, float *v)
682 {
683   int a1 = I->atmToIdx(at);
684 
685   if(a1 < 0)
686     return false;
687 
688   copy3f(I->coordPtr(a1), v);
689   return true;
690 }
691 
692 
693 /*========================================================================*/
CoordSetGetAtomTxfVertex(const CoordSet * I,int at,float * v)694 int CoordSetGetAtomTxfVertex(const CoordSet * I, int at, float *v)
695 {
696   ObjectMolecule *obj = I->Obj;
697   int a1 = I->atmToIdx(at);
698 
699   if(a1 < 0)
700     return false;
701 
702   copy3f(I->coordPtr(a1), v);
703 
704   /* apply state transformation */
705   if (!I->Matrix.empty() && SettingGet<int>(I->G, obj->Setting, I->Setting,
706                                 cSetting_matrix_mode) > 0) {
707     transform44d3f(I->Matrix.data(), v, v);
708   }
709 
710   /* object transformation */
711   if(obj->TTTFlag) {
712     transformTTT44f3f(obj->TTT, v, v);
713   }
714 
715   return true;
716 }
717 
718 
719 /*========================================================================*/
CoordSetSetAtomVertex(CoordSet * I,int at,const float * v)720 int CoordSetSetAtomVertex(CoordSet * I, int at, const float *v)
721 {
722   int a1 = I->atmToIdx(at);
723 
724   if(a1 < 0)
725    return false;
726 
727   copy3f(v, I->coordPtr(a1));
728   return true;
729 }
730 
731 
732 /*========================================================================*/
CoordSetRealToFrac(CoordSet * I,const CCrystal * cryst)733 void CoordSetRealToFrac(CoordSet * I, const CCrystal * cryst)
734 {
735   int a;
736   float* v = I->Coord.data();
737   for(a = 0; a < I->NIndex; a++) {
738     transform33f3f(cryst->RealToFrac, v, v);
739     v += 3;
740   }
741 }
742 
743 
744 /*========================================================================*/
CoordSetTransform44f(CoordSet * I,const float * mat)745 void CoordSetTransform44f(CoordSet * I, const float *mat)
746 {
747   int a;
748   float* v = I->Coord.data();
749   for(a = 0; a < I->NIndex; a++) {
750     transform44f3f(mat, v, v);
751     v += 3;
752   }
753 }
754 
755 
756 /*========================================================================*/
757 
CoordSetTransform33f(CoordSet * I,const float * mat)758 void CoordSetTransform33f(CoordSet * I, const float *mat)
759 {
760   int a;
761   float* v = I->Coord.data();
762   for(a = 0; a < I->NIndex; a++) {
763     transform33f3f(mat, v, v);
764     v += 3;
765   }
766 }
767 
768 
769 /*========================================================================*/
CoordSetGetAverage(const CoordSet * I,float * v0)770 void CoordSetGetAverage(const CoordSet * I, float *v0)
771 {
772   int a;
773   double accum[3];
774   if(I->NIndex) {
775     const float* v = I->Coord.data();
776     accum[0] = *(v++);
777     accum[1] = *(v++);
778     accum[2] = *(v++);
779     for(a = 1; a < I->NIndex; a++) {
780       accum[0] += *(v++);
781       accum[1] += *(v++);
782       accum[2] += *(v++);
783     }
784     v0[0] = (float) (accum[0] / I->NIndex);
785     v0[1] = (float) (accum[1] / I->NIndex);
786     v0[2] = (float) (accum[2] / I->NIndex);
787   }
788 }
789 
790 
791 /*========================================================================*/
CoordSetFracToReal(CoordSet * I,const CCrystal * cryst)792 void CoordSetFracToReal(CoordSet * I, const CCrystal * cryst)
793 {
794   int a;
795   float* v = I->Coord.data();
796   for(a = 0; a < I->NIndex; a++) {
797     transform33f3f(cryst->FracToReal, v, v);
798     v += 3;
799   }
800 }
801 
802 /*
803  * Apply the `sca` (SCALEn) transformation to transform to fractional space,
804  * and then the crystals `FracToReal` transformation to transform back to
805  * cartesian space.
806  *
807  * Don't do anything if pdb_insure_orthogonal=off.
808  *
809  * Don't do anything if SCALEn or CRYST1 look bogus. There is a number of
810  * structures in the PDB which have meaningless values for those.
811  *
812  * Without this, creating symmetry mates might produce wrong results.
813  */
CoordSetInsureOrthogonal(PyMOLGlobals * G,CoordSet * cset,const float * sca,const CCrystal * cryst,bool quiet)814 bool CoordSetInsureOrthogonal(PyMOLGlobals * G,
815     CoordSet * cset,            // coord set to modify
816     const float * sca,          // 4x4 SCALE
817     const CCrystal *cryst,
818     bool quiet)
819 {
820   if (!SettingGetGlobal_b(G, cSetting_pdb_insure_orthogonal))
821     return false;
822 
823   if (!cryst)
824     cryst = &cset->Symmetry->Crystal;
825 
826   const float * r2f = cryst->RealToFrac;
827 
828   // are the matrices sufficiently close to be the same?
829   if (!sca[3] && !sca[7] && !sca[11] &&
830       is_allclosef(3, r2f, 3, sca, 4, R_SMALL4)) {
831     return false;
832   }
833 
834   // is the cell a orthogonal 1x1x1? If so, then it should probably be ignored...
835   // is SCALEn the identity matrix?  If so, then it should probably be ignored...
836   if (is_identityf(3, r2f, R_SMALL4) ||
837       is_identityf(4, sca, R_SMALL4)) {
838     PRINTFB(G, FB_ObjectMolecule, FB_Blather)
839       " ObjectMolReadPDBStr: ignoring SCALEn (identity matrix).\n" ENDFB(G);
840     return false;
841   }
842 
843   // is SCALEn invalid?  If so, then it should definitely be ignored...
844   if (determinant33f(sca, 4) < R_SMALL8 ||
845       determinant33f(r2f, 3) < R_SMALL8) {
846     PRINTFB(G, FB_ObjectMolecule, FB_Blather)
847       " ObjectMolReadPDBStr: ignoring SCALEn (invalid matrix).\n" ENDFB(G);
848     return false;
849   }
850 
851   PRINTFB(G, FB_ObjectMolecule, quiet ? FB_Blather : FB_Actions)
852     " ObjectMolecule: using SCALEn to compute orthogonal coordinates.\n"
853     ENDFB(G);
854 
855   CoordSetTransform44f(cset, sca);
856   CoordSetFracToReal(cset, cryst);
857 
858   return true;
859 }
860 
861 /*========================================================================*/
862 /* Rotates the ANISOU vector
863  *
864  * matrix: flat 4x4, but only rotation (upper left 3x3) is considered
865  * anisou: has 6 elements (of symmetric 3x3) and will be rotated in-place
866  */
RotateU(const double * matrix,float * anisou)867 bool RotateU(const double *matrix, float *anisou)
868 {
869   int i, j, k;
870   float Re[3][3];
871   double e_val[3], e_vec[3][3];
872   double U[3][3] = {
873     { anisou[0], anisou[3], anisou[4] },
874     { anisou[3], anisou[1], anisou[5] },
875     { anisou[4], anisou[5], anisou[2] },
876   };
877 
878   // e_val, e_vec = linalg.eigh(U)
879   if(!xx_matrix_jacobi_solve(*e_vec, e_val, &i, *U, 3))
880     return false;
881 
882   // Re = dot(matrix[:3,:3], e_vec)
883   for (i = 0; i < 3; i++)
884     for (j = 0; j < 3; j++) {
885       Re[i][j] = 0.0;
886       for (k = 0; k < 3; k++)
887         Re[i][j] += matrix[i * 4 + k] * e_vec[k][j];
888     }
889 
890   // U = dot(Re * e_val, Re.T)
891   for (i = 0; i < 3; i++)
892     for (j = 0; j <= i; j++) {
893       U[i][j] = 0.0;
894       for (k = 0; k < 3; k++)
895         U[i][j] += Re[i][k] * e_val[k] * Re[j][k];
896     }
897 
898   anisou[0] = U[0][0];
899   anisou[1] = U[1][1];
900   anisou[2] = U[2][2];
901   anisou[3] = U[1][0];
902   anisou[4] = U[2][0];
903   anisou[5] = U[2][1];
904 
905   return true;
906 }
907 
908 /*========================================================================*/
CoordSetAtomToPDBStrVLA(PyMOLGlobals * G,char ** charVLA,int * c,const AtomInfoType * ai,const float * v,int cnt,const PDBInfoRec * pdb_info,const double * matrix)909 void CoordSetAtomToPDBStrVLA(PyMOLGlobals * G, char **charVLA, int *c,
910                              const AtomInfoType * ai,
911                              const float *v, int cnt,
912                              const PDBInfoRec * pdb_info,
913                              const double *matrix)
914 /*
915  * v: 3x1 vertex in final output space
916  * matrix: 4x4 homogenous transformation matrix from model space to output
917  *         space (view matrix * state matrix). Used for ANISOU.
918  */
919 {
920   char *aType;
921   AtomName name;
922   ResName resn;
923   lexidx_t chain;
924 
925   char formalCharge[4];
926   int ignore_pdb_segi = SettingGetGlobal_b(G, cSetting_ignore_pdb_segi);
927   WordType x, y, z;
928 
929   AtomInfoGetAlignedPDBResidueName(G, ai, resn);
930   AtomInfoGetAlignedPDBAtomName(G, ai, resn, name);
931 
932   formalCharge[0] = 0;
933   if(SettingGetGlobal_b(G, cSetting_pdb_formal_charges)) {
934     if((ai->formalCharge > 0) && (ai->formalCharge < 10)) {
935       sprintf(formalCharge, "%d+", ai->formalCharge);
936     } else if((ai->formalCharge < 0) && (ai->formalCharge > -10)) {
937       sprintf(formalCharge, "%d-", -ai->formalCharge);
938     }
939   }
940 
941   if(ai->hetatm)
942     aType = sHETATM;
943   else
944     aType = sATOM;
945 
946   char inscode = ai->getInscode(true);
947 
948   VLACheck(*charVLA, char, (*c) + 1000);
949 
950   if(SettingGetGlobal_b(G, cSetting_pdb_retain_ids)) {
951     cnt = ai->id - 1;
952   }
953   if(cnt > 99998)
954     cnt = 99998;
955 
956   if((!pdb_info) || (!pdb_info->is_pqr_file())) { /* relying upon short-circuit */
957     short linelen;
958     sprintf(x, "%8.3f", v[0]);
959     x[8] = 0;
960     sprintf(y, "%8.3f", v[1]);
961     y[8] = 0;
962     sprintf(z, "%8.3f", v[2]);
963     z[8] = 0;
964     linelen =
965       sprintf((*charVLA) + (*c),
966               "%6s%5i %-4s%1s%-4s%1.1s%4i%c   %s%s%s%6.2f%6.2f      %-4.4s%2s%2s\n", aType,
967               cnt + 1, name, ai->alt, resn, LexStr(G, ai->chain), ai->resv % 10000, inscode, x, y, z, ai->q, ai->b,
968               ignore_pdb_segi ? "" :
969               LexStr(G, ai->segi), ai->elem, formalCharge);
970     if(ai->anisou) {
971       // Columns 7 - 27 and 73 - 80 are identical to the corresponding ATOM/HETATM record.
972       char *atomline = (*charVLA) + (*c);
973       char *anisoline = atomline + linelen;
974       float anisou[6];
975       std::copy_n(ai->anisou, 6, anisou);
976 
977       if(matrix && !RotateU(matrix, anisou)) {
978         PRINTFB(G, FB_CoordSet, FB_Errors) "RotateU failed\n" ENDFB(G);
979         return;
980       }
981       strncpy(anisoline + 6, atomline + 6, 22);
982       sprintf(anisoline + 28,
983           "%7.0f%7.0f%7.0f%7.0f%7.0f%7.0f",
984           anisou[0] * 1e4, anisou[1] * 1e4, anisou[2] * 1e4,
985           anisou[3] * 1e4, anisou[4] * 1e4, anisou[5] * 1e4);
986       strcpy(anisoline + 70, atomline + 70);
987       strncpy(anisoline, "ANISOU", 6);
988       (*c) += linelen;
989     }
990     (*c) += linelen;
991   } else {
992     Chain alt;
993     if(pdb_info->is_pqr_file() && pdb_info->pqr_workarounds) {
994       inscode = ' ';
995       chain = 0;                /* no chain IDs */
996       alt[0] = 0;               /* not alt conf identifiers */
997     } else {
998       alt[0] = ai->alt[0];
999       alt[1] = 0;
1000       chain = ai->chain;
1001     }
1002     sprintf(x, "%8.3f", v[0]);
1003     if(x[0] != 32)
1004       sprintf(x, " %7.2f", v[0]);
1005     x[8] = 0;
1006     sprintf(y, "%8.3f", v[1]);
1007     y[8] = 0;
1008     if(y[0] != 32)
1009       sprintf(y, " %7.2f", v[1]);
1010     y[8] = 0;
1011     sprintf(z, "%8.3f", v[2]);
1012     if(z[0] != 32)
1013       sprintf(z, " %7.2f", v[2]);
1014     z[8] = 0;
1015 
1016     (*c) += sprintf((*charVLA) + (*c), "%6s%5i %-4s%1s%-4s%1.1s%4i%c   %s%s%s %11.8f %7.3f\n",
1017                     aType, cnt + 1, name, alt, resn,
1018                     LexStr(G, chain), ai->resv, inscode, x, y, z, ai->partialCharge, ai->elec_radius);
1019   }
1020 
1021 }
1022 
1023 
1024 /*========================================================================*/
CoordSetAtomToChemPyAtom(PyMOLGlobals * G,AtomInfoType * ai,const float * v,const float * ref,int index,const double * matrix)1025 PyObject *CoordSetAtomToChemPyAtom(PyMOLGlobals * G, AtomInfoType * ai, const float *v,
1026                                    const float *ref, int index, const double *matrix)
1027 {
1028 #ifdef _PYMOL_NOPY
1029   return NULL;
1030 #else
1031   PyObject *atom = PYOBJECT_CALLMETHOD(P_chempy, "Atom", "");
1032   if(!atom)
1033     ErrMessage(G, "CoordSetAtomToChemPyAtom", "can't create atom");
1034   else {
1035     float tmp_array[6] = { 0.f, 0.f, 0.f, 0.f, 0.f, 0.f };
1036 
1037     if (ai->anisou) {
1038       std::copy_n(ai->anisou, 6, tmp_array);
1039       if (matrix)
1040       RotateU(matrix, tmp_array);
1041     }
1042 
1043     PConvFloat3ToPyObjAttr(atom, "coord", v);
1044     if(ref)
1045       PConvFloat3ToPyObjAttr(atom, "ref_coord", ref);
1046     if (ai->name)
1047     PConvStringToPyObjAttr(atom, "name", LexStr(G, ai->name));
1048     PConvStringToPyObjAttr(atom, "symbol", ai->elem);
1049 // TODO defaults to UNK    if (ai->resn)
1050     PConvStringToPyObjAttr(atom, "resn", LexStr(G, ai->resn));
1051     if (ai->inscode) {
1052       char ins_code[2] = {ai->inscode, 0};
1053       PConvStringToPyObjAttr(atom, "ins_code", ins_code);
1054     }
1055     if (ai->ssType[0])
1056     PConvStringToPyObjAttr(atom, "ss", ai->ssType);
1057 // TODO defaults to 1    if (ai->resv)
1058     PConvIntToPyObjAttr(atom, "resi_number", ai->resv);
1059     if (ai->stereo)
1060     PConvIntToPyObjAttr(atom, "stereo", ai->stereo);
1061     if (ai->chain)
1062     PConvStringToPyObjAttr(atom, "chain", LexStr(G, ai->chain));
1063     if(ai->alt[0])
1064       PConvStringToPyObjAttr(atom, "alt", ai->alt);
1065     if (ai->segi)
1066     PConvStringToPyObjAttr(atom, "segi", LexStr(G, ai->segi));
1067     if (ai->q != 1.)
1068     PConvFloatToPyObjAttr(atom, "q", ai->q);
1069     if (ai->b)
1070     PConvFloatToPyObjAttr(atom, "b", ai->b);
1071     if (ai->anisou)
1072     {
1073       {
1074         PyObject *tmp_obj = PConvFloatArrayToPyList(tmp_array, 6);
1075         if(tmp_obj) {
1076           PyObject_SetAttrString(atom, "u_aniso", tmp_obj);
1077           Py_XDECREF(tmp_obj);
1078         }
1079       }
1080     }
1081     PConvFloatToPyObjAttr(atom, "vdw", ai->vdw);
1082     if (ai->elec_radius)
1083     PConvFloatToPyObjAttr(atom, "elec_radius", ai->elec_radius);
1084     if (ai->partialCharge)
1085     PConvFloatToPyObjAttr(atom, "partial_charge", ai->partialCharge);
1086     if (ai->formalCharge)
1087     PConvIntToPyObjAttr(atom, "formal_charge", ai->formalCharge);
1088 // TODO customType=0 from most files
1089     if(ai->customType != -9999)
1090       PConvIntToPyObjAttr(atom, "numeric_type", ai->customType);
1091     if (ai->textType)
1092     PConvStringToPyObjAttr(atom, "text_type", LexStr(G, ai->textType));
1093     if (ai->custom)
1094     PConvStringToPyObjAttr(atom, "custom", LexStr(G, ai->custom));
1095 
1096     PConvIntToPyObjAttr(atom, "hetatm", ai->hetatm);
1097     PConvIntToPyObjAttr(atom, "flags", ai->flags);
1098     PConvIntToPyObjAttr(atom, "id", ai->id);    /* not necc. unique */
1099     PConvIntToPyObjAttr(atom, "index", index + 1);      /* fragile */
1100 
1101 #ifdef _PYMOL_IP_PROPERTIES
1102 #endif
1103 
1104   }
1105   if(PyErr_Occurred())
1106     PyErr_Print();
1107   return (atom);
1108 #endif
1109 }
1110 
1111 
1112 /*========================================================================*/
invalidateRep(int type,int level)1113 void CoordSet::invalidateRep(int type, int level)
1114 {
1115   CoordSet * I = this;
1116   int a;
1117   if(level >= cRepInvVisib) {
1118     if (I->Obj)
1119       I->Obj->RepVisCacheValid = false;
1120   }
1121   /* graphical representations need redrawing */
1122   if(level == cRepInvVisib) {
1123     /* cartoon_side_chain_helper */
1124     if(SettingGet<bool>(G, I->Setting, I->Obj->Setting,
1125                     cSetting_cartoon_side_chain_helper)) {
1126       if((type == cRepCyl) || (type == cRepLine) || (type == cRepSphere))
1127         invalidateRep(cRepCartoon, cRepInvVisib2);
1128       else if(type == cRepCartoon) {
1129         invalidateRep(cRepLine, cRepInvVisib2);
1130         invalidateRep(cRepCyl, cRepInvVisib2);
1131         invalidateRep(cRepSphere, cRepInvVisib2);
1132       }
1133     }
1134     /* ribbon_side_chain_helper */
1135     if(SettingGet<bool>(G, I->Setting, I->Obj->Setting,
1136                     cSetting_ribbon_side_chain_helper)) {
1137       if((type == cRepCyl) || (type == cRepLine) || (type == cRepSphere))
1138         invalidateRep(cRepRibbon, cRepInvVisib2);
1139       else if(type == cRepRibbon) {
1140         invalidateRep(cRepLine, cRepInvVisib2);
1141         invalidateRep(cRepCyl, cRepInvVisib2);
1142         invalidateRep(cRepSphere, cRepInvVisib2);
1143       }
1144     }
1145     /* line_stick helper  */
1146     if(SettingGet<bool>(G, I->Setting, I->Obj->Setting,
1147                     cSetting_line_stick_helper)) {
1148       if(type == cRepCyl)
1149         invalidateRep(cRepLine, cRepInvVisib2);
1150       else if(type == cRepLine) {
1151         invalidateRep(cRepCyl, cRepInvVisib2);
1152       }
1153     }
1154   }
1155 
1156   if(!I->Spheroid.empty())
1157     if (I->Spheroid.size() !=
1158         I->NAtIndex * GetSpheroidSphereRec(G)->nDot) {
1159       I->Spheroid.clear();
1160       I->SpheroidNormal.clear();
1161     }
1162 
1163   /* invalidate basd on one representation, 'type' */
1164   for (RepIterator iter(G, type); iter.next(); ){
1165     int eff_level = level;
1166     a = iter.rep;
1167     if(level == cRepInvPick) {
1168       switch (a) {
1169       case cRepSurface:
1170       case cRepMesh:
1171       case cRepDot:
1172         /* skip the expensive to recompute, non-pickable
1173            representations */
1174         break;
1175       default:               /* default behavior is to blow away the representation */
1176         eff_level = cRepInvRep;
1177         break;
1178       }
1179     }
1180     if(eff_level >= cRepInvVisib)     /* make active if visibility has changed */
1181       I->Active[a] = true;
1182     if(I->Rep[a]) {
1183       if(I->Rep[a]->fInvalidate && (eff_level < cRepInvPurge))
1184         I->Rep[a]->fInvalidate(I->Rep[a], I, eff_level);
1185       else if(eff_level >= cRepInvExtColor) {
1186         I->Rep[a]->fFree(I->Rep[a]);
1187         I->Rep[a] = NULL;
1188       }
1189     }
1190   }
1191 
1192   if(level >= cRepInvCoord) {   /* if coordinates change, then this map becomes invalid */
1193     MapFree(I->Coord2Idx);
1194     I->Coord2Idx = NULL;
1195     ExecutiveInvalidateSelectionIndicatorsCGO(G);
1196     SceneInvalidatePicking(G);
1197     /* invalidate distances */
1198   }
1199 
1200 #ifndef NO_MMLIBS
1201   if (level >= cRepInvProp) {
1202     if (validMMStereo == MMPYMOLX_PROP_STATE_AUTO)
1203       validMMStereo = MMPYMOLX_PROP_STATE_NULL;
1204     if (validTextType == MMPYMOLX_PROP_STATE_AUTO)
1205       validTextType = MMPYMOLX_PROP_STATE_NULL;
1206   }
1207 #endif
1208 
1209   SceneChanged(G);
1210 }
1211 
1212 
1213 /*========================================================================*/
1214 
1215 #define RepUpdateMacro(I,rep,new_fn,state) {\
1216   if(I->Active[rep]&&(!G->Interrupt)) {\
1217     if(!I->Rep[rep]) {\
1218       I->Rep[rep]=new_fn(I,state);\
1219       if(I->Rep[rep]){ \
1220          I->Rep[rep]->fNew=(struct Rep *(*)(struct CoordSet *,int state))new_fn;\
1221          SceneInvalidatePicking(G);\
1222       } else {  \
1223 	I->Active[rep] = false;			\
1224       }         \
1225     } else {\
1226       if(I->Rep[rep]->fUpdate)\
1227          I->Rep[rep] = I->Rep[rep]->fUpdate(I->Rep[rep],I,state,rep);\
1228     }\
1229   }\
1230 OrthoBusyFast(G,rep,cRepCnt);\
1231 }
1232 
1233 /*========================================================================*/
update(int state)1234 void CoordSet::update(int state)
1235 {
1236   CoordSet * I = this;
1237   int a;
1238   assert(G == I->Obj->G);
1239 
1240   PRINTFB(G, FB_CoordSet, FB_Blather) " CoordSetUpdate-Entered: object %s state %d cset %p\n",
1241     I->Obj->Name, state, (void *) I
1242     ENDFB(G);
1243 
1244   OrthoBusyFast(G, 0, cRepCnt);
1245   RepUpdateMacro(I, cRepLine, RepWireBondNew, state);
1246   RepUpdateMacro(I, cRepCyl, RepCylBondNew, state);
1247   RepUpdateMacro(I, cRepDot, RepDotNew, state);
1248   RepUpdateMacro(I, cRepMesh, RepMeshNew, state);
1249   RepUpdateMacro(I, cRepSphere, RepSphereNew, state);
1250   RepUpdateMacro(I, cRepRibbon, RepRibbonNew, state);
1251   RepUpdateMacro(I, cRepCartoon, RepCartoonNew, state);
1252   RepUpdateMacro(I, cRepSurface, RepSurfaceNew, state);
1253   RepUpdateMacro(I, cRepLabel, RepLabelNew, state);
1254   RepUpdateMacro(I, cRepNonbonded, RepNonbondedNew, state);
1255   RepUpdateMacro(I, cRepNonbondedSphere, RepNonbondedSphereNew, state);
1256   RepUpdateMacro(I, cRepEllipsoid, RepEllipsoidNew, state);
1257 
1258   for(a = 0; a < cRepCnt; a++)
1259     if(!I->Rep[a])
1260       I->Active[a] = false;
1261 
1262   SceneInvalidate(G);
1263   OrthoBusyFast(G, 1, 1);
1264   if(Feedback(G, FB_CoordSet, FB_Blather)) {
1265     printf(" CoordSetUpdate-Leaving: object %s state %d cset %p\n",
1266            I->Obj->Name, state, (void *) I);
1267   }
1268 }
1269 
1270 
1271 /*========================================================================*/
CoordSetUpdateCoord2IdxMap(CoordSet * I,float cutoff)1272 void CoordSetUpdateCoord2IdxMap(CoordSet * I, float cutoff)
1273 {
1274   if(cutoff < R_SMALL4)
1275     cutoff = R_SMALL4;
1276   if(I->NIndex > 10) {
1277     if(I->Coord2Idx) {
1278       if((I->Coord2IdxDiv < cutoff) ||
1279          (((cutoff - I->Coord2IdxReq) / I->Coord2IdxReq) < -0.5F)) {
1280         MapFree(I->Coord2Idx);
1281         I->Coord2Idx = NULL;
1282       }
1283     }
1284     if(I->NIndex && (!I->Coord2Idx)) {  /* NOTE: map based on stored coords */
1285       I->Coord2IdxReq = cutoff;
1286       I->Coord2IdxDiv = cutoff * 1.25F;
1287       I->Coord2Idx = MapNew(I->G, I->Coord2IdxDiv, I->Coord, I->NIndex, NULL);
1288       if(I->Coord2IdxDiv < I->Coord2Idx->Div)
1289         I->Coord2IdxDiv = I->Coord2Idx->Div;
1290     }
1291   }
1292 }
1293 
1294 /*
1295  * RepToTransparencySetting: maps rep to transparency setting
1296  */
RepToTransparencySetting(const int rep_id)1297 static int RepToTransparencySetting(const int rep_id){
1298   switch(rep_id){
1299     // based on the rep,
1300     // transparency is set with different settings
1301   case cRepCyl:
1302     return cSetting_stick_transparency;
1303   case cRepSurface:
1304     return cSetting_transparency;
1305   case cRepSphere:
1306     return cSetting_sphere_transparency;
1307   case cRepEllipsoid:
1308     return cSetting_ellipsoid_transparency;
1309   case cRepCartoon:
1310     return cSetting_cartoon_transparency;
1311   case cRepRibbon:
1312     return cSetting_ribbon_transparency;
1313   }
1314   return 0;
1315 }
1316 
1317 /*========================================================================*/
render(RenderInfo * info)1318 void CoordSet::render(RenderInfo * info)
1319 {
1320   CoordSet * I = this;
1321   PRINTFD(G, FB_CoordSet)
1322     " CoordSetRender: entered (%p).\n", (void *) I ENDFD;
1323 
1324   if(!(info->ray || info->pick) &&
1325      (SettingGet_i(G, I->Setting, I->Obj->Setting,
1326                    cSetting_defer_builds_mode) == 5)) {
1327     if(!info->pass) {
1328       ObjectUseColor((CObject *) I->Obj);
1329       if(I->Active[cRepLine])
1330         RepWireBondRenderImmediate(I, info);
1331       if(I->Active[cRepNonbonded])
1332         RepNonbondedRenderImmediate(I, info);
1333       if(I->Active[cRepSphere])
1334         RepSphereRenderImmediate(I, info);
1335       if(I->Active[cRepCyl])
1336         RepCylBondRenderImmediate(I, info);
1337       if(I->Active[cRepRibbon])
1338         RepRibbonRenderImmediate(I, info);
1339     }
1340   } else {
1341     int pass = info->pass;
1342     CRay *ray = info->ray;
1343     auto pick = info->pick;
1344     int a, aa, abit, aastart = 0, aaend = cRepCnt;
1345     ::Rep *r;
1346     int sculpt_vdw_vis_mode = SettingGet_i(G, I->Setting,
1347 					   I->Obj->Setting,
1348 					   cSetting_sculpt_vdw_vis_mode);
1349     if((!pass) && sculpt_vdw_vis_mode &&
1350        I->SculptCGO && (I->Obj->visRep & cRepCGOBit)) {
1351       if(ray) {
1352         int ok = CGORenderRay(I->SculptCGO, ray, info,
1353 			      ColorGet(G, I->Obj->Color), NULL, I->Setting, I->Obj->Setting);
1354 	if (!ok){
1355 	  CGOFree(I->SculptCGO);
1356 	  CGOFree(I->SculptShaderCGO);
1357 	}
1358       } else if(G->HaveGUI && G->ValidContext) {
1359         if(!pick) {
1360 	  int use_shader = SettingGetGlobal_b(G, cSetting_use_shaders);
1361 	  if (use_shader){
1362 	    if (!I->SculptShaderCGO){
1363 	      CGO *convertcgo = NULL;
1364 	      convertcgo = CGOCombineBeginEnd(I->SculptCGO, 0);
1365 	      if (convertcgo){
1366 		I->SculptShaderCGO = CGOOptimizeToVBONotIndexed(convertcgo, 0);
1367 		I->SculptShaderCGO->use_shader = true;
1368 		CGOFree(convertcgo);
1369 	      }
1370 	    }
1371 	  } else {
1372 	    CGOFree(I->SculptShaderCGO);
1373 	  }
1374 	  if (I->SculptShaderCGO){
1375 	    CGORenderGL(I->SculptShaderCGO, NULL,
1376 			I->Setting, I->Obj->Setting, info, NULL);
1377 	  } else {
1378 	    CGORenderGL(I->SculptCGO, NULL,
1379 			I->Setting, I->Obj->Setting, info, NULL);
1380 	  }
1381         }
1382       }
1383     }
1384 
1385     if (pick){
1386       int pick_labels = SettingGet_i(G, I->Setting,
1387 				     I->Obj->Setting,
1388 				     cSetting_pick_labels);
1389       if (pick_labels == 2){ // only pick labels
1390 	aastart = cRepLabel;
1391 	aaend = aastart + 1;
1392       }
1393     }
1394     for(aa = aastart; aa < aaend; aa++) {
1395       if(aa == cRepSurface) {   /* reorder */
1396         a = cRepCell;
1397       } else if(aa == cRepCell) {
1398         a = cRepSurface;
1399       } else {
1400         a = aa;
1401       }
1402 
1403       abit = (1 << a);
1404 
1405       if(I->Active[a] && I->Rep[a]) {
1406         r = I->Rep[a];
1407         if(!ray) {
1408           ObjectUseColor((CObject *) I->Obj);
1409         } else {
1410           if(I->Obj)
1411             ray->wobble(
1412                          SettingGet_i(G, I->Setting,
1413                                       I->Obj->Setting,
1414                                       cSetting_ray_texture),
1415                          SettingGet_3fv(G, I->Setting,
1416                                         I->Obj->Setting,
1417                                         cSetting_ray_texture_settings));
1418           else
1419             ray->wobble(
1420                          SettingGet_i(G, I->Setting,
1421                                       NULL, cSetting_ray_texture),
1422                          SettingGet_3fv(G, I->Setting, NULL,
1423                                         cSetting_ray_texture_settings));
1424           ray->color3fv(ColorGet(G, I->Obj->Color));
1425         }
1426 
1427         if(r->fRender) {        /* do OpenGL rendering in three passes */
1428           if(ray || pick) {
1429 
1430             /* here we need to iterate through and apply coordinate set matrices */
1431 
1432             r->fRender(r, info);
1433           } else {
1434             bool t_mode_3 = SettingGetGlobal_i(G, cSetting_transparency_mode) == 3;
1435             bool render_both = t_mode_3;  // render both opaque (1) and transparent (-1) pass if t_mode_3
1436             /* here we need to iterate through and apply coordinate set matrices */
1437 
1438             switch (a) {
1439             case cRepVolume:
1440               {
1441                 int t_mode = SettingGetGlobal_i(G, cSetting_transparency_mode);
1442                 if (t_mode == 3 && pass == -1){
1443                   r->fRender(r, info);
1444                 }
1445               }
1446               break;
1447             case cRepLabel:
1448               {
1449                 int t_mode_3 = SettingGetGlobal_i(G, cSetting_transparency_mode) == 3;
1450                 if (pass == -1 || (t_mode_3 && pass == 1 /* TODO_OPENVR 0 */))
1451                   r->fRender(r, info);
1452               }
1453               break;
1454             case cRepDot:
1455             case cRepCGO:
1456             case cRepCallback:
1457               if(pass == 1)
1458                 r->fRender(r, info);
1459               break;
1460             case cRepLine:
1461             case cRepMesh:
1462             case cRepDash:
1463             case cRepCell:
1464             case cRepExtent:
1465               if(!pass)
1466                 r->fRender(r, info);
1467               break;
1468             case cRepNonbonded:
1469             case cRepRibbon:
1470               render_both = false; // cartoon and nonbonded do not have atom-level transparency
1471             case cRepNonbondedSphere:
1472             case cRepSurface:
1473             case cRepEllipsoid:
1474             case cRepCyl:
1475             case cRepCartoon:
1476             case cRepSphere:
1477               {
1478                 if (render_both){
1479                   // for transparency_mode 3, render both opaque and transparent pass
1480                   if (pass != 0){
1481                     r->fRender(r, info);
1482                   }
1483                 } else {
1484                   bool checkAlphaCGO = abit & (cRepSurfaceBit);
1485                   int check_setting = RepToTransparencySetting(a);
1486                   bool cont = true;
1487                   if (checkAlphaCGO){
1488                     if(info->alpha_cgo) {
1489                       if(pass == 1){
1490                         r->fRender(r, info);
1491                       }
1492                       cont = false;
1493                     }
1494                   }
1495                   if (cont){
1496                     if(check_setting && SettingGet_f(G, r->cs->Setting,
1497                                                      r->obj->Setting, check_setting) > 0.0001) {
1498                       /* if object has transparency, only render in -1 pass */
1499                       if(pass == -1)
1500                         r->fRender(r, info);
1501                     } else if(pass == 1){
1502                       r->fRender(r, info);
1503                     }
1504                   }
1505                 }
1506               }
1507               break;
1508             }
1509           }
1510         }
1511         /*          if(ray)
1512            ray->fWobble(ray,0,NULL); */
1513       }
1514     }
1515   }
1516   PRINTFD(G, FB_CoordSet)
1517     " CoordSetRender: leaving...\n" ENDFD;
1518 }
1519 
1520 
1521 /*========================================================================*/
CoordSet(PyMOLGlobals * G)1522 CoordSet::CoordSet(PyMOLGlobals* G)
1523     : CObjectState(G)
1524 {
1525 }
1526 
1527 /*========================================================================*/
CoordSet(const CoordSet & cs)1528 CoordSet::CoordSet(const CoordSet& cs)
1529     : CObjectState(cs)
1530 {
1531   this->Obj = cs.Obj;
1532   this->Coord = cs.Coord;
1533   this->IdxToAtm = cs.IdxToAtm;
1534   this->NIndex = cs.NIndex;
1535   this->NAtIndex = cs.NAtIndex;
1536   this->prevNIndex = cs.prevNIndex;
1537   this->prevNAtIndex = cs.prevNAtIndex;
1538   std::copy(std::begin(cs.Rep), std::end(cs.Rep), std::begin(this->Rep));
1539   std::copy(std::begin(cs.Active), std::end(cs.Active), std::begin(this->Active));
1540   this->NTmpBond = cs.NTmpBond;
1541   this->NTmpLinkBond = cs.NTmpLinkBond;
1542   /* deep copy & return ptr to new symmetry */
1543   if(cs.Symmetry) {
1544     this->Symmetry = pymol::make_unique<CSymmetry>(*cs.Symmetry);
1545   }
1546   if(cs.PeriodicBox) {
1547     this->PeriodicBox = pymol::make_unique<CCrystal>(*cs.PeriodicBox);
1548   }
1549   std::copy(std::begin(cs.Name), std::end(cs.Name), std::begin(this->Name));
1550   this->PeriodicBoxType = cs.PeriodicBoxType;
1551   this->tmp_index = cs.tmp_index;
1552   this->Coord2IdxReq = cs.Coord2IdxReq;
1553   this->Coord2IdxDiv = cs.Coord2IdxDiv;
1554   this->objMolOpInvalidated = cs.objMolOpInvalidated;
1555 
1556   // copy VLAs
1557   this->LabPos     = cs.LabPos;
1558   this->RefPos     = cs.RefPos;
1559   this->AtmToIdx   = cs.AtmToIdx;
1560 
1561   UtilZeroMem(this->Rep, sizeof(::Rep *) * cRepCnt);
1562 
1563 #ifdef _PYMOL_IP_PROPERTIES
1564 #endif
1565 
1566 }
1567 
1568 
1569 /*========================================================================*/
extendIndices(int nAtom)1570 int CoordSet::extendIndices(int nAtom)
1571 {
1572   CoordSet * I = this;
1573   int a, b;
1574   ObjectMolecule *obj = I->Obj;
1575   int ok = true;
1576   if(obj->DiscreteFlag) {
1577     ok = obj->setNDiscrete(nAtom);
1578 
1579     if(I->AtmToIdx) {           /* convert to discrete if necessary */
1580       I->AtmToIdx.freeP();
1581       if (ok){
1582 	for(a = 0; a < I->NIndex; a++) {
1583 	  b = I->IdxToAtm[a];
1584 	  obj->DiscreteAtmToIdx[b] = a;
1585 	  obj->DiscreteCSet[b] = I;
1586 	}
1587       }
1588     }
1589   }
1590   if(ok && I->NAtIndex < nAtom) {
1591     if(I->AtmToIdx) {
1592       VLASize(I->AtmToIdx, int, nAtom);
1593       CHECKOK(ok, I->AtmToIdx);
1594       if(ok && nAtom) {
1595         for(a = I->NAtIndex; a < nAtom; a++)
1596           I->AtmToIdx[a] = -1;
1597       }
1598       I->NAtIndex = nAtom;
1599     } else if(!obj->DiscreteFlag) {
1600       I->AtmToIdx = pymol::vla<int>(nAtom);
1601       CHECKOK(ok, I->AtmToIdx);
1602       if (ok){
1603 	for(a = 0; a < nAtom; a++)
1604 	  I->AtmToIdx[a] = -1;
1605 	I->NAtIndex = nAtom;
1606       }
1607     }
1608   }
1609   return ok;
1610 }
1611 
1612 
1613 /*========================================================================*/
appendIndices(int offset)1614 void CoordSet::appendIndices(int offset)
1615 {
1616   CoordSet * I = this;
1617   int a, b;
1618   ObjectMolecule *obj = I->Obj;
1619 
1620   I->IdxToAtm = pymol::vla<int>(I->NIndex);
1621   if(I->NIndex) {
1622     ErrChkPtr(G, I->IdxToAtm);
1623     for(a = 0; a < I->NIndex; a++)
1624       I->IdxToAtm[a] = a + offset;
1625   }
1626   if(obj->DiscreteFlag) {
1627     VLACheck(obj->DiscreteAtmToIdx, int, I->NIndex + offset);
1628     VLACheck(obj->DiscreteCSet, CoordSet *, I->NIndex + offset);
1629     for(a = 0; a < I->NIndex; a++) {
1630       b = a + offset;
1631       obj->DiscreteAtmToIdx[b] = a;
1632       obj->DiscreteCSet[b] = I;
1633     }
1634   } else {
1635     I->AtmToIdx = pymol::vla<int>(I->NIndex + offset);
1636     if(I->NIndex + offset) {
1637       ErrChkPtr(G, I->AtmToIdx);
1638       for(a = 0; a < offset; a++)
1639         I->AtmToIdx[a] = -1;
1640       for(a = 0; a < I->NIndex; a++)
1641         I->AtmToIdx[a + offset] = a;
1642     }
1643   }
1644   I->NAtIndex = I->NIndex + offset;
1645 }
1646 
1647 
1648 /*========================================================================*/
enumIndices()1649 void CoordSet::enumIndices()
1650 {
1651   CoordSet * I = this;
1652   /* set up for simple case where 1 = 1, etc. */
1653   int a;
1654   I->AtmToIdx = pymol::vla<int>(I->NIndex);
1655   I->IdxToAtm = pymol::vla<int>(I->NIndex);
1656   if(I->NIndex) {
1657     ErrChkPtr(G, I->AtmToIdx);
1658     ErrChkPtr(G, I->IdxToAtm);
1659     for(a = 0; a < I->NIndex; a++) {
1660       I->AtmToIdx[a] = a;
1661       I->IdxToAtm[a] = a;
1662     }
1663   }
1664   I->NAtIndex = I->NIndex;
1665 }
1666 
1667 
1668 /*========================================================================*/
~CoordSet()1669 CoordSet::~CoordSet()
1670 {
1671   CoordSet * I = this;
1672   int a;
1673   ObjectMolecule *obj;
1674   if(I) {
1675 #ifdef _PYMOL_IP_PROPERTIES
1676 #endif
1677 
1678     if (I->has_atom_state_settings){
1679       for(a = 0; a < I->NIndex; a++){
1680         if (I->has_atom_state_settings[a]){
1681           SettingUniqueDetachChain(G, I->atom_state_setting_id[a]);
1682         }
1683       }
1684     }
1685     for(a = 0; a < cRepCnt; a++)
1686       if(I->Rep[a])
1687         I->Rep[a]->fFree(I->Rep[a]);
1688     obj = I->Obj;
1689     if(obj)
1690       if(obj->DiscreteFlag)     /* remove references to the atoms in discrete objects */
1691         for(a = 0; a < I->NIndex; a++) {
1692           obj->DiscreteAtmToIdx[I->IdxToAtm[a]] = -1;
1693           obj->DiscreteCSet[I->IdxToAtm[a]] = NULL;
1694         }
1695     MapFree(I->Coord2Idx);
1696     SettingFreeP(I->Setting);
1697     CGOFree(I->SculptCGO);
1698   }
1699 }
fFree()1700 void CoordSet::fFree()
1701 {
1702   delete this;
1703 }
LabPosTypeCopy(const LabPosType * src,LabPosType * dst)1704 void LabPosTypeCopy(const LabPosType * src, LabPosType * dst){
1705   dst->mode = src->mode;
1706   copy3f(src->pos, dst->pos);
1707   copy3f(src->offset, dst->offset);
1708 }
RefPosTypeCopy(const RefPosType * src,RefPosType * dst)1709 void RefPosTypeCopy(const RefPosType * src, RefPosType * dst){
1710   copy3f(src->coord, dst->coord);
1711   dst->specified = src->specified;
1712 }
1713 
1714 #ifndef _PYMOL_NOPY
CoordSetSetSettingFromPyObject(PyMOLGlobals * G,CoordSet * cs,int at,int setting_id,PyObject * val)1715 int CoordSetSetSettingFromPyObject(PyMOLGlobals * G, CoordSet *cs, int at, int setting_id, PyObject *val){
1716   if (val == Py_None)
1717     val = NULL;
1718 
1719   if (!val) {
1720     if (!cs->has_atom_state_settings ||
1721         !cs->has_atom_state_settings[at])
1722       return true;
1723   }
1724 
1725   CoordSetCheckUniqueID(G, cs, at);
1726   cs->has_atom_state_settings[at] = true;
1727 
1728   return SettingUniqueSetPyObject(G, cs->atom_state_setting_id[at], setting_id, val);
1729 }
1730 #endif
1731 
CoordSetCheckSetting(PyMOLGlobals * G,CoordSet * cs,int at,int setting_id)1732 int CoordSetCheckSetting(PyMOLGlobals * G, CoordSet *cs, int at, int setting_id){
1733   if(!cs->has_atom_state_settings || !cs->has_atom_state_settings[at]) {
1734     return 0;
1735   } else {
1736     if(!SettingUniqueCheck(G, cs->atom_state_setting_id[at], setting_id)) {
1737       return 0;
1738     } else {
1739       return 1;
1740     }
1741   }
1742 }
1743 
SettingGetIfDefinedPyObject(PyMOLGlobals * G,CoordSet * cs,int at,int setting_id)1744 PyObject *SettingGetIfDefinedPyObject(PyMOLGlobals * G, CoordSet *cs, int at, int setting_id){
1745   if(cs->has_atom_state_settings && cs->has_atom_state_settings[at]){
1746     return SettingUniqueGetPyObject(G, cs->atom_state_setting_id[at], setting_id);
1747   }
1748   return NULL;
1749 }
1750 
CoordSetCheckUniqueID(PyMOLGlobals * G,CoordSet * I,int at)1751 int CoordSetCheckUniqueID(PyMOLGlobals * G, CoordSet *I, int at){
1752   if (!I->atom_state_setting_id){
1753     I->atom_state_setting_id = pymol::vla<int>(I->NIndex);
1754   }
1755   if (!I->has_atom_state_settings){
1756     I->has_atom_state_settings = pymol::vla<char>(I->NIndex);
1757   }
1758   if (!I->atom_state_setting_id[at]){
1759     I->atom_state_setting_id[at] = AtomInfoGetNewUniqueID(G);
1760   }
1761   return I->atom_state_setting_id[at];
1762 }
1763 
1764 template <typename V>
AtomStateGetSetting(ATOMSTATEGETSETTINGARGS,V * out)1765 void AtomStateGetSetting(ATOMSTATEGETSETTINGARGS, V * out) {
1766   if (cs->has_atom_state_settings &&
1767       cs->has_atom_state_settings[idx] &&
1768       SettingUniqueGetIfDefined(G, cs->atom_state_setting_id[idx], setting_id, out))
1769     return;
1770 
1771   if (AtomSettingGetIfDefined(G, ai, setting_id, out))
1772     return;
1773 
1774   *out = SettingGet<V>(G, cs->Setting, obj->Setting, setting_id);
1775 }
1776 
1777 template void AtomStateGetSetting(ATOMSTATEGETSETTINGARGS, int * out);
1778 template void AtomStateGetSetting(ATOMSTATEGETSETTINGARGS, float * out);
1779 template void AtomStateGetSetting(ATOMSTATEGETSETTINGARGS, const float ** out);
1780