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

DrawMolItemRibbons.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: DrawMolItemRibbons.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.157 $      $Date: 2021/10/28 21:12:15 $
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 "DrawMolecule.h"
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 #include "WKFUtils.h"
00041 
00042 #define RIBBON_ERR_NOTENOUGH            1
00043 #define RIBBON_ERR_PROTEIN_MESS         2
00044 #define RIBBON_ERR_MISSING_PHOSPHATE    4
00045 #define RIBBON_ERR_MISSING_O1P_O2P      8
00046 #define RIBBON_ERR_BADNUCLEIC          16
00047 #define RIBBON_ERR_BADPYRIMIDINE       32
00048 #define RIBBON_ERR_BADPURINE           64
00049 
00050 #if 0
00051 #define RIBBON_ERR_NOTENOUGH_STR         "Cannot draw a ribbon for a nucleic acid with only one element"
00052 #define RIBBON_ERR_PROTEIN_MESS_STR      "Someone's been messing around with the definition of a protein!\nThings are going to look messy"
00053 #define RIBBON_ERR_MISSING_PHOSPHATE_STR "Cannot find first phosphate of a nucleic acid, so it won't be drawn"
00054 #define RIBBON_ERR_MISSING_O1P_O2P_STR   "Cannot find both O1P and O2P of a nucleic acid, so it won't be drawn"
00055 #define RIBBON_ERR_BADNUCLEIC_STR        "Someone has redefined the nucleic acid atom names!"
00056 #define RIBBON_ERR_BADPYRIMIDINE_STR     "Tried to draw a nucleic acid residue I thought was a pyrimidine, but it doesn't have the right atom names."
00057 #define RIBBON_ERR_BADPURINE_STR         "Tried to draw a nucleic acid residue I thought was a purine, but it doesn't have the right atom names."
00058 #endif
00059 
00060 // draw a spline ribbon.  The math is the same as a tube. 
00061 // However, two splines are calculated, one for coords+perp
00062 // and the other along coords-perp.  Triangles are used to
00063 // connect the curves.  If cylinders are used,  they are
00064 // drawn along the spline paths.  The final result should
00065 // look something like:
00066 //    ooooooooooooooooooo   where the oooo are the (optional) cylinders
00067 //    ! /| /| /| /| /| /!   the ! are a control point, in the ribbon center
00068 //    |/ |/ |/ |/ |/ |/ |   the /| are the triangles  (size per interpolation)
00069 //    ooooooooooooooooooo   the edges go through coords[i] +/- offset[i]
00070 void DrawMolItem::draw_spline_ribbon(int num, float *coords,
00071         float *offset, int *idx, int use_cyl, float b_rad, int b_res) {
00072   ResizeArray<float> pickpointcoords;
00073   ResizeArray<int> pickpointindices;
00074   float q[4][3];
00075   // for the next 2 variables, the information for this residue
00076   // starts at element 2; elements 0 and 1 are copies of the last
00077   // two elements of the previous residue
00078   float pts[2][9][3];  // data for 2 edges, 7 points (+2 overlaps from
00079                        // the previous interpolation) and x,y,z
00080    
00081   //  contains the norms as [edge][interpolation point][x/y/z]
00082   float norms[2][8][3]; 
00083   // contains the summed norms of the sequential triangles
00084   // the first element (# 0) contains the info for what should have
00085   // been the last two triangles of the previous residue.  In other
00086   // words, tri_norm[0][0] contains the summed norm for the point
00087   // located at pts[1][0]
00088   float tri_norm[2][7][3]; 
00089 
00090   int last_loop = -10;
00091   int loop, i, j;
00092   float *tmp_ptr = (float *) malloc(2L*(num+4L) * sizeof(float) * 3L);
00093   if (tmp_ptr == NULL) {
00094     msgErr << "Cannot make a ribbon; not enough memory!" << sendmsg;
00095     return;
00096   }
00097   float *edge[2];
00098 
00099   // XXX disgusting array math here, rewrite!
00100   // copy the coordinates +/- the offsets into the temp ("edge") arrays
00101   edge[0] = tmp_ptr + 2*3L;
00102   memcpy(edge[0]-2*3L, coords-2*3L, (num+4L)*sizeof(float)*3L);
00103   edge[1] = edge[0] + (num+4)*3;
00104   memcpy(edge[1]-2*3L, coords-2*3L, (num+4L)*sizeof(float)*3L);
00105   for (j=-2*3; j<(num+2)*3L-1; j++) {
00106     edge[0][j] += offset[j];
00107     edge[1][j] -= offset[j];
00108   }
00109 
00110   // go through the data points
00111   for (loop=-1; loop<num; loop++) {
00112     int j;
00113 
00114     // If I'm to draw anything....
00115     if ((idx[loop] >= 0 && atomSel->on[idx[loop]]) ||
00116         (idx[loop+1] >= 0 && atomSel->on[idx[loop+1]])) {
00117 
00118       // construct the interpolation points (into the "pts" array)
00119       // remember, these are offset by two to keep some information
00120       // about the previous residue in the first 2 elements
00121       for (i=0; i<=1; i++) {
00122         make_spline_Q_matrix(q, spline_basis, edge[i]+(loop-1)*3);
00123         make_spline_interpolation(pts[i][2], 0.0f/6.0f, q);
00124         make_spline_interpolation(pts[i][3], 1.0f/6.0f, q);
00125         make_spline_interpolation(pts[i][4], 2.0f/6.0f, q);
00126         make_spline_interpolation(pts[i][5], 3.0f/6.0f, q);
00127         make_spline_interpolation(pts[i][6], 4.0f/6.0f, q);
00128         make_spline_interpolation(pts[i][7], 5.0f/6.0f, q);
00129         make_spline_interpolation(pts[i][8], 6.0f/6.0f, q);
00130       }
00131 
00132       // make the normals for each new point.
00133       for (i=2; i<8; i++) {
00134         float diff1[3], diff2[3];
00135         vec_sub(diff1, pts[1][i+1], pts[1][i+0]);
00136         vec_sub(diff2, pts[1][i+1], pts[0][i+0]);
00137         cross_prod(norms[1][i], diff1, diff2);
00138         vec_sub(diff1, pts[0][i+1], pts[0][i+0]);
00139         cross_prod(norms[0][i], diff1, diff2);
00140       }
00141 
00142       // if this wasn't a continuation, I need to set the
00143       // first 2 elements properly so the norms are smooth
00144       if (last_loop != loop-1) {
00145         vec_copy(norms[0][0], norms[0][2]);
00146         vec_copy(norms[1][0], norms[1][2]);
00147         vec_copy(norms[0][1], norms[0][2]);
00148         vec_copy(norms[1][1], norms[1][2]);
00149       }
00150 
00151       // sum up the values of neighboring triangles to make
00152       // a smooth normal
00153       for (j=0; j<=1; j++) {
00154         for (i=0; i<8-1; i++) {
00155           vec_add(tri_norm[j][i], norms[j][i+1],     // "this" triangle
00156                   norms[1-j][i+1]);  // opposite strand
00157           vec_add(tri_norm[j][i], tri_norm[j][i],
00158                   norms[j][i]);      // prev on strand
00159         }
00160       }
00161 
00162       // pre-normalize the normals so we don't need to have
00163       // OpenGL doing this for us repetitively on every frame,
00164       // this allows use to use the GL_RESCALE_NORMAL extension
00165       for (j=0; j<=1; j++) {
00166         for (i=0; i<8-1; i++) {
00167           vec_normalize(tri_norm[j][i]);
00168         } 
00169       }
00170          
00171       // draw what I need for atom 'loop'
00172       if (idx[loop] >= 0 &&          // this is a real atom
00173           atomSel->on[idx[loop]]) {  // and it it turned on
00174 
00175         if (last_loop != loop - 1) {  
00176           // do prev. points exist? if not then I don't know if the color 
00177           // was properly set, so set it here
00178           cmdColorIndex.putdata(atomColor->color[idx[loop]], cmdList);
00179         }
00180 
00181         // draw the cylinders to finish off the last residue, if
00182         // need be, and draw the ones for this half of the residue
00183         // Cylinders are drawn on the top and bottom of the residues
00184         if (use_cyl) { // draw top/bot edge cylinders if need be
00185           // Special case the first cylinder because I want
00186           // it to be a smooth continuation from the previous
00187           // cylinder; assuming it exists
00188           if (last_loop != loop-1) {  // continue from previous?
00189             int ii;
00190             for (ii=0; ii<=1; ii++) {  // doesn't exist, so
00191               make_connection(NULL, pts[ii][2], pts[ii][3],
00192                               pts[ii][4], b_rad, b_res, use_cyl);
00193             }
00194           } else { // there was a previous cylinder, so be smooth
00195             for (i=0; i<=1; i++) {
00196                make_connection(pts[i][0], pts[i][1], pts[i][2], pts[i][3],
00197                                b_rad, b_res, use_cyl);
00198                make_connection(pts[i][1], pts[i][2], pts[i][3], pts[i][4],
00199                                b_rad, b_res, use_cyl);
00200             }
00201           }
00202                
00203           // and draw the rest of the cylinders for this 1/2 residue
00204           for (i=0; i<=1; i++) {
00205             make_connection(pts[i][2], pts[i][3], pts[i][4], pts[i][5],
00206                             b_rad, b_res, use_cyl);
00207             make_connection(pts[i][3], pts[i][4], pts[i][5], pts[i][6],
00208                             b_rad, b_res, use_cyl);
00209           }
00210         } // drew cylinders
00211 
00212         // Draw the triangles that connect the cylinders and make up the 
00213         // ribbon proper.  The funky start condition is so that it starts at
00214         // pts[][1] if I need to finish the previous residue, or
00215         // pts[][2] if I don't
00216         for (i= (last_loop == loop-1 ? 1: 2); i<5; i++) {
00217           cmdTriangle.putdata(pts[1][1+i], pts[1][0+i], pts[0][0+i],
00218                               tri_norm[1][0+i], tri_norm[1][-1+i],
00219                               tri_norm[0][-1+i], cmdList);
00220           cmdTriangle.putdata(pts[0][1+i], pts[1][1+i], pts[0][0+i],
00221                               tri_norm[0][0+i], tri_norm[1][0+i],
00222                               tri_norm[0][-1+i], cmdList);
00223 
00224           // indicate this atom can be picked
00225           // this spot is in the middle of the ribbon, both in length
00226           // and in width
00227           int pidx = loop * 3;
00228           pickpointcoords.append3(&coords[pidx]);
00229           pickpointindices.append(idx[loop]);
00230         }
00231       }  // atom 'loop' finished
00232          
00233       // draw what I can for atom 'loop+1'
00234       if (idx[loop+1] >= 0 && atomSel->on[idx[loop+1]]) {
00235         // If this is on, then I may have to change the color,
00236         // since I'm lazy, I won't check to see if I _have_ to change it
00237         // but assume I do
00238         cmdColorIndex.putdata(atomColor->color[idx[loop+1]], cmdList);
00239         // I can draw the first couple of cylinders, but I need knowledge
00240         // of what goes on next to get a smooth fit.  Thus, I won't
00241         // draw them here.
00242         // I can't draw the last two cylinders.
00243         if (use_cyl) {
00244           for (i=0; i<=1; i++) {
00245             make_connection(pts[i][4], pts[i][5], pts[i][6], pts[i][7],
00246                             b_rad, b_res, use_cyl);
00247             make_connection(pts[i][5], pts[i][6], pts[i][7], pts[i][8],
00248                             b_rad, b_res, use_cyl);
00249           }
00250         }
00251 
00252         // I can draw 3 of the four triangles, but I need
00253         // the normals to do the last one properly
00254         for (i= 5; i<8-1; i++) {
00255 /*
00256           cmdTriangle.putdata(pts[0][0+i], pts[1][0+i], pts[1][1+i],
00257                               tri_norm[0][-1+i], tri_norm[1][-1+i],
00258                               tri_norm[1][0+i], cmdList);
00259           cmdTriangle.putdata(pts[0][0+i], pts[1][1+i], pts[0][1+i],
00260                               tri_norm[0][-1+i], tri_norm[1][0+i],
00261                               tri_norm[0][0+i], cmdList);
00262 */
00263           cmdTriangle.putdata(pts[1][1+i], pts[1][0+i], pts[0][0+i],
00264                               tri_norm[1][0+i], tri_norm[1][-1+i],
00265                               tri_norm[0][-1+i], cmdList);
00266           cmdTriangle.putdata(pts[0][1+i], pts[1][1+i], pts[0][0+i],
00267                               tri_norm[0][0+i], tri_norm[1][0+i],
00268                               tri_norm[0][-1+i], cmdList);
00269         }
00270         last_loop = loop;
00271       } // atom 'loop+1' finished
00272          
00273       // save infor for next loop, for smoothing cylinders and normals
00274       for (i=0; i<=1; i++) {
00275         vec_copy(pts[i][0], pts[i][6]);  // remember, because of the spline,
00276         vec_copy(pts[i][1], pts[i][7]);  // loop pts[][8] is loop+1 pts[][2]
00277         vec_copy(norms[i][0], norms[i][6]);
00278         vec_copy(norms[i][1], norms[i][7]);
00279       }
00280     } 
00281   } // gone down the fragment
00282 
00283   free(tmp_ptr);
00284 
00285   // draw the pickpoints if we have any
00286   if (pickpointindices.num() > 0) {
00287     DispCmdPickPointArray pickPointArray;
00288     pickPointArray.putdata(pickpointindices.num(), &pickpointindices[0],
00289                            &pickpointcoords[0], cmdList);
00290   }
00291 }
00292 
00293 
00294 
00295 
00296 // draw ribbons along the protein backbone
00297 // part of this method taken (with permission) from the 'ribbon.f' in the
00298 // Raster3D package:
00299 //  Merritt, Ethan A. and Murphy, Michael E.P. (1994).
00300 //   "Raster3D Version 2.0, a Program for Photorealistic Molecular Graphics"
00301 //        Acta Cryst. D50, 869-873.
00302 
00303 // That method was based on ideas from  Carson & Bugg, J. Molec. Graphics
00304 //   4,121-122 (1986)
00305 void DrawMolItem::draw_ribbons(float *framepos, float brad, int bres, float thickness) {
00306   sprintf(commentBuffer,"Mol[%d] Rep[%d] Ribbons", mol->id(), repNumber);
00307   cmdCommentX.putdata(commentBuffer, cmdList);
00308 
00309   // find out if I'm using lines or cylinders
00310   int use_cyl = TRUE; // use cylinders by default
00311   int rc = 0;  // clear errors
00312   
00313   if (bres <= 2 || brad < 0.01f) { 
00314     use_cyl = FALSE; // the representation will be as lines
00315   } 
00316 
00317   append(DMATERIALON); // enable shading 
00318 
00319   float ribbon_width = thickness;
00320   if (ribbon_width < 1.0) 
00321     ribbon_width = 1.0;
00322   ribbon_width /= 7.0;
00323 
00324   rc |= draw_protein_ribbons_old(framepos, bres, brad, ribbon_width, use_cyl);
00325   rc |= draw_nucleic_ribbons(framepos, bres, brad, ribbon_width, use_cyl, 0, 0);
00326   rc |= draw_base_sugar_rings(framepos, bres, brad, ribbon_width, use_cyl);
00327 
00328   // XXX put the more specific error messages back in here if we feel
00329   //     that they are helpful to the user.
00330   if (rc != 0) {
00331     if (emitstructwarning())
00332       msgErr << "Warning: ribbons code encountered an unusual structure, geometry may not look as expected." << sendmsg;
00333   }
00334 }
00335 
00336 
00337 
00338 
00339 int DrawMolItem::draw_protein_ribbons_old(float *framepos, int b_res, float b_rad,
00340                                 float ribbon_width, int use_cyl) {
00341   float *real_coords = NULL, *coords;
00342   float *real_perps = NULL, *perps;
00343   int *real_idx = NULL, *idx;
00344   float *capos, *last_capos;
00345   float *opos,  *last_opos;
00346   int onum, canum, res, frag, num;
00347   int rc=0;
00348 
00349   // these are the variables used in the Raster3D package
00350   float a[3], b[3], c[3], d[3], e[3], g[3];  
00351 
00352   // go through each protein and find the CA and O
00353   // from that, construct the offsets to use (called "perps")
00354   // then do two splines and connect them
00355   for (frag=0; frag<mol->pfragList.num(); frag++) {
00356     int loop;
00357 
00358     num = mol->pfragList[frag]->num();  // number of residues in this fragment
00359     if (num < 2) {
00360       rc |= RIBBON_ERR_NOTENOUGH;
00361       continue; // can't draw a ribbon with only one element, so skip
00362     }
00363 
00364     // check that we have a valid structure before continuing
00365       res = (*mol->pfragList[frag])[0];
00366     canum = mol->find_atom_in_residue("CA", res); // get CA 
00367      onum = mol->find_atom_in_residue("O", res);  // get O
00368 
00369     if (canum < 0 || onum < 0) {
00370       continue; // can't find 1st CA or O of the protein, so don't draw
00371     }
00372 
00373     if (real_coords) {
00374       free(real_coords);
00375       free(real_idx);
00376       free(real_perps);
00377       real_coords = NULL;
00378       real_idx = NULL;
00379       real_perps = NULL;
00380     }
00381 
00382     // allocate memory for control points, perps, and index storage
00383     // The arrays have 2 extra elements at the front and end to allow
00384     // duplication of control points necessary for the old ribbons code.
00385     real_coords = (float *) malloc((num+4) * sizeof(float)*3);
00386        real_idx =   (int *) malloc((num+4) * sizeof(int));
00387      real_perps = (float *) malloc((num+4) * sizeof(float)*3);
00388 
00389     coords = real_coords + 2*3L;
00390        idx = real_idx + 2L;
00391      perps = real_perps + 2*3L;
00392 
00393     // initialize CA and O atom pointers
00394     capos = framepos + 3L*canum;
00395     last_capos = capos;
00396     opos = framepos + 3L*onum;
00397     last_opos = opos;
00398 
00399     // duplicate first control point data
00400     vec_copy(coords-6, capos);
00401     vec_copy(coords-3, capos);
00402     idx[-2] = idx[-1] = -1;
00403 
00404     // now go through and set the coordinates and the perps
00405     e[0] = e[1] = e[2] = 0.0;
00406     vec_copy(g, e);
00407     for (loop=0; loop<num; loop++) {
00408       res = (*mol->pfragList[frag])[loop];
00409       canum = mol->find_atom_in_residue("CA", res);
00410       if (canum >= 0) {
00411         capos = framepos + 3L*canum;
00412       }
00413       onum = mol->find_atom_in_residue("O", res);
00414       if (onum < 0) {
00415         onum = mol->find_atom_in_residue("OT1", res);
00416       }
00417 
00418       if (onum >= 0) {
00419         opos = framepos + 3L*onum;
00420       } else {
00421         rc |= RIBBON_ERR_PROTEIN_MESS;
00422         opos = last_opos; // try to cope if we have no oxygen index
00423       }
00424 
00425       vec_copy(coords+loop*3, capos);
00426       idx[loop] = canum;
00427 
00428       // now I need to figure out where the ribbon goes
00429       vec_sub(a, capos, last_capos);     // A=(pos(CA(res)) - pos(CA(res-1)))
00430       vec_sub(b, last_opos, last_capos); // B=(pos(O(res-1)) - pos(CA(res-1)))
00431       cross_prod(c, a, b);               // C=AxB, define res-1's peptide plane
00432       cross_prod(d, c, a);               // D=CxA, normal to plane and backbone
00433       // if normal is > 90 deg from  previous one, negate the new one
00434       if (dot_prod(d, g) < 0) {
00435         vec_negate(b, d);
00436       } else {
00437         vec_copy(b, d);
00438       }
00439       vec_add(e, g, b);   // average to the sum of the previous vectors
00440       vec_normalize(e);
00441       vec_scale(&perps[3L*loop], ribbon_width, e); // compute perps from normal
00442       vec_copy(g, e);     // make a cumulative sum; cute, IMHO
00443       last_capos = capos;
00444       last_opos = opos;
00445     }
00446 
00447     // duplicate the last control point into two extra copies
00448     vec_copy(coords+3L*num, capos);
00449     vec_copy(coords+3L*(num+1), capos);
00450     idx[num] = idx[num+1] = -1;
00451 
00452     // copy the second perp to the first since the first one didn't have
00453     // enough data to construct a good perpendicular.
00454     vec_copy(perps, perps+3);
00455 
00456     // duplicate the first and last perps into the extra control points
00457     vec_copy(perps-3, perps);
00458     vec_copy(perps-6, perps);
00459     vec_copy(perps+3L*num, perps+3L*num-3);
00460     vec_copy(perps+3L*num+3, perps+3L*num);
00461 
00462     draw_spline_ribbon(num, coords, perps, idx, use_cyl, b_rad, b_res);
00463 
00464     if (real_coords) {
00465       free(real_coords);
00466       free(real_idx);
00467       free(real_perps);
00468       real_coords = NULL;
00469       real_idx = NULL;
00470       real_perps = NULL;
00471     }
00472   } // drew the protein fragment ribbon
00473 
00474   return rc;
00475 }
00476 
00477 
00478 int DrawMolItem::draw_protein_ribbons_new(float *framepos, int b_res, float b_rad,
00479                                 float ribbon_width, int use_cyl) {
00480   float *coords = NULL;
00481   float *perps = NULL;
00482   int *idx = NULL;
00483   float *capos, *last_capos, *opos, *last_opos;
00484   int onum, canum, frag, num, res;
00485   const char *modulatefield = NULL; // XXX this needs to become a parameter
00486   const float *modulatedata = NULL; // data field to use for width modulation
00487   float *modulate = NULL;           // per-control point width values
00488   int rc=0;
00489 
00490   // these are the variables used in the Raster3D package
00491   float a[3], b[3], c[3], d[3], e[3], g[3];  
00492 
00493   // Lookup atom typecodes ahead of time so we can use the most 
00494   // efficient variant of find_atom_in_residue() in the main loop.
00495   // If we can't find the atom types we need, bail out immediately
00496   int CAtypecode  = mol->atomNames.typecode("CA");
00497   int Otypecode   = mol->atomNames.typecode("O");
00498   int OT1typecode = mol->atomNames.typecode("OT1");
00499 
00500   if (CAtypecode < 0 || ((Otypecode < 0) && (OT1typecode < 0))) {
00501     return rc; // can't draw a ribbon without CA and O atoms for guidance
00502   }
00503 
00504   // allocate for the maximum possible per-residue control points, perps, 
00505   // and indices so we don't have to reallocate for every fragment
00506   coords = (float *) malloc((mol->nResidues) * sizeof(float)*3);
00507      idx =   (int *) malloc((mol->nResidues) * sizeof(int));
00508    perps = (float *) malloc((mol->nResidues) * sizeof(float)*3);
00509 
00510 #if 1
00511   // XXX hack to let users try various stuff
00512   modulatefield = getenv("VMDMODULATERIBBON");
00513 #endif
00514   if (modulatefield != NULL) {
00515     if (!strcmp(modulatefield, "user")) {
00516       modulatedata = mol->current()->user;
00517       // XXX user field can be NULL on some timesteps
00518     } else {
00519       modulatedata = mol->extraflt.data(modulatefield);
00520     } 
00521 
00522     // allocate for the maximum possible per-residue modulation values
00523     // so we don't have to reallocate for every fragment
00524     // XXX the modulate array is allocated cleared to zeros so that
00525     // in the case we get a NULL modulatedata pointer (user field for
00526     // example, we'll just end up with no modulation.
00527     modulate = (float *) calloc(1, mol->nResidues * sizeof(float));
00528   }
00529 
00530   // go through each protein and find the CA and O
00531   // from that, construct the offsets to use (called "perps")
00532   int numpfrag = mol->pfragList.num();
00533   for (frag=0; frag<numpfrag; frag++) {
00534     int cyclic=mol->pfragCyclic[frag];  // whether fragment is cyclic
00535     num = mol->pfragList[frag]->num();  // number of residues in this fragment
00536     if (num < 2) {
00537       rc |= RIBBON_ERR_NOTENOUGH;
00538       continue; // can't draw a ribbon with only one element, so skip
00539     }
00540 
00541     // check that we have a valid structure before continuing
00542       res = (*mol->pfragList[frag])[0];
00543     canum = mol->find_atom_in_residue(CAtypecode, res);
00544      onum = mol->find_atom_in_residue(Otypecode, res);
00545 
00546     if (canum < 0 || onum < 0) {
00547       continue; // can't find 1st CA or O of the protein, so don't draw
00548     }
00549 
00550     // initialize CA and O atom pointers
00551     capos = framepos + 3L*canum;
00552     last_capos = capos;
00553     opos = framepos + 3L*onum;
00554     last_opos = opos;
00555 
00556     // now go through and set the coordinates and the perps
00557     e[0] = e[1] = e[2] = g[0] = g[1] = g[2] = 0.0;
00558 
00559     // for a cyclic structure, we use the positions of the last residue
00560     // to seed the initial direction vectors for the ribbon
00561     if (cyclic) {
00562       int lastres = (*mol->pfragList[frag])[num-1];
00563       int lastcanum = mol->find_atom_in_residue(CAtypecode, lastres);
00564       last_capos = framepos + 3L*lastcanum;
00565 
00566       int lastonum = mol->find_atom_in_residue(Otypecode, lastres);
00567       if (lastonum < 0 && OT1typecode >= 0) {
00568         lastonum = mol->find_atom_in_residue(OT1typecode, lastres);
00569       }
00570       last_opos = framepos + 3L*lastonum;
00571 
00572       // now I need to figure out where the ribbon goes
00573       vec_sub(a, capos, last_capos);     // A=(pos(CA(res)) - pos(CA(res-1)))
00574       vec_sub(b, last_opos, last_capos); // B=(pos(O(res-1)) - pos(CA(res-1)))
00575       cross_prod(c, a, b);               // C=AxB, define res-1's peptide plane
00576       cross_prod(d, c, a);               // D=CxA, normal to plane and backbone
00577       // if normal is > 90 deg from  previous one, negate the new one
00578       if (dot_prod(d, g) < 0) {
00579         vec_negate(b, d);
00580       } else {
00581         vec_copy(b, d);
00582       }
00583       vec_add(e, g, b);            // average to the sum of the previous vectors
00584       vec_normalize(e);
00585       vec_copy(g, e);              // make a cumulative sum; cute, IMHO
00586     }
00587 
00588     int loop;
00589     for (loop=0; loop<num; loop++) {
00590       res = (*mol->pfragList[frag])[loop];
00591 
00592       canum = mol->find_atom_in_residue(CAtypecode, res);
00593       if (canum >= 0) {
00594         capos = framepos + 3L*canum;
00595       }
00596 
00597       onum = mol->find_atom_in_residue(Otypecode, res);
00598       if (onum < 0 && OT1typecode >= 0) {
00599         onum = mol->find_atom_in_residue(OT1typecode, res);
00600       }
00601       if (onum >= 0) {
00602         opos = framepos + 3L*onum;
00603       } else {
00604         rc |= RIBBON_ERR_PROTEIN_MESS;
00605         opos = last_opos; // try to cope if we have no oxygen index
00606       }
00607 
00608       // copy the CA coordinate into the control point array
00609       vec_copy(coords+loop*3L, capos);
00610 
00611       // modulate the ribbon width by user-specified per-atom data
00612       if (modulatedata != NULL)
00613         modulate[loop] = modulatedata[canum];
00614  
00615       idx[loop] = canum;
00616 
00617       // now I need to figure out where the ribbon goes
00618       vec_sub(a, capos, last_capos);     // A=(pos(CA(res)) - pos(CA(res-1)))
00619       vec_sub(b, last_opos, last_capos); // B=(pos(O(res-1)) - pos(CA(res-1)))
00620       cross_prod(c, a, b);               // C=AxB, define res-1's peptide plane
00621       cross_prod(d, c, a);               // D=CxA, normal to plane and backbone
00622       // if normal is > 90 deg from  previous one, negate the new one
00623       if (dot_prod(d, g) < 0) {
00624         vec_negate(b, d);
00625       } else {
00626         vec_copy(b, d);
00627       }
00628       vec_add(e, g, b);            // average to the sum of the previous vectors
00629       vec_normalize(e);
00630       vec_copy(&perps[3L*loop], e); // compute perps from the normal
00631       vec_copy(g, e);              // make a cumulative sum; cute, IMHO
00632       last_capos = capos;
00633       last_opos = opos;
00634     }
00635 
00636     if (!cyclic) {
00637       // copy the second perp to the first since the first one didn't have
00638       // enough data to construct a good perpendicular.
00639       vec_copy(perps, perps+3);
00640     } 
00641 
00642     if (modulate != NULL) {
00643       // modulate ribbon width by user-specified per-atom field value
00644       float *widths = (float *) malloc(num * sizeof(float));
00645       float *heights = (float *) malloc(num * sizeof(float));
00646       float m_fac;
00647       int i;
00648       for (i=0; i<num; i++) {
00649 #if 1
00650         // only allow modulation values > 0. otherwise fall back on old scheme.
00651         m_fac = modulate[i];
00652         if (m_fac <= 0.0f)
00653            m_fac = 1.0f;
00654 
00655         // modulate both width and height
00656         widths[i] = 7.0f * ribbon_width * b_rad * m_fac;
00657         heights[i] = b_rad*m_fac;
00658 #else
00659         // modulate only width, and only additively
00660         widths[i] = 7.0f * ribbon_width * b_rad + m_fac;
00661         heights[i] = b_rad;
00662 #endif
00663       }
00664       draw_spline_new(num, coords, perps, idx, widths, heights, num, b_res, cyclic);
00665       free(widths);
00666       free(heights);
00667     } else {
00668       // draw normal unmodulated ribbons
00669       float widths = 7.0f * ribbon_width * b_rad;
00670       float heights = b_rad;
00671       draw_spline_new(num, coords, perps, idx, &widths, &heights, 1, b_res, cyclic);
00672     }
00673   } // drew the protein fragment ribbon
00674 
00675   if (coords) {
00676     free(coords);
00677     free(idx);
00678     free(perps);
00679   }
00680 
00681   if (modulate) {
00682     free(modulate);
00683   }
00684 
00685   return rc;
00686 }
00687 
00688 
00689 
00690 int DrawMolItem::draw_nucleic_ribbons(float *framepos, int b_res, float b_rad,
00691                                 float ribbon_width, int use_cyl, int use_new, 
00692                                 int use_carb) {
00694   float *real_coords = NULL, *coords;
00695   float *real_perps = NULL, *perps;
00696   int *real_idx = NULL, *idx;
00697   float *ppos, *last_ppos;
00698   float opos[3], last_opos[3];
00699   int cpnum;         // index of spline control point atom
00700   int o1num, o2num;  // indices of atoms used to construct perps
00701   int frag, num;
00702   int rc=0;
00703 
00704   // these are the variables used in the Raster3D package
00705   float a[3], b[3], c[3], d[3], e[3], g[3];  
00706 
00707   // go through each nucleic acid and find the phospate
00708   // then find the O1P/O2P, or OP1/OP2 (new nomenclature). 
00709   // From those construct the perps then do two splines and connect them
00710   for (frag=0; frag<mol->nfragList.num(); frag++) {
00711     int loop;
00712 
00713     num = mol->nfragList[frag]->num(); // number of residues in this fragment
00714     if (num < 2) {
00715       rc |= RIBBON_ERR_NOTENOUGH;
00716       continue;
00717     }
00718     if (real_coords) {
00719       free(real_coords);
00720       free(real_idx);
00721       free(real_perps);
00722       real_coords = NULL;
00723       real_idx = NULL;
00724       real_perps = NULL;
00725     }
00726     real_coords = (float *) malloc((num+4)*sizeof(float)*3);
00727     real_idx = (int *) malloc((num+4) * sizeof(int));
00728     real_perps = (float *) malloc((num+4) * sizeof(float)*3);
00729     coords = real_coords + 2*3L;
00730        idx = real_idx + 2L;
00731      perps = real_perps + 2*3L;
00732          
00733     // okay, I've got space for the coordinates, the index translations,
00734     // and the perpendiculars, now initialize everything
00735     int res = (*mol->nfragList[frag])[0];
00736 
00737     // get atoms to use as the spline control points
00738     if (use_carb) {
00739       // use the carbons
00740       cpnum = mol->find_atom_in_residue("C3'", res);
00741       if (cpnum < 0)
00742         cpnum = mol->find_atom_in_residue("C3*", res);
00743       // use the phosphates if no carbon control points found
00744       if (cpnum < 0)
00745         cpnum = mol->find_atom_in_residue("P", res);
00746     } else { 
00747       // use the phosphates 
00748       cpnum = mol->find_atom_in_residue("P", res);
00749     }
00750 
00751     // if no P found, check the terminal atom names
00752     if (cpnum < 0) {
00753       cpnum = mol->find_atom_in_residue("H5T", res);
00754     }
00755     if (cpnum < 0) {
00756       cpnum = mol->find_atom_in_residue("H3T", res);
00757     }
00758 
00759     if (cpnum < 0) {
00760       rc |= RIBBON_ERR_MISSING_PHOSPHATE;
00761       continue;
00762     }
00763    
00764     o1num = mol->find_atom_in_residue("O1P", res);   //  and an oxygen
00765     if (o1num < 0)
00766       o1num = mol->find_atom_in_residue("OP1", res); //  and an oxygen
00767 
00768     o2num = mol->find_atom_in_residue("O2P", res);   //  and an oxygen
00769     if (o2num < 0)
00770       o2num = mol->find_atom_in_residue("OP2", res); //  and an oxygen
00771 
00772     // if we failed to find these on the terminal residue, try the next one..
00773     if (o1num  < 0 || o2num < 0) {
00774       int nextres = (*mol->nfragList[frag])[1];
00775       o1num = mol->find_atom_in_residue("O1P", nextres);   //  and an oxygen
00776       if (o1num < 0) 
00777         o1num = mol->find_atom_in_residue("OP1", nextres); //  and an oxygen
00778 
00779       o2num = mol->find_atom_in_residue("O2P", nextres);   //  and an oxygen
00780       if (o2num < 0)
00781         o2num = mol->find_atom_in_residue("OP2", nextres); //  and an oxygen
00782 
00783       if (o1num  < 0 || o2num < 0) {
00784         rc |= RIBBON_ERR_MISSING_O1P_O2P;
00785         continue;
00786       }
00787     }
00788 
00789     ppos = framepos + 3L*cpnum;
00790     vec_add(opos, framepos + 3L*o1num, framepos + 3L*o2num);  // along the bisector
00791     vec_copy(last_opos, opos);
00792     last_ppos = ppos;
00793 
00794     vec_copy(coords-6, ppos);
00795     vec_copy(coords-3, ppos);
00796     idx[-2] = idx[-1] = -1;
00797 
00798     // now go through and set the coordinates and the perps
00799     e[0] = e[1] = e[2] = 0.0;
00800     vec_copy(g, e);
00801     int abortfrag=0;
00802     for (loop=0; (loop<num && abortfrag==0); loop++) {
00803       res = (*mol->nfragList[frag])[loop];
00804 
00805       // get atoms to use as the spline control points
00806       if (use_carb) {
00807         // use the carbons
00808         cpnum = mol->find_atom_in_residue("C3'", res);
00809         if (cpnum < 0)
00810           cpnum = mol->find_atom_in_residue("C3*", res);
00811         // use the phosphates if no carbon control points found
00812         if (cpnum < 0)
00813           cpnum = mol->find_atom_in_residue("P", res);
00814       } else { 
00815         // use the phosphates 
00816         cpnum = mol->find_atom_in_residue("P", res);
00817       }
00818 
00819       // if no P found, check the terminal atom names
00820       if (cpnum < 0) {
00821         cpnum = mol->find_atom_in_residue("H5T", res);
00822       }
00823       if (cpnum < 0) {
00824         cpnum = mol->find_atom_in_residue("H3T", res);
00825       }
00826 
00827       // cpnum must be set to a valid atom or we'll crash
00828       if (cpnum >= 0) {
00829         ppos = framepos + 3L*cpnum;
00830         idx[loop] = cpnum;
00831       } else {
00832         rc |= RIBBON_ERR_MISSING_PHOSPHATE;
00833         abortfrag = 1;
00834         break;
00835       }
00836         
00837       o1num = mol->find_atom_in_residue("O1P", res);   //  and an oxygen
00838       if (o1num < 0)
00839         o1num = mol->find_atom_in_residue("OP1", res); //  and an oxygen
00840 
00841       o2num = mol->find_atom_in_residue("O2P", res);   //  and an oxygen
00842       if (o2num < 0)
00843         o2num = mol->find_atom_in_residue("OP2", res); //  and an oxygen
00844 
00845       if (o1num < 0 || o2num < 0) {
00846         rc |= RIBBON_ERR_PROTEIN_MESS;
00847         vec_copy(opos, last_opos);
00848       } else {
00849         float tmp[3];
00850         vec_sub(tmp, framepos + 3L*o1num, ppos);
00851         vec_sub(opos, framepos + 3L*o2num, ppos);
00852         vec_add(opos, tmp, opos);  // along the bisector
00853       }
00854       vec_copy(coords+loop*3, ppos);
00855  
00856       // now I need to figure out where the ribbon goes
00857       vec_sub(a, ppos, last_ppos);      // A=(pos(P(res)) - pos(P(res-1)))
00858 //    vec_sub(b, last_opos, last_ppos); // B=(pos(Obisector(res-1)) - pos(P(res-1)))
00859 //    cross_prod(c, a, b); // C=AxB defines res-1's peptide plane
00860       vec_copy(c, opos);   // already have the normal to the ribbon
00861       cross_prod(d, c, a); // D=CxA normal to plane and backbone
00862 
00863       // if normal is > 90 deg from  previous one, invert the new one
00864       if (dot_prod(d, g) < 0) { 
00865         vec_negate(b, d);
00866       } else {
00867         vec_copy(b, d);
00868       }
00869       vec_add(e, g, b);   // average to the sum of the previous vectors
00870       vec_normalize(e);
00871       vec_scale(&perps[3L*loop], ribbon_width, e); // compute perps from normal
00872       vec_copy(g, e);     // make a cumulative sum; cute, IMHO
00873       last_ppos = ppos;
00874       vec_copy(last_opos, opos);
00875     }
00876 
00877     // abort drawing the entire fragment if unfixable problems occur
00878     if (abortfrag)
00879       continue;
00880 
00881     // and set the final points to the last element
00882     vec_copy(coords+3L*num, ppos);
00883     vec_copy(coords+3L*(num+1), ppos);
00884     idx[num] = idx[num+1] = -1;
00885 
00886     // copy the second perp to the first since the first one didn't have
00887     // enough data to construct a good perpendicular.
00888     vec_copy(perps, perps+3);
00889 
00890     // now set the first and last perps correctly
00891     vec_copy(perps-3, perps);
00892     vec_copy(perps-6, perps);
00893     vec_copy(perps+3L*num, perps+3L*num-3);
00894     vec_copy(perps+3L*num+3, perps+3L*num);
00895 
00896     // draw the nucleic acid fragment ribbon
00897     if (use_new) {
00898       float widths = 7.0f * ribbon_width * b_rad;
00899       float heights = b_rad; 
00900       draw_spline_new(num, coords, perps, idx, &widths, &heights, 1, b_res, 0);
00901     } else {
00902       draw_spline_ribbon(num, coords, perps, idx, use_cyl, b_rad, b_res);
00903     }
00904   }
00905 
00906   if (real_coords) {
00907     free(real_coords);
00908     free(real_idx);
00909     free(real_perps);
00910     real_coords = NULL;
00911     real_idx = NULL;
00912     real_perps = NULL;
00913   }
00914 
00915   return rc;
00916 }
00917 
00918 
00919 
00920 
00921 int DrawMolItem::draw_base_sugar_rings(float *framepos, int b_res, float b_rad,
00922                               float ribbon_width, int use_cyl) {
00923   float *real_coords = NULL;
00924   int frag, num;
00925   int rc=0;
00926 
00928   // ribose 
00929   //   O4',C1',C2',C3'C4' or
00930   //   O4*,C1*,C2*,C3*C4*
00931   // purines (CYT,THY,URA) 
00932   //    N1,C2,N3,C4,C5,C6
00933   // pyrimidines (ADE,GUA) 
00934   //      N1,C2,N3,C4,C5,C6,N7,N8,N9
00935   // O4',C1' and C4' are define the ribose plane and 
00936   // C2' and C3' then define the pucker of the ring
00937   // sugar -- base bonds
00938   //    pyrimidine      C1' to N9 
00939   //    purine          C1' to N1  
00940   float *o4ppos=NULL, *c1ppos=NULL, *c2ppos=NULL; 
00941   float *c3ppos=NULL, *c4ppos=NULL;
00942   int o4pnum, c1pnum, c2pnum, c3pnum, c4pnum;
00943   float *n1pos,*c2pos,*n3pos,*c4pos,*c5pos,*c6pos,*n7pos,*c8pos,*n9pos;
00944   int n1num,c2num,n3num,c4num,c5num,c6num,n7num,c8num,n9num;
00945   float rescentra[3], rescentrb[3];
00946   float midptc1pc4p[3];
00947       
00948   for (frag=0; frag<mol->nfragList.num(); frag++) {
00949     int loop;
00950 
00951     num = mol->nfragList[frag]->num(); // num of residues in this fragment
00952     if (real_coords) {
00953       free(real_coords);
00954       real_coords = NULL;
00955     }
00956 
00957     // 5atoms for the ribose but only 4 triangles
00958     // 9atoms max for a base
00959     real_coords = (float *) malloc( 14L * num * sizeof(float)*3L);
00960          
00961     // okay, I've got space for the coordinates now go
00962     for (loop=0; loop<num; loop++) {
00963       // the furanose
00964       int res = (*mol->nfragList[frag])[loop];
00965 
00966       c1pnum = mol->find_atom_in_residue("C1'", res);
00967       if (c1pnum < 0) {
00968         c1pnum = mol->find_atom_in_residue("C1*", res);
00969         if (c1pnum < 0) {
00970           rc |= RIBBON_ERR_BADNUCLEIC;
00971           continue;
00972         }
00973       } 
00974       c1ppos = framepos + 3L*c1pnum;
00975 
00976       if (atomSel->on[c1pnum]) { // switch drawing by C1' atom
00977         o4pnum = mol->find_atom_in_residue("O4'", res);
00978         if (o4pnum < 0) {
00979           o4pnum = mol->find_atom_in_residue("O4*", res);
00980         }
00981         c2pnum = mol->find_atom_in_residue("C2'", res); 
00982         if (c2pnum < 0) {
00983           c2pnum = mol->find_atom_in_residue("C2*", res);
00984         }
00985         c3pnum = mol->find_atom_in_residue("C3'", res); 
00986         if (c3pnum < 0) {
00987           c3pnum = mol->find_atom_in_residue("C3*", res);
00988         }
00989         c4pnum = mol->find_atom_in_residue("C4'", res);
00990         if (c4pnum < 0) {
00991           c4pnum = mol->find_atom_in_residue("C4*", res);
00992         }
00993 
00994         if (o4pnum < 0 || c2pnum < 0 || c3pnum < 0 || c4pnum < 0) {
00995           rc |= RIBBON_ERR_BADNUCLEIC;
00996           continue;
00997         }
00998 
00999         o4ppos = framepos + 3L*o4pnum;
01000         c2ppos = framepos + 3L*c2pnum;
01001         c3ppos = framepos + 3L*c3pnum;
01002         c4ppos = framepos + 3L*c4pnum;
01003          
01004         midpoint(midptc1pc4p, c1ppos, c4ppos);
01005          
01006         // now display triangles 
01007         cmdColorIndex.putdata(atomColor->color[c1pnum], cmdList);
01008         
01009         cmdTriangle.putdata(c4ppos,c1ppos,o4ppos,cmdList);
01010         cmdTriangle.putdata(c3ppos,midptc1pc4p,c4ppos,cmdList);
01011         cmdTriangle.putdata(c2ppos,midptc1pc4p,c3ppos,cmdList);
01012         cmdTriangle.putdata(c1ppos,midptc1pc4p,c2ppos,cmdList);
01013       } 
01014 
01015       // begin bases
01016       rescentra[0]=rescentra[1]=rescentra[2]=0.0;
01017       rescentrb[0]=rescentrb[1]=rescentrb[2]=0.0;
01018  
01019       // check for purine and pyrimidine specific atoms
01020       n9num = mol->find_atom_in_residue("N9", res);     
01021       n1num = mol->find_atom_in_residue("N1", res);
01022 
01023       // if there is a N9, then this is a pyrimidine
01024       if ((n9num >= 0) && (atomSel->on[n9num])) {
01025         c8num = mol->find_atom_in_residue("C8", res); 
01026         n7num = mol->find_atom_in_residue("N7", res); 
01027         c6num = mol->find_atom_in_residue("C6", res); 
01028         c5num = mol->find_atom_in_residue("C5", res);
01029         c4num = mol->find_atom_in_residue("C4", res);
01030         n3num = mol->find_atom_in_residue("N3", res);
01031         c2num = mol->find_atom_in_residue("C2", res);
01032         n1num = mol->find_atom_in_residue("N1", res);
01033 
01034         if (c8num < 0 || n7num < 0 || c6num < 0 || c5num < 0 ||
01035             c4num < 0 || n3num < 0 || c2num < 0 || n1num < 0) {
01036           rc |= RIBBON_ERR_BADPYRIMIDINE;
01037           continue;
01038         }
01039 
01040         n9pos = framepos + 3L*n9num;            
01041         vec_add(rescentra,rescentra,n9pos);
01042         c8pos = framepos + 3L*c8num; 
01043         vec_add(rescentra,rescentra,c8pos);
01044         n7pos = framepos + 3L*n7num; 
01045         vec_add(rescentra,rescentra,n7pos);
01046 
01047         c5pos = framepos + 3L*c5num;
01048         vec_add(rescentra,rescentra,c5pos);
01049         vec_add(rescentrb,rescentrb,c5pos);
01050         c4pos = framepos + 3L*c4num;
01051         vec_add(rescentra,rescentra,c5pos);
01052         vec_add(rescentrb,rescentrb,c5pos);
01053 
01054         c6pos = framepos + 3L*c6num; 
01055         vec_add(rescentrb,rescentrb,c6pos);
01056         n3pos = framepos + 3L*n3num;
01057         vec_add(rescentrb,rescentrb,n3pos);
01058         c2pos = framepos + 3L*c2num;
01059         vec_add(rescentrb,rescentrb,c2pos);
01060         n1pos = framepos + 3L*n1num;
01061         vec_add(rescentrb,rescentrb,n1pos);
01062                 
01063         rescentrb[0] = rescentrb[0]/6.0f;
01064         rescentrb[1] = rescentrb[1]/6.0f;
01065         rescentrb[2] = rescentrb[2]/6.0f;
01066 
01067         rescentra[0] = rescentra[0]/5.0f;
01068         rescentra[1] = rescentra[1]/5.0f;
01069         rescentra[2] = rescentra[2]/5.0f;
01070 
01071         // draw bond from ribose to base
01072         cmdCylinder.putdata(c1ppos, n9pos, b_rad, b_res, 0, cmdList);
01073         // now display triangles
01074         cmdColorIndex.putdata(atomColor->color[n9num], cmdList);
01075 
01076         cmdTriangle.putdata(n1pos,rescentrb,c2pos,cmdList);
01077         cmdTriangle.putdata(c2pos,rescentrb,n3pos,cmdList);
01078         cmdTriangle.putdata(n3pos,rescentrb,c4pos,cmdList);
01079         cmdTriangle.putdata(c4pos,rescentrb,c5pos,cmdList);
01080         cmdTriangle.putdata(c5pos,rescentrb,c6pos,cmdList);
01081         cmdTriangle.putdata(c6pos,rescentrb,n1pos,cmdList);
01082 
01083         cmdTriangle.putdata(n9pos,rescentra,c8pos,cmdList);
01084         cmdTriangle.putdata(c8pos,rescentra,n7pos,cmdList);
01085         cmdTriangle.putdata(n7pos,rescentra,c5pos,cmdList);
01086         cmdTriangle.putdata(c5pos,rescentra,c4pos,cmdList);
01087         cmdTriangle.putdata(c5pos,rescentra,c4pos,cmdList);
01088         cmdTriangle.putdata(c4pos,rescentra,n9pos,cmdList);
01089       }     
01090       else if (( n1num >= 0) && (atomSel->on[n1num])){
01091         // residue is purine and turned on
01092         c6num = mol->find_atom_in_residue("C6", res); 
01093         c5num = mol->find_atom_in_residue("C5", res);
01094         c4num = mol->find_atom_in_residue("C4", res);
01095         n3num = mol->find_atom_in_residue("N3", res);
01096         c2num = mol->find_atom_in_residue("C2", res);
01097 
01098         if (c6num < 0 || c5num < 0 || c4num < 0 || n3num < 0 || c2num < 0) {
01099           rc |= RIBBON_ERR_BADPURINE; 
01100           continue;
01101         }
01102 
01103         c6pos = framepos + 3L*c6num; 
01104         vec_add(rescentrb,rescentrb,c6pos);
01105         c5pos = framepos + 3L*c5num;
01106         vec_add(rescentrb,rescentrb,c5pos);
01107         c4pos = framepos + 3L*c4num;
01108         vec_add(rescentrb,rescentrb,c5pos);
01109         n3pos = framepos + 3L*n3num;
01110         vec_add(rescentrb,rescentrb,n3pos);
01111         c2pos = framepos + 3L*c2num;
01112         vec_add(rescentrb,rescentrb,c2pos);
01113         n1pos = framepos + 3L*n1num;
01114         vec_add(rescentrb,rescentrb,n1pos);
01115 
01116         rescentrb[0] = rescentrb[0]/6.0f;
01117         rescentrb[1] = rescentrb[1]/6.0f;
01118         rescentrb[2] = rescentrb[2]/6.0f;
01119 
01120         // draw bond from ribose to base
01121         cmdCylinder.putdata(c1ppos, n1pos, b_rad, b_res, 0, cmdList);
01122         cmdColorIndex.putdata(atomColor->color[n1num], cmdList);
01123 
01124         cmdTriangle.putdata(n1pos,rescentrb,c2pos,cmdList);
01125         cmdTriangle.putdata(c2pos,rescentrb,n3pos,cmdList);
01126         cmdTriangle.putdata(n3pos,rescentrb,c4pos,cmdList);
01127         cmdTriangle.putdata(c4pos,rescentrb,c5pos,cmdList);
01128         cmdTriangle.putdata(c5pos,rescentrb,c6pos,cmdList);
01129         cmdTriangle.putdata(c6pos,rescentrb,n1pos,cmdList);
01130       }
01131     }
01132   }
01133 
01134   return rc;
01135 }
01136 
01137 int DrawMolItem::draw_nucleotide_cylinders(float *framepos, int b_res, float b_rad, float ribbon_width, int use_cyl) {
01138   int frag, num;
01139   int rc=0;
01140   int lastcolor = -1; // only emit color changes when necessary
01141 
01142   b_rad *= 1.5; // hack to make it look a little better by default
01143 
01144   // 
01145   // XXX
01146   // Match residue names so we know which atom names to look for when we
01147   // generate nucleotide cyliners.  This approach requires no pre-analysis
01148   // of the molecule, but has the problem that performance is lost to the
01149   // redundant structure analysis work that goes on here.  Compared to the
01150   // rendering workload it is usually less significant, but it would be
01151   // better to build an acceleration structure that is updated on-demand,
01152   // e.g. when residue names have been changed by scripts or by newly 
01153   // loaded structures.  We should not have to grind through this process
01154   // again for every trajectory frame update.  This would allow us to use
01155   // more sophisticated approaches as well.
01156   //
01157 
01158   //
01159   // PDB/CHARMM residue names
01160   //
01161   int nuctypes[10];
01162   nuctypes[0] = mol->resNames.typecode((char *) "A");
01163   nuctypes[1] = mol->resNames.typecode((char *) "ADE");
01164   nuctypes[2] = mol->resNames.typecode((char *) "C");
01165   nuctypes[3] = mol->resNames.typecode((char *) "CYT");
01166   nuctypes[4] = mol->resNames.typecode((char *) "G");
01167   nuctypes[5] = mol->resNames.typecode((char *) "GUA");
01168   nuctypes[6] = mol->resNames.typecode((char *) "T");
01169   nuctypes[7] = mol->resNames.typecode((char *) "THY");
01170   nuctypes[8] = mol->resNames.typecode((char *) "U");
01171   nuctypes[9] = mol->resNames.typecode((char *) "URA");
01172 
01173   // 
01174   // AMBER uses a different nucleic residue naming convention which
01175   // is a combination one of 'D' or 'R' for DNA/RNA respectively,
01176   // followed by one of A, C, G, T, U, possibly followed by '3' or '5',
01177   // as described here: http://ffamber.cnsm.csulb.edu/
01178   // 
01179   int ambernuctypes[4L*2*3];
01180   const char *dnarnastr = "DR";
01181   const char *nucstr = "ACGTU";
01182   const char *nuctermstr = "35";
01183 
01184   int amberind=0;
01185   int foundambertypes=0;
01186   for (int nucind=0; nucind<5; nucind++) {
01187     char resstr[32];
01188     const char nuc = nucstr[nucind];  
01189     resstr[1] = nuc;
01190 
01191     for (int dnarna=0; dnarna<=1; dnarna++) {
01192       char dr = dnarnastr[dnarna];
01193       resstr[0] = dr;
01194 
01195       if (nucind == 3 && dnarna == 1)
01196         continue;
01197       if (nucind == 4 && dnarna == 0)
01198         continue;
01199 
01200       for (int termind=0; termind<3; termind++) {
01201         if (termind < 2) 
01202           resstr[2] = nuctermstr[termind];
01203         else 
01204           resstr[2] = '\0';
01205         resstr[3] = '\0';
01206 
01207         int resnameind = mol->resNames.typecode(resstr);
01208         ambernuctypes[amberind] = resnameind;
01209         if (resnameind >= 0) {
01210           foundambertypes++;
01211         }
01212         amberind++;
01213       }
01214     } 
01215   }
01216 
01217   ResizeArray<float> centers;
01218   ResizeArray<float> radii;
01219   ResizeArray<float> colors;
01220   ResizeArray<float> pickpointcoords;
01221   ResizeArray<int>   pickpointindices;
01222 
01223   // draw all fragments
01224   for (frag=0; frag<mol->nfragList.num(); frag++) {
01225     int loop;
01226 
01227     num = mol->nfragList[frag]->num(); // number of residues in this fragment
01228 
01229     // loop over all of the residues drawing nucleotide cylinders where we can
01230     for (loop=0; loop<num; loop++) {
01231       int istart = -1;  // init to invalid atom index
01232       int iend = -1;    // init to invalid atom index
01233       int res = (*mol->nfragList[frag])[loop];
01234       int myatom = mol->residueList[res]->atoms[0];
01235       int resnameindex = mol->atom(myatom)->resnameindex;
01236 
01237       if (istart < 0)
01238         istart = mol->find_atom_in_residue("C3'", res);
01239       if (istart < 0)
01240         istart = mol->find_atom_in_residue("C3*", res);
01241       if (istart < 0)
01242         istart = mol->find_atom_in_residue("C1'", res);
01243       if (istart < 0)
01244         istart = mol->find_atom_in_residue("C1*", res);
01245 
01246       // catch single nucleotides of interest
01247       if (istart < 0)
01248         istart = mol->find_atom_in_residue("P", res);
01249 
01250       // skip this nucleotide if a starting atom can't be found
01251       if (istart < 0)
01252         continue;
01253 
01254       // XXX now that we're using the C3 atoms as the backbone,
01255       // we'll only use the end atoms to select the nucleotide
01256       // cylinder, otherwise one can't use "backbone" and "not backbone"
01257       // selections to color the ribbon and nucleotide cylinders separately.  
01258 #if 0
01259       // skip this nucleotide if the starting atom is turned off
01260       if (!(atomSel->on[istart]))
01261         continue;
01262 #endif
01263 
01264       //
01265       // If nucleotide end atoms unassigned, try PDB/CHARMM residue names
01266       //
01267       if (iend < 0) {
01268         // assign end atoms using PDB/CHARMM residue names
01269         if (resnameindex == nuctypes[0] || resnameindex == nuctypes[1]) {
01270           // ADE 
01271           int n1num = mol->find_atom_in_residue("N1", res);
01272           if (n1num < 0 || !(atomSel->on[n1num]))
01273             continue;
01274           iend = n1num;
01275         } else if (resnameindex == nuctypes[2] || resnameindex == nuctypes[3]) {
01276           // CYT
01277           int n3num = mol->find_atom_in_residue("N3", res);
01278           if (n3num < 0 || !(atomSel->on[n3num]))
01279             continue;
01280           iend = n3num;
01281         } else if (resnameindex == nuctypes[4] || resnameindex == nuctypes[5]) {
01282           // GUA
01283           int n1num = mol->find_atom_in_residue("N1", res);
01284           if (n1num < 0 || !(atomSel->on[n1num]))
01285             continue;
01286           iend = n1num;
01287         } else if (resnameindex == nuctypes[6] || resnameindex == nuctypes[7] ||
01288                    resnameindex == nuctypes[8] || resnameindex == nuctypes[9]) {
01289           // THY or URA
01290           int o4num = mol->find_atom_in_residue("O4", res);
01291           if (o4num < 0 || !(atomSel->on[o4num]))
01292             continue;
01293           iend = o4num;
01294         } 
01295       }
01296 
01297       //
01298       // If nucleotide end atoms unassigned, try AMBER residue names
01299       //
01300       if (iend < 0 && foundambertypes) {
01301         // assign end atoms using AMBER residue names
01302         if (resnameindex == ambernuctypes[0] ||
01303             resnameindex == ambernuctypes[1] ||
01304             resnameindex == ambernuctypes[2] ||
01305             resnameindex == ambernuctypes[3] ||
01306             resnameindex == ambernuctypes[4] ||
01307             resnameindex == ambernuctypes[5]) {
01308           // ADE 
01309           int n1num = mol->find_atom_in_residue("N1", res);
01310           if (n1num < 0 || !(atomSel->on[n1num]))
01311             continue;
01312           iend = n1num;
01313         } else if (resnameindex == ambernuctypes[6] ||
01314                    resnameindex == ambernuctypes[7] ||
01315                    resnameindex == ambernuctypes[8] ||
01316                    resnameindex == ambernuctypes[9] ||
01317                    resnameindex == ambernuctypes[10] ||
01318                    resnameindex == ambernuctypes[11]) {
01319           // CYT
01320           int n3num = mol->find_atom_in_residue("N3", res);
01321           if (n3num < 0 || !(atomSel->on[n3num]))
01322             continue;
01323           iend = n3num;
01324         } else if (resnameindex == ambernuctypes[12] ||
01325                    resnameindex == ambernuctypes[13] ||
01326                    resnameindex == ambernuctypes[14] ||
01327                    resnameindex == ambernuctypes[15] ||
01328                    resnameindex == ambernuctypes[16] ||
01329                    resnameindex == ambernuctypes[17]) {
01330           // GUA
01331           int n1num = mol->find_atom_in_residue("N1", res);
01332           if (n1num < 0 || !(atomSel->on[n1num]))
01333             continue;
01334           iend = n1num;
01335         } else if (resnameindex == ambernuctypes[18] ||
01336                    resnameindex == ambernuctypes[19] ||
01337                    resnameindex == ambernuctypes[20] ||
01338                    resnameindex == ambernuctypes[21] ||
01339                    resnameindex == ambernuctypes[22] ||
01340                    resnameindex == ambernuctypes[23]) {
01341           // THY or URA
01342           int o4num = mol->find_atom_in_residue("O4", res);
01343           if (o4num < 0 || !(atomSel->on[o4num]))
01344             continue;
01345           iend = o4num;
01346         }
01347       }
01348 
01349       //
01350       // If we get to this point and no endpoint was defined, we try
01351       // to catch single nucleotides of interest.
01352       // This is a desperation move to identify modified bases that
01353       // use unusual residue names that can't be recognized by VMD.
01354       //
01355       if (iend < 0) {
01356         int o4num = mol->find_atom_in_residue("O4", res);
01357         if (o4num < 0 || !(atomSel->on[o4num]))
01358           continue;
01359         iend = o4num;
01360       }
01361  
01362       // add pick points for nucleotides 
01363       int pidx = 3L * istart;
01364       pickpointcoords.append3(&framepos[pidx]);
01365 
01366       pidx = 3L * iend;
01367       pickpointcoords.append3(&framepos[pidx]);
01368 
01369       pickpointindices.append(istart);
01370       pickpointindices.append(iend);
01371 
01372       // only emit color command if it has changed
01373       if (lastcolor != atomColor->color[istart]) {
01374         lastcolor = atomColor->color[istart]; 
01375         cmdColorIndex.putdata(lastcolor, cmdList);
01376       }
01377 
01378       // get coordinates of start and end atom
01379       float *cstart = framepos + 3L*istart;
01380       float *cend = framepos + 3L*iend;
01381 
01382       // draw the nucleotide cylinder
01383       cmdCylinder.putdata(cstart, cend, b_rad, b_res, 0, cmdList); 
01384 
01385       // cap ends with spheres
01386       const float *cp = scene->color_value(lastcolor);
01387 
01388 #if 0
01389       // only draw the ribbon-end sphere the starting atom is turned off
01390       // (thus no ribbon is being drawn..)  If the ribbon is being drawn,
01391       // we'll leave that end of the cylinder totally open.
01392       if (!(atomSel->on[istart])) {
01393         centers.append3(&cstart[0]);
01394 
01395         radii.append(b_rad);
01396 
01397         colors.append3(&cp[0]);
01398       }
01399 #endif
01400 
01401       centers.append3(&cend[0]);
01402       radii.append(b_rad);
01403       colors.append3(&cp[0]);
01404     } 
01405   } 
01406 
01407   // draw the spheres if we have any
01408   if (radii.num() > 0) {
01409     cmdSphereArray.putdata((float *) &centers[0],
01410                            (float *) &radii[0],
01411                            (float *) &colors[0],
01412                            radii.num(),
01413                            b_res,
01414                            cmdList);
01415   }
01416 
01417   // draw the pickpoints if we have any
01418   if (pickpointindices.num() > 0) {
01419     DispCmdPickPointArray pickPointArray;
01420     pickPointArray.putdata(pickpointindices.num(), &pickpointindices[0],
01421                            &pickpointcoords[0], cmdList);
01422   }
01423 
01424   return rc;
01425 }
01426 
01427 
01428 // draw ribbons along the protein backbone
01429 void DrawMolItem::draw_ribbons_new(float *framepos, float brad, int bres, int use_bspline, float thickness) {
01430   sprintf(commentBuffer,"Mol[%d] Rep[%d] NewRibbons", mol->id(), repNumber);
01431   cmdCommentX.putdata(commentBuffer, cmdList);
01432 
01433   // find out if I'm using lines or cylinders
01434   int use_cyl = TRUE; // use cylinders by default
01435   int rc = 0;  // clear errors
01436   
01437   if (bres <= 2 || brad < 0.01f) {
01438     use_cyl = FALSE; // the representation will be as lines
01439   } 
01440 
01441   if (use_bspline) {
01442     create_Bspline_basis(spline_basis);
01443   }
01444 
01445   append(DMATERIALON); // enable shading 
01446 
01447   float ribbon_width = thickness;
01448   if (ribbon_width < 1.0) 
01449     ribbon_width = 1.0;
01450   ribbon_width /= 7.0;
01451 
01452   rc |= draw_protein_ribbons_new(framepos, bres, brad, ribbon_width, use_cyl);
01453   rc |= draw_nucleic_ribbons(framepos, bres, brad, ribbon_width, use_cyl, 1, 1);
01454   rc |= draw_base_sugar_rings(framepos, bres, brad, ribbon_width, use_cyl);
01455 
01456   if (rc != 0) {
01457     if (emitstructwarning())
01458       msgErr << "Warning: ribbons code encountered an unusual structure, geometry may not look as expected." << sendmsg;
01459   }
01460 
01461   // set the spline basis back to the default CR so that if we change rep
01462   // styles the other reps don't get messed up.  Ahh, the joys of one big
01463   // monolithic class...
01464   if (use_bspline) {
01465     create_modified_CR_spline_basis(spline_basis, 1.25f);
01466   }
01467 }
01468 
01469 
01470 
01471 //
01472 // The new ribbons code more clearly separates the calculation of
01473 // the ribbon spline itself and the construction of the extrusion along
01474 // the spline.  Future implementations may be able to perform the 
01475 // extrusion work on graphics accelerator hardware rather than on the host
01476 // CPU.
01477 //
01478 void DrawMolItem::draw_spline_new(int num, const float *coords,
01479                                   const float *offset, const int *idx, 
01480                                   const float *cpwidths, const float *cpheights,
01481                                   int cpscalefactors, int b_res, int cyclic) {
01482   float q[4][3] = { {0, 0, 0},
01483                     {0, 0, 0},
01484                     {0, 0, 0},
01485                     {0, 0, 0},
01486                    }; // spline Q matrix
01487   float *pts = NULL;
01488   float *prps = NULL;
01489   int   *cols = NULL;
01490   float *interpwidths = NULL;
01491   float *interpheights = NULL;
01492   const float *widths = NULL;
01493   const float *heights = NULL;
01494   int i, j;
01495   int state, spanlen, numscalefactors;
01496   ResizeArray<float> pickpointcoords;
01497   ResizeArray<int> pickpointindices;
01498 
01499   // number of spline interpolations between control points
01500   int splinedivs = b_res + 2; 
01501   float invsplinedivs = 1.0f / ((float) splinedivs);
01502   int hdivs = splinedivs / 2;
01503 
01504   // allocations for interpolated spline points, perps, and color indices 
01505   // there can be a max of splinedivs * (num + 1) sections rendered, so
01506   // we allocate the max we'll ever need and go from there.
01507   pts  = (float *) malloc(sizeof(float) * splinedivs * 3L * (num + 1L));
01508   prps = (float *) malloc(sizeof(float) * splinedivs * 3L * (num + 1L));
01509   cols = (int *)   malloc(sizeof(int)   * splinedivs * (num + 1L));
01510 
01511   // allocate memory for per-section interpolated width/height values
01512   if (cpscalefactors == num) {
01513     numscalefactors = splinedivs * (num + 1);
01514     interpwidths  = (float *) malloc(sizeof(float) * numscalefactors);
01515     interpheights = (float *) malloc(sizeof(float) * numscalefactors);
01516     widths = interpwidths;
01517     heights = interpheights;
01518   } else if (cpscalefactors == 1) {
01519     numscalefactors = 1;
01520     widths  = cpwidths;
01521     heights = cpheights;
01522   } else {
01523     return; // XXX error, this should never happen
01524   }
01525 
01526   if ((pts != NULL) && (prps != NULL) && (cols != NULL)) { 
01527     // place pick points for CA atoms
01528     for (i=0; i<num; i++) {
01529       if (atomSel->on[idx[i]]) {
01530         pickpointcoords.append3(&coords[3L * i]);
01531         pickpointindices.append(idx[i]);
01532       }
01533     }
01534 
01535     // draw spans of connected 'on' segments 
01536     state=0;     // no 'on' segments yet
01537     spanlen = 0; // reset span length, no segs yet (# of splinedivs) 
01538     for (i=0; i<num; i++) {
01539       if (atomSel->on[idx[i]]) {
01540         int cindex;
01541         int iminus1min, iminus2min, iplus1max;
01542         if (!cyclic) {
01543           // duplicate start/end control points for non-cyclic fragments
01544           iminus1min = ((i - 1) >= 0) ? (i - 1) : 0;
01545           iminus2min = ((i - 2) >= 0) ? (i - 2) : 0;
01546           iplus1max  = (i + 1 < num) ? (i + 1) : (num-1);
01547         } else {
01548           // cyclic structures use modular control point indexing
01549           iminus1min = (num + i - 1) % num;
01550           iminus2min = (num + i - 2) % num;
01551           iplus1max  = (num + i + 1) % num;
01552         }
01553 
01554         // reset and start a new span
01555         if (state == 0) {
01556           state = 1;   // a new span is in progress
01557           spanlen = 0; // reset span length, no segs yet (# of splinedivs) 
01558         }
01559 
01560         // construct the interpolation points (into the "pts" array)
01561         // remember, these are offset by two to keep some information
01562         // about the previous residue in the first 2 elements
01563         make_spline_Q_matrix_noncontig(q, spline_basis, 
01564                           &coords[iminus2min * 3],
01565                           &coords[iminus1min * 3],
01566                           &coords[i          * 3],
01567                           &coords[iplus1max  * 3]);
01568 
01569         // calculated interpolated spline points
01570         for (j=0; j<splinedivs; j++) {    
01571           make_spline_interpolation(&pts[(spanlen + j) * 3], 
01572                                     j * invsplinedivs, q);
01573         }
01574 
01575         // range-clamp control point index since the first few are copies
01576         int cpind = iminus1min * 3;
01577 
01578         // interpolate perps
01579         for (j=0; j<splinedivs; j++) {  
01580           int ind = (spanlen + j)*3;
01581           float v = j * invsplinedivs;
01582           float vv = (1.0f - v); 
01583           prps[ind    ] = vv * offset[cpind    ] + v * offset[cpind + 3];
01584           prps[ind + 1] = vv * offset[cpind + 1] + v * offset[cpind + 4];
01585           prps[ind + 2] = vv * offset[cpind + 2] + v * offset[cpind + 5];
01586           vec_normalize(&prps[ind]);
01587         }
01588 
01589         // interpolate width/height values
01590         if (cpscalefactors == num) {
01591           float wminus1 = cpwidths[iminus1min];
01592           float w = cpwidths[i];
01593           float hminus1 = cpheights[iminus1min];
01594           float h = cpheights[i];
01595 
01596           if (w >= 0 && wminus1 >= 0) {
01597             // if both widths are positive, interpolate between them, 
01598             // otherwise take the previous width value rather than 
01599             // interpolating, used for drawing beta arrows.
01600             for (j=0; j<splinedivs; j++) {  
01601               int ind = spanlen + j;
01602               float v = j * invsplinedivs;
01603               float vv = (1.0f - v); 
01604               interpwidths[ind]  = vv * wminus1 + v * w;
01605               interpheights[ind] = vv * hminus1 + v * h;
01606             }
01607           } else {
01608             float wplus1 = cpwidths[iplus1max];
01609             if (wplus1 < 0) 
01610               wplus1 = -wplus1; // undo width negation
01611 
01612             if (w < 0) {
01613               w = -w; // undo width negation
01614 
01615               for (j=0; j<hdivs; j++) {  
01616                 int ind = spanlen + j;
01617                 float v = j * invsplinedivs;
01618                 interpwidths[ind]  = wminus1;
01619                 interpheights[ind] = (1.0f - v) * hminus1 + v * h;
01620               }
01621               for (j=hdivs; j<splinedivs; j++) {  
01622                 int ind = spanlen + j;
01623                 float v = j * invsplinedivs;
01624                 float nv = (j-hdivs) * invsplinedivs;
01625                 interpwidths[ind]  = (1.0f - nv) * w + nv * wplus1;
01626                 interpheights[ind] = (1.0f - v) * hminus1 + v * h;
01627               }
01628             } else {
01629               wminus1 = -wminus1; // undo width negation
01630 
01631               for (j=0; j<hdivs; j++) {  
01632                 int ind = spanlen + j;
01633                 float v = j * invsplinedivs;
01634                 float nv = (j + (splinedivs - hdivs)) * invsplinedivs;
01635                 interpwidths[ind]  = (1.0f - nv) * wminus1 + nv * w;
01636                 interpheights[ind] = (1.0f - v) * hminus1 + v * h;
01637               }
01638               for (j=hdivs; j<splinedivs; j++) {  
01639                 int ind = spanlen + j;
01640                 float v = j * invsplinedivs;
01641                 interpwidths[ind]  = w;
01642                 interpheights[ind] = (1.0f - v) * hminus1 + v * h;
01643               }
01644             }
01645           }
01646         }
01647 
01648         // lookup atom color indices and assign to the interpolated points
01649         cindex = atomColor->color[idx[iminus1min]]; 
01650         for (j=0; j<hdivs; j++) {    
01651           cols[spanlen + j] = cindex; // set color index for each point
01652         } 
01653 
01654         cindex = atomColor->color[idx[i]]; 
01655         for (j=hdivs; j<splinedivs; j++) {    
01656           cols[spanlen + j] = cindex; // set color index for each point
01657         }
01658   
01659         spanlen += splinedivs; // number of segments rendered (cp * splinediv)
01660       } else {
01661         // render a finished span
01662         if (state == 1) {
01663           state = 0; // span is finished, ready to start a new one
01664 
01665           // draw last section to connect with other reps...
01666           int iminus1min, iminus2min, iplus1max;
01667           if (!cyclic) {
01668             // duplicate start/end control points for non-cyclic fragments
01669             iminus1min = ((i - 1) >= 0) ? (i - 1) : 0;
01670             iminus2min = ((i - 2) >= 0) ? (i - 2) : 0;
01671             iplus1max  = (i + 1 < num) ? (i + 1) : (num-1);
01672           } else {
01673             // cyclic structures use modular control point indexing
01674             iminus1min = (num + i - 1) % num;
01675             iminus2min = (num + i - 2) % num;
01676             iplus1max  = (num + i + 1) % num;
01677           }
01678 
01679           int ind = spanlen * 3;
01680           make_spline_interpolation(&pts[ind], 1.0, q);
01681           int cpind = iminus2min * 3;
01682           prps[ind    ] = offset[cpind + 3]; 
01683           prps[ind + 1] = offset[cpind + 4]; 
01684           prps[ind + 2] = offset[cpind + 5]; 
01685           vec_normalize(&prps[ind]);
01686 
01687           if (cpscalefactors == num) {
01688             float w = cpwidths[i];
01689             float wminus1 = cpwidths[iminus1min];
01690             // float h = cpheights[i];
01691             float hminus1 = cpheights[iminus1min];
01692 
01693             if (w >= 0 && wminus1 >= 0) {
01694               interpwidths[spanlen] = wminus1;
01695               interpheights[spanlen] = hminus1; 
01696             } else {
01697               float wplus1 = cpwidths[iplus1max];
01698               if (wplus1 < 0) 
01699                 wplus1 = -wplus1; // undo width negation
01700 
01701               if (w < 0) {
01702                 w = -w; // undo width negation
01703 
01704                 float nv = (splinedivs - hdivs) * invsplinedivs;
01705                 interpwidths[spanlen]  = (1.0f - nv) * w + nv * wplus1;
01706                 interpheights[spanlen] = wminus1;
01707               } else {
01708                 if (w < 0) 
01709                   w = -w; // undo width negation
01710                 wminus1 = -wminus1; // undo width negation
01711                 interpwidths[spanlen]  = w;
01712                 interpheights[spanlen] = hminus1;
01713               }
01714             }
01715           }
01716 
01717           cols[spanlen] = atomColor->color[idx[iminus1min]];
01718           spanlen++;
01719 
01720           // draw a ribbon with the spline data
01721           int scalenum = (cpscalefactors == num) ? spanlen : 1;
01722           draw_ribbon_from_points(spanlen, pts, prps, cols, 
01723                                   b_res, widths, heights, scalenum);
01724         }
01725       }
01726     } 
01727 
01728     if (state == 1) {
01729       // draw last section to connect with other reps...
01730       int iminus1min, iminus2min, iplus1max;
01731       if (!cyclic) {
01732         // duplicate start/end control points for non-cyclic fragments
01733         iminus1min = ((i - 1) >= 0) ? (i - 1) : 0;
01734         iminus2min = ((i - 2) >= 0) ? (i - 2) : 0;
01735         iplus1max  = (i + 1 < num) ? (i + 1) : (num-1);
01736       } else {
01737         // cyclic structures use modular control point indexing
01738         iminus1min = (num + i - 1) % num;
01739         iminus2min = (num + i - 2) % num;
01740         iplus1max  = (num + i + 1) % num;
01741       }
01742 
01743       int ind = spanlen * 3;
01744       make_spline_interpolation(&pts[ind], 1.0, q);
01745       int cpind = iminus2min * 3;
01746       prps[ind    ] = offset[cpind + 3]; 
01747       prps[ind + 1] = offset[cpind + 4]; 
01748       prps[ind + 2] = offset[cpind + 5]; 
01749       vec_normalize(&prps[ind]);
01750 
01751       if (cpscalefactors == num) {
01752         float w = cpwidths[iminus1min];
01753         float wminus1 = cpwidths[iminus1min];
01754         // float h = cpheights[iminus1min];
01755         float hminus1 = cpheights[iminus1min];
01756 
01757         if (w >= 0 && wminus1 >= 0) {
01758           interpwidths[spanlen] = wminus1;
01759           interpheights[spanlen] = hminus1; 
01760         } else {
01761           float wplus1 = cpwidths[iplus1max];
01762           if (wplus1 < 0) 
01763             wplus1 = -wplus1; // undo width negation
01764 
01765           if (w < 0) {
01766             w = -w; // undo width negation
01767 
01768             float nv = (splinedivs - hdivs) * invsplinedivs;
01769             float nvv = (1.0f - nv); 
01770             interpwidths[spanlen]  = nvv * w + nv * wplus1;
01771             interpheights[spanlen] = wminus1;
01772           } else {
01773             if (w < 0) 
01774               w = -w; // undo width negation
01775             wminus1 = -wminus1; // undo width negation
01776             interpwidths[spanlen]  = w;
01777             interpheights[spanlen] = hminus1;
01778           }
01779         }
01780       }
01781 
01782       cols[spanlen] = atomColor->color[idx[iminus1min]];
01783       spanlen++;
01784 
01785       // draw a ribbon with the spline data
01786       int scalenum = (cpscalefactors == num) ? spanlen : 1;
01787       draw_ribbon_from_points(spanlen, pts, prps, cols, 
01788                               b_res, widths, heights, scalenum);
01789     }
01790   }
01791 
01792   if (pts != NULL)
01793     free(pts);
01794 
01795   if (prps != NULL) 
01796     free(prps);
01797 
01798   if (cols != NULL)
01799     free(cols);
01800 
01801   // deallocate memory for per-section interpolated width/height values
01802   if (interpwidths != NULL)
01803     free(interpwidths);
01804 
01805   if (interpheights != NULL)
01806     free(interpheights);
01807 
01808   // draw the pickpoints if we have any
01809   if (pickpointindices.num() > 0) {
01810     DispCmdPickPointArray pickPointArray;
01811     pickPointArray.putdata(pickpointindices.num(), &pickpointindices[0],
01812                            &pickpointcoords[0], cmdList);
01813   }
01814 }
01815 
01816 
01817 // 
01818 // Routine to render a ribbon or tube based on an array of points that
01819 // define the backbone spline on which an extrusion is built.  The 
01820 // cross-section of the extrusion can be any oval shape described by
01821 // height/width parameters that are used to scale a circular cross 
01822 // section.
01823 //
01824 void DrawMolItem::draw_ribbon_from_points(int numpoints, const float *points, 
01825                   const float *perps, const int *cols, int numpanels, 
01826                   const float *widths, const float *heights, 
01827                   int numscalefactors) {
01828   int numverts, numsections;
01829   int point, section, panel, panelindex;
01830   int index; // array index of the point coordinate we're working with
01831   float *vertexarray, *normalarray, *colorarray;
01832   float *panelshape, *panelverts, *panelnorms;
01833   float curdir[3], perpdir[3], updir[3];
01834   float width, height, lastwidth, lastheight;
01835 
01836   // we must have at least two points and two panels per section
01837   if (numpoints < 2 || numpanels < 2)
01838     return; // skip drawing if basic requirements are not met
01839  
01840   numsections = numpoints - 1;
01841   numverts = numpoints * numpanels;  
01842 
01843   // storage for final vertex array data
01844   vertexarray = (float *) malloc(numverts * 3L * sizeof(float));
01845   normalarray = (float *) malloc(numverts * 3L * sizeof(float));
01846    colorarray = (float *) malloc(numverts * 3L * sizeof(float));
01847 
01848   // storage for panel cross-section data
01849   panelverts = (float *) malloc(numpanels * 2L * sizeof(float));
01850   panelnorms = (float *) malloc(numpanels * 2L * sizeof(float));
01851   panelshape = (float *) malloc(numpanels * 2L * sizeof(float));
01852 
01853   // bail out if any of the memory allocations fail
01854   if (!vertexarray || !normalarray || !colorarray ||
01855       !panelverts  || !panelnorms  || !panelshape) {
01856     if (vertexarray) free(vertexarray);
01857     if (normalarray) free(normalarray);
01858     if (colorarray)  free(colorarray);
01859     if (panelverts)  free(panelverts);
01860     if (panelnorms)  free(panelnorms);
01861     if (panelshape)  free(panelshape);
01862     return;
01863   }
01864 
01865   // Pre-calculate cross-section template shape.  The current template shape 
01866   // is a unit circle, but it could be any shape that can be scaled efficiently
01867   // by the code that generates section offsets and normals.
01868   float fracangle = ((float) VMD_TWOPI) / ((float) numpanels);
01869   for (panel=0; panel<numpanels; panel++) {
01870     int pidx = panel * 2;
01871     float radangle = panel * fracangle;
01872     sincosf(radangle, &panelshape[pidx], &panelshape[pidx + 1]);
01873   }
01874 
01875   // initialize with invalid values so we recalculate the section 
01876   // vertex offsets and normals the first time through.
01877   lastwidth  = -999;
01878   lastheight = -999;
01879 
01880   // setup first width/height values
01881   width  =  widths[0];
01882   height = heights[0];
01883 
01884   // calculate each cross section from the spline points and "perp" vectors
01885   // that are given to us by caller, combined with the vertex offsets and
01886   // normals we've already calculated.
01887   for (point=0, index=0; point<numpoints; point++, index+=3) {
01888 
01889     // If we are provided with width/height values for every section, we
01890     // have to update the cross-sections each time the values change.
01891     // Otherwise we just re-use the same cross-section vertex offsets and
01892     // normals over and over.
01893     if (numscalefactors == numpoints) {
01894       width  = widths[point];
01895       height = heights[point];
01896     }
01897 
01898     // Calculate cross-section vertex offsets and normals. For Cartoon 
01899     // reps these have to be recalculated for each section, but for
01900     // a plain ribbon rep, the width/height are constant and the data can
01901     // be re-used over and over again for all sections.   We test to see
01902     // if width or height have changed, and only recalculate when necessary.
01903     if (width != lastwidth || height != lastheight) {
01904       float invwidth, invheight;
01905 
01906       lastwidth = width;
01907       lastheight = width;
01908 
01909       // calculate inverse scaling coefficients which are applied to normals
01910       invwidth = 1.0f / width;
01911       invheight = 1.0f / height;
01912 
01913 #if 1
01914       // Standard elliptical cross-section created by scaling a unit circle
01915       // and applying inverse scaling to the normals.  Simple to calculate,
01916       // generates decent looking results.
01917       for (panel=0; panel<numpanels; panel++) {
01918         int pidx = panel * 2;
01919         int pidy = pidx + 1;
01920         float xn, yn, invlen;
01921 
01922         // calculate vertices for a cross-section with a given width/height
01923         panelverts[pidx] = width  * panelshape[pidx];
01924         panelverts[pidy] = height * panelshape[pidy];
01925 
01926         // calculate normals for a cross-section with a given width/height
01927         xn = invwidth  * panelshape[pidx];
01928         yn = invheight * panelshape[pidy];
01929         invlen = 1.0f / sqrtf(xn*xn + yn*yn);
01930         panelnorms[pidx] = xn * invlen;
01931         panelnorms[pidy] = yn * invlen;
01932       }
01933 #else
01934       // Alternative cross-section shape built from a completely flat ribbon 
01935       // capped by two semi-circles at the edges.  Looks similar to the
01936       // elliptical ribbon, but has a flat face, so vertices are concentrated
01937       // at the edges, giving better looking shading results for the same number
01938       // of drawn vertices in many cases.  This might be worth trying out 
01939       // in place of the elliptical implementation, though it requires a 
01940       // few more more calculations.
01941   
01942       // calculate vertices for semi-circle cap on the right hand side
01943       for (panel=0; panel<(numpanels/2); panel++) {
01944         int pidx = panel * 2;
01945         int pidy = pidx + 1;
01946 
01947         // calculate vertices for a cross-section with a given width/height
01948         panelverts[pidx] = (width  / 2.0f) + panelshape[pidx];
01949         panelverts[pidy] = height * panelshape[pidy];
01950       }
01951 
01952       // calculate vertices for semi-circle cap on the left hand side
01953       for (panel=(numpanels/2); panel<numpanels; panel++) {
01954         int pidx = panel * 2;
01955         int pidy = pidx + 1;
01956 
01957         // calculate vertices for a cross-section with a given width/height
01958         panelverts[pidx] = (-width  / 2.0f) + panelshape[pidx];
01959         panelverts[pidy] = height * panelshape[pidy];
01960       }
01961 
01962       // calculate normals for both sides
01963       for (panel=0; panel<numpanels; panel++) {
01964         int pidx = panel * 2;
01965         int pidy = pidx + 1;
01966         float xn, yn, invlen;
01967     
01968         // calculate normals for a cross-section with a given width/height
01969         xn = invwidth  * panelshape[pidx];
01970         yn = invheight * panelshape[pidy];
01971         invlen = 1.0f / sqrtf(xn*xn + yn*yn);
01972         panelnorms[pidx] = xn * invlen;
01973         panelnorms[pidy] = yn * invlen;
01974       }
01975 #endif
01976     } // end of vertex offset and normal calculations
01977 
01978 
01979     // calculate and normalize the direction vector, but avoid going out
01980     // of bounds on the last point/index we calculate.
01981     if (point != (numpoints - 1)) {
01982       vec_sub(curdir, &points[index], &points[index + 3]); 
01983     } else {
01984       // handle the endpoint case where we re-use previous direction
01985       // this code requires that numpoints > 1.
01986       vec_sub(curdir, &points[index-3], &points[index]); 
01987     }
01988     vec_normalize(curdir);                 
01989 
01990     // calculate "up" from direction and perp vectors
01991     vec_copy(perpdir, &perps[index]);    // copy perpdir
01992     cross_prod(updir, curdir, perpdir);  // find section "up" direction   
01993     vec_normalize(updir);                // normalize the up direction
01994 
01995     // create vertices for the cross-section
01996     for (panel=0; panel<numpanels; panel++) {
01997       panelindex = (point * numpanels + panel) * 3;
01998       int pidx = panel * 2;
01999       float xv = panelverts[pidx    ];
02000       float yv = panelverts[pidx + 1];
02001 
02002       // use the unnormalized vertex direction vectors to construct the 
02003       // vertex position by offsetting from the original spline point.
02004       vertexarray[panelindex    ] = 
02005         (xv * perpdir[0]) + (yv * updir[0]) + points[index    ]; 
02006       vertexarray[panelindex + 1] = 
02007         (xv * perpdir[1]) + (yv * updir[1]) + points[index + 1]; 
02008       vertexarray[panelindex + 2] = 
02009         (xv * perpdir[2]) + (yv * updir[2]) + points[index + 2]; 
02010     }
02011 
02012     // create the normals for the cross-section
02013     for (panel=0; panel<numpanels; panel++) {
02014       panelindex = (point * numpanels + panel) * 3;
02015       int pidx = panel * 2;
02016       float xn = panelnorms[pidx    ];
02017       float yn = panelnorms[pidx + 1];
02018 
02019       // basis vectors are all pre-normalized so no re-normalization of the
02020       // resulting calculations needs to be done here inside this loop.
02021       normalarray[panelindex    ] = (xn * perpdir[0]) + (yn * updir[0]); 
02022       normalarray[panelindex + 1] = (xn * perpdir[1]) + (yn * updir[1]); 
02023       normalarray[panelindex + 2] = (xn * perpdir[2]) + (yn * updir[2]); 
02024     }
02025 
02026     // assign colors to the vertices in this cross-section
02027     const float *fp = scene->color_value(cols[point]);
02028     for (panel=0; panel<numpanels; panel++) {
02029       panelindex = ((point * numpanels) + panel) * 3;
02030       colorarray[panelindex    ] = fp[0]; // Red
02031       colorarray[panelindex + 1] = fp[1]; // Green
02032       colorarray[panelindex + 2] = fp[2]; // Blue
02033     }
02034   } 
02035 
02036   // generate facet lists from the vertex, normal, and color arrays
02037   // using the triangle strip primitive in VMD.
02038   int numstripverts = numsections * ((numpanels * 2) + 2);
02039   unsigned int * faces = new unsigned int[numstripverts];
02040   int * vertsperstrip = new int[numsections];
02041   int numstrips = numsections;
02042 
02043   int l=0;  
02044   for (section=0; section<numsections; section++) {
02045     vertsperstrip[section] = (numpanels * 2) + 2;
02046 
02047     int panel;
02048     for (panel=0; panel<numpanels; panel++) {
02049       // create a 2 triangles for each panel 
02050       int index = ((section * numpanels) + panel);
02051       faces[l    ] = index + numpanels;
02052       faces[l + 1] = index;
02053       l+=2;
02054     }
02055 
02056     // create a 2 triangles for each panel 
02057     int index  = section * numpanels;
02058     faces[l    ] = index + numpanels;
02059     faces[l + 1] = index;
02060 
02061     l+=2;
02062   }
02063 
02064   // Draw the ribbon!
02065   append(DMATERIALON);         // enable materials and lighting
02066 
02067   // draw triangle strips using single-sided lighting for best speed
02068   DispCmdTriStrips cmdTriStrips;  
02069   cmdTriStrips.putdata(vertexarray, normalarray, colorarray, numverts, 
02070                        vertsperstrip, numstrips, faces, numstripverts, 
02071                        0, cmdList);
02072 
02073   delete [] faces;
02074   delete [] vertsperstrip;
02075 
02076   free(vertexarray);
02077   free(normalarray);
02078   free(colorarray);
02079   free(panelshape);
02080   free(panelverts);
02081   free(panelnorms);
02082 }
02083 
02084 int DrawMolItem::draw_cartoon_ribbons(float *framepos, int b_res, float b_rad,
02085                              float ribbon_width, int use_cyl, int use_bspline) {
02086   float *coords = NULL;
02087   float *perps = NULL;
02088   int *idx = NULL;
02089   float *capos, *last_capos, *opos, *last_opos;
02090   int onum, canum, frag, num, res;
02091   float *widths, *heights;          // per-control point widths/heights
02092   const char *modulatefield = NULL; // XXX this needs to become a parameter
02093   const float *modulatedata = NULL; // data field to use for width modulation
02094   float *modulate = NULL;           // per-control point width values
02095   int rc=0;
02096 
02097   sprintf (commentBuffer,"Mol[%d] Rep[%d] NewCartoon", mol->id(), repNumber);
02098   cmdCommentX.putdata(commentBuffer, cmdList);
02099 
02100 #if defined(VMDFASTRIBBONS)
02101   wkf_timerhandle tm = wkf_timer_create();
02102   wkf_timerhandle tm2 = wkf_timer_create();
02103   wkf_timerhandle tm3 = wkf_timer_create();
02104   wkf_timerhandle tm4 = wkf_timer_create();
02105   wkf_timer_start(tm);
02106   wkf_timer_start(tm3);
02107   double foo=0.0;
02108 #endif
02109 
02110   // these are the variables used in the Raster3D package
02111   float a[3], b[3], c[3], d[3], e[3], g[3];  
02112 
02113   if (use_bspline) {
02114     create_Bspline_basis(spline_basis);
02115   }
02116 
02117   // draw nucleic acid ribbons and nucleotide cylinders
02118   rc |= draw_nucleic_ribbons(framepos, b_res, b_rad, ribbon_width / 7.0f, use_cyl, 1, 1);
02119   rc |= draw_nucleotide_cylinders(framepos, b_res, b_rad, ribbon_width / 7.0f, use_cyl);
02120 
02121   // Lookup atom typecodes ahead of time so we can use the most 
02122   // efficient variant of find_atom_in_residue() in the main loop.
02123   // If we can't find the atom types we need, bail out immediately
02124   int CAtypecode  = mol->atomNames.typecode("CA");
02125   int Otypecode   = mol->atomNames.typecode("O");
02126   int OT1typecode = mol->atomNames.typecode("OT1");
02127 
02128   if (CAtypecode < 0 || ((Otypecode < 0) && (OT1typecode < 0))) {
02129     return rc; // can't draw a ribbon without CA and O atoms for guidance
02130   }
02131 #if defined(VMDFASTRIBBONS)
02132   wkf_timer_stop(tm3);
02133   msgInfo << "Cartoon nucleotide time: " << wkf_timer_time(tm3) << sendmsg;
02134 #endif
02135 
02136   // Indicate that I need secondary structure information
02137   mol->need_secondary_structure(1);  // calculate 2ndary structure if need be
02138 
02139 #if defined(VMDFASTRIBBONS)
02140   wkf_timer_start(tm2);
02141 #endif
02142 
02143   // allocate for the maximum possible per-residue control points, perps, 
02144   // and indices so we don't have to reallocate for every fragment
02145   coords = (float *) malloc((mol->nResidues) * sizeof(float)*3);
02146      idx =   (int *) malloc((mol->nResidues) * sizeof(int));
02147    perps = (float *) malloc((mol->nResidues) * sizeof(float)*3);
02148 
02149   // cross section widths and heights
02150   widths = (float *) malloc(mol->nResidues * sizeof(float));
02151  heights = (float *) malloc(mol->nResidues * sizeof(float));
02152 
02153 #if 1
02154   // XXX hack to let users try various stuff
02155   modulatefield = getenv("VMDMODULATENEWCARTOON");
02156 #endif
02157   if (modulatefield != NULL) {
02158     if (!strcmp(modulatefield, "user")) {
02159       modulatedata = mol->current()->user;
02160       // XXX user field can be NULL on some timesteps
02161     } else {
02162       modulatedata = mol->extraflt.data(modulatefield);
02163     }
02164 
02165     // allocate for the maximum possible per-residue modulation values
02166     // so we don't have to reallocate for every fragment
02167     // XXX the modulate array is allocated cleared to zeros so that
02168     // in the case we get a NULL modulatedata pointer (user field for
02169     // example, we'll just end up with no modulation.
02170     modulate = (float *) calloc(1, mol->nResidues * sizeof(float));
02171   }
02172 
02173   // go through each protein and find the CA and O
02174   // from that, construct the offsets to use (called "perps")
02175   int numpfrag = mol->pfragList.num();
02176   for (frag=0; frag<numpfrag; frag++) {
02177     int cyclic=mol->pfragCyclic[frag];  // whether fragment is cyclic
02178     num = mol->pfragList[frag]->num();  // number of residues in this fragment
02179     if (num < 2) {
02180       rc |= RIBBON_ERR_NOTENOUGH;
02181       continue; // can't draw a ribbon with only one element, so skip
02182     }
02183 
02184     // check that we have a valid structure before continuing
02185       res = (*mol->pfragList[frag])[0];
02186     canum = mol->find_atom_in_residue(CAtypecode, res);
02187      onum = mol->find_atom_in_residue(Otypecode, res);
02188 
02189     if (canum < 0 || onum < 0) {
02190       continue; // can't find 1st CA or O of the protein, so don't draw
02191     }
02192 
02193     // initialize last 4 CA indices for use by helix rendering code
02194     int ca2, ca3, ca4;
02195     ca2=ca3=ca4=canum;
02196 
02197     int starthelix=-1;
02198     int endhelix=-1;    
02199     float starthelixperp[3];
02200 
02201     // initialize CA and O atom pointers
02202     capos = framepos + 3L*canum;
02203     last_capos = capos;
02204     opos = framepos + 3L*onum;
02205     last_opos = opos;
02206       
02207     // now go through and set the coordinates and the perps
02208     e[0] = e[1] = e[2] = 0.0;
02209     vec_copy(g, e);
02210 
02211     // for a cyclic structure, we use the positions of the last residue
02212     // to seed the initial direction vectors for the ribbon
02213     if (cyclic) {
02214       int lastres = (*mol->pfragList[frag])[num-1];
02215       int lastcanum = mol->find_atom_in_residue(CAtypecode, lastres);
02216       last_capos = framepos + 3L*lastcanum;
02217 
02218       int lastonum = mol->find_atom_in_residue(Otypecode, lastres);
02219       if (lastonum < 0 && OT1typecode >= 0) {
02220         lastonum = mol->find_atom_in_residue(OT1typecode, lastres);
02221       }
02222       last_opos = framepos + 3L*lastonum;
02223 
02224       // now I need to figure out where the ribbon goes
02225       vec_sub(a, capos, last_capos);     // A=(pos(CA(res)) - pos(CA(res-1)))
02226       vec_sub(b, last_opos, last_capos); // B=(pos(O(res-1)) - pos(CA(res-1)))
02227       cross_prod(c, a, b);               // C=AxB, define res-1's peptide plane
02228       cross_prod(d, c, a);               // D=CxA, normal to plane and backbone
02229       // if normal is > 90 deg from  previous one, negate the new one
02230       if (dot_prod(d, g) < 0) {
02231         vec_negate(b, d);
02232       } else {
02233         vec_copy(b, d);
02234       }
02235       vec_add(e, g, b);            // average to the sum of the previous vectors
02236       vec_normalize(e);
02237       vec_copy(g, e);              // make a cumulative sum; cute, IMHO
02238     }
02239 
02240     int loop;
02241     for (loop=0; loop<num; loop++) {
02242       res = (*mol->pfragList[frag])[loop];
02243       const int ss = mol->residue(res)->sstruct;
02244       float helixpos[3]; // storage for modified control point position
02245 
02246       int newcanum = mol->find_atom_in_residue(CAtypecode, res);
02247       if (newcanum >= 0) {
02248         ca4 = ca3;
02249         ca3 = ca2;
02250         ca2 = canum;
02251         canum = newcanum; 
02252         capos = framepos + 3L*canum;
02253       }
02254 
02255       onum = mol->find_atom_in_residue(Otypecode, res);
02256       if (onum < 0 && OT1typecode >= 0) {
02257         onum = mol->find_atom_in_residue(OT1typecode, res);
02258       }
02259       if (onum >= 0) {
02260         opos = framepos + 3L*onum;
02261       } else {
02262         rc |= RIBBON_ERR_PROTEIN_MESS;
02263         opos = last_opos; // try to cope if we have no oxygen index
02264       }
02265 
02266       // modulate the ribbon width by user-specified per-atom data
02267       float m_fac = 1.0f;
02268       if (modulatedata != NULL) {
02269         modulate[loop] = modulatedata[canum];
02270 
02271         // only allow modulation values > 0. otherwise fall back on old scheme.
02272         m_fac = modulate[loop];
02273         if (m_fac <= 0.0f)
02274           m_fac = 1.0f;
02275       }
02276 
02277       idx[loop] = canum;
02278 
02279       // calculate ribbon cross section 
02280       int drawhelixwithrods = 0;
02281       int drawbetawithribbons = 0;
02282 #if 0
02283       drawhelixwithrods = (getenv("VMDHELIXRODS") != NULL);
02284       drawbetawithribbons = (getenv("VMDBETARIBBONS") != NULL);
02285 #endif
02286 
02287       switch (ss) {
02288         case SS_HELIX_ALPHA:
02289         case SS_HELIX_3_10:
02290         case SS_HELIX_PI:
02291           if (drawhelixwithrods) { 
02292             // if we just entered the helix, find the end of the helix so 
02293             // we can generate usable perps by interpolation.
02294             if (starthelix == -1) {
02295               int helind;
02296               starthelix = loop;
02297               for (helind=loop; helind<num; helind++) {
02298                 int hres = (*mol->pfragList[frag])[helind];
02299                 const int hss = mol->residue(hres)->sstruct;
02300                 if (hss == SS_HELIX_ALPHA ||
02301                     hss == SS_HELIX_3_10 || 
02302                     hss == SS_HELIX_PI) {
02303                   endhelix = helind;
02304                 } else {
02305                   break; // stop when this helix ends
02306                 }
02307               }
02308 
02309               // save starting perp and ending perp, for interpolation later
02310             }
02311 
02312             widths[loop] = 4 * b_rad * m_fac;
02313             heights[loop] = 4 * b_rad * m_fac;
02314             vec_copy(helixpos, framepos + 3L * canum); 
02315             vec_add(helixpos, helixpos, framepos + 3L * ca2); 
02316             vec_add(helixpos, helixpos, framepos + 3L * ca3); 
02317             vec_add(helixpos, helixpos, framepos + 3L * ca4); 
02318             vec_scale(helixpos, 0.25, helixpos);
02319 
02320             // copy the helix coordinate into the control point array
02321             vec_copy(coords+loop*3, helixpos);
02322             capos = helixpos; // perps are calculated from modified position 
02323           } else {
02324             widths[loop] = ribbon_width * b_rad * m_fac;
02325             heights[loop] = b_rad * m_fac;
02326             // copy the CA coordinate into the control point array
02327             vec_copy(coords+loop*3, capos);
02328           }
02329           break;
02330 
02331         case SS_TURN:
02332         case SS_COIL:
02333         case SS_BRIDGE:
02334         default:
02335           widths[loop] = b_rad * m_fac;
02336           heights[loop] = b_rad * m_fac;
02337           // copy the CA coordinate into the control point array
02338           vec_copy(coords+loop*3, capos);
02339           break;
02340 
02341         case SS_BETA:
02342           widths[loop] = ribbon_width * b_rad * m_fac;
02343           heights[loop] = b_rad * m_fac;
02344 
02345           if (drawbetawithribbons) {
02346             // copy the CA coordinate into the control point array
02347             vec_copy(coords+loop*3, capos);
02348           } else {
02349             float betapos[3]; 
02350             int drawarrowhead = 0;
02351 
02352             // filter the control point positions, giving equal weight
02353             // to the current CA atom and it's two neighbors, so they
02354             // ideally cancel out the wiggles in the beta sheet, without
02355             // having to switch spline basis etc.
02356             int caplus1 = -1; // make invalid initially
02357 
02358             if ((loop+1) < num) {
02359               int nextres = (*mol->pfragList[frag])[loop+1];
02360               caplus1 = mol->find_atom_in_residue(CAtypecode, nextres);
02361 
02362               // draw directionality arrow if we're at the end
02363               if (mol->residue(nextres)->sstruct != SS_BETA) 
02364                 drawarrowhead = 1;
02365             } else {
02366               // draw directionality arrow if we're at the end
02367               drawarrowhead = 1;
02368             }
02369 
02370             // mark this width to avoid interpolating it, and make
02371             // the arrowhead wider than the beta sheet body
02372             // XXX non-interpolated control points are negated
02373             if (drawarrowhead) 
02374               widths[loop] = -ribbon_width * b_rad * 1.75f;
02375 
02376             if (caplus1 < 0)
02377               caplus1 = canum;
02378 
02379             vec_copy(betapos, framepos + 3L * canum); 
02380             vec_scale(betapos, 2.0f, betapos);
02381             vec_add(betapos, betapos, framepos + 3L * caplus1); 
02382             vec_add(betapos, betapos, framepos + 3L * ca2); 
02383             vec_scale(betapos, 0.25f, betapos);
02384 
02385             // copy the beta coordinate into the control point array
02386             vec_copy(coords+loop*3, betapos);
02387 
02388             // perps will be generated using the original CA position
02389           }
02390           break;
02391       }
02392 
02393 
02394       // now I need to figure out where the ribbon goes
02395       vec_sub(a, capos, last_capos);     // A=(pos(CA(res)) - pos(CA(res-1)))
02396       vec_sub(b, last_opos, last_capos); // B=(pos(O(res-1)) - pos(CA(res-1)))
02397       cross_prod(c, a, b);               // C=AxB, define res-1's peptide plane
02398       cross_prod(d, c, a);               // D=CxA, normal to plane and backbone
02399 
02400       // if normal is > 90 deg from  previous one, negate the new one
02401       if (dot_prod(d, g) < 0) {
02402         vec_negate(b, d);
02403       } else {
02404         vec_copy(b, d);
02405       }
02406       vec_add(e, g, b);            // average to the sum of the previous vectors
02407       vec_normalize(e);
02408       vec_copy(&perps[3L*loop], e); // compute perps from the normal
02409       vec_copy(g, e);              // make a cumulative sum; cute, IMHO
02410       last_capos = capos;
02411       last_opos = opos;
02412 
02413       if (drawhelixwithrods) {
02414         if (loop == starthelix) {
02415           // save starting helix perp
02416           vec_copy(starthelixperp, e); 
02417         } else  if (loop > starthelix && loop < endhelix) {
02418           // should interpolate rather than copy, but this is just for testing
02419           vec_copy(&perps[3L*loop], starthelixperp);
02420         } else if (loop == endhelix) {
02421           // reset helix start and end if we've reached the end
02422           starthelix = -1; 
02423           endhelix = -1; 
02424         }
02425       }
02426                      
02427     }
02428 
02429     if (!cyclic) {
02430       // copy the second perp to the first since the first one didn't have
02431       // enough data to construct a good perpendicular.
02432       vec_copy(perps, perps+3);
02433     }
02434 
02435 #if defined(VMDFASTRIBBONS)
02436     wkf_timer_start(tm4);
02437 #endif
02438     // draw the cartoon 
02439     draw_spline_new(num, coords, perps, idx, widths, heights, num, b_res, cyclic);
02440 #if defined(VMDFASTRIBBONS)
02441     wkf_timer_stop(tm4);
02442     foo+=wkf_timer_time(tm4);
02443 #endif
02444   } // drew the protein fragment ribbon
02445 
02446   if (coords) {
02447     free(coords);
02448     free(idx);
02449     free(perps);
02450     free(widths);
02451     free(heights);
02452   }
02453 
02454   if (modulate) {
02455     free(modulate);
02456   }
02457 
02458   // XXX set the spline basis back to the default CR so that if we change rep
02459   // styles the other reps don't get messed up.  Ahh, the joys of one big
02460   // monolithic class...
02461   if (use_bspline) {
02462     create_modified_CR_spline_basis(spline_basis, 1.25f);
02463   }
02464 
02465 #if defined(VMDFASTRIBBONS)
02466   wkf_timer_stop(tm2);
02467   msgInfo << "Cartoon spline time: " << wkf_timer_time(tm2) << sendmsg;
02468   msgInfo << "     subspline time: " << foo << sendmsg;
02469 #endif
02470 #if defined(VMDFASTRIBBONS)
02471   wkf_timer_stop(tm);
02472   msgInfo << "Cartoon regen time: " << wkf_timer_time(tm) << sendmsg;
02473   wkf_timer_destroy(tm);
02474   wkf_timer_destroy(tm2);
02475   wkf_timer_destroy(tm3);
02476 #endif
02477 
02478   return rc;
02479 }
02480 

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