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

DrawRingsUtils.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: DrawRingsUtils.C,v $
00013  *  $Author: johns $  $Locker:  $    $State: Exp $
00014  *  $Revision: 1.42 $  $Date: 2019/10/28 21:02:44 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *
00019  * Ulities for calculating ring axes, ring puckering and displacement of
00020  * atoms from the mean ring plane.
00021  *
00022  ***************************************************************************/
00023 
00024 #include <math.h>
00025 #include "utilities.h"
00026 #include "AtomColor.h"
00027 #include "Scene.h"
00028 #include "DrawRingsUtils.h"
00029 
00030 /*
00031  * Hill-Reilly pucker-paramter based coloring
00032  */
00033 void lerp_color_range(float *newcol, float val, float rmin, float rmax,
00034                       float *mincol, float *maxcol) {
00035   float range = rmax-rmin;
00036   float lerpval = (val - rmin) / range;
00037   vec_lerp(newcol, mincol, maxcol, lerpval);
00038 }
00039 
00040 
00041 void hotcold_gradient_lerp(float pucker_sum, float *rgb) {
00042   vec_zero(rgb); // set default color to black
00043 
00044   // hot to cold color map
00045   // Red (1, 0, 0) -> Yellow (1, 1, 0) -> Green (0, 1, 0) -> Cyan (0, 1, 1) -> blue (0, 0, 1)
00046   float     red[3] = {1.0f, 0.0f, 0.0f};
00047   float  yellow[3] = {1.0f, 1.0f, 0.0f};
00048 //  float yellow2[3] = {0.8f, 1.0f, 0.0f};
00049   float   green[3] = {0.0f, 1.0f, 0.0f};
00050   float  green2[3] = {0.6f, 1.0f, 0.0f};
00051   float    cyan[3] = {0.0f, 1.0f, 1.0f};
00052   float   cyan2[3] = {0.0f, 1.0f, 0.8f};
00053   float    blue[3] = {0.0f, 0.0f, 1.0f};
00054 
00055   if (pucker_sum < 0.25f) {
00056     lerp_color_range(rgb, pucker_sum, 0.00f, 0.25f, red, yellow);
00057   } else if (pucker_sum < 0.45f) {
00058     vec_copy(rgb, yellow);
00059   } else if (pucker_sum < 0.55f) {
00060     lerp_color_range(rgb, pucker_sum, 0.45f, 0.55f, yellow, green2);
00061   } else if (pucker_sum < 0.75f) {
00062     lerp_color_range(rgb, pucker_sum, 0.55f, 0.75f, green, cyan2);
00063   } else {
00064     lerp_color_range(rgb, pucker_sum, 0.75f, 1.00f, cyan, blue);
00065   }
00066 
00067   clamp_color(rgb); // clamp color values to legal range
00068 }
00069 
00070 
00071 void hotcold_gradient(float pucker_sum, float *rgb) {
00072   vec_zero(rgb); // set default color to black
00073 
00074   // hot to cold color map
00075   // Red (1, 0, 0) -> Yellow (1, 1, 0) -> Green (0, 1, 0) -> Cyan (0, 1, 1) -> blue (0, 0, 1) -> magenta (1, 0, 1)
00076   if (pucker_sum < 0.40f) {  //MK - envelopes here
00077     rgb[0] = 1.0f;  // red
00078     rgb[1] = pucker_sum * 2.5f;  // MK from red increasing green -> yellow -  adjusted multiplier for large range
00079     rgb[2] = 0.0f;
00080   } else if (pucker_sum < 0.56f) {
00081     rgb[0] = 1.0f - (pucker_sum - 0.40f) * 6.25f; // from Yellow, decrease red -> green adjusted multiplier for small range
00082     rgb[1] = 1.0f;
00083     rgb[2] = 0.0f;
00084   } else if (pucker_sum < 0.64f) {
00085     rgb[0] = 0.0f;
00086     rgb[1] = 1.0f; //green
00087     rgb[2] = (pucker_sum - 0.56f) * 12.5f; // from green, increasing blue ->  cyan,  adjusted multiplier for small range
00088   } else if (pucker_sum < 0.76f) {
00089     rgb[0] = 0.0f;
00090     rgb[1] = 1.0f - (pucker_sum - 0.64f) * 5.0f; // from cyan, decrease green -> blue, adjusted multiplier for small range
00091     rgb[2] = 1.0f;
00092   } else {
00093     rgb[0] = (pucker_sum - 0.76f) * 0.8f; // from blue, increase red to get magenta, adjusted multiplier for very large range
00094     rgb[1] = 0.0f;
00095     rgb[2] = 1.0f;
00096   }
00097 
00098   clamp_color(rgb); // clamp color values to legal range
00099 }
00100 
00101 
00102 float hill_reilly_ring_pucker(SmallRing &ring, float *framepos) {
00103   int N = ring.num(); // the number of atoms in the current ring
00104 
00105 #if 0
00106   // return the default color if this isn't a 5 or 6 ring atom
00107   if (N != 5 && N != 6)
00108     return 0.0;
00109     //MK added
00110     if (N==6) {
00111       //MK do Hill-Reilly for 6-membered rings
00112       int NP = N-3; // number of puckering parameters
00113       float *X = new float[N*3L]; // atom co-ordinates
00114       float *r = new float[N*3L]; // bond vectors
00115       float *a = new float[NP*3L]; // puckering axes
00116       float *q = new float[NP*3L]; // normalized puckering vectors
00117       float *n = new float[3L]; // normal to reference plane
00118       float *p = new float[3L]; // a flap normal
00119       float *theta = new float[NP]; // puckering parameters
00120       float pucker_sum;
00121       float max_pucker_sum;
00122       float *atompos;
00123       int curatomid, i, j, k, l;
00124     
00125       // load ring co-ordinates
00126       for (i=0; i<N; i++) {
00127         curatomid = ring[i];
00128         atompos = framepos + 3L*curatomid;
00129         X[3L*i  ] = atompos[0];
00130         X[3L*i+1] = atompos[1];
00131         X[3L*i+2] = atompos[2];
00132       }     
00133     
00134       // calculate bond vectors
00135       for (i=0; i<N; i++) {
00136         j = (i+1) % N;
00137         vec_sub(r+3L*i, X+3L*j, X+3L*i);
00138       }
00139     
00140       // calculate puckering axes, flap normals and puckering vectors
00141       for (i=0; i<NP; i++) {
00142         k = (L2*(i+1)) % N;
00143         j = (L2*i) % N;
00144         l = (L2*i+1) % N;
00145         vec_sub(a+3L*i, X+3L*k, X+3L*j);
00146         cross_prod(p, r+3L*j, r+3L*l);
00147         cross_prod(q+3L*i, a+3L*i, p);
00148         vec_normalize(q+3L*i);
00149       }
00150     
00151       // reference normal
00152       cross_prod(n, a+3L*0, a+3L*1);
00153       vec_normalize(n);
00154     
00155       // calculate the puckering parameters
00156       pucker_sum = 0.0;
00157     
00158       for (i=0; i<NP; i++) {
00159         theta[i] = (float(VMD_PI)/2.0f) - acosf(dot_prod(q+3L*i, n));
00160         pucker_sum += theta[i];
00161       }
00162     
00163     
00164       // 0.6154 radians (35.26 degrees) has significance for perfect tetrahedral bond geometry (see Hill paper)
00165       max_pucker_sum = NP * 0.6154f;
00166       float pucker_scaled = pucker_sum/max_pucker_sum;
00167       pucker_sum = fabsf((pucker_scaled < 1.0f) ? pucker_scaled : 1.0f);
00168       pucker_sum = (pucker_sum < 1.0f) ? pucker_sum : 1.0f;
00169     
00170       delete [] X;
00171       delete [] r;
00172       delete [] a;
00173       delete [] q;
00174       delete [] n;
00175       delete [] p;
00176       delete [] theta;
00177       return pucker_sum;
00178     }  //end MK if N==6
00179     else {  //N==5 
00180 #endif
00181     float *xring = new float[N];
00182     float *yring = new float[N];
00183     float *zring = new float[N];
00184     float *displ = new float[N];
00185     float *q = new float[N];
00186     float *phi = new float[N];
00187     float Q;
00188     int m;
00189     float *atompos;
00190     int curatomid;
00191     
00192     for (int i=0; i<N; i++) {
00193       curatomid = ring[i];
00194       atompos = framepos + 3L*curatomid; // pointer arithmetic is evil :)
00195       xring[i] = atompos[0];
00196       yring[i] = atompos[1];
00197       zring[i] = atompos[2];
00198     }     
00199     
00200     atom_displ_from_mean_plane(xring, yring, zring, displ, N);
00201 
00202     delete [] xring;
00203     delete [] yring;
00204     delete [] zring;
00205 
00206     if (cremer_pople_params(N, displ, q, phi, m, Q)) {
00207       // Q is the puckering amplitude - i.e. the intensity of the pucker.
00208       Q = (Q < 2.0f) ? Q : 2.0f;  //truncate amplitude at 2
00209       delete [] displ;
00210       delete [] q;
00211       delete [] phi;
00212       return Q;
00213     } else {
00214       delete [] displ;
00215       delete [] q;
00216       delete [] phi;
00217       return 0.0;
00218     }    
00219 #if 0
00220   }
00221 #endif
00222 }
00223 
00224 
00225 // Calculate Hill-Reilly Pucker Parameters and convert these to a ring colour
00226 void hill_reilly_ring_color(SmallRing &ring, float *framepos, float *rgb) {
00227   float pucker_sum = hill_reilly_ring_pucker(ring, framepos); 
00228   hotcold_gradient(pucker_sum, rgb);
00229   //MK added and now removed   hotcold_gradient(pucker_sum/0.8, rgb);  //scale, assuming amplitude value not bigger than 0.8
00230 }
00231 
00232 void hill_reilly_ring_colorscale(SmallRing &ring, float *framepos, 
00233                                  float vmin, float vmax,
00234                                  const Scene *scene, float *rgb) {
00235   float pucker_sum = hill_reilly_ring_pucker(ring, framepos); 
00236 
00237   // map data min/max to range 0->1
00238   // values must be clamped before use, since user-specified
00239   // min/max can cause out-of-range color indices to be generated
00240   float vscale;
00241   float vrange = vmax - vmin;
00242   if (fabsf(vrange) < 0.00001f)
00243     vscale = 0.0f;
00244   else
00245     vscale = 1.00001f / vrange;
00246 
00247   float level = (pucker_sum - vmin) * vscale;
00248   int colindex = (int)(level * MAPCLRS-1);
00249   if (colindex < 0)
00250     colindex = 0;
00251   else if (colindex >= MAPCLRS)
00252     colindex = MAPCLRS-1;
00253  
00254   const float *scrgb = scene->color_value(MAPCOLOR(colindex));
00255   
00256   rgb[0] = scrgb[0];
00257   rgb[1] = scrgb[1];
00258   rgb[2] = scrgb[2];
00259 }
00260 
00261 
00262 
00263 /*
00264  * Cremer-Pople pucker-parameter based coloring
00265  */
00266 
00267 // return sum + (x,y,z)
00268 void vec_incr(float sum[3], const float x, const float y, const float z) {
00269   sum[0] += x;
00270   sum[1] += y;
00271   sum[2] += z;
00272 }
00273 
00274 // Calculates cartesian axes based on coords of nuclei in ring
00275 // using cremer-pople algorithm. It is assumed that the
00276 // centre of geometry is the centre of the ring.
00277 void ring_axes(const float * X, const float * Y, const float * Z, int N, 
00278              float x[3], float y[3], float z[3]) {
00279   float Rp[3] = {0.0, 0.0, 0.0}; float Rpp[3] = {0.0, 0.0, 0.0};
00280   int j;
00281   for (j=0; j<N; j++) {
00282     float ze_angle = 2.0f * float(VMD_PI) * float(j-1) / float(N);
00283     float ze_sin, ze_cos;
00284     sincosf(ze_angle, &ze_sin, &ze_cos);
00285     vec_incr(Rp,  X[j]*ze_sin, Y[j]*ze_sin, Z[j]*ze_sin);
00286     vec_incr(Rpp, X[j]*ze_cos, Y[j]*ze_cos, Z[j]*ze_cos);
00287   }
00288 
00289   cross_prod(z, Rp, Rpp);
00290   vec_normalize(z);
00291 
00292   /* 
00293    * OK, now we have z, the norm to the central plane, we need
00294    * to calculate y as the projection of Rp onto the plane
00295    * and x as y x z
00296    */
00297   float lambda = dot_prod(z, Rp);
00298   y[0] = Rp[0] - z[0]*lambda;
00299   y[1] = Rp[1] - z[1]*lambda;
00300   y[2] = Rp[2] - z[2]*lambda;
00301   vec_normalize(y);
00302   cross_prod(x, y, z);    // voila !
00303 }
00304 
00305 
00306 // Calculate distances of atoms from the mean ring plane
00307 void atom_displ_from_mean_plane(float * X, float * Y, float * Z, 
00308                                 float * displ, int N) {
00309   float cog[3] = {0.0, 0.0, 0.0};
00310   float x_axis[3], y_axis[3], z_axis[3];
00311   int i;
00312 
00313   // calculate centre of geometry
00314   for (i=0; i<N; i++) {
00315       cog[0] += X[i]; cog[1] += Y[i]; cog[2] += Z[i];
00316   }
00317   cog[0] /= float(N); 
00318   cog[1] /= float(N); 
00319   cog[2] /= float(N); 
00320 
00321   // centre the ring
00322   for (i=0; i<N; i++) {
00323     X[i] -= cog[0]; 
00324     Y[i] -= cog[1]; 
00325     Z[i] -= cog[2];
00326   }
00327 
00328   ring_axes( X, Y, Z, N, x_axis, y_axis, z_axis );
00329 
00330   // calculate displacement from mean plane
00331   for (i=0; i<N; i++) {
00332     displ[i] = X[i]*z_axis[0] + Y[i]*z_axis[1] + Z[i]*z_axis[2];
00333   }
00334 }
00335 
00336 
00337 // Calculate Cremer-Pople puckering parameters
00338 int cremer_pople_params(int N_ring_atoms, float * displ, float * q, 
00339                         float * phi, int  & m , float & Q) {
00340   int i, j, k;
00341   if (N_ring_atoms<3)
00342     return -1;
00343 
00344   float N = float(N_ring_atoms);
00345   phi[0]=0;  q[0]=0;  //no puckering parameters for m=1
00346 
00347   // if even no ring atoms, first calculate unpaired q puck parameter
00348   if (fmod(N, 2.0f)==0) {
00349     float sum =0;
00350     m = N_ring_atoms/2 -1;
00351 
00352     for (i=0; i<N_ring_atoms; i++)         
00353       sum += displ[i]*cosf(i*float(VMD_PI));
00354 
00355     q[int(N_ring_atoms/2)-1]=sqrtf(1.0f/N)*sum;
00356   } else {
00357     m = int(N_ring_atoms-1)/2;
00358   }
00359 
00360   // calculate paired puckering parameters
00361   for (i=1; i<m; i++) {
00362     float q_cosphi=0, q_sinphi=0;
00363     for (j=0; j<N_ring_atoms; j++) {
00364       float ang = 2.0f*float(VMD_PI)*float((i+1)*j)/N;
00365       float sa, ca;
00366       sincosf(ang, &sa, &ca);
00367       q_cosphi += displ[j]*ca;
00368       q_sinphi += displ[j]*sa;
00369     }
00370 
00371     q_cosphi *=  sqrtf(2.0f/N);
00372     q_sinphi *= -sqrtf(2.0f/N);
00373     phi[i]=atanf(q_sinphi/q_cosphi);
00374 
00375     if (q_cosphi < 0)
00376       phi[i]+=float(VMD_PI);
00377     else if (q_sinphi < 0) 
00378       phi[i]+=2.0f*float(VMD_PI);
00379 
00380     q[i]=q_cosphi/phi[i];
00381     //phi[i]*=180.0f/VMD_PI;  //convert to degrees
00382 
00383     //calculate puckering amplitude
00384     Q=0.0f;
00385     for (k=0; k<N_ring_atoms; k++)
00386       Q+=displ[k]*displ[k];
00387     Q=sqrtf(Q);  
00388   }
00389 
00390   return 1;
00391 }
00392 
00393 
00394 // Calculate Cremer-Pople Pucker Parameters and convert these to a ring colour
00395 void cremer_pople_ring_color(SmallRing &ring, float *framepos, float *rgb) {
00396   int N = ring.num(); //the number of atoms in the current ring
00397   float *xring = new float[N];
00398   float *yring = new float[N];
00399   float *zring = new float[N];
00400   float *displ = new float[N];
00401   float *q = new float[N];
00402   float *phi = new float[N];
00403   float Q=0;
00404   int m;
00405   float *atompos;
00406   int curatomid;
00407 
00408   vec_zero(rgb); // set default color to black
00409 
00410   for (int i=0; i<N; i++) {
00411     curatomid = ring[i];
00412     atompos = framepos + 3L*curatomid; // pointer arithmetic is evil :)
00413     xring[i] = atompos[0];
00414     yring[i] = atompos[1];
00415     zring[i] = atompos[2];
00416   }     
00417 
00418   atom_displ_from_mean_plane(xring, yring, zring, displ, N);
00419          
00420   if (N==6) { //special case - pyranose rings
00421     if (cremer_pople_params(N, displ, q, phi, m, Q)) {
00422       float cosTheta = q[2]/Q;
00423       float theta = acosf(cosTheta); 
00424       float sinTheta = sinf(theta);
00425 
00426       // Q is the puckering amplitude - i.e. the intensity of the pucker.
00427       // multiply by Q to show intensity, particularly for rings with 
00428       // little pucker (black)
00429       // NOTE -using abs - polar positions therefore equivalent
00430       float intensity = Q;
00431       
00432       rgb[0] = fabsf(sinTheta)*intensity;
00433       rgb[1] = fabsf(cosTheta)*intensity;
00434       rgb[2] = fabsf(sinf(3.0f*phi[1])*sinTheta)*intensity;
00435     }
00436   } else if (N==5) { //special case - furanose rings
00437     if (cremer_pople_params(N, displ, q, phi, m, Q)) {
00438       rgb[0] = 0;
00439       rgb[1] = 0;
00440       rgb[2] = Q;
00441     }
00442   }
00443 
00444   // clamp color values to legal range
00445   clamp_color(rgb);
00446  
00447   delete [] xring;
00448   delete [] yring;
00449   delete [] zring;
00450   delete [] displ;
00451   delete [] q;
00452   delete [] phi;
00453 }
00454 
00455 
00456 /*
00457  * Ribbon spline handling
00458  */
00459 // Calculates the position at point t along the spline with co-efficients
00460 // A, B, C and D.
00461 // spline(t) = ((A * t + B) * t + C) * t + D
00462 void ribbon_spline(float *pos, const float * const A, const float * const B,
00463                    const float * const C, const float * const D, const float t) {
00464   vec_copy(pos,D);
00465   vec_scaled_add(pos,t,C);
00466   vec_scaled_add(pos,t*t,B);
00467   vec_scaled_add(pos,t*t*t,A);
00468 }
00469 
00470 
00471 
00472 

Generated on Tue Oct 15 02:44:12 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002