00001 /*************************************************************************** 00002 *cr 00003 *cr (C) Copyright 1995-2019 The Board of Trustees of the 00004 *cr University of Illinois 00005 *cr All Rights Reserved 00006 *cr 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * RCS INFORMATION: 00011 * 00012 * $RCSfile: GeometryMol.C,v $ 00013 * $Author: johns $ $Locker: $ $State: Exp $ 00014 * $Revision: 1.59 $ $Date: 2020年07月08日 04:40:23 $ 00015 * 00016 *************************************************************************** 00017 * DESCRIPTION: 00018 * 00019 * A base class for all Geometry objects which measure information about 00020 * atoms in a molecule. A molecule Geometry monitor is assumed to operate 00021 * on N atoms, and be able to calculate a single floating-point value for 00022 * those atoms. (i.e. the angle formed by three atoms in space) 00023 * 00024 ***************************************************************************/ 00025 00026 #include "GeometryMol.h" 00027 #include "MoleculeList.h" 00028 #include "Molecule.h" 00029 #include "Displayable.h" 00030 #include "DispCmds.h" 00031 #include "utilities.h" 00032 #include "VMDApp.h" 00033 #include "TextEvent.h" 00034 #include "CommandQueue.h" 00035 #include "PeriodicTable.h" 00036 00038 class GeometryMonitor : public DrawMoleculeMonitor { 00039 private: 00040 GeometryMol *mygeo; 00041 00042 public: 00043 GeometryMonitor(GeometryMol *g, int id) : mygeo(g), molid(id) { } 00044 void notify(int id) { 00045 mygeo->calculate(); 00046 mygeo->update(); 00047 } 00048 const int molid; 00049 }; 00050 00052 GeometryMol::GeometryMol(int n, int *mols, int *atms, const int *cells, 00053 MoleculeList *mlist, CommandQueue *cq, Displayable *d) 00054 : Displayable(d) { 00055 00056 objIndex = new int[n]; 00057 comIndex = new int[n]; 00058 cellIndex = new int[3L*n]; 00059 if (cells) { 00060 memcpy(cellIndex, cells, 3L*n*sizeof(int)); 00061 } else { 00062 memset(cellIndex, 0, 3L*n*sizeof(int)); 00063 } 00064 geomValue = 0.0; 00065 numItems = n; 00066 hasValue = TRUE; 00067 00068 my_color = 8; 00069 my_text_size = 1.0f; 00070 my_text_thickness = 1.0f; 00071 my_text_offset[0] = my_text_offset[1] = 0; 00072 my_text_format = "%R%d:%a"; 00073 00074 molList = mlist; 00075 cmdqueue = cq; 00076 gmName = NULL; 00077 uniquegmName = NULL; 00078 00079 for(int i=0; i < numItems; i++) { 00080 objIndex[i] = mols[i]; 00081 comIndex[i] = atms[i]; 00082 00083 // make sure an atom is not repeated in this list 00084 if(i > 0 && objIndex[i-1]==objIndex[i] && comIndex[i-1]==comIndex[i] && 00085 !memcmp(cellIndex+3L*i, cellIndex+3L*(i-1), 3L*sizeof(int))) { 00086 // set a bogus value for the first atom index, to make the 00087 // check_mol routine fail 00088 comIndex[0] = (-1); 00089 } 00090 00091 Molecule *m = molList->mol_from_id(mols[i]); 00092 if (m) { 00093 GeometryMonitor *mon = new GeometryMonitor(this, objIndex[i]); 00094 m->register_monitor(mon); 00095 monitors.append(mon); 00096 } 00097 } 00098 00099 // sort the items properly 00100 sort_items(); 00101 00102 // create the name for this object 00103 geom_set_name(); 00104 } 00105 00106 00108 GeometryMol::~GeometryMol(void) { 00109 // delete name if necessary 00110 if(gmName) 00111 delete [] gmName; 00112 if(uniquegmName) 00113 delete [] uniquegmName; 00114 delete [] objIndex; 00115 delete [] comIndex; 00116 delete [] cellIndex; 00117 for (int i=0; i<monitors.num(); i++) { 00118 GeometryMonitor *mon = monitors[i]; 00119 Molecule *m = molList->mol_from_id(mon->molid); 00120 if (m) m->unregister_monitor(mon); 00121 delete mon; 00122 } 00123 } 00124 00125 00127 00128 void GeometryMol::atom_full_name(char *buf, Molecule *mol, int ind) { 00129 sprintf(buf, "%-d/%-d", mol->id(), ind); 00130 } 00131 00132 void GeometryMol::atom_short_name(char *buf, Molecule *mol, int ind) { 00133 MolAtom *nameAtom = mol->atom(ind); 00134 sprintf(buf, "%s%d:%s", 00135 mol->resNames.name(nameAtom->resnameindex), 00136 nameAtom->resid, 00137 mol->atomNames.name(nameAtom->nameindex)); 00138 } 00139 00140 void GeometryMol::atom_formatted_name(JString &str, Molecule *mol, int ind) { 00141 MolAtom *atom = mol->atom(ind); 00142 str = my_text_format; 00143 char buf[1024]; 00144 JString resname = mol->resNames.name(atom->resnameindex); 00145 Timestep *ts = mol->current(); 00146 float totalforce = 0.0f; 00147 double physical_time = 0.0f; 00148 if (ts != NULL) { 00149 if (ts->force != NULL) { 00150 float ufx = ts->force[ind*3 ]; 00151 float ufy = ts->force[ind*3 + 1]; 00152 float ufz = ts->force[ind*3 + 2]; 00153 totalforce = sqrtf(ufx*ufx + ufy*ufy + ufz*ufz); 00154 } 00155 physical_time = ts->physical_time; 00156 } 00157 00158 // '%a' atom name 00159 str.gsub("%a", mol->atomNames.name(atom->nameindex)); 00160 00161 // '%d' resid 00162 sprintf(buf, "%d", atom->resid); 00163 str.gsub("%d", buf); 00164 00165 // '%i' atom index (0-based) 00166 sprintf(buf, "%d", ind); 00167 str.gsub("%i", buf); 00168 00169 // '%1i' atom index (1-based) 00170 sprintf(buf, "%d", ind+1); 00171 str.gsub("%1i", buf); 00172 00173 // '%e' atomic element 00174 str.gsub("%e", get_pte_label(atom->atomicnumber)); 00175 00176 // '%b' beta 00177 sprintf(buf, "%4.3f", mol->beta()[ind]); 00178 str.gsub("%b", buf); 00179 00180 // '%c' chain 00181 str.gsub("%c", mol->chainNames.name(atom->chainindex)); 00182 00183 // '%C' conformation 00184 str.gsub("%C", mol->altlocNames.name(atom->altlocindex)); 00185 00186 // '%f' user-applied force 00187 sprintf(buf, "%g", totalforce); 00188 str.gsub("%f", buf); 00189 00190 // '%F' current trajectory frame 00191 sprintf(buf, "%d", mol->frame()); 00192 str.gsub("%F", buf); 00193 00194 // '%m' mass 00195 sprintf(buf, "%4.3f", mol->mass()[ind]); 00196 str.gsub("%m", buf); 00197 00198 // '%n' molecule index 00199 sprintf(buf, "%d", mol->id()); 00200 str.gsub("%n", buf); 00201 00202 // '%N' molecule name 00203 str.gsub("%N", mol->molname()); 00204 00205 // '%o' occupancy 00206 sprintf(buf, "%4.3f", mol->occupancy()[ind]); 00207 str.gsub("%o", buf); 00208 00209 // '%p' atom periodic element number 00210 sprintf(buf, "%d", atom->atomicnumber); 00211 str.gsub("%p", buf); 00212 00213 // '%q' atom charge 00214 sprintf(buf, "%4.3f", mol->charge()[ind]); 00215 str.gsub("%q", buf); 00216 00217 // '%R' resname in upper-case 00218 str.gsub("%R", (const char *)resname); 00219 00220 // '%1R' 1-char resname in upper-case 00221 sprintf(buf, "%c", resname[0]); 00222 str.gsub("%1R", buf); 00223 00224 // '%r' resname in camel case 00225 resname.to_camel(); 00226 str.gsub("%r", (const char *)resname); 00227 00228 // '%s' segname 00229 str.gsub("%s", mol->segNames.name(atom->segnameindex)); 00230 00231 // '%t' atom type 00232 str.gsub("%t", mol->atomTypes.name(atom->typeindex)); 00233 00234 // '%T' physical time 00235 sprintf(buf, "%g", physical_time); 00236 str.gsub("%T", buf); 00237 00238 // '%x', '%y', '%z' 00239 if ((ts != NULL) && (ts->pos != NULL)) { 00240 sprintf(buf, "%g", ts->pos[ind*3 ]); 00241 str.gsub("%x", buf); 00242 sprintf(buf, "%g", ts->pos[ind*3 + 1]); 00243 str.gsub("%y", buf); 00244 sprintf(buf, "%g", ts->pos[ind*3 + 2]); 00245 str.gsub("%z", buf); 00246 } 00247 } 00248 00249 // set the name of this item 00250 void GeometryMol::geom_set_name(void) { 00251 char namebuf[256] = { 0 }; 00252 char namebuf2[256] = { 0 }; 00253 char cellbuf[256] = { 0 }; 00254 Molecule *mol; 00255 int i; 00256 00257 if(items() < 1) 00258 return; 00259 00260 // put first name in string 00261 if((mol = check_mol(objIndex[0], comIndex[0]))) { 00262 atom_full_name(namebuf, mol, comIndex[0]); 00263 atom_short_name(namebuf2, mol, comIndex[0]); 00264 sprintf(cellbuf, "-%d-%d-%d", cellIndex[0], cellIndex[1], cellIndex[2]); 00265 strcat(namebuf, cellbuf); 00266 } else 00267 return; 00268 00269 // put rest of names in string in format: n1/n2/n3/.../nN 00270 for(i = 1; i < items(); i++) { 00271 if((mol = check_mol(objIndex[i], comIndex[i]))) { 00272 strcat(namebuf, "/"); 00273 atom_full_name(namebuf+strlen(namebuf), mol, comIndex[i]); 00274 strcat(namebuf2, "/"); 00275 atom_short_name(namebuf2+strlen(namebuf2), mol, comIndex[i]); 00276 sprintf(cellbuf, "-%d-%d-%d", cellIndex[3L*i+0], cellIndex[3L*i+1], 00277 cellIndex[3L*i+2]); 00278 strcat(namebuf, cellbuf); 00279 } else { 00280 return; 00281 } 00282 } 00283 00284 // now make a copy of this name 00285 if(gmName) 00286 delete [] gmName; 00287 if(uniquegmName) 00288 delete [] uniquegmName; 00289 00290 uniquegmName = stringdup(namebuf); 00291 gmName = stringdup(namebuf2); 00292 } 00293 00294 00295 // sort the elements in the list, so that the lowest atom index is first 00296 // (but preserve the relative order, i.e. a-b-c or c-b-a) 00297 void GeometryMol::sort_items(void) { 00298 int i,j; 00299 00300 // swap order if first component index > last component index 00301 if( (comIndex[0] > comIndex[items()- 1]) || 00302 (comIndex[0] == comIndex[items()-1] && 00303 objIndex[0] > objIndex[items()-1]) ) { 00304 for(i=0, j=(items() - 1); i < j; i++, j--) { 00305 int tmpindex = comIndex[i]; 00306 comIndex[i] = comIndex[j]; 00307 comIndex[j] = tmpindex; 00308 tmpindex = objIndex[i]; 00309 objIndex[i] = objIndex[j]; 00310 objIndex[j] = tmpindex; 00311 int celltmp[3]; 00312 memcpy(celltmp, cellIndex+3L*i, 3L*sizeof(int)); 00313 memcpy(cellIndex+3L*i, cellIndex+3L*j, 3L*sizeof(int)); 00314 memcpy(cellIndex+3L*j, celltmp, 3L*sizeof(int)); 00315 } 00316 } 00317 } 00318 00319 00320 // check whether the given molecule & atom index is OK 00321 // if OK, return Molecule pointer; otherwise, return NULL 00322 Molecule *GeometryMol::check_mol(int m, int a) { 00323 00324 Molecule *mol = molList->mol_from_id(m); 00325 00326 if(!mol || a < 0 || a >= mol->nAtoms) 00327 mol = NULL; 00328 00329 return mol; 00330 } 00331 00332 00333 // for the given Molecule, find the TRANSFORMED coords for the given atom 00334 // return Molecule pointer if successful, NULL otherwise. 00335 Molecule *GeometryMol::transformed_atom_coord(int ind, float *pos) { 00336 Molecule *mol; 00337 00338 // only return value if molecule is legal and atom is displayed 00339 int a=comIndex[ind]; 00340 if((mol = normal_atom_coord(ind, pos)) && mol->atom_displayed(a)) { 00341 00342 // now multiply it by the molecule's tranformation matrix 00343 (mol->tm).multpoint3d(pos, pos); 00344 00345 // calculation was successful; return the molecule pointer 00346 return mol; 00347 } 00348 00349 // if here, error (i.e. atom not displayed, or not proper mol id) 00350 return NULL; 00351 } 00352 00353 00354 // for the given Molecule, find the UNTRANSFORMED coords for the given atom 00355 // return Molecule pointer if successful, NULL otherwise. 00356 Molecule *GeometryMol::normal_atom_coord(int ind, float *pos) { 00357 Timestep *now; 00358 Molecule *mol; 00359 00360 int m=objIndex[ind]; 00361 int a=comIndex[ind]; 00362 const int *cell = cellIndex+3L*ind; 00363 00364 // get the molecule pointer, and get the coords for the current timestep 00365 if((mol = check_mol(m, a)) && (now = mol->current())) { 00366 memcpy((void *)pos, (void *)(now->pos + 3L*a), 3L*sizeof(float)); 00367 00368 // Apply periodic image transformation before returning 00369 Matrix4 mat; 00370 now->get_transform_from_cell(cell, mat); 00371 mat.multpoint3d(pos, pos); 00372 00373 return mol; 00374 } 00375 00376 // if here, error (i.e. atom not displayed, or not proper mol id) 00377 return NULL; 00378 } 00379 00380 00381 // draws a line between the two given points 00382 void GeometryMol::display_line(float *p1, float *p2, VMDDisplayList *d) { 00383 DispCmdLine cmdLineGeom; 00384 cmdLineGeom.putdata(p1, p2, d); 00385 } 00386 00387 00388 // print given text at current valuePos position 00389 void GeometryMol::display_string(const char *str, VMDDisplayList *d) { 00390 DispCmdText cmdTextGeom; 00391 cmdTextGeom.putdata(valuePos, str, my_text_thickness, my_text_size, 00392 my_text_offset[0], my_text_offset[1], d); 00393 } 00394 00396 00397 // return the name of this geometry marker; by default, just blank 00398 const char *GeometryMol::name(void) { 00399 return (gmName ? gmName : ""); 00400 } 00401 00402 00403 // return 'unique' name of the marker, which should be different than 00404 // other names for different markers of this same type 00405 const char *GeometryMol::unique_name(void) { 00406 return (uniquegmName ? uniquegmName : name()); 00407 } 00408 00409 00410 // check whether the geometry value can still be calculated 00411 int GeometryMol::ok(void) { 00412 int i; 00413 00414 for(i=0; i < numItems; i++) 00415 if(!check_mol(objIndex[i], comIndex[i])) 00416 return FALSE; 00417 00418 return TRUE; 00419 } 00420 00421 00422 // calculate a whole list of items, if this object can do so. Return success. 00423 // We must loop over the maximum number of frames for the set of molecules 00424 // involved in the label in order for multi-molecule labels 00425 // to be correcty updated. 00426 int GeometryMol::calculate_all(ResizeArray<float> &valArray) { 00427 int i, j, n=items(); 00428 00429 if (!has_value()) 00430 return FALSE; 00431 00432 // find the highest number of frames in any molecule in the label, and 00433 // cache the current frame so we can restore it later 00434 int maxframes=0; 00435 ResizeArray<int> curframe(n); 00436 for (i=0; i<n; i++) { 00437 Molecule * mol = molList->mol_from_id(objIndex[i]); 00438 if (!mol) 00439 return FALSE; // WHAT?? should never happen 00440 curframe[i]=mol->frame(); 00441 if (curframe[i]<0) { 00442 // then this molecule has no frames, so we can't calculate anything 00443 return FALSE; 00444 } 00445 if (mol->numframes() > maxframes) 00446 maxframes=mol->numframes(); 00447 } 00448 00449 // go through all the frames, calculating values 00450 for (j=0; j<maxframes; j++) { 00451 /* override frame */ 00452 for (i=0; i<n; i++) { 00453 Molecule * mol = molList->mol_from_id(objIndex[i]); 00454 if (mol) 00455 mol->override_current_frame(j); // paranoia 00456 } 00457 00458 /* calculate value */ 00459 valArray.append(calculate()); 00460 00461 /* restore frame */ 00462 for (i=0; i<n; i++) { 00463 Molecule * mol = molList->mol_from_id(objIndex[i]); 00464 if (mol) 00465 mol->override_current_frame(curframe[i]); // paranoia 00466 } 00467 } 00468 00469 return TRUE; 00470 } 00471 00472 00473 // save the text of a selection to the interpreter variable "vmd_pick_selection" 00474 // The format is of the form "index %d %d %d .... %d" for each atom picked. 00475 // (If nothing is picked, the text is "none") 00476 void GeometryMol::set_pick_selection(int , int num, int *atoms) { 00477 cmdqueue->runcommand(new PickSelectionEvent(num, atoms)); 00478 } 00479 00480 void GeometryMol::set_pick_selection() { 00481 cmdqueue->runcommand(new PickSelectionEvent(0, 0)); 00482 } 00483 00484 void GeometryMol::set_pick_value(double newval) { 00485 cmdqueue->runcommand(new PickValueEvent(newval)); 00486 } 00487