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

TclMDFF.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: TclMDFF.C,v $
00013  *      $Author: johns $        $Locker:  $             $State: Exp $
00014  *      $Revision: 1.34 $      $Date: 2020/10/15 17:44:27 $
00015  *
00016  ***************************************************************************
00017  * DESCRIPTION:
00018  *   Tcl bindings for MDFF functions
00019  *
00020  ***************************************************************************/
00021 #include <stdio.h>
00022 #include <stdlib.h>
00023 #include <string.h>
00024 #include <float.h> // FLT_MAX etc
00025 
00026 #include "Inform.h"
00027 #include "utilities.h"
00028 
00029 #include "AtomSel.h"
00030 #include "VMDApp.h"
00031 #include "MoleculeList.h"
00032 #include "VolumetricData.h"
00033 #include "VolMapCreate.h" // volmap_write_dx_file()
00034 
00035 #include "CUDAMDFF.h"
00036 #include "MDFF.h"
00037 #include <math.h>
00038 #include <tcl.h>
00039 #include "TclCommands.h"
00040 #include "TclMDFF.h"
00041 #include "Voltool.h"
00042 
00043 
00044 int mdff_sim(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
00045   int verbose = 0;
00046   if (argc < 3) {
00047      // "     options: --allframes (average over all frames)\n"
00048     Tcl_SetResult(interp, (char *) "usage: voltool "
00049       "sim: <selection> [options]\n"
00050 //      "              --weight (weight density with atomic numbers)\n"
00051       "              -o <output map> \n"
00052       "              -res <target resolution in Angstroms> (default 10.0)\n"
00053       "              -spacing <grid spacing in Angstroms> (default based on resolution)\n",
00054       TCL_STATIC);
00055     return TCL_ERROR;
00056   }
00057 
00058   //atom selection
00059   AtomSel *sel = NULL;
00060   sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
00061   if (!sel) {
00062     Tcl_AppendResult(interp, "volmap: no atom selection.", NULL);
00063     return TCL_ERROR;
00064   }
00065   if (!sel->selected) {
00066     Tcl_AppendResult(interp, "volmap: no atoms selected.", NULL);
00067     return TCL_ERROR;
00068   }
00069   if (!app->molecule_valid_id(sel->molid())) {
00070     Tcl_AppendResult(interp, "invalide mol id.", NULL);
00071     return TCL_ERROR;
00072   }
00073 
00074 //  int ret_val=0;
00075   float radscale;
00076   double gspacing = 0;
00077   double resolution = 10.0;
00078 //  bool use_all_frames = false;
00079 //  bool useweight = false;
00080 //  VolumetricData *volmap = NULL;
00081   MoleculeList *mlist = app->moleculeList;
00082   Molecule *mymol = mlist->mol_from_id(sel->molid());
00083 //  const char *outputmap = Tcl_GetStringFromObj(objv[2], NULL);
00084   const char *outputmap = NULL;
00085 
00086   for (int i=0; i < argc; i++) {
00087     char *opt = Tcl_GetStringFromObj(objv[i], NULL);
00088 //    if (!strcmp(opt, "--weight")) {useweight = true;}
00089 //    if (!strcmp(opt, "--allframes")) {use_all_frames = true;}
00090     if (!strcmp(opt, "-res")) {
00091       if (i == argc-1) {
00092         Tcl_AppendResult(interp, "No resolution specified",NULL);
00093         return TCL_ERROR;
00094       } else if (Tcl_GetDoubleFromObj(interp, objv[i+1], &resolution) != TCL_OK) { 
00095         Tcl_AppendResult(interp, "\nResolution incorrectly specified",NULL);
00096         return TCL_ERROR;
00097       }
00098     }
00099     if (!strcmp(opt, "-o")) {
00100       if (i == argc-1) {
00101         Tcl_AppendResult(interp, "No output file specified",NULL);
00102         return TCL_ERROR;
00103       } else {
00104         outputmap = Tcl_GetStringFromObj(objv[i+1], NULL);
00105       }
00106     }
00107     if (!strcmp(opt, "-spacing")) {
00108       if (i == argc-1) {
00109         Tcl_AppendResult(interp, "No spacing specified",NULL);
00110         return TCL_ERROR;
00111       } else if ( Tcl_GetDoubleFromObj(interp, objv[i+1], &gspacing) != TCL_OK) {
00112         Tcl_AppendResult(interp, "\nspacing incorrectly specified",NULL);
00113         return TCL_ERROR;
00114       }
00115     }
00116 
00117     if (!strcmp(opt, "-verbose") || (getenv("VMDMDFFVERBOSE") != NULL)) {
00118       verbose = 1;
00119     }
00120   }
00121 
00122   // use quicksurf to compute simulated density map
00123   const float *framepos = sel->coordinates(app->moleculeList);
00124   const float *radii = mymol->radius();
00125   radscale = 0.2f * float(resolution);
00126 
00127   if (gspacing == 0) {
00128     gspacing = 1.5*radscale;
00129   }
00130 
00131   int quality = 0;
00132   if (resolution >= 9)
00133     quality = 0;
00134   else
00135     quality = 3;
00136 
00137   if (verbose)
00138     printf("MDFF dens: radscale %f gspacing %f\n", radscale, gspacing);
00139 
00140   int cuda_err = -1;
00141 #if defined(VMDCUDA)
00142   VolumetricData *synthvol=NULL;
00143   if (getenv("VMDNOCUDA") == NULL) {
00144     cuda_err = vmd_cuda_calc_density(sel, app->moleculeList, quality, radscale, gspacing, &synthvol, NULL, NULL, verbose);
00145     if(cuda_err != -1) {
00146       int newvolmolid = init_new_volume_molecule(app, synthvol, "sim_map");
00147       Molecule *newvolmol = mlist->mol_from_id(newvolmolid);
00148       if (outputmap != NULL) {
00149         if (!write_file(app, newvolmol, 0, outputmap)){
00150           Tcl_AppendResult(interp, "\n no target volume correctly specified",NULL);
00151       return TCL_ERROR;
00152         }
00153       }
00154     }
00155   }
00156 #endif
00157 
00158   // If CUDA failed, we use CPU fallback, and we have to prevent QuickSurf
00159   // from using the GPU either...
00160   if (cuda_err == -1) {
00161     const int force_cpu_calc=1;
00162     QuickSurf *qs = new QuickSurf(force_cpu_calc);
00163     VolumetricData *volmap = NULL;
00164     volmap = qs->calc_density_map(sel, mymol, framepos, radii,
00165                                   quality, (float)radscale, (float)gspacing);
00166    int newvolmolid = init_new_volume_molecule(app, volmap, "sim_map");
00167   Molecule *newvolmol = mlist->mol_from_id(newvolmolid);
00168   if (outputmap != NULL) {
00169     if (!write_file(app, newvolmol, 0, outputmap)){
00170       Tcl_AppendResult(interp, "\n no target volume correctly specified",NULL);
00171       return TCL_ERROR;
00172     }
00173   }
00174     delete qs;
00175   }
00176 
00177   return TCL_OK;
00178 }
00179 
00180 
00181 
00182 int mdff_cc(VMDApp *app, int argc, Tcl_Obj * const objv[], Tcl_Interp *interp) {
00183   if (argc < 3) {
00184     Tcl_SetResult(interp,
00185       (char *) "usage: voltool "
00186       "cc: <selection> -res <target resolution in A> [options]\n"
00187       "     options: --allframes (average over all frames)\n"
00188     //  "              --weight (weight simulated density with atomic numbers)\n"
00189       "              -i <input map> specifies new target density filename to load.\n"
00190       "              -mol <molid> specifies an already loaded density's molid for use as target\n"
00191       "              -thresholddensity <x> (ignores voxels with values below x threshold)\n"
00192       "              -vol <volume id> specifies an already loaded density's volume id for use as target. Defaults to 0.\n"
00193       "              -spacing <grid spacing in Angstroms> (default based on resolution)\n"
00194       "              -calcsynthmap append the simulated density map to the molecule.\n" 
00195       "              -calcdiffmap append a difference map (between the simulated and target densities) to the molecule.(CUDA only)\n" 
00196       "              -calcspatialccmap append to the molecule a spatial map of the correlations between 8x8 voxel regions of the simulated and target densities.(CUDA only)\n" 
00197       "              -savesynthmap <output file> save the simulated density to a file.\n"
00198       "              -savespatialccmap <output file> save the spatial correlaiton map to a file.(CUDA only)\n"
00199       ,
00200       TCL_STATIC);
00201     return TCL_ERROR;
00202   }
00203 
00204   //atom selection
00205   AtomSel *sel = NULL;
00206   sel = tcl_commands_get_sel(interp, Tcl_GetStringFromObj(objv[1],NULL));
00207   if (!sel) {
00208     Tcl_AppendResult(interp, "volmap: no atom selection.", NULL);
00209     return TCL_ERROR;
00210   }
00211   if (!sel->selected) {
00212     Tcl_AppendResult(interp, "volmap: no atoms selected.", NULL);
00213     return TCL_ERROR;
00214   }
00215   if (!app->molecule_valid_id(sel->molid())) {
00216     Tcl_AppendResult(interp, "invalide mol id.", NULL);
00217     return TCL_ERROR;
00218   }
00219 
00220 //  int ret_val=0;
00221 //  bool useweight = false;
00222   int verbose = 0;
00223   float radscale;
00224   int calcsynthmap = 0;
00225   const char *savesynthmap = NULL;
00226   int calcdiffmap = 0;
00227   int calcspatialccmap = 0;
00228   const char *savespatialccmap = NULL;
00229   double gspacing = 0;
00230   double thresholddensity = -FLT_MAX;
00231 
00232   int molid = -1;
00233   int volid = 0;
00234   double resolution;
00235   const char *input_map = NULL;
00236   bool use_all_frames = false;
00237 //  VolumetricData *volmap = NULL;
00238   MoleculeList *mlist = app->moleculeList;
00239   Molecule *mymol = mlist->mol_from_id(sel->molid());
00240 
00241 #if 0
00242   if ( Tcl_GetDoubleFromObj(interp, objv[2], &resolution) != TCL_OK){
00243     Tcl_AppendResult(interp, "\nResolution incorrectly specified",NULL);
00244     return TCL_ERROR;
00245   }
00246 #endif
00247 
00248   for (int i=0; i < argc; i++) {
00249     char *opt = Tcl_GetStringFromObj(objv[i], NULL);
00250 //    if (!strcmp(opt, "--weight")) {useweight = true;}
00251     if (!strcmp(opt, "--allframes")) {use_all_frames = true;}
00252     if (!strcmp(opt, "-i")) {
00253       if (i == argc-1) {
00254         Tcl_AppendResult(interp, "No input map specified",NULL);
00255         return TCL_ERROR;
00256       }
00257 
00258       FileSpec spec;
00259       spec.waitfor = FileSpec::WAIT_BACK; // shouldn't this be waiting for all?
00260       input_map = Tcl_GetStringFromObj(objv[1+i], NULL);
00261       molid = app->molecule_new(input_map,0,1);
00262       //sel->molid()
00263       int ret_val = app->molecule_load(molid, input_map,app->guess_filetype(input_map),&spec);
00264       if (ret_val < 0) return TCL_ERROR;
00265     }
00266 
00267     if (!strcmp(opt, "-res")) {
00268       if (i == argc-1){
00269         Tcl_AppendResult(interp, "No resolution specified",NULL);
00270         return TCL_ERROR;
00271       } else if (Tcl_GetDoubleFromObj(interp, objv[i+1], &resolution) != TCL_OK){ 
00272         Tcl_AppendResult(interp, "\n resolution incorrectly specified",NULL);
00273         return TCL_ERROR;
00274       }
00275     }
00276 
00277     if (!strcmp(opt, "-spacing")) {
00278       if (i == argc-1) {
00279         Tcl_AppendResult(interp, "No spacing specified",NULL);
00280         return TCL_ERROR;
00281       } else if (Tcl_GetDoubleFromObj(interp, objv[i+1], &gspacing) != TCL_OK) {
00282           Tcl_AppendResult(interp, "\n spacing incorrectly specified",NULL);
00283           return TCL_ERROR;
00284       }
00285     }
00286 
00287     if (!strcmp(opt, "-mol")) {
00288       if (i == argc-1) {
00289         Tcl_AppendResult(interp, "No molid specified",NULL);
00290         return TCL_ERROR;
00291       } else if ( Tcl_GetIntFromObj(interp, objv[i+1], &molid) != TCL_OK) {
00292         Tcl_AppendResult(interp, "\n molid incorrectly specified",NULL);
00293         return TCL_ERROR;
00294       }
00295     }
00296 
00297     if (!strcmp(opt, "-vol")) {
00298       if (i == argc-1){
00299         Tcl_AppendResult(interp, "No volume id specified",NULL);
00300         return TCL_ERROR;
00301       } else if ( Tcl_GetIntFromObj(interp, objv[i+1], &volid) != TCL_OK) {
00302         Tcl_AppendResult(interp, "\n volume id incorrectly specified",NULL);
00303         return TCL_ERROR;
00304       }
00305     }
00306 
00307     if (!strcmp(opt, "-thresholddensity")) {
00308       if (i == argc-1) {
00309         Tcl_AppendResult(interp, "No threshold specified",NULL);
00310         return TCL_ERROR;
00311       } else if (Tcl_GetDoubleFromObj(interp, objv[i+1], &thresholddensity) != TCL_OK) { 
00312         Tcl_AppendResult(interp, "\nthreshold incorrectly specified",NULL);
00313         return TCL_ERROR;
00314       }
00315     }
00316 
00317     // calculate (and/or save) simulated density map
00318     if (!strcmp(opt, "-calcsynthmap")) {
00319       calcsynthmap = 1;
00320     }
00321     if (!strcmp(opt, "-savesynthmap")) {
00322       calcsynthmap = 1;
00323       if (i == argc-1){
00324         Tcl_AppendResult(interp, "No output dx file specified", NULL);
00325         return TCL_ERROR;
00326       } else {
00327         savesynthmap = Tcl_GetStringFromObj(objv[i+1], NULL);
00328 #if 0
00329         Tcl_AppendResult(interp, "\n volume id incorrectly specified",NULL);
00330         return TCL_ERROR;
00331 #endif
00332       }
00333     }
00334 
00335 
00336     if (!strcmp(opt, "-calcdiffmap")) {
00337       calcdiffmap = 1;
00338     }
00339 
00340 
00341     // calculate (and/or save) spatialcc map
00342     if (!strcmp(opt, "-calcspatialccmap")) {
00343       calcspatialccmap = 1;
00344     }
00345     if (!strcmp(opt, "-savespatialccmap")) {
00346       calcspatialccmap = 1;
00347       if (i == argc-1){
00348         Tcl_AppendResult(interp, "No output dx file specified", NULL);
00349         return TCL_ERROR;
00350       } else {
00351         savespatialccmap = Tcl_GetStringFromObj(objv[i+1], NULL);
00352 #if 0
00353         Tcl_AppendResult(interp, "\n volume id incorrectly specified",NULL);
00354         return TCL_ERROR;
00355 #endif
00356       }
00357     }
00358 
00359 
00360     if (!strcmp(opt, "-verbose") || (getenv("VMDMDFFVERBOSE") != NULL)) {
00361       verbose = 1;
00362     }
00363   }
00364 
00365 
00366   const VolumetricData *volmapA = NULL;
00367   if (molid > -1) {
00368     Molecule *volmol = mlist->mol_from_id(molid);
00369     if (volmol == NULL) {
00370       Tcl_AppendResult(interp, "\n invalid molecule specified",NULL);
00371       return TCL_ERROR;
00372     }
00373 
00374     if (volmapA == NULL) 
00375       volmapA = volmol->get_volume_data(volid);
00376   } else {
00377     Tcl_AppendResult(interp, "\n no target volume specified",NULL);
00378     return TCL_ERROR;
00379   }
00380   if (volmapA == NULL) {
00381     Tcl_AppendResult(interp, "\n no target volume correctly specified",NULL);
00382     return TCL_ERROR;
00383   }
00384 
00385   float return_cc = 0;
00386   int start, end;
00387   if (use_all_frames) {
00388     start = 0;
00389     end = mymol->numframes()-1;
00390   } else {
00391     // start = mymol->frame();
00392     // end = mymol->frame();
00393     start = sel->which_frame;
00394     end = sel->which_frame;
00395   }
00396 
00397   PROFILE_PUSH_RANGE("MDFF cross correlation", 5);
00398 
00399   // use quicksurf density map algorithm
00400   for (int frame = start; frame <= end; frame++) {
00401     const float *framepos = new float [sel->num_atoms*3L];
00402  
00403     // frame is -2 or -1, meaning last or current but 
00404     // either way not allframes, so current sel frame is fine
00405     if (frame < 0)
00406       framepos = sel->coordinates(app->moleculeList);
00407     else
00408       framepos = (mymol->get_frame(frame))->pos;
00409     const float *radii = mymol->radius();
00410 
00411     radscale = 0.2f * float(resolution);
00412  //   if(gspacing == 0){gspacing = 1.5*radscale;}
00413     int quality = 0;
00414     if (resolution >=9){quality = 0;}
00415     else{quality = 3;}
00416 
00417     VolumetricData *synthvol=NULL;
00418     VolumetricData *diffvol=NULL;
00419 
00420 #if 0
00421     double cc=0.0;
00422     if (threshold != 0.0) {
00423       printf("Calculating threshold...\n");
00424       int size = synthvol->xsize * synthvol->ysize * synthvol->zsize;
00425       double mean = 0.;
00426       int i;
00427       for (i=0; i<size; i++)
00428         mean += synthvol->data[i];
00429       mean /= size;
00430 
00431       double sigma = 0.;
00432       for (i=0; i<size; i++)
00433         sigma += (synthvol->data[i] - mean)*(synthvol->data[i] - mean);
00434       sigma /= size;
00435       sigma = sqrt(sigma);
00436 
00437       threshold = mean + threshold*sigma;
00438       printf("Threshold: %f\n", threshold);
00439     }
00440 #endif
00441 
00442     VolumetricData *synthmap = NULL;
00443     VolumetricData *diffmap = NULL;
00444     VolumetricData *spatialccmap = NULL;
00445 
00446     if (verbose) {
00447       if (calcsynthmap)
00448         printf("MDFF calculating simulated map\n");
00449       if (calcdiffmap)
00450         printf("MDFF calculating difference map\n");
00451       if (calcspatialccmap)
00452         printf("MDFF calculating spatial CC map\n");
00453     }
00454     
00455 
00456     // 
00457     // Try computing with CUDA first
00458     //
00459     int cuda_err = -1;
00460 #if defined(VMDCUDA)
00461     VolumetricData **synthpp = NULL;
00462     VolumetricData **diffpp = NULL;
00463     VolumetricData **spatialccpp = NULL;
00464 
00465     if (calcsynthmap) {
00466       synthpp = &synthmap;
00467     }
00468 
00469     if (calcdiffmap) {
00470       diffpp = &diffmap;
00471     }
00472 
00473     if (calcspatialccmap) {
00474       spatialccpp = &spatialccmap;
00475     }
00476 
00477     if (gspacing == 0 && (getenv("VMDNOCUDA") == NULL)) {
00478       if (verbose)
00479         printf("Computing CC on GPU...\n");
00480 
00481 #if 0
00482       if (verbose) {
00483         printf("TclMDFF: prep for vmd_cuda_compare_sel_refmap():\n");
00484         printf("  refmap xaxis: %lf %lf %lf\n",
00485                volmapA->xaxis[0], volmapA->xaxis[1], volmapA->xaxis[2]);
00486         printf("  refmap size: %d %d %d\n",
00487                volmapA->xsize, volmapA->ysize, volmapA->zsize);
00488         printf("  gridspacing (orig): %f\n", gspacing);
00489       }
00490 #endif
00491 
00492       cuda_err = vmd_cuda_compare_sel_refmap(sel, app->moleculeList, quality, 
00493                                     radscale, gspacing, 
00494                                     volmapA, synthpp, diffpp, spatialccpp, 
00495                                     &return_cc, thresholddensity, verbose);
00496     }
00497 #endif
00498 
00499     // If CUDA failed, we use CPU fallback, and we have to prevent QuickSurf
00500     // from using the GPU either...
00501     if (cuda_err == -1) {
00502       const int force_cpu_calc=1;
00503       if (verbose)
00504         printf("Computing CC on CPUs...\n");
00505 
00506       if (gspacing == 0) {
00507         gspacing = 1.5*radscale;
00508       }
00509 
00510       QuickSurf *qs = new QuickSurf(force_cpu_calc);
00511       VolumetricData *volmap = NULL;
00512       volmap = qs->calc_density_map(sel, mymol, framepos, radii,
00513                                     quality, (float)radscale, (float)gspacing);
00514       double cc = 0.0;
00515 
00516 #if 0
00517       // this is 'old style' sigma threshold 
00518       if (thresholddensity != -FLT_MAX) {
00519         printf("Calculating threshold...\n"); 
00520         int size = volmap->xsize*volmap->ysize*volmap->zsize;
00521         double mean = 0.;
00522         int i;
00523         for (i=0; i<size; i++)
00524           mean += volmap->data[i];
00525         mean /= size;
00526  
00527         double sigma = 0.;
00528         for (i=0; i<size; i++)
00529           sigma += (volmap->data[i] - mean)*(volmap->data[i] - mean);
00530         sigma /= size;
00531         sigma = sqrt(sigma);
00532  
00533         thresholddensity = mean + thresholddensity*sigma;
00534         printf("Threshold: %f\n", thresholddensity);
00535       }
00536 #endif
00537   
00538       cc_threaded(volmap, volmapA, &cc, thresholddensity);
00539       return_cc += float(cc);
00540       delete qs;
00541   
00542       if (start != end) {
00543         return_cc /= mymol->numframes();
00544       }
00545       Tcl_SetObjResult(interp, Tcl_NewDoubleObj(return_cc));
00546   
00547       VolumetricData *vol = NULL;
00548       if (calcsynthmap) {
00549         vol=volmap;
00550         if (savesynthmap) {
00551           printf("Writing simulated map to '%s'\n", savesynthmap);
00552           volmap_write_dx_file(vol, savesynthmap);
00553         } else {
00554           printf("Adding simulated map to molecule\n");
00555           app->molecule_add_volumetric(molid, vol->name, vol->origin, 
00556                                        vol->xaxis, vol->yaxis, vol->zaxis,
00557                                        vol->xsize, vol->ysize, vol->zsize, 
00558                                        vol->data);
00559           vol->data = NULL; // prevent destruction of density array;
00560         }
00561       }
00562 
00563       delete volmap;
00564 
00565       PROFILE_POP_RANGE(); // first return point
00566       return TCL_OK;
00567     }
00568 
00569     VolumetricData *v = NULL;
00570     if (calcsynthmap) {
00571       v=synthmap;
00572       if (savesynthmap) {
00573         printf("Writing simulated map to '%s'\n", savesynthmap);
00574         volmap_write_dx_file(v, savesynthmap);
00575       } else {
00576         printf("Adding simulated map to molecule\n");
00577         app->molecule_add_volumetric(molid, v->name, v->origin, 
00578                                      v->xaxis, v->yaxis, v->zaxis,
00579                                      v->xsize, v->ysize, v->zsize, v->data);
00580         v->data = NULL; // prevent destruction of density array;
00581       }
00582       delete v;
00583     }
00584  
00585     if (calcdiffmap) {
00586       printf("Adding difference map to molecule\n");
00587       v=diffmap;
00588       app->molecule_add_volumetric(molid, v->name, v->origin, 
00589                                    v->xaxis, v->yaxis, v->zaxis,
00590                                    v->xsize, v->ysize, v->zsize, v->data);
00591       v->data = NULL; // prevent destruction of density array;
00592       delete v;
00593     }
00594 
00595     if (calcspatialccmap) {
00596       v=spatialccmap;
00597       if (savespatialccmap) {
00598         printf("Writing spatialcc map to '%s'\n", savespatialccmap);
00599         volmap_write_dx_file(v, savespatialccmap);
00600       } else {
00601         printf("Adding spatial CC map to molecule\n");
00602         app->molecule_add_volumetric(molid, v->name, v->origin, 
00603                                      v->xaxis, v->yaxis, v->zaxis,
00604                                      v->xsize, v->ysize, v->zsize, v->data);
00605         v->data = NULL; // prevent destruction of density array;
00606       }
00607       delete v;
00608     }
00609 
00610     delete synthvol;
00611     delete diffvol;
00612   }
00613 
00614   if (start != end)
00615     return_cc /= mymol->numframes(); 
00616 
00617   PROFILE_POP_RANGE(); // second return point
00618 
00619   Tcl_SetObjResult(interp, Tcl_NewDoubleObj(return_cc));
00620   return TCL_OK;
00621 }
00622 
00623 
00624 
00625 
00626 
00627 int obj_mdff_cc(ClientData cd, Tcl_Interp *interp, int argc,
00628                             Tcl_Obj * const objv[]){
00629   if (argc < 2) {
00630     Tcl_SetResult(interp,
00631     (char *) "Usage: mdffi <command> [args...]\n"
00632       "Commands:\n"
00633       "cc   -- calculates the cross-correlation coefficient\n"
00634       "sim  -- creates a simulated map from an atomic structure\n"
00635       ,
00636       TCL_STATIC);
00637     return TCL_ERROR;
00638   }
00639   char *argv1 = Tcl_GetStringFromObj(objv[1],NULL);
00640 
00641   VMDApp *app = (VMDApp *)cd;
00642   if (!strupncmp(argv1, "cc", CMDLEN))
00643     return mdff_cc(app, argc-1, objv+1, interp);
00644   if (!strupncmp(argv1, "sim", CMDLEN))
00645     return mdff_sim(app, argc-1, objv+1, interp);
00646 
00647   Tcl_SetResult(interp, (char *) "Type 'mdffi' for summary of usage\n", TCL_VOLATILE);
00648   return TCL_OK;
00649 }
00650 
00651 

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