Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

DrawMolItem.C

Go to the documentation of this file.
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: DrawMolItem.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.375 $      $Date: 2024/03/01 01:33:49 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *
00019  * Child Displayable component of a molecule; this is responsible for doing
00020  * the actual drawing of a molecule.  It contains an atom color, atom
00021  * selection, and atom representation object to specify how this component
00022  * should look.
00023  *
00024  ***************************************************************************/
00025 
00026 #include <stdio.h>
00027 #include <stdlib.h>
00028 #include <sys/types.h>
00029 #include <math.h>
00030 #include <ctype.h>              // for isdigit()
00031 
00032 #include "DrawMolItem.h"
00033 #include "Molecule.h"  // for file_in_progress()
00034 #include "DispCmds.h"
00035 #include "Inform.h"
00036 #include "Scene.h"
00037 #include "TextEvent.h"
00038 #include "BondSearch.h"
00039 #include "DisplayDevice.h"
00040 #ifdef VMDMSMS
00041 #include "MSMSInterface.h" // Interface to MSMS surface generation program
00042 #endif
00043 #ifdef VMDNANOSHPER
00044 #include "NanoShaperInterface.h" // Interface to NanoShaper surface program
00045 #endif
00046 #ifdef VMDSURF
00047 #include "Surf.h"          // this is an interface to the SURF program
00048 #endif
00049 #include "VMDApp.h"        // for vmd_alloc/vmd_dealloc
00050 #include "VolumetricData.h"
00051 
00053 DrawMolItem::DrawMolItem(const char *nm, DrawMolecule *dm, AtomColor *ac, 
00054         AtomRep *ar, AtomSel *as) 
00055         : Displayable(dm) {
00056 
00057   // save data and pointers to drawing method objects
00058   mol = dm;
00059   avg = new float[3L*mol->nAtoms];
00060   avgsize = 0;
00061   atomColor = ac;
00062   atomRep = ar;
00063   atomSel = as;
00064   structwarningcount = 0;
00065 
00066   name = stringdup(nm);
00067   framesel = stringdup("now");
00068   tubearray = NULL;
00069 
00070   colorlookups = new ColorLookup[MAXCOLORS]; // for color sorted line drawing
00071 
00072   // use non-standard Catmull-Rom spline, with slope of 1.25 instead of 2.0 
00073   create_modified_CR_spline_basis(spline_basis, 1.25f);
00074 
00075   // set orbital and volumetric data objects sentinel values or to NULL
00076   waveftype = -1;
00077   wavefspin = -1;
00078   wavefexcitation = -1;
00079   gridorbid = -1;
00080   orbgridspacing = -1.0f; 
00081   orbgridisdensity = -1;
00082   orbvol = NULL;
00083 
00084   // initialize volume texture data to invalid values
00085   voltexVolid = -1;
00086   voltexColorMethod = -1;
00087   voltexDataMin = 0;
00088   voltexDataMax = 0;
00089 
00090   // signal that we need a complete refresh next prepare
00091   needRegenerate = MOL_REGEN | COL_REGEN | SEL_REGEN;
00092   update_pbc = 0;
00093   update_ss = 0;
00094   update_ts = 0;
00095   update_traj = 0;
00096   update_instances = 0;
00097   isOn = TRUE;     // newly created reps default to on
00098 }
00099 
00101 DrawMolItem::~DrawMolItem(void) {
00102   if (tubearray) {
00103     for (int i=0; i<tubearray->num(); i++) 
00104       delete (*tubearray)[i];
00105     delete tubearray;
00106     tubearray = NULL;
00107   }
00108 
00109   delete atomColor;
00110   delete atomRep;
00111   delete atomSel;
00112   delete [] colorlookups; // color/atomid lookup for color sorted line drawing
00113   delete [] name;
00114   delete [] avg;
00115   delete [] framesel;
00116   delete orbvol;
00117 }
00118 
00119 int DrawMolItem::emitstructwarning(void) {
00120   if (structwarningcount < 30) {
00121     structwarningcount++;
00122     return 1;
00123   } 
00124 
00125   if (structwarningcount == 30) {
00126     msgErr << "Maximum structure display warnings reached, further warnings will be suppressed" << sendmsg; 
00127     structwarningcount++;
00128     return 0;
00129   }
00130 
00131   structwarningcount++;
00132   return 0;
00133 }
00134 
00135 void DrawMolItem::update_lookups(AtomColor *ac, AtomSel *sel, 
00136                                     ColorLookup *lookups) {
00137   int i;
00138   for (i=0; i<MAXCOLORS; i++) {
00139     lookups[i].num = 0;
00140   }
00141   for (i=sel->firstsel; i<=sel->lastsel; i++) {
00142     if (sel->on[i]) {
00143       int color = ac->color[i];
00144       lookups[color].append(i);
00145     }
00146   }
00147 }
00148   
00149 void DrawMolItem::set_pbc(int pbc) {
00150   cmdList->pbc = pbc;
00151   change_pbc(); // tell rep to update PBC transformation matrices next prepare
00152 
00153   // tell the molecule to notify its geometry monitors that the number of
00154   // periodic images has changed so that they can be turned off or on if
00155   // necessary.  XXX This may not be currently implemented in GeometryMol.
00156   mol->notify();
00157 }
00158 
00159 int DrawMolItem::get_pbc() const {
00160   return cmdList->pbc;
00161 }
00162 
00163 void DrawMolItem::set_pbc_images(int n) {
00164   if (n < 1) return;
00165   cmdList->npbc = n;
00166   need_matrix_recalc();
00167   mol->notify();
00168 }
00169 int DrawMolItem::get_pbc_images() const {
00170   return cmdList->npbc;
00171 }
00172 
00173 
00174 void DrawMolItem::set_instances(int inst) {
00175   cmdList->instanceset = inst;
00176   change_instances(); // tell rep to update instance matrices next prepare
00177 
00178   // tell the molecule to notify its geometry monitors that the number of
00179   // instance images has changed so that they can be turned off or on if
00180   // necessary.  XXX This may not be currently implemented in GeometryMol.
00181   mol->notify();
00182 }
00183 
00184 int DrawMolItem::get_instances() const {
00185   return cmdList->instanceset;
00186 }
00187 
00188 
00189 static void parse_frames(const char *beg, int len, int maxframe, 
00190     ResizeArray<int>& frames) {
00191   if (!len) {
00192     msgErr << "parse_frames got zero-length string!" << sendmsg;
00193     return;
00194   }
00195   const char *firstsep = NULL, *secondsep = NULL;
00196   int i;
00197   for (i=0; i<len; i++) {
00198     if (beg[i] == ':') {
00199       firstsep = beg+i;
00200       break;
00201     }
00202   }
00203   for (++i; i < len; i++) {
00204     if (beg[i] == ':') {
00205       secondsep = beg+i;
00206       break;
00207     }
00208   }
00209   char *endptr;
00210   long first, second, third;
00211   first = strtol(beg, &endptr, 0);
00212   if (endptr == beg) {
00213     msgErr << "frame element is invalid: " << beg << sendmsg;
00214     return;
00215   }
00216   second = third = first; // quiet compiler warnings
00217   if (firstsep) {
00218     firstsep++;
00219     // if no text after separator, then assume last frame.
00220     if (!(*firstsep)) {
00221       second = maxframe;
00222     } else {
00223       second = strtol(firstsep, &endptr, 0);
00224       if (endptr == firstsep) {
00225         msgErr << "frame element is invalid: " << beg << sendmsg;
00226         return;
00227       }
00228     }
00229   }
00230   if (secondsep) {
00231     secondsep++;
00232     if (!(*secondsep)) {
00233       third = maxframe;
00234     } else {
00235       third = strtol(secondsep, &endptr, 0);
00236       if (endptr == secondsep) {
00237         msgErr << "frame element is invalid: " << beg << sendmsg;
00238         return;
00239       }
00240     }
00241   }
00242   // append frames based on what we got:
00243   if (!firstsep) {
00244     //printf("adding frame %d\n", (int)first);
00245     frames.append((int)first);
00246   } else if (!secondsep) {
00247     for (long i=first; i <= second; i++) {
00248       //printf("adding frame %d\n", (int)i);
00249       frames.append((int)i);
00250     }
00251   } else {
00252     // prevent zero, or negative frame step sizes
00253     if (second < 1) {
00254       msgErr << "zero or negative step size is invalid: " << second << sendmsg;
00255       return;
00256     }
00257     // printf("first: %ld third: %ld second: %ld\n", first, third, second);
00258 
00259     for (long i=first; i <= third; i += second) {
00260       //printf("adding frame %d\n", (int)i);
00261       frames.append((int)i);
00262     } 
00263   }
00264 }
00265 
00266 void DrawMolItem::set_drawframes(const char *frames) {
00267   if (!frames) return;
00268   while (isspace(*frames)) frames++;
00269   if (!strncmp(frames, "now", 3) && !strncmp(framesel, "now", 3)) return;
00270   delete [] framesel;
00271   framesel = stringdup(frames);
00272   needRegenerate |= MOL_REGEN;
00273 }
00274 
00276 
00277 void DrawMolItem::create_cmdlist() {
00278   int i;
00279 
00280   // do we need to recreate everything?
00281   if (needRegenerate) {
00282     // Identify the representation number for output comments
00283     repNumber = -1;   // Nonsense default number
00284     for (i = 0; i < mol->components(); i++) {
00285       if (this == mol->component(i))
00286         repNumber = i;
00287     }
00288 
00289     if (needRegenerate & COL_REGEN)
00290       atomColor->find(mol); // execute any pending color updates
00291 
00292     if (needRegenerate & SEL_REGEN)
00293       atomSel->change(NULL, mol); // update atom selection onoff flags
00294 
00295     if (needRegenerate & REP_REGEN) {
00296       // dup cmdStr to prevent overlapping strcpy() buffer pointers
00297       char *newcmdstr = stringdup(atomRep->cmdStr);
00298       atomRep->change(newcmdstr); // update the representation text
00299       delete [] newcmdstr;
00300     }
00301  
00302     reset_disp_list(); // regenerate both data block and display commands
00303 
00304     //
00305     // Record the start of a new representation geometry group
00306     // This can only be called once per rep, as the generated 
00307     // group names need to be unique.
00308     char repbuf[2048];
00309     sprintf(repbuf, "vmd_mol%d_rep%d", mol->id(), repNumber);
00310     cmdBeginRepGeomGroup.putdata(repbuf, cmdList);
00311 
00312     // If coloring by a volume texture, setup the graphics state
00313     //
00314     // We want to generate the fewest number of volume textures possible,
00315     // so we actually put this code outside all per-timestep and
00316     // per-pbc image drawing loops.  That way we'll minimize the number of
00317     // OpenGL texture downloads.  I might be able to make low level renderer
00318     // smart enough to recognize textures by some unique serial number, but
00319     // that code doesn't exist yet, so the next best thing is to be extra 
00320     // careful in generating them in the first place.
00321     //
00322     // XXX One gotcha WRT volume textures, is that since we clamp them at
00323     //     the edge of the texture, there is no provision to make PBC images
00324     //     of volumes render nicely.  Since the PBC box isn't necessarily the
00325     //     same coordinate system as the volume set, not much we can do there.
00326     //     We'd have to implement a separate PBC transform feature for the
00327     //     volumetric data in order to correctly solve this.  For now, we'll
00328     //     just live without it.
00329     if (atomColor->method() == AtomColor::VOLUME ||
00330         atomRep->method() == AtomRep::VOLSLICE) {
00331 
00332       updateVolumeTexture();
00333       // Pass a pointer to the texture map, retained in shared memory
00334       float v0[3], v1[3], v2[3], v3[3];
00335       volumeTexture.calculateTexgenPlanes(v0, v1, v2, v3);
00336       DispCmdVolumeTexture cmdVolTexture;
00337       cmdVolTexture.putdata(
00338           volumeTexture.getTextureID(),
00339           volumeTexture.getTextureSize(),
00340           volumeTexture.getTextureMap(),
00341           v0, v1, v2, v3,
00342           cmdList);
00343     }
00344 
00345     if (strcmp(framesel, "now")) {
00346       // parse the frame selection into a list of frames to process 
00347       ResizeArray<int> frames;
00348       const char *endptr = framesel;
00349       int maxframe = mol->numframes() - 1;
00350       do {
00351         const char *begptr = endptr;
00352         while (isspace(*begptr)) begptr++;
00353         endptr = begptr;
00354         while (*endptr && !isspace(*endptr)) endptr++;
00355         parse_frames(begptr, endptr-begptr, maxframe, frames);
00356       } while (*endptr);
00357 
00358       int curframe = mol->frame(); // save current frame for later
00359 
00360       // loop over the selected frames and render them all
00361       for (i=0; i<frames.num(); i++) {
00362         int frame = frames[i];
00363         if (frame < 0 || frame > maxframe) 
00364           continue;
00365 
00366         mol->override_current_frame(frame);  // override current frame
00367 
00368         // force timestep color to update, redraw or recolor affected geometry
00369         // XXX for the POS[XYZ] coloring methods, the user will need to 
00370         // select "update color every frame" in order to get time-varying color
00371         switch (atomColor->method()) {
00372           case AtomColor::USER:
00373           case AtomColor::USER2:
00374           case AtomColor::USER3:
00375           case AtomColor::USER4:
00376           case AtomColor::PHYSICALTIME:
00377           case AtomColor::TIMESTEP:
00378           case AtomColor::VELOCITY:
00379             atomColor->find(mol);  // recalc atom colors and redraw
00380             break;
00381         }   
00382         do_create_cmdlist();     // draw selected timestep(s)
00383       }
00384       mol->override_current_frame(curframe); // restore previous frame
00385     } else {
00386       do_create_cmdlist();       // draw the current timestep only
00387     }
00388 
00389     // If we have a volume texture enabled, we must turn it back off here
00390     if (atomColor->method() == AtomColor::VOLUME ||
00391         atomRep->method() == AtomRep::VOLSLICE) {
00392       append(DVOLTEXOFF);
00393     }
00394 
00395     cacheskip(1);
00396   } else {
00397     cacheskip(0);
00398   }
00399 
00400   // XXX work around a bug with texturing and display list caching
00401   //     The OpenGL texture object is not stored in the display list 
00402   //     unless the rep was created while display list caching was enabled.
00403   //     If the rep was created beforehand, then the display list will not
00404   //     contain the texture data, and thus things will be drawn incorrectly.
00405   //     For now we'll take the easy way out and disable caching of geometry
00406   //     that is textured.  The right thing to do would be to add a new 
00407   //     flag to tell the lowest level render() routine to regen the textures
00408   //     at the time display lists are created.  This is a reasonable short
00409   //     term solution however.
00410   if (atomColor->method() == AtomColor::VOLUME ||
00411       atomRep->method() == AtomRep::VOLSLICE) {
00412     cacheskip(1); // don't allow this rep to be cached
00413   }
00414 
00415   needRegenerate = NO_REGEN; // done regenerating
00416 }
00417 
00418 
00419 // regenerate the command list
00420 void DrawMolItem::do_create_cmdlist(void) {
00421   float *framepos;
00422   int i;
00423 
00424   // only render frame-oriented reps if there is a current frame
00425   if (mol->current()) {
00426     if (atomSel->do_update)
00427       atomSel->change(NULL, mol);     // update atom selection if necessary
00428 
00429     if (atomColor->do_update)
00430       atomColor->find(mol);           // update colors if necessary
00431 
00432     framepos = (mol->current())->pos; // get atom coordinates for this frame
00433 
00434     // average atom coordinates using a window of size ((2*avgsize) + 1)
00435     if (avgsize && atomSel->selected > 0) {
00436       const int curframe = mol->frame();
00437       const int begframe = curframe - avgsize;
00438       const int endframe = curframe + avgsize;
00439       const int lastframe = mol->numframes() - 1;
00440 #if 0
00441       // Only smooth coordinates of selected atoms
00442       // XXX This would be a big performance gain in most cases, but
00443       //     we can't use it if we're drawing a rep that references
00444       //     coordinates for atoms outside the selection...
00445       //     In order to use this, we'll have to validate that the
00446       //     active representation doesn't touch non-selected atom coords
00447       const int atomloopfirst = atomSel->firstsel;
00448       const int atomlooplast  = 3L*(atomSel->lastsel+1);
00449 #else
00450       // Smooth all atom coordinates
00451       const int atomloopfirst = 0;
00452       const int atomlooplast  = 3L*mol->nAtoms;
00453 #endif
00454       const float rescale = 1.0f/(2.0f*avgsize+1.0f);
00455 
00456 #define VMDPBCSMOOTH 1
00457 #ifdef VMDPBCSMOOTH
00458       // for pbc aware averaging
00459       int isortho,usepbc;
00460       float a,b,c,alpha,beta,gamma;
00461 
00462       // set some arbitrary defaults
00463       isortho=usepbc=1;
00464       a=b=c=9999999.0f;
00465       alpha=beta=gamma=90.0f;
00466 
00467       // reference timestep for pbc wrap
00468       const float *ref = mol->get_frame(curframe)->pos;
00469 #endif      
00470 
00471       memset(avg, 0, 3L*mol->nAtoms*sizeof(float)); // clear average array
00472 
00473       for (i=begframe; i<=endframe; i++) {
00474         int ind = i;
00475         if (ind < 0) 
00476           ind = 0;
00477         else if (ind > lastframe) 
00478           ind = lastframe;
00479         const float *ts = mol->get_frame(ind)->pos;
00480 
00481 #ifdef VMDPBCSMOOTH
00482         // get periodic cell information for current frame
00483         a = mol->get_frame(ind)->a_length;
00484         b = mol->get_frame(ind)->b_length;
00485         c = mol->get_frame(ind)->c_length;
00486         alpha = mol->get_frame(ind)->alpha;
00487         beta = mol->get_frame(ind)->beta;
00488         gamma = mol->get_frame(ind)->gamma;
00489 
00490         // check validity of PBC cell side lengths
00491         if (fabsf(a*b*c) < 0.0001)
00492           usepbc=0;
00493 
00494         // check PBC unit cell shape to select proper low level algorithm.
00495         if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0))
00496           isortho=0;
00497         
00498         // until we can handle non-orthogonal periodic cells turn off pbc handling
00499         if (!isortho)
00500           usepbc=0;
00501 
00502         if (usepbc) {
00503           const float ahalf=a*0.5f;
00504           const float bhalf=b*0.5f;
00505           const float chalf=c*0.5f;
00506 
00507           // smooth by offsetting with computed ifs using the current frame as reference.
00508           // (in the hope that the compiler can optimize this better than regular nested ifs).
00509           // only works for orthogonal cells.
00510           for (int j=atomloopfirst; j<atomlooplast; j += 3) {
00511             float adiff=ts[j  ]-ref[j  ];
00512             avg[j  ] += ts[j  ] + ((adiff < -ahalf) ? a : ( (adiff > ahalf) ? -a : 0));
00513             float bdiff=ts[j+1]-ref[j+1];
00514             avg[j+1] += ts[j+1] + ((bdiff < -bhalf) ? b : ( (bdiff > bhalf) ? -b : 0));
00515             float cdiff=ts[j+2]-ref[j+2];
00516             avg[j+2] += ts[j+2] + ((cdiff < -chalf) ? c : ( (cdiff > chalf) ? -c : 0));
00517           }
00518         } else {                // normal smoothing
00519           for (int j=atomloopfirst; j<atomlooplast; j++) {
00520             avg[j] += ts[j];
00521           }
00522         }
00523 #else 
00524         for (int j=atomloopfirst; j<atomlooplast; j++) {
00525           avg[j] += ts[j];
00526         }
00527 #endif
00528   
00529       }
00530 
00531       // scale results appropriately
00532       for (int j=atomloopfirst; j<atomlooplast; j++)
00533         avg[j] *= rescale;
00534 
00535       framepos = avg; // framepos points to the new average values
00536     }
00537 
00538     // now put in drawing commands, which index into the above block
00539     switch (atomRep->method()) {
00540       case AtomRep::LINES:        
00541         draw_lines(framepos, 
00542           (int)atomRep->get_data(AtomRep::LINETHICKNESS),
00543           atomRep->get_data(AtomRep::ISOLINETHICKNESS));
00544         place_picks(framepos);
00545         break;
00546 
00547       case AtomRep::BONDS:
00548         draw_bonds(framepos, 
00549           atomRep->get_data(AtomRep::BONDRAD),
00550           (int)atomRep->get_data(AtomRep::BONDRES),
00551           atomRep->get_data(AtomRep::ISOLINETHICKNESS));
00552         place_picks(framepos);
00553         break;
00554 
00555       case AtomRep::DYNAMICBONDS: 
00556         draw_dynamic_bonds(framepos, 
00557           atomRep->get_data(AtomRep::BONDRAD),
00558           (int)atomRep->get_data(AtomRep::BONDRES),
00559           atomRep->get_data(AtomRep::SPHERERAD)); 
00560         break;
00561 
00562       case AtomRep::HBONDS:
00563         draw_hbonds(framepos,
00564           atomRep->get_data(AtomRep::SPHERERAD),
00565           (int)atomRep->get_data(AtomRep::LINETHICKNESS),
00566           atomRep->get_data(AtomRep::BONDRAD)); 
00567         break;
00568 
00569       case AtomRep::POINTS:       
00570         draw_points(framepos, atomRep->get_data(AtomRep::LINETHICKNESS)); 
00571         place_picks(framepos);
00572         break;
00573 
00574       case AtomRep::VDW:          
00575         draw_solid_spheres(framepos,
00576           (int)atomRep->get_data(AtomRep::SPHERERES),
00577           atomRep->get_data(AtomRep::SPHERERAD),
00578           0.0);
00579         place_picks(framepos);
00580         break;
00581 
00582       case AtomRep::CPK:          
00583         draw_cpk_licorice(framepos, 1, 
00584           atomRep->get_data(AtomRep::BONDRAD),           // bond rad
00585           (int)atomRep->get_data(AtomRep::BONDRES),      // bond res
00586           atomRep->get_data(AtomRep::SPHERERAD),         // scaled VDW rad
00587           (int)atomRep->get_data(AtomRep::SPHERERES),    // sphere res
00588           (int)atomRep->get_data(AtomRep::LINETHICKNESS),// line thickness
00589           atomRep->get_data(AtomRep::ISOLINETHICKNESS)); // bonds cutoff
00590         place_picks(framepos);
00591         break;
00592 
00593       case AtomRep::LICORICE:     
00594         draw_cpk_licorice(framepos, 0, 
00595           atomRep->get_data(AtomRep::BONDRAD),            // bond rad
00596           (int)atomRep->get_data(AtomRep::BONDRES),       // bond res
00597           atomRep->get_data(AtomRep::SPHERERAD),          // scaled VDW rad
00598           (int)atomRep->get_data(AtomRep::SPHERERES),     // sphere res
00599           (int)atomRep->get_data(AtomRep::LINETHICKNESS), // line thickness
00600           atomRep->get_data(AtomRep::ISOLINETHICKNESS));  // bonds cutoff
00601         place_picks(framepos);
00602         break;
00603 
00604 #ifdef VMDPOLYHEDRA
00605       case AtomRep::POLYHEDRA:     
00606         draw_polyhedra(framepos, atomRep->get_data(AtomRep::SPHERERAD)); 
00607         break;
00608 #endif
00609 
00610       case AtomRep::TRACE:
00611         draw_trace(framepos,
00612           atomRep->get_data(AtomRep::BONDRAD),
00613           (int)atomRep->get_data(AtomRep::BONDRES),
00614           (int)atomRep->get_data(AtomRep::LINETHICKNESS));
00615         place_picks(framepos);
00616         break;
00617 
00618       case AtomRep::TUBE:
00619         draw_tube(framepos, 
00620           atomRep->get_data(AtomRep::BONDRAD),
00621           (int)atomRep->get_data(AtomRep::BONDRES));
00622         break;
00623 
00624       case AtomRep::RIBBONS:
00625         draw_ribbons(framepos, 
00626           atomRep->get_data(AtomRep::BONDRAD) / 3.0f,
00627           (int)atomRep->get_data(AtomRep::BONDRES),
00628           atomRep->get_data(AtomRep::LINETHICKNESS));
00629         break;
00630 
00631       case AtomRep::NEWRIBBONS:
00632         draw_ribbons_new(framepos,
00633           atomRep->get_data(AtomRep::BONDRAD) / 3.0f,
00634           (int)atomRep->get_data(AtomRep::BONDRES),
00635           (int)atomRep->get_data(AtomRep::SPHERERAD), // use bspline or not
00636           atomRep->get_data(AtomRep::LINETHICKNESS));
00637         break;
00638 
00639 #ifdef VMDWITHCARBS
00640       case AtomRep::RINGS_PAPERCHAIN:
00641         draw_rings_paperchain(framepos,
00642                               atomRep->get_data(AtomRep::LINETHICKNESS), // bipyramid_height
00643                               (int)atomRep->get_data(AtomRep::ISOSTEPSIZE) // maximum ring size
00644                               ); 
00645         place_picks(framepos); // XXX add better pick points in the ring traversal code
00646         break;
00647 
00648       case AtomRep::RINGS_TWISTER:
00649         draw_rings_twister(framepos,
00650                            (int)atomRep->get_data(AtomRep::LINETHICKNESS), // start_end_centroid
00651                            (int)atomRep->get_data(AtomRep::BONDRES), // hide_shared_links
00652                            (int)atomRep->get_data(AtomRep::SPHERERAD), // rib_steps
00653                            atomRep->get_data(AtomRep::BONDRAD), // rib_width
00654                            atomRep->get_data(AtomRep::SPHERERES), // rib_height
00655                            (int)atomRep->get_data(AtomRep::ISOSTEPSIZE), // maximum ring size
00656                            (int)atomRep->get_data(AtomRep::ISOLINETHICKNESS) // maximum link length
00657                            );
00658         place_picks(framepos); // XXX add better pick points in the ring traversal code
00659         break;
00660 #endif
00661 
00662       case AtomRep::NEWCARTOON:   
00663         draw_cartoon_ribbons(framepos, 
00664                              (int)atomRep->get_data(AtomRep::BONDRES),
00665                              atomRep->get_data(AtomRep::BONDRAD), 
00666                              (float) atomRep->get_data(AtomRep::LINETHICKNESS),
00667                              1,
00668                              (int)atomRep->get_data(AtomRep::SPHERERAD));
00669         break;
00670 
00671       case AtomRep::STRUCTURE:
00672         draw_structure(framepos,
00673           atomRep->get_data(AtomRep::BONDRAD),
00674           (int)atomRep->get_data(AtomRep::BONDRES),
00675           (int)atomRep->get_data(AtomRep::LINETHICKNESS));
00676         break;
00677 
00678 #if defined(VMDNANOSHAPER)
00679       case AtomRep::NANOSHAPER:
00680         draw_nanoshaper(framepos, 
00681           (int)atomRep->get_data(AtomRep::LINETHICKNESS), // surftype
00682           (int)atomRep->get_data(AtomRep::BONDRES),       // draw wireframe
00683           atomRep->get_data(AtomRep::GRIDSPACING),        // grid spacing
00684           atomRep->get_data(AtomRep::SPHERERAD),          // probe radius
00685           atomRep->get_data(AtomRep::SPHERERES),          // skin parm
00686           atomRep->get_data(AtomRep::BONDRAD));           // blob parm
00687         break;
00688 #endif
00689 #ifdef VMDMSMS
00690       case AtomRep::MSMS:
00691         draw_msms(framepos, 
00692           (int)atomRep->get_data(AtomRep::BONDRES),    // draw wireframe
00693           (atomRep->get_data(AtomRep::LINETHICKNESS) < 0.5), // all / selected
00694           atomRep->get_data(AtomRep::SPHERERAD),  // probe radius
00695           atomRep->get_data(AtomRep::SPHERERES)); // point density 
00696         break;
00697 #endif
00698 #ifdef VMDSURF   
00699       case AtomRep::SURF:
00700         draw_surface(framepos,
00701           (int)atomRep->get_data(AtomRep::BONDRES),    // draw wireframe
00702           atomRep->get_data(AtomRep::SPHERERAD)); // probe radius
00703         break;
00704 #endif
00705 #ifdef VMDQUICKSURF   
00706       case AtomRep::QUICKSURF:
00707         draw_quicksurf(framepos,
00708           int(atomRep->get_data(AtomRep::BONDRES)), // quality level
00709           atomRep->get_data(AtomRep::SPHERERAD),    // sphere radius scale
00710           atomRep->get_data(AtomRep::BONDRAD),      // density isovalue
00711           atomRep->get_data(AtomRep::GRIDSPACING)); // grid spacing
00712         break;
00713 #endif
00714 
00715       case AtomRep::VOLSLICE:
00716         draw_volslice((int)atomRep->get_data(AtomRep::SPHERERES),
00717           atomRep->get_data(AtomRep::SPHERERAD),
00718           (int)atomRep->get_data(AtomRep::LINETHICKNESS),
00719           (int)atomRep->get_data(AtomRep::BONDRES));
00720         break;
00721 
00722       case AtomRep::FIELDLINES:
00723         // XXX use environment variables to control new parameters until
00724         // we have a chance to update the GUI code...
00725         draw_volume_field_lines((int)atomRep->get_data(AtomRep::SPHERERES),
00726           int(atomRep->get_data(AtomRep::FIELDLINESEEDUSEGRID)),
00727           (getenv("VMDFIELDLINEMAXSEEDS")) ?  atoi(getenv("VMDFIELDLINEMAXSEEDS")) : 50000,
00728           atomRep->get_data(AtomRep::SPHERERAD),
00729           atomRep->get_data(AtomRep::FIELDLINEDELTA),
00730           atomRep->get_data(AtomRep::BONDRAD),
00731           atomRep->get_data(AtomRep::BONDRES),
00732           int(atomRep->get_data(AtomRep::FIELDLINESTYLE)),
00733           (getenv("VMDFIELDLINETUBERES")) ? atoi(getenv("VMDFIELDLINETUBERES")) : 12,
00734           (float) atomRep->get_data(AtomRep::LINETHICKNESS));
00735         break;
00736 
00737       case AtomRep::ISOSURFACE:   
00738         draw_isosurface((int)atomRep->get_data(AtomRep::SPHERERES),
00739           atomRep->get_data(AtomRep::SPHERERAD),
00740           (int)atomRep->get_data(AtomRep::LINETHICKNESS),
00741           (int)atomRep->get_data(AtomRep::BONDRES),
00742           (int)atomRep->get_data(AtomRep::ISOSTEPSIZE),
00743           (int)atomRep->get_data(AtomRep::ISOLINETHICKNESS));
00744         break;
00745 
00746       case AtomRep::ORBITAL:   
00747         draw_orbital(
00748 #if 1
00749                      (getenv("VMDMODENSITY") != NULL),
00750 #else
00751                      0,
00752 #endif
00753                      (int)atomRep->get_data(AtomRep::WAVEFNCTYPE),
00754                      (int)atomRep->get_data(AtomRep::WAVEFNCSPIN),
00755                      (int)atomRep->get_data(AtomRep::WAVEFNCEXCITATION),
00756                      (int)atomRep->get_data(AtomRep::SPHERERES),
00757                      atomRep->get_data(AtomRep::SPHERERAD),
00758                      (int)atomRep->get_data(AtomRep::LINETHICKNESS),
00759                      (int)atomRep->get_data(AtomRep::BONDRES),
00760                      atomRep->get_data(AtomRep::GRIDSPACING),
00761                      (int)atomRep->get_data(AtomRep::ISOSTEPSIZE),
00762                      (int)atomRep->get_data(AtomRep::ISOLINETHICKNESS));
00763         break;
00764 
00765       case AtomRep::BEADS:
00766         draw_residue_beads(framepos,
00767           (int)atomRep->get_data(AtomRep::SPHERERES),
00768           atomRep->get_data(AtomRep::SPHERERAD));
00769         break;
00770 
00771       case AtomRep::DOTTED:       
00772         draw_dotted_spheres(framepos, 
00773           atomRep->get_data(AtomRep::SPHERERAD),
00774           (int)atomRep->get_data(AtomRep::SPHERERES)); 
00775         place_picks(framepos);
00776         break;
00777 
00778       case AtomRep::SOLVENT:
00779         draw_dot_surface(framepos,
00780           atomRep->get_data(AtomRep::SPHERERAD),
00781           (int)atomRep->get_data(AtomRep::SPHERERES),
00782           (int)atomRep->get_data(AtomRep::LINETHICKNESS) - 1); // method
00783         place_picks(framepos);
00784         break;
00785 
00786 #ifdef VMDLATTICECUBES
00787       case AtomRep::LATTICECUBES:          
00788         draw_solid_cubes(framepos, atomRep->get_data(AtomRep::SPHERERAD));
00789         place_picks(framepos);
00790         break;
00791 #endif
00792 
00793       default:
00794         msgErr << "Illegal atom representation in DrawMolecule." << sendmsg;
00795     }
00796   } else {
00797     // do reps that don't require a current frame
00798     switch (atomRep->method()) {
00799       case AtomRep::VOLSLICE:
00800         draw_volslice((int)atomRep->get_data(AtomRep::SPHERERES),
00801           atomRep->get_data(AtomRep::SPHERERAD),
00802           (int)atomRep->get_data(AtomRep::LINETHICKNESS),
00803           (int)atomRep->get_data(AtomRep::BONDRES));
00804         break;
00805 
00806       case AtomRep::ISOSURFACE:   
00807         draw_isosurface((int)atomRep->get_data(AtomRep::SPHERERES),
00808           atomRep->get_data(AtomRep::SPHERERAD),
00809           (int)atomRep->get_data(AtomRep::LINETHICKNESS),
00810           (int)atomRep->get_data(AtomRep::BONDRES),
00811           (int)atomRep->get_data(AtomRep::ISOSTEPSIZE),
00812           (int)atomRep->get_data(AtomRep::ISOLINETHICKNESS));
00813         break;
00814 
00815       case AtomRep::FIELDLINES:
00816         // XXX use environment variables to control new parameters until
00817         // we have a chance to update the GUI code...
00818         draw_volume_field_lines((int)atomRep->get_data(AtomRep::SPHERERES),
00819           int(atomRep->get_data(AtomRep::FIELDLINESEEDUSEGRID)),
00820           (getenv("VMDFIELDLINEMAXSEEDS")) ?  atoi(getenv("VMDFIELDLINEMAXSEEDS")) : 50000,
00821           atomRep->get_data(AtomRep::SPHERERAD),
00822           atomRep->get_data(AtomRep::FIELDLINEDELTA),
00823           atomRep->get_data(AtomRep::BONDRAD),
00824           atomRep->get_data(AtomRep::BONDRES),
00825           int(atomRep->get_data(AtomRep::FIELDLINESTYLE)),
00826           (getenv("VMDFIELDLINETUBERES")) ? atoi(getenv("VMDFIELDLINETUBERES")) : 12,
00827           (float) atomRep->get_data(AtomRep::LINETHICKNESS));
00828         break;
00829     }
00830   }
00831 }
00832 
00833 
00834 // enable pick points for all atoms that are turned on
00835 void DrawMolItem::place_picks(float *pos) {
00836   DispCmdPickPointArray cmdPickPointArray;
00837 
00838   // optimize pick point generation by taking advantage of
00839   // the firstsel/lastsel indices for the active selection
00840   int selseglen = atomSel->lastsel-atomSel->firstsel+1;
00841   cmdPickPointArray.putdata(selseglen, atomSel->selected, atomSel->firstsel,
00842                             &atomSel->on[atomSel->firstsel], 
00843                             pos + 3L*atomSel->firstsel, cmdList);
00844 }
00845 
00846 
00848 
00849 void DrawMolItem::do_color_changed(int ccat) {
00850   // check Display and Axes categories, since volumetric and other
00851   // representations that uses these colors will need regeneration
00852   if ((ccat == scene->category_index("Display")) ||
00853       (ccat == scene->category_index("Axes"))) {
00854     needRegenerate |= COL_REGEN; // recalc colors and redraw
00855   }
00856 
00857   // check the atom category ... see if it's for the current coloring method
00858   if (atomColor && atomColor->current_color_use(ccat)) {
00859     change_color(atomColor);    // recalc atom colors and redraw
00860   }
00861 }
00862 
00863 void DrawMolItem::do_color_rgb_changed(int color) {
00864   int ccat;
00865 
00866   // check Display and Axes categories, since volumetric and other
00867   // representations that uses these colors will need regeneration
00868   ccat = scene->category_index("Display");
00869   if ((color == scene->category_item_value(ccat, "Background")) ||
00870       (color == scene->category_item_value(ccat, "Foreground")))
00871     needRegenerate |= COL_REGEN; // recalc colors and redraw
00872 
00873   ccat = scene->category_index("Axes");
00874   if ((color == scene->category_item_value(ccat, "X")) ||
00875       (color == scene->category_item_value(ccat, "Y")) ||
00876       (color == scene->category_item_value(ccat, "Z")) ||
00877       (color == scene->category_item_value(ccat, "Origin")) ||
00878       (color == scene->category_item_value(ccat, "Labels")))
00879     needRegenerate |= COL_REGEN; // recalc colors and redraw
00880 
00881   // check all atom colors 
00882   const int *colors = atomColor->color;
00883   for (int i=0; i<mol->nAtoms; i++) {
00884     if (colors[i] == color) {
00885       change_color(atomColor);   // recalc atom colors and redraw
00886       break;
00887     }
00888   }
00889 }
00890 
00891 void DrawMolItem::do_color_scale_changed() {
00892   if (atomColor->uses_colorscale()) {
00893     atomColor->find(mol);
00894     needRegenerate |= COL_REGEN; // recalc atom colors and redraw
00895   }
00896 }
00897   
00899 
00900 // change the atom coloring method.  Return success.
00901 int DrawMolItem::change_color(AtomColor *ac) {
00902   if (ac) {
00903     *atomColor = *ac;
00904     needRegenerate |= COL_REGEN;
00905     return TRUE;
00906   } else 
00907     return FALSE;
00908 }
00909 
00910 // return which representation this is, i.e. an index from 0 ... # reps - 1
00911 int DrawMolItem::representation_index(void) {
00912   int totalreps = mol->components();
00913   for (int i=0; i < totalreps; i++) 
00914     if (mol->component(i) == this)
00915       return i;
00916 
00917   msgErr << "Unknown molecular representation found." << sendmsg;
00918   return 0; // couldn't find anything that made sense.
00919 }
00920 
00921 // prepare for drawing ... do any updates needed right before draw.
00922 void DrawMolItem::prepare() {
00923   // if the rep is off, do nothing; leave the state flags alone and
00924   // do all required actions when we're turned back on again.
00925   if (!displayed()) return;
00926 
00927   // force throb atom colors to update, recoloring affected geometry 
00928   if (atomColor->method() == AtomColor::THROB) { 
00929     needRegenerate |= COL_REGEN; // recalc atom colors and redraw    
00930   }   
00931 
00932   // update periodic image display
00933   if ((update_pbc || update_ts) && cmdList->pbc) {
00934     update_pbc_transformations();
00935     update_pbc = 0;
00936     need_matrix_recalc(); // sets the _needUpdate flag
00937   }
00938 
00939   // update instance image display
00940   if ((update_instances || update_ts) && cmdList->instanceset) {
00941     update_instance_transformations();
00942     update_instances = 0;
00943     need_matrix_recalc(); // sets the _needUpdate flag
00944   }
00945 
00946   // update secondary structure
00947   if (update_ss) {
00948     if (atomRep->method() == AtomRep::STRUCTURE ||
00949         atomRep->method() == AtomRep::NEWCARTOON)
00950       needRegenerate |= MOL_REGEN;
00951 
00952     if (atomColor->method() == AtomColor::STRUCTURE)
00953       needRegenerate |= COL_REGEN;
00954 
00955     update_ss = 0;
00956   }
00957 
00958   // update timestep
00959   if (update_ts) {
00960     // If we're drawing the current frame, treat the same as MOL_REGEN.
00961     // Otherwise, we're just using cached geometry and there's nothing to do.
00962     if (!strcmp(framesel, "now")) {
00963       needRegenerate |= MOL_REGEN;
00964 
00965       // force timestep color to update, redraw or recolor affected geometry
00966       // XXX for the POS[XYZ] coloring methods, the user will need to 
00967       // select "update color every frame" in order to get time-varying color
00968       switch (atomColor->method()) {
00969         case AtomColor::USER:
00970         case AtomColor::USER2:
00971         case AtomColor::USER3:
00972         case AtomColor::USER4:
00973         case AtomColor::PHYSICALTIME:
00974         case AtomColor::TIMESTEP:
00975         case AtomColor::VELOCITY:
00976           needRegenerate |= COL_REGEN; // recalc atom colors and redraw
00977           break;
00978       }   
00979     }
00980 
00981     update_ts = 0; // update complete
00982   }
00983 
00984   // new trajectory frames, update reps that draw multiple frames
00985   if (update_traj) {
00986     // If we're drawing multiple frames, wait for any molecule I/O to finish
00987     // before recalculating the entire set of geometry.
00988     if (strcmp(framesel, "now")) {  // we have a frame selection...
00989       // XXX ick - cast to Molecule so that we can access Molecule method
00990       Molecule *m = (Molecule *)mol; 
00991       if (!m->file_in_progress()) { // ...and file I/O is done.
00992         needRegenerate |= MOL_REGEN;
00993         update_traj = 0; // update complete
00994       }
00995     }
00996   }
00997 
00998   create_cmdlist();
00999 }
01000 
01001 void DrawMolItem::update_pbc_transformations() {
01002   const Timestep *ts = mol->current();
01003   if (!ts) return;
01004 
01005   cmdList->transX.identity();
01006   cmdList->transY.identity();
01007   cmdList->transZ.identity();
01008   ts->get_transforms(cmdList->transX, cmdList->transY, cmdList->transZ);
01009   cmdList->transXinv = cmdList->transX;
01010   cmdList->transYinv = cmdList->transY;
01011   cmdList->transZinv = cmdList->transZ;
01012   cmdList->transXinv.inverse();
01013   cmdList->transYinv.inverse();
01014   cmdList->transZinv.inverse();
01015 }
01016 
01017 
01018 void DrawMolItem::update_instance_transformations() {
01019   cmdList->instances.clear();
01020   int txcnt = mol->instances.num();
01021 #if 0
01022   printf("drawmolitem::update_instance_trans(): cnt %d\n", txcnt);
01023 #endif
01024   int i;
01025   for (i=0; i<txcnt; i++) {
01026     cmdList->instances.append(mol->instances[i]);
01027   }
01028 }
01029 
01030 
01031 // change the atom rep method.  Return success.
01032 int DrawMolItem::change_rep(AtomRep *ar) {
01033   if (!ar) return FALSE;  // XX does this actually occur?
01034 
01035 #ifdef VMDMSMS
01036   // If we're changing from MSMS to something else, free MSMS memory
01037   if (atomRep->method() == AtomRep::MSMS && ar->method() != AtomRep::MSMS)
01038     msms.clear();
01039 #endif
01040 
01041 #ifdef VMDNANOSHAPER
01042   // If we're changing from NanoShaper to something else, free NanoShaper memory
01043   if (atomRep->method() == AtomRep::NANOSHAPER && ar->method() != AtomRep::NANOSHAPER)
01044     nanoshaper.clear();
01045 #endif
01046 
01047 #ifdef VMDSURF
01048   // If we're changing from Surf to something else, free Surf memory
01049   if (atomRep->method() == AtomRep::SURF && ar->method() != AtomRep::SURF)
01050     surf.clear();
01051 #endif
01052 
01053   // If we're changing from Orbital to something else, free the grid
01054   if (atomRep->method() == AtomRep::ORBITAL && ar->method() != AtomRep::ORBITAL) {
01055     waveftype = -1;
01056     wavefspin = -1;
01057     wavefexcitation = -1;
01058     gridorbid=-1;
01059     orbgridspacing = -1.0f; 
01060     delete orbvol;
01061     orbvol=NULL;
01062   }   
01063 
01064   // Free tubearray if changing from tube or cartoon to something else.
01065   // We need to check for cartoon because cartoon uses draw_tube to draw
01066   // its coils. 
01067   if ((atomRep->method() == AtomRep::TUBE ||
01068        atomRep->method() == AtomRep::STRUCTURE) && 
01069       (ar->method() != AtomRep::TUBE &&
01070        ar->method() != AtomRep::STRUCTURE)) {
01071     if (tubearray) {
01072       for (int i=0; i<tubearray->num(); i++) 
01073         delete (*tubearray)[i];
01074       delete tubearray;
01075       tubearray = NULL;
01076     }
01077   }
01078     
01079   *atomRep = *ar;
01080   needRegenerate |= REP_REGEN;
01081   return TRUE;
01082 }
01083 
01084 
01085 // change the atom selection method.  Return success.
01086 int DrawMolItem::change_sel(const char *cmdStr) {
01087   if (!cmdStr) return TRUE; // no problem
01088   if (atomSel->change(cmdStr, mol) == AtomSel::NO_PARSE) return FALSE;
01089   mol->notify(); // _after_ changing atomSel, tell molecule to update
01090                  // geometry monitors since they need to be switched off
01091                  // if they point to atoms that are no longer on.
01092   needRegenerate |= SEL_REGEN;
01093   return TRUE;
01094 }
01095 
01096 // force a recalculation
01097 void DrawMolItem::force_recalc(int r) {
01098   needRegenerate |= r;
01099 }
01100 
01101 // return whether the Nth atom is displayed in this representation
01102 int DrawMolItem::atom_displayed(int n) {
01103   return (atomSel != NULL && mol != NULL && 
01104           n >= 0 && n < mol->nAtoms && atomSel->on[n]);
01105 }
01106 
01108 void DrawMolItem::draw_lines(float *framepos, int thickness, float cutoff) {
01109   update_lookups(atomColor, atomSel, colorlookups); // update line color table
01110   int *nbonds = NULL;
01111   int *bondlists = NULL;
01112   if (cutoff > 0) {
01113     nbonds = new int[mol->nAtoms];
01114     memset(nbonds, 0, mol->nAtoms*sizeof(int));
01115     bondlists = new int[MAXATOMBONDS * mol->nAtoms];
01116     GridSearchPair *pairlist = vmd_gridsearch1(framepos, mol->nAtoms, atomSel->on, 
01117         cutoff, 0, mol->nAtoms * 27L);
01118     GridSearchPair *p, *tmp;
01119     for (p=pairlist; p != NULL; p=tmp) {
01120       MolAtom *atom1 = mol->atom(p->ind1);
01121       MolAtom *atom2 = mol->atom(p->ind2);
01122   
01123       // don't bond atoms that aren't part of the same conformation    
01124       // or that aren't in the all-conformations part of the structure 
01125       if ((atom1->altlocindex != atom2->altlocindex) &&
01126           ((mol->altlocNames.name(atom1->altlocindex)[0] != '\0') &&
01127           (mol->altlocNames.name(atom2->altlocindex)[0] != '\0'))) {
01128         tmp = p->next;
01129         free(p);
01130         continue;
01131       }
01132       // Prevent hydrogens from bonding with each other.
01133       // Use atomType info derived during initial molecule analysis for speed.
01134       if (atom1->atomType == ATOMHYDROGEN &&
01135           atom2->atomType == ATOMHYDROGEN) {
01136         tmp = p->next;
01137         free(p);
01138         continue;
01139       }
01140       bondlists[p->ind1 * MAXATOMBONDS + nbonds[p->ind1]] = p->ind2;
01141       bondlists[p->ind2 * MAXATOMBONDS + nbonds[p->ind2]] = p->ind1;
01142       nbonds[p->ind1]++;
01143       nbonds[p->ind2]++;
01144       tmp = p->next;
01145       free(p);
01146       continue;
01147     }
01148   }
01149 
01150   sprintf(commentBuffer, "Mol[%d] Rep[%d] Lines", mol->id(), repNumber);
01151   cmdCommentX.putdata(commentBuffer, cmdList);
01152 
01153   // turn off material characteristics, set line style
01154   append(DMATERIALOFF);
01155   cmdLineType.putdata(SOLIDLINE, cmdList);
01156   cmdLineWidth.putdata(thickness, cmdList);
01157 
01158   for (int i=0; i<MAXCOLORS; i++) {
01159     const ColorLookup &cl = colorlookups[i];
01160 
01161     if (cl.num == 0) continue; // skip if no bonds of this color
01162     
01163     cmdColorIndex.putdata(i, cmdList); // set color for these half-bonds
01164   
01165     ResizeArray<float> verts;
01166 
01167     // maintain total vert count for max OpenGL buffer size management
01168     int totalverts=0;
01169 
01170     // loop over all half-bonds to be drawn in this color
01171     for (int j=0; j<cl.num; j++) { 
01172       const int id = cl.idlist[j];
01173       float *fp1 = framepos + 3L*id;
01174       const MolAtom *a1 = mol->atom(id);
01175       int bondsdrawn = 0;
01176       
01177       // draw half-bond to each bonded, displayed partner
01178       // NOTE: bond mid-points are calculated TWICE, once for each atom.
01179       // Why? So each atom is drawn entirely, including all it's half-bonds,
01180       // after a single set-color command, instead of having a set-color
01181       // command for each bond.  This reduces total graphics state changes
01182       // considerably, at the cost of extra calculation here.
01183       int n = cutoff > 0 ? nbonds[id] : a1->bonds;
01184       for (int k=0; k < n; k++) {
01185         int a2n = cutoff > 0 ? bondlists[MAXATOMBONDS*id + k] : a1->bondTo[k];
01186         if (atomSel->on[a2n]) {    // bonded atom displayed?
01187           float *fp2 = framepos + 3L*a2n;
01188           if (atomColor->color[a2n] == i) {
01189             // same color, so just draw a whole bond, but only if the atom
01190             // id of the other atom is higher to avoid drawing it twice.
01191             if (a2n > id) {
01192               verts.append3(&fp1[0]);
01193               verts.append3(&fp2[0]);
01194               totalverts+=2;
01195             }
01196           } else {
01197             float mid[3];
01198             mid[0] = 0.5f * (fp1[0] + fp2[0]); 
01199             mid[1] = 0.5f * (fp1[1] + fp2[1]); 
01200             mid[2] = 0.5f * (fp1[2] + fp2[2]); 
01201             verts.append3(&fp1[0]);
01202             verts.append3(&mid[0]);
01203             totalverts+=2;
01204           }
01205           bondsdrawn++;                // increment counter of bonds to atom i
01206         }
01207       }
01208      
01209       // if atom is all alone (no bonds drawn) draw a point 
01210       if (!bondsdrawn) {  
01211         DispCmdPoint cmdPoint;
01212         cmdPoint.putdata(fp1, cmdList);
01213       }
01214 
01215       // Ensure we don't generate a vertex buffer so large that 
01216       // we'd crash from internal integer overflows in renderer 
01217       // subclasses or back-end renderers
01218       if (totalverts > VMDMAXVERTEXBUFSZ) {
01219         DispCmdLineArray cmdLineArray;
01220         cmdLineArray.putdata(&verts[0], verts.num()/6, cmdList);
01221         verts.clear(); // clear vertex buffer
01222         totalverts=0;  // reset counter
01223       }
01224     }
01225 
01226     if (verts.num()) {
01227       DispCmdLineArray cmdLineArray;
01228       cmdLineArray.putdata(&verts[0], verts.num()/6, cmdList);
01229     }
01230   }
01231   if (cutoff > 0) {
01232     delete [] nbonds;
01233     delete [] bondlists;
01234   }
01235 }
01236 
01237 
01238 // draw all lattice site particles as cubes
01239 // radius is normally set to half of lattice site cube side length dimension
01240 void DrawMolItem::draw_solid_cubes(float * framepos, float radscale) {
01241   int i;
01242 
01243   // set cube type and shading characteristics
01244   sprintf(commentBuffer,"Mol[%d] Rep[%d] LatticeCubes", mol->id(), repNumber);
01245   cmdCommentX.putdata(commentBuffer, cmdList);
01246   append(DMATERIALON); // turn on lighting
01247 
01248   // maintain total vert count for max OpenGL buffer size management
01249   int totalverts=0;
01250 
01251   // draw cubes using cube array primitive
01252   if (radscale > 0) {           // don't draw zero scaled cubes
01253     long ind = 0;
01254     ResizeArray<float> centers;
01255     ResizeArray<float> radii;
01256     ResizeArray<float> colors;
01257     const float *radius = mol->radius();
01258 
01259     ind = atomSel->firstsel * 3L;    
01260     for (i=atomSel->firstsel; i <= atomSel->lastsel; i++) {
01261       // draw a cube for each selected lattice site particle
01262       if (atomSel->on[i]) {
01263         totalverts++;
01264         float *fp = framepos + ind;
01265         const float *cp; 
01266 
01267         centers.append3(&fp[0]);
01268         radii.append(radius[i]*radscale);
01269      
01270         cp = scene->color_value(atomColor->color[i]);
01271         colors.append3(&cp[0]);
01272       }
01273       ind += 3;
01274 
01275       // Ensure we don't generate a vertex buffer so large that 
01276       // we'd crash from internal integer overflows in renderer 
01277       // subclasses or back-end renderers
01278       if (totalverts > VMDMAXVERTEXBUFSZ) {
01279         cmdCubeArray.putdata((float *) &centers[0], 
01280                              (float *) &radii[0], 
01281                              (float *) &colors[0], 
01282                              radii.num(), 
01283                              cmdList);
01284 
01285         centers.clear(); // clear vertices
01286         radii.clear();   // clear vertices
01287         colors.clear();  // clear vertices
01288         totalverts=0;    // reset counter
01289       }
01290     }
01291    
01292     if (radii.num() > 0) {
01293       cmdCubeArray.putdata((float *) &centers[0], 
01294                            (float *) &radii[0], 
01295                            (float *) &colors[0], 
01296                            radii.num(), 
01297                            cmdList);
01298     }
01299   }
01300 }
01301 
01302 
01303 // draw all atoms as spheres 
01304 // radius is set to vdw radius * radscale + fixrad
01305 // This allows same code to be used for VDW, CPK, Licorice, etc.
01306 void DrawMolItem::draw_solid_spheres(float * framepos, int res, 
01307                                      float radscale, float fixrad) {
01308   int i;
01309 
01310   // set sphere type, resolution, and shading characteristics
01311   sprintf(commentBuffer,"Mol[%d] Rep[%d] VDW", mol->id(), repNumber);
01312   cmdCommentX.putdata(commentBuffer, cmdList);
01313   append(DMATERIALON); // turn on lighting
01314   cmdSphres.putdata(res, cmdList);
01315   cmdSphtype.putdata(SOLIDSPHERE, cmdList);
01316 
01317   const char *modulatefield = NULL; // XXX this needs to become a parameter
01318   const float *modulatedata = NULL; // data field to use for width modulation
01319   int  modulateoffs  = 0;    // data offset
01320   long modulatemult  = 1;    // data multiples
01321 
01322 #if 1
01323   // XXX hack to let users try various stuff
01324   modulatefield = getenv("VMDMODULATERADIUS");
01325 #endif
01326   if (modulatefield != NULL) {
01327     if (!strcmp(modulatefield, "user")) {
01328       modulatedata = mol->current()->user;
01329       // XXX user field can be NULL on some timesteps
01330     } else if (!strcmp(modulatefield, "user2")) {
01331       modulatedata = mol->current()->user2;
01332       // XXX user2 field can be NULL on some timesteps
01333     } else if (!strcmp(modulatefield, "user3")) {
01334       modulatedata = mol->current()->user3;
01335       // XXX user3 field can be NULL on some timesteps
01336     } else if (!strcmp(modulatefield, "user4")) {
01337       modulatedata = mol->current()->user4;
01338       // XXX user1 field can be NULL on some timesteps
01339     } else if (!strcmp(modulatefield, "vx")) {
01340       modulatedata = mol->current()->vel;
01341       // XXX vx field can be NULL on some timesteps
01342       modulateoffs = 0; modulatemult = 3;
01343     } else if (!strcmp(modulatefield, "vy")) {
01344       modulatedata = mol->current()->vel;
01345       // XXX vy field can be NULL on some timesteps
01346       modulateoffs = 1; modulatemult = 3;
01347     } else if (!strcmp(modulatefield, "vz")) {
01348       modulatedata = mol->current()->vel;
01349       // XXX vz field can be NULL on some timesteps
01350       modulateoffs = 2; modulatemult = 3;
01351     } else {
01352       modulatedata = mol->extraflt.data(modulatefield);
01353     } 
01354   }
01355 
01356   // maintain total vert count for max OpenGL buffer size management
01357   int totalverts=0;
01358 
01359   // draw spheres using new sphere array primitive
01360   if ((radscale + fixrad) > 0) {           // don't draw zero scaled spheres
01361     long ind = 0;
01362     ResizeArray<float> centers;
01363     ResizeArray<float> radii;
01364     ResizeArray<float> colors;
01365     const float *radius = mol->radius();
01366 
01367     ind = atomSel->firstsel * 3L;    
01368     for (i=atomSel->firstsel; i <= atomSel->lastsel; i++) {
01369       // draw a sphere for each selected atom
01370       if (atomSel->on[i]) {
01371         totalverts++;
01372         centers.append3(framepos + ind);
01373     
01374         float my_rad = 1.0f;
01375         if (modulatedata != NULL) {
01376           my_rad = modulatedata[i*modulatemult + modulateoffs];
01377           if (my_rad <= 0.0f)
01378             my_rad = 1.0f;
01379         }
01380         radii.append(radius[i]*radscale*my_rad + fixrad);
01381         colors.append3(scene->color_value(atomColor->color[i]));
01382       }
01383       ind += 3;
01384 
01385       // Ensure we don't generate a vertex buffer so large that 
01386       // we'd crash from internal integer overflows in renderer 
01387       // subclasses or back-end renderers
01388       if (totalverts > VMDMAXVERTEXBUFSZ) {
01389         cmdSphereArray.putdata((float *) &centers[0], 
01390                                (float *) &radii[0], 
01391                                (float *) &colors[0], 
01392                                radii.num(), 
01393                                res,
01394                                cmdList);
01395 
01396         centers.clear(); // clear vertices
01397         radii.clear();   // clear vertices
01398         colors.clear();  // clear vertices
01399         totalverts=0;    // reset counter
01400       }
01401     }
01402    
01403     if (radii.num() > 0) {
01404       cmdSphereArray.putdata((float *) &centers[0], 
01405                              (float *) &radii[0], 
01406                              (float *) &colors[0], 
01407                              radii.num(), 
01408                              res,
01409                              cmdList);
01410     }
01411   }
01412 }
01413 
01414 
01415 // draw residues as beads
01416 // radius is set to vdw radius * radscale + fixrad
01417 // This allows same code to be used for VDW, CPK, Licorice, etc.
01418 void DrawMolItem::draw_residue_beads(float * framepos, int sres, float radscale) {
01419   int i, resid, numres;
01420 
01421   // set sphere type, resolution, and shading characteristics
01422   sprintf(commentBuffer,"Mol[%d] Rep[%d] Beads", mol->id(), repNumber);
01423   cmdCommentX.putdata(commentBuffer, cmdList);
01424   append(DMATERIALON); // turn on lighting
01425   cmdSphres.putdata(sres, cmdList);
01426   cmdSphtype.putdata(SOLIDSPHERE, cmdList);
01427 
01428   // draw spheres using new sphere array primitive
01429   ResizeArray<float> centers;
01430   ResizeArray<float> radii;
01431   ResizeArray<float> colors;
01432   ResizeArray<int> pickpointindices;
01433 
01434   // draw a bead for each residue
01435   numres = mol->residueList.num();
01436   for (resid=0; resid<numres; resid++) {
01437     float com[3] = {0.0, 0.0, 0.0};
01438     const ResizeArray<int> &atoms = mol->residueList[resid]->atoms;
01439     int numatoms = atoms.num();
01440     int oncount = 0;
01441     int pickindex=-1;
01442     
01443     // find COM for residue
01444     for (i=0; i<numatoms; i++) {
01445       int idx = atoms[i];
01446       if (atomSel->on[idx]) {
01447         oncount++;
01448         vec_add(com, com, framepos + 3L*idx);
01449       }
01450     }
01451 
01452     if (oncount < 1)
01453       continue; // exit if there weren't any atoms
01454 
01455     vec_scale(com, 1.0f / (float) oncount, com);
01456 
01457     // find radius of bounding sphere and save last atom index for color
01458     int atomcolorindex=0; // initialize, to please compilers
01459     float boundradsq = 0.0f;
01460 #ifdef BEADELLIPSOID
01461     float avdist = 0.0f;
01462     float majoraxis[3] = {0.0, 0.0, 0.0};
01463 #endif
01464     for (i=0; i<numatoms; i++) {
01465       int idx = atoms[i];
01466       if (atomSel->on[idx]) {
01467         float tmpdist[3];
01468         atomcolorindex = idx;
01469         vec_sub(tmpdist, com, framepos + 3L*idx);
01470         float distsq = dot_prod(tmpdist, tmpdist);
01471         if (distsq > boundradsq) {
01472 #ifdef BEADELLIPSOID
01473           // XXX if we want to draw an ellipsoid rather than a sphere,
01474           // we could keep track of the direction from the COM to the 
01475           // most distant atom so far, and use that as the major axis of
01476           // the ellipsoid, and use the average distance as the minor axis
01477           // for a rotationally symmetric ellipsoid.  
01478           vec_copy(majoraxis, tmpdist);
01479           avdist += sqrtf(distsq);
01480 #endif
01481           boundradsq = distsq;
01482         }
01483       }
01484     }
01485 
01486 #ifdef BEADELLIPSOID
01487     avdist /= (float) oncount;
01488     float cep1[3], cep2[3];
01489     vec_copy(cep1, majoraxis);
01490     vec_scale(cep2, -1, majoraxis);
01491     vec_add(cep1, cep1, com);
01492     vec_add(cep2, cep2, com);
01493 
01494     cmdColorIndex.putdata(atomColor->color[atomcolorindex], cmdList);
01495 
01496     // we'd draw a real ellipsoid here rather than a cylinder, but
01497     // drawing a cylinder is a pretty representative test of the idea.
01498     cmdCylinder.putdata(cep1, cep2, avdist, 8, 0, cmdList);
01499 #endif
01500 
01501     centers.append3(&com[0]);
01502     radii.append(radscale * (sqrtf(boundradsq) + 1.0f));
01503 
01504     const float *cp = scene->color_value(atomColor->color[atomcolorindex]);
01505     colors.append3(&cp[0]);
01506 
01507     // Add a pick point
01508     // for proteins use the CA atom
01509     if (pickindex < 0)
01510       pickindex = mol->find_atom_in_residue("CA", resid); 
01511 
01512     // for nucleic acids, use the C3', C3*, or P atom
01513     if (pickindex < 0)
01514       pickindex = mol->find_atom_in_residue("C3'", resid); 
01515     if (pickindex < 0)
01516       pickindex = mol->find_atom_in_residue("C3*", resid); 
01517     if (pickindex < 0)
01518       pickindex = mol->find_atom_in_residue("P", resid);
01519 
01520     // if nothing better is found, use the first atom
01521     if (pickindex < 0)
01522       pickindex = atoms[0];
01523 
01524     pickpointindices.append(pickindex);
01525   }
01526 
01527 #if !defined(BEADELLIPSOID)
01528   if (radii.num() > 0) {
01529     cmdSphereArray.putdata((float *) &centers[0],
01530                            (float *) &radii[0],
01531                            (float *) &colors[0],
01532                            radii.num(),
01533                            sres,
01534                            cmdList);
01535 
01536     DispCmdPickPointArray pickPointArray;
01537     pickPointArray.putdata(pickpointindices.num(), &pickpointindices[0],
01538                            &centers[0], cmdList);
01539   }
01540 #endif
01541 }
01542 
01543 
01544 
01545 // draw all atoms as dotted spheres of vdw radius
01546 void DrawMolItem::draw_dotted_spheres(float * framepos, float srad, int sres) {
01547   float radscale;
01548   int i;
01549   const float *radius = mol->radius();
01550 
01551   // set sphere type, resolution, and shading characteristics
01552   sprintf(commentBuffer,"Mol[%d] Rep[%d] Dotted", mol->id(), repNumber);
01553   cmdCommentX.putdata(commentBuffer, cmdList);
01554   append(DMATERIALOFF); // turn off lighting
01555   cmdSphres.putdata(sres, cmdList);
01556   cmdSphtype.putdata(POINTSPHERE, cmdList);
01557 
01558   // XXX we should teach this code to sort by color and by radius, 
01559   //     which would allow us to avoid dynamic scaling of the spheres at
01560   //     display time.
01561   radscale = srad;
01562   if (radscale > 0) {                  // don't draw zero scaled spheres
01563     int lastcolor = -1;
01564     for (i=atomSel->firstsel; i <= atomSel->lastsel; i++) {
01565       // draw a sphere for each selected atom
01566       if (atomSel->on[i]) {
01567         if (lastcolor != atomColor->color[i]) {   // set color
01568           lastcolor = atomColor->color[i];
01569           cmdColorIndex.putdata(lastcolor, cmdList);
01570         }
01571       
01572         cmdSphere.putdata(framepos+3L*i, radius[i]*radscale,cmdList);
01573       }
01574     }
01575   }
01576 }
01577 
01578 
01579 // draws each atom as a dot
01580 void DrawMolItem::draw_points(float *framepos, float pointsize) {
01581   sprintf(commentBuffer, "Mol[%d] Rep[%d] Points", mol->id(), repNumber);
01582   cmdCommentX.putdata(commentBuffer, cmdList);
01583   append(DMATERIALOFF);
01584 
01585   if (atomSel->selected > 0) {
01586     // optimize pick point generation by taking advantage of
01587     // the firstsel/lastsel indices for the active selection
01588     int selseglen = atomSel->lastsel-atomSel->firstsel+1;
01589 
01590     // draw spheres and pick points
01591     cmdPointArray.putdata(framepos + 3L*atomSel->firstsel,
01592                           atomColor->color + atomSel->firstsel,
01593                           scene,
01594                           pointsize,
01595                           selseglen,
01596                           atomSel->on + atomSel->firstsel,
01597                           atomSel->selected,
01598                           cmdList);
01599   }
01600   append(DMATERIALON);
01601 
01602 }
01603 
01604 
01605 void DrawMolItem::draw_cpk_licorice(float *framepos, int cpk, float brad, int bres, float srad, int sres, int linethickness, float cutoff) {
01606   MolAtom *a1;
01607   int i, j, a2n;
01608   float radscale, fixrad; 
01609   int lastcolor = -1;
01610   int use_bonds = TRUE;
01611 
01612   if (cpk) { 
01613         brad *= 0.25f;
01614     radscale = srad * 0.25f; // scaled VDW rad
01615       fixrad = 0.0; // no fixed radius
01616   } else {
01617     sprintf(commentBuffer, "Mol[%d] Rep[%d] Licorice", mol->id(), repNumber);
01618     cmdCommentX.putdata(commentBuffer, cmdList);
01619 
01620     radscale = 0.0;  // no VDW scaling
01621       fixrad = brad; // used fixed radius spheres only
01622 
01623     // special case -- if there is no thickness, call the wireframe
01624     if (brad == 0) {
01625       draw_lines(framepos, linethickness, cutoff);
01626       return;
01627     }
01628   }
01629 
01630   if (bres <= 2 || brad < 0.01) { 
01631     use_bonds = FALSE; // draw bonds using lines
01632   }
01633 
01634   int *nbonds = NULL;
01635   int *bondlists = NULL;
01636   if (cutoff > 0) {
01637     nbonds = new int[mol->nAtoms];
01638     memset(nbonds, 0, mol->nAtoms*sizeof(int));
01639     bondlists = new int[MAXATOMBONDS * mol->nAtoms];
01640     GridSearchPair *pairlist = vmd_gridsearch1(framepos, mol->nAtoms, atomSel->on, 
01641         cutoff, 0, mol->nAtoms * 27L);
01642     GridSearchPair *p, *tmp;
01643     for (p=pairlist; p != NULL; p=tmp) {
01644       MolAtom *atom1 = mol->atom(p->ind1);
01645       MolAtom *atom2 = mol->atom(p->ind2);
01646   
01647       // don't bond atoms that aren't part of the same conformation    
01648       // or that aren't in the all-conformations part of the structure 
01649       if ((atom1->altlocindex != atom2->altlocindex) &&
01650           ((mol->altlocNames.name(atom1->altlocindex)[0] != '\0') &&
01651           (mol->altlocNames.name(atom2->altlocindex)[0] != '\0'))) {
01652         tmp = p->next;
01653         free(p);
01654         continue;
01655       }
01656       // Prevent hydrogens from bonding with each other.
01657       // Use atomType info derived during initial molecule analysis for speed.
01658       if (atom1->atomType == ATOMHYDROGEN &&
01659           atom2->atomType == ATOMHYDROGEN) {
01660         tmp = p->next;
01661         free(p);
01662         continue;
01663       }
01664       bondlists[p->ind1 * MAXATOMBONDS + nbonds[p->ind1]] = p->ind2;
01665       bondlists[p->ind2 * MAXATOMBONDS + nbonds[p->ind2]] = p->ind1;
01666       nbonds[p->ind1]++;
01667       nbonds[p->ind2]++;
01668       tmp = p->next;
01669       free(p);
01670       continue;
01671     }
01672   }
01673   sprintf(commentBuffer,"Mol[%d] Rep[%d] CPK", mol->id(), repNumber);
01674   cmdCommentX.putdata(commentBuffer, cmdList);
01675 
01676   append(DMATERIALON);
01677  
01678   // only draw spheres if either radscale or fixrad is nonzero 
01679   if ((radscale + fixrad) > 0) {   
01680     draw_solid_spheres(framepos, sres, radscale, fixrad);
01681   }
01682 
01683   if (use_bonds) {
01684 #if defined(VMDOPTIXRTRT)
01685     ResizeArray <float> cylpoints;
01686     ResizeArray <float> cylradii;
01687     ResizeArray <float> cylcolors;
01688 #endif
01689 
01690     // draw bonds between atoms with cylinders
01691     for (i=atomSel->firstsel; i <= atomSel->lastsel; i++) {
01692       // for each atom, draw half-bond to other displayed atoms
01693       if (atomSel->on[i]) {
01694         float mid[3], *fp1, *fp2;
01695         fp1 = framepos + 3L*i; // position of atom 'i'
01696         a1 = mol->atom(i);
01697 
01698 #if !defined(VMDOPTIXRTRT)
01699         if (lastcolor != atomColor->color[i]) {
01700           lastcolor = atomColor->color[i];
01701           cmdColorIndex.putdata(lastcolor, cmdList);
01702         }
01703 #endif
01704 
01705         // draw half-bond to each bonded, displayed partner
01706         int n = cutoff > 0 ? nbonds[i] : a1->bonds;
01707         for (j=0; j < n; j++) {
01708           a2n = cutoff > 0 ? bondlists[MAXATOMBONDS*i + j] : a1->bondTo[j];
01709           if (atomSel->on[a2n]) {      // bonded atom 'a2n' displayed?
01710             fp2 = framepos + 3L*a2n;   // position of atom 'a2n'
01711             // find the bond midpoint 'mid' between atoms 'i' and 'a2n'
01712             mid[0] = 0.5f * (fp1[0] + fp2[0]);
01713             mid[1] = 0.5f * (fp1[1] + fp2[1]);
01714             mid[2] = 0.5f * (fp1[2] + fp2[2]);
01715 
01716 #if defined(VMDOPTIXRTRT)
01717             cylpoints.append3(fp1);
01718             cylpoints.append3(mid);
01719             cylradii.append(brad);
01720             cylcolors.append3(scene->color_value(atomColor->color[i]));
01721 #else
01722             // XXX here we should be aggregating into a cyl array primitive
01723             cmdCylinder.putdata(fp1, mid, brad, bres, 0, cmdList);
01724 #endif
01725           }
01726         }
01727       }
01728     }
01729 
01730 #if defined(VMDOPTIXRTRT)
01731     DispCmdCylinderArray cmdCylinderArray;
01732     cmdCylinderArray.putdata(&cylpoints[0], &cylradii[0], &cylcolors[0], 
01733                              cylradii.num(), bres, 0, cmdList);
01734 #endif
01735   }
01736 
01737   if (cutoff > 0) {
01738     delete [] nbonds;
01739     delete [] bondlists;
01740   }
01741 }
01742 
01743 
01744 
01745 // This draws a cylindrical bond from atom 'i' to atom 'k', but
01746 // only if both are selected.  Each end of the cylinder is colored
01747 // appropriately.  If there are bonds connected to the ends of this
01748 // i-k bond, then the i-k bond is extended so that the bonds have
01749 // no gap between them
01750 void DrawMolItem::draw_bonds(float *framepos, float brad, int bres, float cutoff) {
01751   MolAtom *a1, *a2;
01752   int g=0, h=0, i=0, j=0, k=0, l=0, m=0;
01753   int lastcolor = -1;
01754   int use_cyl;
01755   const float *bondorders = mol->bondorders();
01756 
01757   sprintf(commentBuffer,"Mol[%d] Rep[%d] Bonds", mol->id(), repNumber);
01758   cmdCommentX.putdata(commentBuffer, cmdList);
01759 
01760   if (bres <= 2 || brad < 0.01 ) {
01761     // draw bonds with lines
01762     use_cyl = FALSE;
01763     append(DMATERIALOFF);
01764     cmdLineType.putdata(SOLIDLINE, cmdList);
01765     cmdLineWidth.putdata(2, cmdList);
01766   } else {
01767     use_cyl = TRUE;
01768     // set up general drawing characteristics for this drawing method
01769     // turn on material characteristics
01770     append(DMATERIALON);
01771   }
01772 
01773   for (i=atomSel->firstsel; i <= atomSel->lastsel; i++) {
01774     // Only bonds to 'on' atoms are considered.
01775     if (atomSel->on[i]) {
01776       float *p2 = framepos + 3L*i;   // set p2 to atom 'i'
01777       a1 = mol->atom(i);            // find a selected atom
01778       float idouble[3], kdouble[3]; 
01779       float itriple[3], ktriple[3]; 
01780       const float *bondorderi = bondorders + (i * MAXATOMBONDS);
01781 
01782       for (j=0; j<a1->bonds; j++) {
01783         k = a1->bondTo[j];          // find a bonded atom that is turned on
01784         if (k > i && atomSel->on[k]) {
01785           float *p1, p3[3], *p4, *p5;
01786           a2 = mol->atom(k);
01787           p4 = framepos + 3L*k;      // set p4 to atom 'k'
01788 
01789           // find the bond midpoint 'p3' between atoms 'i' and 'k' 
01790           p3[0] = 0.5f * (p2[0] + p4[0]);
01791           p3[1] = 0.5f * (p2[1] + p4[1]);
01792           p3[2] = 0.5f * (p2[2] + p4[2]);
01793 
01794           // We extend the ends of the bond to touch neighboring bonds
01795           // Look for atom 'm' bonded to 'k' that isn't 'i', to extend to
01796           p5 = NULL;                // only used if not NULL
01797           for (l=a2->bonds-1; l>=0; l--) {
01798             m = a2->bondTo[l];
01799             if (m != i && atomSel->on[m]) {
01800               p5 = framepos + 3L*m; // set p5 to atom 'm'
01801               break;
01802             }
01803           }
01804   
01805           // Look for atom 'g' bonded to 'i' that isn't 'k', to extend to
01806           p1 = NULL;                // only used if not NULL
01807           for (h=a1->bonds-1; h>=0; h--) {
01808             g = a1->bondTo[h];
01809             if (g != k && atomSel->on[g]) {
01810               p1 = framepos + 3L*g; // set p1 to atom 'g'
01811               break;
01812             }
01813           }
01814 
01815           // determine the bond order (-1 (unset), 1, 2, or 3)
01816           float order = 0;
01817           if (bondorders != NULL) {
01818             order = bondorderi[j];
01819             if (order > 1) {
01820               int lv;
01821               for (lv=0; lv<3; lv++) {
01822                 idouble[lv] = p2[lv] + 0.333f * (p4[lv] - p2[lv]);
01823                 kdouble[lv] = p2[lv] + 0.666f * (p4[lv] - p2[lv]);
01824                 itriple[lv] = p2[lv] + 0.450f * (p4[lv] - p2[lv]);
01825                 ktriple[lv] = p2[lv] + 0.550f * (p4[lv] - p2[lv]);
01826               } 
01827             }
01828           }
01829 
01830           // now I can draw the bonds, setting the color for each
01831           if (lastcolor != atomColor->color[i]) {
01832             lastcolor = atomColor->color[i];
01833             cmdColorIndex.putdata(lastcolor, cmdList);
01834           }
01835           make_connection(p1, p2, p3, p4, brad, bres, use_cyl);
01836           if (order > 1) {
01837             make_connection(NULL, idouble, p3, kdouble, brad * 1.5f, bres, use_cyl);
01838             if (order > 2)
01839               make_connection(NULL, itriple, p3, ktriple, brad * 2.0f, bres, use_cyl);
01840           }
01841 
01842           if (lastcolor != atomColor->color[k]) {
01843             lastcolor = atomColor->color[k];
01844             cmdColorIndex.putdata(lastcolor, cmdList);
01845           }
01846           make_connection(p2, p3, p4, p5, brad, bres, use_cyl);
01847           if (order > 1) {
01848             make_connection(idouble, p3, kdouble, NULL, brad * 1.5f, bres, use_cyl);
01849             if (order > 2)
01850               make_connection(itriple, p3, ktriple, NULL, brad * 2.0f, bres, use_cyl);
01851           }
01852         } // found k atom
01853       } // searching along i's bonds
01854     } // found i
01855   } // searching each atom
01856 }
01857 
01858 
01859 // given the current cylinder and the previous and next cylinders,
01860 // find and draw the cylinder (or line) that best fits without making
01861 // an overlap.  prev can be NULL, in which case it acts as if the
01862 // previous cylinder was linear with the one to draw.  Ditto (but
01863 // opposite) for next.
01864 //
01865 // The problem is that if two cylinders touch head to tail,
01866 // but at an angle, there is a gap between them.  This routine
01867 // eliminates the gap by extending the cylinder endpoints.
01868 void DrawMolItem::make_connection(float *prev, float *start, float *end,
01869                                   float *next, float rad, int res, int use_cyl)
01870 {
01871   if (!start || !end) {
01872     msgErr << "Trying to make an extended cylinder with NULL end point(s)"
01873            << sendmsg;
01874     return;
01875   }
01876 
01877   if (!use_cyl) {
01878     cmdLine.putdata(start, end, cmdList);
01879     return;
01880   }
01881 
01882   float new_start[3], new_end[3];
01883   float cthis[3];
01884   float cnext[3];
01885   float cprev[3];
01886   float length;
01887   float costheta;
01888 
01889   vec_sub(cthis, end, start);    // cthis = end - start, cylinder (b)
01890   vec_normalize(cthis);
01891 
01892   length = 0.0f;                 // initialize to zero length by default
01893   if (prev != NULL) {            // previous cylinder (a)
01894     vec_sub(cprev, start, prev); // cprev = start - prev
01895     vec_normalize(cprev);
01896 
01897     // Compute length for aligning this cylinder with previous
01898     costheta = dot_prod(cprev, cthis);
01899 
01900     // does this exceed threshhold to start extension (might speed things up)
01901     if ((costheta > 0.0f) && (costheta < 0.9999f)) 
01902       length = rad * sqrtf((1.0f-costheta) / (1.0f+costheta));
01903   }
01904   new_start[0] = start[0] - length * cthis[0];
01905   new_start[1] = start[1] - length * cthis[1];  
01906   new_start[2] = start[2] - length * cthis[2];
01907 
01908   length = 0.0f;                 // initialize to zero length by default
01909   if (next != NULL) {            // along the next cylinder (a)
01910     vec_sub(cnext, next, end);   // cnext = next - end
01911     vec_normalize(cnext);
01912 
01913     // Now move the ending point to line up with the next cylinder.
01914     costheta = dot_prod(cthis,cnext);
01915 
01916     // does this exceed threshhold to start extension (might speed things up)
01917     if ((costheta > 0.0f) && (costheta < 0.9999f)) 
01918       length = rad * sqrtf((1.0f-costheta) / (1.0f+costheta));
01919   }
01920   new_end[0] = end[0] + length * cthis[0];
01921   new_end[1] = end[1] + length * cthis[1];
01922   new_end[2] = end[2] + length * cthis[2];
01923 
01924   // XXX here we should be aggregating into a cyl array primitive
01925   cmdCylinder.putdata(new_start, new_end, rad, res, 0, cmdList);
01926 }
01927 
01928 
01929 // Make a spline curve from a set of coordinates.  There are 'num' coords.
01930 // Coords goes from -2 to num+1 and has 3 floats per coord
01931 //    idx goes from -2 to num+1 and has either the atom number or -1
01932 void DrawMolItem::draw_spline_curve(int num, float *coords, int *idx,
01933                                     int use_cyl, float b_rad, int b_res) {
01934   float q[4][3];       // spline Q matrix
01935   float prev[2][3];    // the last two points for the previous spline
01936   float final[7][3];   // the points of the current approximation
01937   int last_loop = -10; // what is this value for?
01938   int loop;
01939 
01940   ResizeArray<float> pickpointcoords;
01941   ResizeArray<int> pickpointindices;
01942 
01943 #if 0
01944   const char *modulatefield = NULL; // XXX this needs to become a parameter
01945   const float *modulatedata = NULL; // data field to use for width modulation
01946   float *modulate = NULL;           // per-control point width values
01947 
01948   // XXX hack to let users try various stuff
01949   modulatefield = getenv("VMDMODULATERIBBON");
01950   if (modulatefield != NULL) {
01951     if (!strcmp(modulatefield, "user")) {
01952       modulatedata = mol->current()->user;
01953       // XXX user field can be NULL on some timesteps
01954     } else {
01955       modulatedata = mol->extraflt.data(modulatefield);
01956     } 
01957 
01958     // allocate for the maximum possible per-residue modulation values
01959     // so we don't have to reallocate for every fragment
01960     // XXX the modulate array is allocated cleared to zeros so that
01961     // in the case we get a NULL modulatedata pointer (user field for
01962     // example, we'll just end up with no modulation.
01963     modulate = (float *) calloc(1, mol->nResidues * sizeof(float));
01964   }
01965 #endif
01966 
01967   for (loop=-1; loop<num; loop++) { // go through the array
01968     // check if we need to do any computations, makes code faster
01969     // but its a bit more of a headache later on
01970     if ((idx[loop  ] >=0 && atomSel->on[idx[loop  ]]) ||
01971         (idx[loop+1] >=0 && atomSel->on[idx[loop+1]])) {
01972       make_spline_Q_matrix(q, spline_basis, coords+(loop-1)*3);
01973          
01974       // evaluate the interpolation between atom "loop" and
01975       // "loop+1".  I'll make this in 6 pieces, 7 coords.
01976       make_spline_interpolation(final[0], 0.0f/6.0f, q);
01977       make_spline_interpolation(final[1], 1.0f/6.0f, q);
01978       make_spline_interpolation(final[2], 2.0f/6.0f, q);
01979       make_spline_interpolation(final[3], 3.0f/6.0f, q);
01980       make_spline_interpolation(final[4], 4.0f/6.0f, q);
01981       make_spline_interpolation(final[5], 5.0f/6.0f, q);
01982       make_spline_interpolation(final[6], 6.0f/6.0f, q);
01983          
01984       // the possible values of 'on' are 0 (off), 1 (on),
01985       // 2 (draw 1st half), and 3 (draw 2nd half)
01986       // However, if the 1st half is drawn, the prev atom must
01987       // be on, and if the 2nd half is drawn, the new one must
01988       // be on.
01989 
01990       // draw what I need for atom 'loop'
01991       if (idx[loop] >= 0 &&          // this is a real atom
01992           (atomSel->on[idx[loop  ]] == 1 ||    // and it is turned on
01993           (atomSel->on[idx[loop  ]] == 3 && idx[loop+1] >= 0 &&
01994            atomSel->on[idx[loop+1]]))) {
01995     
01996         // here is the trickiness I was talking about
01997         // I need to see if there was a computation for
01998         // the first part of this atom
01999         if (last_loop != loop - 1) {  // no there wasn't
02000            cmdColorIndex.putdata(atomColor->color[idx[loop]], cmdList);
02001            make_connection(NULL, final[0], final[1], final[2],
02002                            b_rad, b_res, use_cyl);
02003         } else {
02004           // finish up drawing the previous part of this residue (I couldn't
02005           // do this before since I didn't know how to angle the final cylinder
02006           // to get a nice match up to the first of the current cylinders)
02007           make_connection(prev[0], prev[1], final[0], final[1],
02008                           b_rad, b_res, use_cyl);
02009           make_connection(prev[1], final[0], final[1], final[2],
02010                           b_rad, b_res, use_cyl);
02011         }
02012         make_connection(final[0], final[1], final[2], final[3],
02013                         b_rad, b_res, use_cyl);
02014         make_connection(final[1], final[2], final[3], final[4],
02015                         b_rad, b_res, use_cyl);
02016 
02017         // indicate this atom can be picked
02018         int pidx = 3L * loop;
02019         pickpointcoords.append3(&coords[pidx]);
02020         pickpointindices.append(idx[loop]);
02021       }
02022          
02023       // draw what I can for atom 'loop+1'
02024       if (idx[loop+1] >= 0 &&
02025           (atomSel->on[idx[loop+1]] == 1 ||
02026           (atomSel->on[idx[loop+1]] == 2 && idx[loop] >= 0 &&
02027            atomSel->on[idx[loop]]))) {
02028         cmdColorIndex.putdata(atomColor->color[idx[loop+1]], cmdList);
02029         make_connection(final[2], final[3], final[4], final[5],
02030                         b_rad, b_res, use_cyl);
02031         make_connection(final[3], final[4], final[5], final[6],
02032                         b_rad, b_res, use_cyl);
02033         last_loop = loop;
02034       }
02035       vec_copy(prev[0], final[4]);  // save for the
02036       vec_copy(prev[1], final[5]);  // next interation
02037     } 
02038   } // gone down the fragment
02039 
02040   // draw the pickpoints if we have any
02041   if (pickpointindices.num() > 0) {
02042     DispCmdPickPointArray pickPointArray;
02043     pickPointArray.putdata(pickpointindices.num(), &pickpointindices[0], 
02044                            &pickpointcoords[0], cmdList);
02045   }
02046 }
02047 
02048 
02049 void DrawMolItem::generate_tubearray() {
02050   if (tubearray) {
02051     for (int i=0; i<tubearray->num(); i++) delete (*tubearray)[i];
02052     delete tubearray;
02053   }
02054   tubearray = new ResizeArray<TubeIndexList *>;
02055 
02056   int frag;
02057   // start with protein
02058   for (frag = 0; frag < mol->pfragList.num(); frag++) {
02059     int num = mol->pfragList[frag]->num(); // number of residues
02060     if (num < 2) continue;
02061     TubeIndexList *idx = new TubeIndexList;
02062     // I'll go ahead and pad the beginning and end with two -1's for now.
02063     // It shouldn't be necessary but that's how draw_spline_curve works.
02064     idx->append2(-1, -1);
02065     for (int loop = 0; loop < num; loop++) {
02066       int res = (*mol->pfragList[frag])[loop];
02067       // Find the CA in this residue.  If not found, use the first atom
02068       // in the residue, since there has to be at least one atom.
02069       int atomnum = mol->find_atom_in_residue("CA", res);
02070       if (atomnum < 0) {
02071         atomnum = mol->atom_residue(res)->atoms[0];
02072       }
02073       idx->append(atomnum);
02074     }
02075     idx->append2(-1, -1);
02076     tubearray->append(idx);
02077   }
02078 
02079   // If there were no pfrags, check for a molecule with only CA's, 
02080   // and draw tubes connecting CA's with consecutive resid's.
02081   if (mol->pfragList.num() == 0) {
02082     int num = mol->nAtoms;
02083     int ca_num = mol->atomNames.typecode("CA"); // for fast lookup
02084     int last_resid = -10000;
02085     int resid;
02086     MolAtom *atm = NULL;
02087     
02088     TubeIndexList *idx = NULL;
02089     for (int i=0; i<=num; i++) {
02090       if (i != num) {
02091         atm = mol->atom(i);
02092         resid = atm->resid;
02093       } else {
02094         resid = -1000; // want it be out of range, but not equal to -10000.
02095       }
02096       // find a sequential CA
02097       if (atm->nameindex == ca_num && resid == last_resid + 1) {
02098         if (idx == NULL) { 
02099           msgErr << "Internal error in draw_tube for CA atoms: last_resid = "
02100             << last_resid << " but no first atom was added." << sendmsg;
02101           msgErr << "Tubes may be incomplete." << sendmsg;
02102           return;
02103         }
02104         idx->append(i);
02105       } else {
02106         if (idx) { // were there any points?
02107           idx->append2(-1, -1);
02108           tubearray->append(idx);
02109           idx = NULL;
02110         }
02111         if (atm->nameindex == ca_num && i != num) {
02112           idx = new TubeIndexList;
02113           idx->append3(-1, -1, i);
02114         }
02115       }
02116       last_resid = resid;
02117     }
02118   }
02119 
02120   // Now for nucleic acid backbone
02121   // This is just like we did for protein, except that instead of looking
02122   // for CA we look for P, or C5', or C5*.
02123   for (frag = 0; frag < mol->nfragList.num(); frag++) {
02124     int num = mol->nfragList[frag]->num(); // number of residues
02125     if (num < 2) continue;
02126     TubeIndexList *idx = new TubeIndexList;
02127     // I'll go ahead and pad the beginning and end with two -1's for now.
02128     // It shouldn't be necessary.
02129     idx->append2(-1, -1);
02130     for (int loop = 0; loop < num; loop++) {
02131       int res = (*mol->nfragList[frag])[loop];
02132       // Find the P, C5', or C5* in this residue.  
02133       // If not found, use the first atom
02134       // in the residue, since there has to be at least one atom.
02135       int atomnum = mol->find_atom_in_residue("P", res);
02136       if (atomnum < 0) {
02137         atomnum = mol->find_atom_in_residue("C5'", res);
02138         if (atomnum < 0) {
02139           atomnum = mol->find_atom_in_residue("C5*", res);
02140           if (atomnum < 0) { 
02141             atomnum = mol->atom_residue(res)->atoms[0];
02142           }
02143         }
02144       }
02145       idx->append(atomnum);
02146     }
02147     idx->append2(-1, -1);
02148     tubearray->append(idx);
02149   }
02150 
02151   // If there were no nfrags, check for a molecule with only P's, 
02152   // and draw tubes connecting P's with consecutive resid's.
02153   // This is duplicated from what we did for CA's - the only line that's
02154   // changed is the definition of ca_num.
02155 
02156   if (mol->nfragList.num() == 0) {
02157     int num = mol->nAtoms;
02158     int ca_num = mol->atomNames.typecode("P"); // for fast lookup
02159     int last_resid = -10000;
02160     int resid;
02161     MolAtom *atm = NULL;
02162     
02163     TubeIndexList *idx = NULL;
02164     for (int i=0; i<=num; i++) {
02165       if (i != num) {
02166         atm = mol->atom(i);
02167         resid = atm->resid;
02168       } else {
02169         resid = -1000; // want it be out of range, but not equal to -10000.
02170       }
02171       // find a sequential P
02172       if (atm->nameindex == ca_num && resid == last_resid + 1) {
02173         if (idx == NULL) { 
02174           msgErr << "Internal error in draw_tube for P atoms: last_resid = "
02175             << last_resid << " but no first atom was added." << sendmsg;
02176           msgErr << "Tubes may be incomplete." << sendmsg;
02177           return;
02178         }
02179         idx->append(i);
02180       } else {
02181         if (idx) { // were there any points?
02182           idx->append2(-1, -1);
02183           tubearray->append(idx);
02184           idx = NULL;
02185         }
02186         if (atm->nameindex == ca_num && i != num) {
02187           idx = new TubeIndexList;
02188           idx->append3(-1, -1, i);
02189         }
02190       }
02191       last_resid = resid;
02192     }
02193   }
02194 }
02195 
02196 
02197 // draw a cubic spline (in this case, using a modified Catmull-Rom basis)
02198 // through the C alpha of the protein and the P of the nucleic acids
02199 void DrawMolItem::draw_tube(float *framepos, float b_rad, int b_res) {
02200   // If we don't have tubearray yet, we need to generate it.  We don't
02201   // have to worry about selection or color because that's handled by
02202   // draw_spline_curve().
02203   if (!tubearray) {
02204     generate_tubearray();
02205   }
02206 
02207   sprintf(commentBuffer,"Mol[%d] Rep[%d] Tube", mol->id(), repNumber);
02208   cmdCommentX.putdata(commentBuffer, cmdList);
02209 
02210   // find out if I'm using lines or cylinders
02211   int use_cyl = FALSE;
02212   if (b_res <= 2 || b_rad < 0.01) { // then going to do lines
02213     append(DMATERIALOFF);
02214     cmdLineType.putdata(SOLIDLINE, cmdList);
02215     cmdLineWidth.putdata(2, cmdList);
02216   } else {
02217     use_cyl = TRUE;
02218     append(DMATERIALON);
02219   }
02220 
02221   for (int i=0; i<tubearray->num(); i++) {
02222     // XXX copy the current coordinates into a temporary array.
02223     TubeIndexList &indxlist = *(*tubearray)[i];
02224     int num = indxlist.num();
02225     float *coords = new float[3L*num];
02226     // First two and last two points are always just the third and third to
02227     // last real control points.
02228     int firstind = indxlist[2];
02229     int lastind  = indxlist[num-3];
02230     memcpy(coords,   framepos+3L*firstind, 3L*sizeof(float));
02231     memcpy(coords+3, framepos+3L*firstind, 3L*sizeof(float));
02232     for (int i=2; i<num-2; i++) {
02233       memcpy(coords+3L*i, framepos+3L*indxlist[i], 3L*sizeof(float));
02234     }
02235     memcpy(coords+3L*(num-2), framepos+3L*lastind, 3L*sizeof(float));
02236     memcpy(coords+3L*(num-1), framepos+3L*lastind, 3L*sizeof(float));
02237     draw_spline_curve(num-4, coords+6, &(indxlist[0])+2, use_cyl, b_rad, b_res);
02238     delete [] coords;
02239   }
02240 }
02241 
02242 
02243 //  This draws a hydrogen bond from atom 'i' to the hydrogen of atom 'k', only
02244 //  if both are selected and meet the hbond distance and/or angle criteria
02245 //  distance is measured from heavy atom to heavy atom
02246 //  the angle is from heavy atom to hydrogen to heavy atom
02247 //  coloring is by the donor atom color specification
02248 void DrawMolItem::draw_hbonds(float *framepos, float maxangle, int thickness, float maxdist) {
02249   int i, k;
02250   float  donortoH[3],Htoacceptor[3];
02251   float cosmaxangle2 = cosf(maxangle); cosmaxangle2 *= cosmaxangle2; 
02252   ResizeArray<float> pickpointcoords;
02253   ResizeArray<int> pickpointindices;
02254 
02255   // protect against ridiculous values which might crash us or
02256   // run the system out of memory.
02257   if (maxdist <= 0.1 || atomSel->selected == 0) {
02258     return;
02259   }
02260    
02261   int *onlist = (int *) calloc(1, mol->nAtoms * sizeof(int)); // clear to zeros
02262   for (i=atomSel->firstsel; i <= atomSel->lastsel; i++) {
02263     if (atomSel->on[i] && mol->atom(i)->atomType != ATOMHYDROGEN)
02264       onlist[i] = 1;
02265   } 
02266 
02267   sprintf(commentBuffer,"Mol[%d] Rep[%d] HBonds", mol->id(), repNumber);
02268   cmdCommentX.putdata(commentBuffer, cmdList);
02269    
02270   // set up general drawing characteristics for this drawing method
02271   append(DMATERIALOFF);
02272   cmdLineType.putdata(DASHEDLINE, cmdList);
02273   cmdLineWidth.putdata(thickness, cmdList);
02274 
02275   // XXX the actual code for measuring hbonds doesn't belong here, it should
02276   //     be moved into Measure.[Ch] where it really belongs.  This file
02277   //     only implements the rendering interface, and should not be doing the
02278   //     hard core math, particularly since we want to expose the same
02279   //     feature via scripting interfaces etc.  Also, having a single
02280   //     implementation avoids having different bugs in the
02281   //     long-term.  Too late to do anything about this now, but should be
02282   //     addressed for the next major version when time allows.
02283 
02284   // XXX This code is similar, but no identical to the code in TclMeasure.C
02285   //     It would be better written if rather than emitting lines and 
02286   //     pickpoints piecemeal, if it instead got a complete list of HBonds
02287   //     and then built a single vertex array instead. 
02288  
02289   // loop over all SELECTED atoms, inspite of the fact that non-hbonding
02290   // atoms may actually be selected, searching for 
02291   // any pair that are closer than distance, not bonded,
02292   // and have a hydrogen that forms an acceptable angle.. all said 
02293   // and done it works pretty well only catching C-H---O hbonds 
02294   // only when the distance and angle are set rather unrealistic
02295   GridSearchPair *pairlist = vmd_gridsearch1(framepos, mol->nAtoms, onlist, maxdist, 0, mol->nAtoms * 27);
02296   GridSearchPair *p, *tmp;
02297   for (p=pairlist; p != NULL; p=tmp) {
02298     MolAtom *a1 = mol->atom(p->ind1);
02299     MolAtom *a2 = mol->atom(p->ind2);
02300 
02301     // ignore if bonded
02302     if (!a2->bonded(p->ind1)) {
02303       int b1 = a1->bonds;
02304       int b2 = a2->bonds;
02305       float *coor1 = framepos + 3L*p->ind1; 
02306       float *coor2 = framepos + 3L*p->ind2; 
02307   
02308       for (k=0; k < b2; k++) {
02309         if (mol->atom(a2->bondTo[k])->atomType == ATOMHYDROGEN) {
02310           float *hydrogen = framepos + 3L*a2->bondTo[k];
02311           vec_sub(donortoH,hydrogen,coor2);
02312           vec_sub(Htoacceptor,coor1,hydrogen);
02313           if (angle(donortoH, Htoacceptor)  < maxangle ) {
02314             cmdColorIndex.putdata(atomColor->color[p->ind2], cmdList);
02315             cmdLine.putdata(coor1,hydrogen, cmdList); // draw line
02316 
02317             // indicate the bonded atoms can be picked
02318             int pidx = 3L * a2->bondTo[k];
02319             pickpointcoords.append3(&framepos[pidx]);
02320 
02321             pidx = 3L * p->ind1;
02322             pickpointcoords.append3(&framepos[pidx]);
02323             pickpointindices.append2(a2->bondTo[k], p->ind1);
02324           }
02325         }
02326       }
02327       for (k=0; k < b1; k++){
02328         if (mol->atom(a1->bondTo[k])->atomType == ATOMHYDROGEN) {
02329           float *hydrogen = framepos + 3L*a1->bondTo[k];
02330           vec_sub(donortoH,hydrogen,coor1);
02331           vec_sub(Htoacceptor,coor2,hydrogen);
02332           if (angle(donortoH, Htoacceptor)  < maxangle ) {
02333             cmdColorIndex.putdata(atomColor->color[p->ind1], cmdList);
02334             cmdLine.putdata(hydrogen,coor2, cmdList); // draw line
02335 
02336             // indicate the bonded atoms can be picked
02337             int pidx = 3L * a1->bondTo[k];
02338             pickpointcoords.append3(&framepos[pidx]);
02339 
02340             pidx = 3L * p->ind2;
02341             pickpointcoords.append3(&framepos[pidx]);
02342             pickpointindices.append2(a1->bondTo[k], p->ind2);
02343           }
02344         }
02345       }
02346     }
02347     tmp = p->next;
02348     free(p); 
02349   }
02350   free(onlist);
02351 
02352   // Revert line state to solid lines just in case other code has failed to 
02353   // ensure the correct line stipple state is in effect prior to drawing
02354   cmdLineType.putdata(SOLIDLINE, cmdList);
02355 
02356   // draw the pickpoints if we have any
02357   if (pickpointindices.num() > 0) {
02358     DispCmdPickPointArray pickPointArray;
02359     pickPointArray.putdata(pickpointindices.num(), &pickpointindices[0],
02360                            &pickpointcoords[0], cmdList);
02361   }
02362 }
02363 
02364 
02365 void DrawMolItem::draw_dynamic_bonds(float *framepos, float brad, int bres, float maxdist) {
02366   int lastcolor = -1;
02367   ResizeArray<float> pickpointcoords;
02368   ResizeArray<int> pickpointindices;
02369 
02370   // protect against ridiculous values which might crash us or
02371   // run the system out of memory.
02372   if (maxdist <= 0.1 || atomSel->selected == 0) {
02373     return;
02374   }
02375    
02376   sprintf(commentBuffer,"Mol[%d] Rep[%d] DynamicBonds", mol->id(), repNumber);
02377   cmdCommentX.putdata(commentBuffer, cmdList);
02378    
02379   // set up general drawing characteristics for this drawing method
02380   append(DMATERIALON);
02381 #if defined(VMDOPTIXRTRT)
02382   ResizeArray <float> cylpoints;
02383   ResizeArray <float> cylradii;
02384   ResizeArray <float> cylcolors;
02385 #endif
02386  
02387   // loop over all SELECTED atoms, in spite of the fact that non-bonded
02388   // atoms may actually be selected, searching for 
02389   // any pair that are closer than distance, not bonded,
02390   GridSearchPair *pairlist = vmd_gridsearch1(framepos, mol->nAtoms, atomSel->on, maxdist, 0, mol->nAtoms * 27);
02391   GridSearchPair *p, *tmp;
02392   for (p=pairlist; p != NULL; p=tmp) {
02393     MolAtom *atom1 = mol->atom(p->ind1);
02394     MolAtom *atom2 = mol->atom(p->ind2);
02395 
02396     // don't bond atoms that aren't part of the same conformation    
02397     // or that aren't in the all-conformations part of the structure 
02398     if ((atom1->altlocindex != atom2->altlocindex) &&
02399         ((mol->altlocNames.name(atom1->altlocindex)[0] != '\0') &&
02400          (mol->altlocNames.name(atom2->altlocindex)[0] != '\0'))) {
02401       tmp = p->next;
02402       free(p);
02403       continue;
02404     }
02405 
02406     // Prevent hydrogens from bonding with each other.
02407     // Use atomType info derived during initial molecule analysis for speed.
02408     if (!(atom1->atomType == ATOMHYDROGEN) ||
02409         !(atom2->atomType == ATOMHYDROGEN)) {
02410       float *coor1 = framepos + 3L*p->ind1; 
02411       float *coor2 = framepos + 3L*p->ind2; 
02412       float mid[3];
02413 #if 0
02414       if (cutoff < 0) { // Do atom-specific distance check
02415         float d2 = distance2(coor1, coor2);
02416         float r1 = atom1->extra[ATOMRAD];
02417         float r2 = atom2->extra[ATOMRAD];
02418         float cut = 0.6f * (r1 + r2);
02419         if (d2 < cut*cut)
02420           // mol->add_bond(p->ind1, p->ind2);
02421       } else
02422 #endif
02423       // mol->add_bond(p->ind1, p->ind2);
02424 
02425       // draw half-bond to each bonded, displayed partner
02426       // find the bond midpoint 'mid' between atoms 'i' and 'a2n'
02427       mid[0] = 0.5f * (coor1[0] + coor2[0]);
02428       mid[1] = 0.5f * (coor1[1] + coor2[1]);
02429       mid[2] = 0.5f * (coor1[2] + coor2[2]);
02430 
02431 #if defined(VMDOPTIXRTRT)
02432       cylpoints.append3(coor1);
02433       cylpoints.append3(mid);
02434       cylradii.append(brad);
02435       cylcolors.append3(scene->color_value(atomColor->color[p->ind1]));
02436       cylpoints.append3(mid);
02437       cylpoints.append3(coor2);
02438       cylradii.append(brad);
02439       cylcolors.append3(scene->color_value(atomColor->color[p->ind2]));
02440 #else
02441       if (lastcolor != atomColor->color[p->ind1]) {
02442         lastcolor = atomColor->color[p->ind1];
02443         cmdColorIndex.putdata(lastcolor, cmdList);
02444       }
02445       // XXX here we should be aggregating into a cyl array primitive
02446       cmdCylinder.putdata(coor1, mid, brad, bres, 0, cmdList);
02447 
02448       if (lastcolor != atomColor->color[p->ind2]) {
02449         lastcolor = atomColor->color[p->ind2];
02450         cmdColorIndex.putdata(lastcolor, cmdList);
02451       }
02452       // XXX here we should be aggregating into a cyl array primitive
02453       cmdCylinder.putdata(mid, coor2, brad, bres, 0, cmdList);
02454 #endif
02455 
02456       // indicate the bonded atoms can be picked
02457       int pidx = 3L * p->ind1;
02458       pickpointcoords.append3(&framepos[pidx]);
02459 
02460       pidx = 3L * p->ind2;
02461       pickpointcoords.append3(&framepos[pidx]);
02462       pickpointindices.append2(p->ind1, p->ind2);
02463     }
02464   
02465     tmp = p->next;
02466     free(p); 
02467   }
02468 
02469 #if defined(VMDOPTIXRTRT)
02470   DispCmdCylinderArray cmdCylinderArray;
02471   cmdCylinderArray.putdata(&cylpoints[0], &cylradii[0], &cylcolors[0], 
02472                            cylradii.num(), bres, 0, cmdList);
02473 #endif
02474 
02475   // draw the pickpoints if we have any
02476   if (pickpointindices.num() > 0) {
02477     DispCmdPickPointArray pickPointArray;
02478     pickPointArray.putdata(pickpointindices.num(), &pickpointindices[0],
02479                            &pickpointcoords[0], cmdList);
02480   }
02481 }
02482 
02483 
02484 #ifdef VMDPOLYHEDRA
02485 // Based on original code by Dr. Francois-Xavier Coudert<f.coudert@ucl.ac.uk>
02486 #define PLYMAXNB 16
02487 void DrawMolItem::draw_polyhedra(float *framepos, float maxdist) {
02488   int i;
02489   int lastcolor = -1;
02490 
02491   // protect against ridiculous values which might crash us or
02492   // run the system out of memory.
02493   if (maxdist <= 0.1) {
02494     return;
02495   }
02496 
02497   // turn on all atoms for grid search (not using a real selection presently)
02498   int natoms = mol->nAtoms;
02499   int *onlist = new int[natoms];
02500   for (i=0; i<natoms; i++)
02501     onlist[i] = 1;
02502 
02503   // allocate array for selected neighbor list
02504   int *nblist = (int *) malloc(natoms * sizeof(int) * PLYMAXNB);
02505   if (nblist == NULL)
02506     return;
02507   memset(nblist, 0, natoms * sizeof(int) * PLYMAXNB);
02508 
02509   sprintf(commentBuffer,"Mol[%d] Rep[%d] Polyhedra", mol->id(), repNumber);
02510   cmdCommentX.putdata(commentBuffer, cmdList);
02511 
02512   // set up general drawing characteristics for this drawing method
02513   append(DMATERIALON);
02514 
02515   // loop over all SELECTED atoms, inspite of the fact that non-bonded
02516   // atoms may actually be selected, searching for 
02517   // any pair that are closer than distance, not bonded,
02518   GridSearchPair *pairlist = vmd_gridsearch1(framepos, natoms, onlist, maxdist, 0, natoms * (PLYMAXNB - 1));
02519   delete [] onlist;
02520 
02521   GridSearchPair *p, *tmp;
02522 
02523   for (p=pairlist; p != NULL; p=tmp) {
02524     MolAtom *atom1 = mol->atom(p->ind1);
02525     MolAtom *atom2 = mol->atom(p->ind2);
02526 
02527     // delete pairs where neither atom is selected
02528     if (!(atomSel->on[p->ind1] || atomSel->on[p->ind2])) {
02529       tmp = p->next;
02530       free(p);
02531       continue;
02532     }
02533 
02534     // don't bond atoms that aren't part of the same conformation    
02535     // or that aren't in the all-conformations part of the structure 
02536     if ((atom1->altlocindex != atom2->altlocindex) &&
02537         ((mol->altlocNames.name(atom1->altlocindex)[0] != '\0') &&
02538          (mol->altlocNames.name(atom2->altlocindex)[0] != '\0'))) {
02539       tmp = p->next;
02540       free(p);
02541       continue;
02542     }
02543 
02544     // record selected neighbors in the array
02545     // increment neighbor counts, clamping at PLYMAXNB-1 neighbors each
02546     int idx1 = p->ind1 * PLYMAXNB;
02547     int idx2 = p->ind2 * PLYMAXNB;
02548     if (nblist[idx1] < (PLYMAXNB-1) && nblist[idx2] < (PLYMAXNB-1)) {
02549       // update neighbor count 
02550       nblist[idx1]++; 
02551       nblist[idx2]++;
02552 
02553       // record new neighbors
02554       nblist[idx1 + nblist[idx1]] = p->ind2;
02555       nblist[idx2 + nblist[idx2]] = p->ind1;
02556     }
02557 
02558     tmp = p->next;
02559     free(p);
02560     continue;
02561   }
02562 
02563   // draw polyhedra composed of triangles for each combination
02564   // of three neighbors, using the neighbor atom coords as the vertices
02565   // XXX we should not be generating a massive list of independent 
02566   //     triangles as this is very wasteful of memory.  Instead, we
02567   //     should be building an indexed list of vertices, and generate
02568   //     one or more vertex array tokens for each of the colors that
02569   //     we'll ultimately draw.  This would save a massive amount of 
02570   //     display list memory that is wasted on needlessly 
02571   //     replicated vertex coordinates in the current implementation.
02572   for (i=0; i < natoms; i++) {
02573     int idx = i*PLYMAXNB;
02574 
02575     // don't draw a polyhedron if the atom isn't selected or
02576     // there are less than three neighbors within the user-selected 
02577     // cutoff distance
02578     if (!atomSel->on[i] || nblist[idx] < 4) 
02579       continue;
02580 
02581     // emit a color command only if necessary
02582     if (lastcolor != atomColor->color[i]) {
02583       lastcolor = atomColor->color[i];
02584       cmdColorIndex.putdata(lastcolor, cmdList);
02585     }
02586 
02587     int index=nblist[idx];
02588     for (int i1=1; i1 <= index - 2; i1++)
02589       for (int i2=i1; i2 <= index - 1; i2++)
02590         for (int i3=i2; i3 <= index; i3++)
02591           cmdTriangle.putdata(framepos + 3L*nblist[idx+i1],
02592                               framepos + 3L*nblist[idx+i2],
02593                               framepos + 3L*nblist[idx+i3], cmdList);
02594   }
02595 
02596   free(nblist);
02597 }
02598 #endif
02599 
02600 
02601 // find x = a*i + b where i = 0..n-1
02602 static void least_squares(int n, const float *x, float *a, float *b) {
02603   float sum = 0;
02604   int i;
02605   for (i=0; i<n; i++) {    // find the sum of x
02606     sum += x[i];
02607   }
02608   float d = (float(n)-1.0f) / 2.0f;
02609   float t, sum_t2 = 0.0f;
02610   *a = 0.0f;
02611   for (i=0; i<n; i++) {
02612     t = (i - d);
02613     sum_t2 += t*t;
02614     *a += t*x[i];
02615   }
02616   *a /= sum_t2;
02617   *b = (sum/float(n) - d*(*a));
02618 }
02619 
02620 // given the info about the helix (length, color, etc.), draw it
02621 // This computes a best fit line long the given points.  For each "on"
02622 // atom, a cylinder is drawn from 1/2 a residue before to 1/2 a residue
02623 // beyond, in the right color.
02624 void DrawMolItem::draw_alpha_helix_cylinders(ResizeArray<float> &x,
02625        ResizeArray<float> &y, ResizeArray<float> &z, ResizeArray<int> &atom_on,
02626        int *color, float bond_rad, int bond_res,
02627        float *start_coord, float *end_coord)
02628 {
02629   // compute the best fit line for each direction
02630   float a[3], b[3];
02631   ResizeArray<float> pickpointcoords;
02632   ResizeArray<int> pickpointindices;
02633 
02634   int num = x.num();
02635   least_squares(num, &(x[0]), a+0, b+0);
02636   least_squares(num, &(y[0]), a+1, b+1);
02637   least_squares(num, &(z[0]), a+2, b+2);
02638   
02639   // draw the cylinder(s)
02640   float start[3], end[3];
02641   
02642   start[0] = a[0] * (-0.5f) + b[0];
02643   start[1] = a[1] * (-0.5f) + b[1];
02644   start[2] = a[2] * (-0.5f) + b[2];
02645   vec_copy(start_coord, start);
02646   for (int i=0; i<x.num(); i++) {
02647     end[0] = a[0] * (float(i)+0.5f) + b[0];
02648     end[1] = a[1] * (float(i)+0.5f) + b[1];
02649     end[2] = a[2] * (float(i)+0.5f) + b[2];
02650     if (atom_on[i] >= 0) {
02651       // indicate this atom can be picked
02652       pickpointcoords.append3(x[i], y[i], z[i]);
02653       pickpointindices.append(atom_on[i]);
02654       cmdColorIndex.putdata(color[atom_on[i]], cmdList);
02655       int caps = 0;
02656       if (i == 0           || (i>0         && atom_on[i-1] < 0)) 
02657         caps |= CYLINDER_TRAILINGCAP;
02658       if (i == x.num() - 1 || (i<x.num()-1 && atom_on[i+1] < 0))
02659         caps |= CYLINDER_LEADINGCAP;
02660       if (bond_res <= 2 || bond_rad < 0.01) {
02661         // the representation will be as lines
02662         cmdLine.putdata(start, end, cmdList);
02663       } else {
02664         // draw the cylinder with closed ends
02665         cmdCylinder.putdata(start, end, bond_rad, bond_res, caps, cmdList);
02666       }
02667     }
02668     memcpy(start, end, 3L*sizeof(float));
02669   }
02670   vec_copy(end_coord, start);
02671 
02672   // draw the pickpoints if we have any
02673   if (pickpointindices.num() > 0) {
02674     DispCmdPickPointArray pickPointArray;
02675     pickPointArray.putdata(pickpointindices.num(), &pickpointindices[0],
02676                            &pickpointcoords[0], cmdList);
02677   }
02678 
02679   // reset for the next round
02680   x.clear();
02681   y.clear();
02682   z.clear();
02683   atom_on.clear();
02684 }
02685 
02686 #define SCALEADD(vec, offset, byscale) {   \
02687     vec[0] += (*(offset  ))*byscale;       \
02688     vec[1] += (*(offset+1))*byscale;       \
02689     vec[2] += (*(offset+2))*byscale;       \
02690 }
02691 #define SCALESUM(vec, term1, term2, byscale) {   \
02692     vec[0] = (*(term1  )) + (*(term2  ))*byscale;       \
02693     vec[1] = (*(term1+1)) + (*(term2+1))*byscale;       \
02694     vec[2] = (*(term1+2)) + (*(term2+2))*byscale;       \
02695 }
02696 #define BETASCALE (0.2f * ribbon_width)
02697 
02698 
02699 
02700 
02701 void DrawMolItem::draw_beta_sheet(ResizeArray<float> &x,
02702        ResizeArray<float> &y, ResizeArray<float> &z, ResizeArray<int> &atom_on,
02703        int *color, float ribbon_width, float *start_coord, float *end_coord)
02704 {
02705   ResizeArray<float> pickpointcoords;
02706   ResizeArray<int> pickpointindices;
02707 
02708   // compute the points between the atoms, and do a linear extrapolation
02709   // to get the first and last point.  (I should base it on something better
02710   // than a linear fit ... ).
02711   //   I am gauranteed at least three residues for initial data
02712   int num = x.num();
02713   int i;
02714   float *centers = new float[3L*(num+1)];
02715   for (i=1; i<num; i++) {   // compute the centers
02716     centers[3L*i+0] = (x[i-1]+x[i])/2;
02717     centers[3L*i+1] = (y[i-1]+y[i])/2;
02718     centers[3L*i+2] = (z[i-1]+z[i])/2;
02719   }
02720   // and the linear extrapolation
02721   centers[0] = 2*centers[3] - centers[6];
02722   centers[1] = 2*centers[4] - centers[7];
02723   centers[2] = 2*centers[5] - centers[8];
02724 
02725   centers[3L*num+0] = 2*centers[3L*num-3] - centers[3L*num-6];
02726   centers[3L*num+1] = 2*centers[3L*num-2] - centers[3L*num-5];
02727   centers[3L*num+2] = 2*centers[3L*num-1] - centers[3L*num-4];
02728 
02729   // now do the normals
02730   float *norms = new float[3L*(num+1)];  // along the width
02731   float *perps = new float[3L*(num+1)];  // along the height
02732   float d[2][3]; // deltas between successive points
02733   d[0][0] = x[1]-x[0];  d[0][1] = y[1]-y[0];  d[0][2] = z[1]-z[0];
02734   for (i=1; i<num-1; i++) {
02735     d[1][0] = x[i+1]-x[i];
02736     d[1][1] = y[i+1]-y[i];
02737     d[1][2] = z[i+1]-z[i];
02738     cross_prod(norms+3L*i, d[0], d[1]);
02739     vec_normalize(norms+3L*i);
02740     if (i%2) { // flip every other normal so they are aligned
02741       norms[3L*i+0] = -norms[3L*i+0];
02742       norms[3L*i+1] = -norms[3L*i+1];
02743       norms[3L*i+2] = -norms[3L*i+2];
02744     }
02745     vec_copy(d[0], d[1]);
02746   }
02747   // for the first one and last two normals, copy the end norms.
02748   vec_copy(norms, norms+3);
02749   vec_copy(norms+3L*num-3, norms+3L*num-6);
02750   vec_copy(norms+3L*num  , norms+3L*num-6);
02751 
02752   // and the perpendiculars
02753   for (i=0; i<num; i++) {
02754     vec_sub(d[0], centers+3L*i, centers+3L*i+3);
02755     cross_prod(perps+3L*i, d[0], norms+3L*i);
02756     vec_normalize(perps+3L*i);
02757   }
02758   vec_copy(perps+3L*num, perps+3L*num-3);
02759 
02760   // Draw rectangular blocks for each segment
02761   //  upper/lower, left/right corners;  this is the beginning of the block
02762   float ul[2][3], ur[2][3], ll[2][3], lr[2][3];
02763   SCALESUM(ul[0], centers, norms, ribbon_width);
02764   vec_copy    (ll[0], ul[0]);
02765   SCALEADD(ul[0], perps, BETASCALE);
02766   SCALEADD(ll[0], perps, -BETASCALE);
02767 
02768   SCALESUM(ur[0] , centers, norms, -ribbon_width);
02769   vec_copy    (lr[0], ur[0]);
02770   SCALEADD(ur[0], perps, BETASCALE);
02771   SCALEADD(lr[0], perps, -BETASCALE);
02772 
02773   // save the initial and final coordinate
02774   vec_copy(start_coord, centers);
02775   vec_copy(end_coord, centers+3L*num);
02776 
02777   int prev_on = 0;
02778   float dot=0.0f;
02779   for (i=0; i<num-1; i++) {  // go down the list of residues
02780     dot = dot_prod(perps+3L*i, perps+3L*i+3);
02781     // check the amount of rotation
02782     if (dot < 0) { // ribbons switched orientation
02783       // swap normals
02784       perps[3L*i+3] = -perps[3L*i+3];
02785       perps[3L*i+4] = -perps[3L*i+4];
02786       perps[3L*i+5] = -perps[3L*i+5];
02787       norms[3L*i+3] = -norms[3L*i+3];
02788       norms[3L*i+4] = -norms[3L*i+4];
02789       norms[3L*i+5] = -norms[3L*i+5];
02790     }
02791 
02792     // calculate the corners of the end of the blocks
02793     SCALESUM(ul[1], centers+3L*i+3, norms+3L*i+3, ribbon_width);
02794     vec_copy    (ll[1], ul[1]);
02795     SCALEADD(ul[1], perps+3L*i+3, BETASCALE);
02796     SCALEADD(ll[1], perps+3L*i+3, -BETASCALE);
02797     
02798     SCALESUM(ur[1] , centers+3L*i+3, norms+3L*i+3, -ribbon_width);
02799     vec_copy    (lr[1], ur[1]);
02800     SCALEADD(ur[1], perps+3L*i+3, BETASCALE);
02801     SCALEADD(lr[1], perps+3L*i+3, -BETASCALE);
02802 
02803 
02804     // draw this section, if it is on
02805     if (atom_on[i] >= 0) {
02806       // indicate this atom can be picked
02807       pickpointcoords.append3(x[i], y[i], z[i]);
02808       pickpointindices.append(atom_on[i]);
02809 
02810       cmdColorIndex.putdata(color[atom_on[i]], cmdList);
02811 
02812       if (ribbon_width != 0) {
02813         // draw 1 or 2 segments based on how much curvature there was
02814         if (i == 0 || dot > .9 || dot < -.9) {
02815           // small flip, so I just need one segment
02816           // draw the sides
02817           cmdTriangle.putdata(ur[0], ul[0], ul[1],  // top
02818                               perps+3L*i, perps+3L*i, perps+3L*i+3, cmdList);
02819           cmdTriangle.putdata(ur[0], ul[1], ur[1],
02820                               perps+3L*i, perps+3L*i+3, perps+3L*i+3, cmdList);
02821           cmdTriangle.putdata(lr[0], ll[0], ll[1],  // bottom
02822                               perps+3L*i, perps+3L*i, perps+3L*i+3, cmdList);
02823           cmdTriangle.putdata(lr[0], ll[1], lr[1],
02824                               perps+3L*i, perps+3L*i+3, perps+3L*i+3, cmdList);
02825           
02826           // and the other sides
02827           cmdTriangle.putdata(ul[0], ll[0], ll[1], 
02828                               norms+3L*i, norms+3L*i, norms+3L*i+3, cmdList);
02829           cmdTriangle.putdata(ul[0], ll[1], ul[1], 
02830                               norms+3L*i, norms+3L*i+3, norms+3L*i+3, cmdList);
02831           cmdTriangle.putdata(ur[0], lr[0], lr[1], 
02832                               norms+3L*i, norms+3L*i, norms+3L*i+3, cmdList);
02833           cmdTriangle.putdata(ur[0], lr[1], ur[1], 
02834                               norms+3L*i, norms+3L*i+3, norms+3L*i+3, cmdList);
02835         } else {
02836           // big turn, so make more lines
02837           // given the construction, I know the first and last residues
02838           // don't twist (the perps are just copies) so I just do a simple
02839           // spline along the centers
02840           float centers_q[4][3];
02841           make_spline_Q_matrix(centers_q, spline_basis, centers + 3L*i-3);
02842           
02843           // make the spline for the perps and norms
02844           float perps_q[4][3];
02845           make_spline_Q_matrix(perps_q, spline_basis, perps + 3L*i-3);
02846           float norms_q[4][3];
02847           make_spline_Q_matrix(norms_q, spline_basis, norms + 3L*i-3);
02848           
02849           // break the section into two parts, so compute the new middle
02850           float new_center[3];
02851           float new_perp[3];
02852           float new_norm[3];
02853           make_spline_interpolation(new_center, 0.5, centers_q);
02854           make_spline_interpolation(new_perp  , 0.5, perps_q);
02855           vec_normalize(new_perp);  // rescale correctly
02856           make_spline_interpolation(new_norm  , 0.5, norms_q);
02857           vec_normalize(new_norm);  // rescale correctly
02858           
02859           // and make the new corners
02860           float ul2[3], ur2[3], ll2[3], lr2[3];
02861           SCALESUM(ul2 , new_center, new_norm, ribbon_width);
02862           vec_copy    (ll2, ul2);
02863           SCALEADD(ul2, new_perp, BETASCALE);
02864           SCALEADD(ll2, new_perp, -BETASCALE);
02865           
02866           SCALESUM(ur2 , new_center, new_norm, -ribbon_width);
02867           vec_copy    (lr2, ur2);
02868           SCALEADD(ur2, new_perp, BETASCALE);
02869           SCALEADD(lr2, new_perp, -BETASCALE);
02870           
02871           // draw the new, intermediate blocks
02872           cmdTriangle.putdata(ur[0], ul[0], ul2,  // top
02873                               perps+3L*i, perps+3L*i, new_perp, cmdList);
02874           cmdTriangle.putdata(ur[0], ul2, ur2,
02875                               perps+3L*i, new_perp, new_perp, cmdList);
02876           cmdTriangle.putdata(lr[0], ll[0], ll2,  // bottom
02877                               perps+3L*i, perps+3L*i, new_perp, cmdList);
02878           cmdTriangle.putdata(lr[0], ll2, lr2,
02879                               perps+3L*i, new_perp, new_perp, cmdList);
02880 
02881         // and the other sides
02882         cmdTriangle.putdata(ul[0], ll[0], ll2, 
02883                             norms+3L*i, norms+3L*i, new_norm, cmdList);
02884         cmdTriangle.putdata(ul[0], ll2, ul2, 
02885                             norms+3L*i, new_norm, new_norm, cmdList);
02886         cmdTriangle.putdata(ur[0], lr[0], lr2, 
02887                             norms+3L*i, norms+3L*i, new_norm, cmdList);
02888         cmdTriangle.putdata(ur[0], lr2, ur2, 
02889                             norms+3L*i, new_norm, new_norm, cmdList);
02890         
02891         // get the 1/2 half
02892         cmdTriangle.putdata(ur2, ul2, ul[1],  // top
02893                             new_perp, new_perp, perps+3L*i+3, cmdList);
02894         cmdTriangle.putdata(ur2, ul[1], ur[1],
02895                             new_perp, perps+3L*i+3, perps+3L*i+3, cmdList);
02896         cmdTriangle.putdata(lr2, ll2, ll[1],  // bottom
02897                             new_perp, new_perp, perps+3L*i+3, cmdList);
02898         cmdTriangle.putdata(lr2, ll[1], lr[1],
02899                             new_perp, perps+3L*i+3, perps+3L*i+3, cmdList);
02900         
02901         // and the other sides
02902         cmdTriangle.putdata(ul2, ll2, ll[1], 
02903                             new_norm, new_norm, norms+3L*i+3, cmdList);
02904         cmdTriangle.putdata(ul2, ll[1], ul[1], 
02905                             new_norm, norms+3L*i+3, norms+3L*i+3, cmdList);
02906         cmdTriangle.putdata(ur2, lr2, lr[1], 
02907                             new_norm, new_norm, norms+3L*i+3, cmdList);
02908         cmdTriangle.putdata(ur2, lr[1], ur[1], 
02909                             new_norm, norms+3L*i+3, norms+3L*i+3, cmdList);
02910         }
02911         
02912         if (!prev_on) {
02913           // draw the base
02914           cmdTriangle.putdata(ul[0], ll[0], ur[0], cmdList);
02915           cmdTriangle.putdata(ll[0], lr[0], ur[0], cmdList);
02916           prev_on = 1;
02917         }
02918       } else {  // ribbon_width == 0
02919         cmdLine.putdata(centers+3L*i, centers+3L*i+3, cmdList);
02920       }
02921     } else { // atom_on[i] >= 0
02922       // this isn't on.  Was prev?  If so, draw its base
02923       if (prev_on && ribbon_width != 0) {
02924         cmdTriangle.putdata(ul[0], ll[0], ur[0], cmdList);
02925         cmdTriangle.putdata(ll[0], lr[0], ur[0], cmdList);
02926         prev_on = 0;
02927       }
02928     }
02929 
02930     vec_copy(ur[0], ur[1]);
02931     vec_copy(ul[0], ul[1]);
02932     vec_copy(lr[0], lr[1]);
02933     vec_copy(ll[0], ll[1]);
02934   } // end of loop though centers
02935   
02936   // but if the second to last one was on, then its end won't be visible
02937   if (prev_on && ribbon_width != 0) {
02938     cmdTriangle.putdata(ul[0], ll[0], ur[0], cmdList);
02939     cmdTriangle.putdata(ll[0], lr[0], ur[0], cmdList);
02940   }
02941 
02942   // swap the normals back, if they were swapped
02943   if (dot < 0) {
02944     perps[3L*i+0] = -perps[3L*i+0];
02945     perps[3L*i+1] = -perps[3L*i+1];
02946     perps[3L*i+2] = -perps[3L*i+2];
02947     norms[3L*i+0] = -norms[3L*i+0];
02948     norms[3L*i+1] = -norms[3L*i+1];
02949     norms[3L*i+2] = -norms[3L*i+2];
02950   }
02951 
02952   // what remains is the last one, which will be drawn as an arrow head
02953   if (atom_on[num-1] >= 0) {
02954     // indicate this atom can be picked
02955     pickpointcoords.append3(x[num-1], y[num-1], z[num-1]);
02956     pickpointindices.append(atom_on[num-1]);
02957 
02958     // recompute the base and tip information correctly
02959     // the 'norms' direction is 50% longer in each direction
02960     norms[3L*num-3] *= 1.5;
02961     norms[3L*num-2] *= 1.5;
02962     norms[3L*num-1] *= 1.5;
02963 
02964     SCALESUM(ul[0], centers+3L*num-3, norms+3L*num-3, ribbon_width);
02965     vec_copy    (ll[0], ul[0]);
02966     SCALEADD(ul[0], perps+3L*num-3, BETASCALE);
02967     SCALEADD(ll[0], perps+3L*num-3, -BETASCALE);
02968 
02969     SCALESUM(ur[0], centers+3L*num-3, norms+3L*num-3, -ribbon_width);
02970     vec_copy    (lr[0], ur[0]);
02971     SCALEADD(ur[0], perps+3L*num-3, BETASCALE);
02972     SCALEADD(lr[0], perps+3L*num-3, -BETASCALE);
02973 
02974     // and the tip has no norms component
02975     vec_copy    (ur[1], centers+3L*num);
02976     vec_copy    (lr[1], ur[1]);
02977     SCALEADD(ur[1], perps+3L*num, BETASCALE);
02978     SCALEADD(lr[1], perps+3L*num, -BETASCALE);
02979 
02980     // draw the arrow
02981     cmdColorIndex.putdata(color[atom_on[num-1]], cmdList);
02982 
02983     if (ribbon_width > 0) {
02984       // the normals are the same (from a copy) so I don't worry about them
02985       cmdTriangle.putdata(ul[0], ur[1], ur[0], cmdList);    // top
02986       // but that means here I _do_ have to switch the normal direction
02987       cmdTriangle.putdata(ll[0], lr[0], lr[1], cmdList);    // bottom
02988     
02989       cmdTriangle.putdata(ul[0], ll[0], lr[1], cmdList);  // sides
02990       cmdTriangle.putdata(ul[0], lr[1], ur[1], cmdList);
02991       cmdTriangle.putdata(ur[0], lr[0], lr[1], cmdList);
02992       cmdTriangle.putdata(ur[0], lr[1], ur[1], cmdList);
02993 
02994       cmdTriangle.putdata(ul[0], ll[0], ur[0], cmdList);  // and base
02995       cmdTriangle.putdata(ll[0], lr[0], ur[0], cmdList);
02996     } else {// ribbon_width == 0
02997       cmdLine.putdata(centers+3L*num-3, centers+3L*num, cmdList);
02998     } // special case for ribbon_width == 0
02999   }
03000   delete [] perps;
03001   delete [] norms;
03002   delete [] centers;
03003 
03004   // draw the pickpoints if we have any
03005   if (pickpointindices.num() > 0) {
03006     DispCmdPickPointArray pickPointArray;
03007     pickPointArray.putdata(pickpointindices.num(), &pickpointindices[0],
03008                            &pickpointcoords[0], cmdList);
03009   }
03010 
03011   // reset for the next round
03012   x.clear();
03013   y.clear();
03014   z.clear();
03015   atom_on.clear();
03016 }
03017 
03018 // draw based on secondary structure elements
03019 void DrawMolItem::draw_structure(float *framepos, float brad, int bres, int linethickness) {
03020   sprintf (commentBuffer,"Mol[%d] Rep[%d] Cartoon", mol->id(), repNumber);
03021   cmdCommentX.putdata(commentBuffer, cmdList);
03022 
03023   // Indicate that I need secondary structure information
03024   mol->need_secondary_structure(1);  // calculate 2ndary structure if need be
03025 
03026   int h_start, b_start; // atom indices that start beta sheets and helices
03027   int atom=0;
03028   int frag, res;
03029 
03030   ResizeArray<float> resx;      // coord of start/end of structure
03031   ResizeArray<float> resy;
03032   ResizeArray<float> resz;
03033   ResizeArray<int> CA_num;      // CA of start/end structure
03034   ResizeArray<int> extra_resid; // CA of very short residues
03035 
03036   ResizeArray<float> x;         // coords (x,y,z) of the residues
03037   ResizeArray<float> y;
03038   ResizeArray<float> z;
03039   ResizeArray<int> atom_on;     // turn on/off unselected parts of the helix
03040 
03041 
03042   //
03043   // draw protein alpha helices as cylinders
03044   //
03045   if (bres <= 2 || brad < 0.01) {
03046     cmdLineType.putdata(SOLIDLINE, cmdList);
03047     cmdLineWidth.putdata(2, cmdList);
03048     append(DMATERIALOFF);
03049   } else {
03050     append(DMATERIALON);
03051   }
03052   
03053   h_start = -1; // reset the helix starting atom index
03054   for (frag=0; frag<mol->pfragList.num(); frag++) {
03055     int num = mol->pfragList[frag]->num();
03056     for (int resindex=0; resindex < num; resindex++) {
03057       res = (*mol->pfragList[frag])[resindex];
03058       const int ss = mol->residue(res)->sstruct;
03059       int residue_is_helix = 
03060         (ss == SS_HELIX_ALPHA || ss == SS_HELIX_3_10 ||
03061          ss == SS_HELIX_PI);
03062 
03063       if (residue_is_helix) {
03064         atom = mol->find_atom_in_residue("CA", res);
03065         if (atom >= 0) {
03066           if (h_start == -1) {
03067             h_start = atom;             // just started a helix
03068           }
03069 
03070           x.append(framepos[3L*atom+0]); // add CA atom to the coordinate arrays
03071           y.append(framepos[3L*atom+1]);
03072           z.append(framepos[3L*atom+2]);
03073 
03074           if (atomSel->on[atom]) {      // atom_on contains either
03075             atom_on.append(atom);       // the atom index (if on)
03076           } else {
03077             atom_on.append(-1);         // or -1 (if off)
03078           }
03079         } else {
03080           msgErr << "Missing a CA in a protein residue!!" << sendmsg;
03081         }
03082       }
03083 
03084       // If we got something that wasn't helix, or if we're on the last
03085       // residue of the fragment, draw what we have.
03086       if (!residue_is_helix || resindex == num-1) {
03087         if (x.num() <= 1) {
03088           // msgWarn << "Cartoon will not draw a helix of length 1" << sendmsg;
03089           if (x.num() == 1) {
03090             extra_resid.append(atom);  // keep track of these for the coil
03091           }
03092           x.clear();                   // clear the coordinate arrays
03093           y.clear();
03094           z.clear();
03095         } else {
03096           float ends[2][3];
03097           CA_num.append(h_start);      // save the residues which start
03098           CA_num.append(atom);         // and end the helix structure
03099           draw_alpha_helix_cylinders(x, y, z, atom_on, atomColor->color, 
03100                                      brad, bres, ends[0], ends[1]);
03101 
03102           resx.append(ends[0][0]);     // save helix start/end coordinates
03103           resy.append(ends[0][1]);
03104           resz.append(ends[0][2]);
03105           resx.append(ends[1][0]);
03106           resy.append(ends[1][1]);
03107           resz.append(ends[1][2]);
03108         }
03109 
03110         h_start = -1; // reset the helix starting atom index
03111       } // drew helix
03112     } // went through residues
03113   } // went through pfrag
03114 
03115   //
03116   // now do beta sheets via another pass through the pfrags
03117   //
03118   if (linethickness == 0) {
03119     cmdLineType.putdata(SOLIDLINE, cmdList);
03120     cmdLineWidth.putdata(2, cmdList);
03121     append(DMATERIALOFF);
03122   } else {
03123     append(DMATERIALON);
03124   }
03125 
03126   b_start = -1; // reset the beta sheet starting atom index
03127   for (frag=0; frag<mol->pfragList.num(); frag++) {
03128     int num = mol->pfragList[frag]->num();
03129     for (int resindex=0; resindex < num; resindex++) {
03130       int res = (*mol->pfragList[frag])[resindex];
03131       const int ss = mol->residue(res)->sstruct;
03132       int is_beta = (ss == SS_BETA || ss == SS_BRIDGE);
03133       if (is_beta) {
03134         atom = mol->find_atom_in_residue("CA", res);
03135         if (atom >= 0) {
03136           if (b_start == -1) {
03137             b_start = atom;             // just started a sheet
03138           }
03139           x.append(framepos[3L*atom+0]); // add CA atom to the coordinate arrays
03140           y.append(framepos[3L*atom+1]);
03141           z.append(framepos[3L*atom+2]);
03142           if (atomSel->on[atom]) {      // atom_on contains either
03143             atom_on.append(atom);       // the atom index (if on)
03144           } else {
03145             atom_on.append(-1);         // or -1 (if off)
03146           }
03147         } else {
03148           msgErr << "Missing a CA in a protein residue!!" << sendmsg;
03149         }
03150       }
03151 
03152       // Draw what we have if we got to a non-beta residue, or if this was
03153       // the last residue.
03154       if (b_start != -1 && (!is_beta || resindex == num-1)) {
03155         if (x.num() <= 2) {
03156           // msgWarn << "Cartoon will not draw a sheet with less than "
03157           //         << "3 residues" << sendmsg;
03158 
03159           extra_resid.append(b_start); // keep track of these for the coil
03160           if (x.num() == 2) {
03161             extra_resid.append(atom);  // keep track of these for the coil
03162           }
03163           x.clear();                   // clear the coordinate arrays
03164           y.clear();
03165           z.clear();
03166           atom_on.clear();             // clear the 'on' array
03167         } else {
03168           float ends[2][3];
03169           CA_num.append(b_start);      // save the residues which start
03170           CA_num.append(atom);         // and end the beta sheet structure
03171           draw_beta_sheet(x, y, z, atom_on, atomColor->color,
03172                linethickness / 5.0f,
03173                ends[0], ends[1]);
03174           resx.append(ends[0][0]);     // save beta sheet start/end coordinates
03175           resy.append(ends[0][1]);
03176           resz.append(ends[0][2]);
03177           resx.append(ends[1][0]);
03178           resy.append(ends[1][1]);
03179           resz.append(ends[1][2]);
03180         }
03181 
03182         b_start = -1; // reset the beta sheet starting atom index
03183       }
03184     }  // went through the residues
03185   } // went through pfrag
03186 
03187 
03188   // Finally, everything left I can draw as a spline.  The tricky part
03189   // is I need to get the start/end coordinates of the cylinders and
03190   // sheets correct, for otherwise  the tube doesn't meet with them nicely.
03191   // Luckily, I kept track of the start/end coordinates as I went along, so
03192   // I swap the coords of that array with the frampos data, replace atomSel'
03193   // 'on' with a new 'on', and call draw_tube.  Code reuse, right?  
03194   // Anyway, when I return I swap everything back
03195   {
03196     // finding the atoms is easy; search for 'protein not (sheet or helix)'
03197     // and not those atoms in the same residue as the extra_resid list.
03198     int nres, i;
03199 
03200     // allocate and clear a new 'on' list
03201     int *on = new int[mol->nAtoms];
03202     memset(on, 0, mol->nAtoms*sizeof(int));
03203 
03204     // find "protein and not (sheet or helix)"
03205     nres=mol->residueList.num();
03206     for (i=0; i<nres; i++) {
03207       const Residue *res = mol->residueList[i];
03208       if (res->residueType == RESPROTEIN) {
03209         const int ss = res->sstruct;
03210         if (ss == SS_TURN || ss == SS_COIL) {
03211           int numatoms = res->atoms.num();
03212           for (int j=0; j<numatoms; j++) {
03213             on[res->atoms[j]] = 1;
03214           }
03215         }
03216       }
03217     }
03218     
03219     // collect the leftover residues in "short" sheets or cylinders
03220     // and draw them as tubes as well
03221     nres=extra_resid.num();
03222     for (i=0; i<nres; i++) {
03223       int res = mol->atom(extra_resid[i])->uniq_resid;
03224       Residue *r = mol->residue(res);
03225       int numatoms = r->atoms.num();
03226       for (int j=0; j<numatoms; j++) {
03227         on[r->atoms[j]] = 1;
03228       }
03229     }
03230 
03231     // swap the real CA atom coordinates with calculated 
03232     // control point coordinates so that draw_tube generates
03233     // the right shape.
03234     float tmp_coords[3];
03235     int atom; 
03236     int numcaatoms = CA_num.num();
03237     for (i=0; i<numcaatoms; i++) {
03238       atom = 3L*CA_num[i];
03239       vec_copy(tmp_coords, framepos+atom);
03240       framepos[atom+0] = resx[i];
03241       framepos[atom+1] = resy[i];
03242       framepos[atom+2] = resz[i];
03243       resx[i] = tmp_coords[0];
03244       resy[i] = tmp_coords[1];
03245       resz[i] = tmp_coords[2];
03246       // turn CA of the start/end on as well
03247       on[CA_num[i]] = 2+i%2;    
03248       // the "2" is a special case for draw_spline_curve which will only
03249       // draw the first 1/2 if it is a 2 or second half if it is a 3
03250     }
03251 
03252     // replace the list of "on" atoms with the new list
03253     for (i=0; i<mol->nAtoms; i++) {
03254       if (!atomSel->on[i]) { // mask the coil against the selected atoms
03255         on[i] = 0;
03256       }
03257     }
03258     int *temp_on = atomSel->on;
03259     atomSel->on = on;
03260 
03261     // make the tube
03262     draw_tube(framepos, 
03263               brad / 7.0f,
03264               int(float(bres) * 2.0f/3.0f + 0.5f));
03265 
03266     // and revert to the original coordinate and 'on' arrays
03267     atomSel->on = temp_on;
03268     for (i=0; i<numcaatoms; i++) {
03269       atom = 3L*CA_num[i];
03270       framepos[atom+0] = resx[i];
03271       framepos[atom+1] = resy[i];
03272       framepos[atom+2] = resz[i];
03273     }
03274 
03275     delete [] on; // delete the temporary 'on' array
03276   }
03277 
03278   append(DMATERIALOFF);
03279 }
03280 
03281 int DrawMolItem::pickable_on() {
03282   return displayed() && mol && mol->displayed();
03283 }
03284 

Generated on Fri Mar 29 02:45:11 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002