VolMapCreate.C

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 * RCS INFORMATION:
00010 *
00011 * $RCSfile: VolMapCreate.C,v $
00012 * $Author: johns $ $Locker: $ $State: Exp $
00013 * $Revision: 1.127 $ $Date: 2024年03月01日 02:39:38 $
00014 *
00015 ***************************************************************************/
00016 
00017 /* Functions for creating useful volumetric maps based on the 3-D molecular structure */
00018 
00019 // Todo List: 
00020 // - Document!
00021 // - Don't just output to a DX file... give user more control (use plugins)
00022 // - Allow a framerange param, don't just use all available frames (and get rid
00023 // of the VMDApp dependency). Also, VMD needs a general FrameRange object...
00024 
00025 // Note:
00026 // All functions in here loop over x fastest and z slowest and also create
00027 // maps with that order of data. This matches the order used in
00028 // VolumetricData which is ultimately dictated by the way graphics hardware 
00029 // works.
00030 
00031 // XXX VolMapCreate provides its own dx file writer.
00032 // Actualy the maps should be generated in memory only and only dumped
00033 // to files using the molfile_plugin interface. This is currently
00034 // not yet possible but will hopefully be enabled soon.
00035 
00036 #include <math.h>
00037 #include <stdlib.h>
00038 #include <stdio.h>
00039 #include <string.h>
00040 #include <errno.h> // only for write_dx_file()
00041 #include <tcl.h>
00042 
00043 #include "VMDApp.h"
00044 #include "MoleculeList.h"
00045 #include "Molecule.h"
00046 #include "Timestep.h"
00047 #include "Measure.h"
00048 #include "SpatialSearch.h"
00049 #include "VolCPotential.h"
00050 #include "VolumetricData.h"
00051 #include "VolMapCreate.h"
00052 #include "utilities.h"
00053 #include "ResizeArray.h"
00054 #include "Inform.h"
00055 #include "WKFUtils.h"
00056 #include "TclCommands.h"
00057 
00058 #if defined(VMDUSEMSMPOT)
00059 #include "msmpot.h"
00060 #endif
00061 
00062 // Conversion factor between raw units (e^2 / A) and kT/e
00063 #define POT_CONV 560.47254
00064 
00065 #define MIN(X,Y) (((X)<(Y))? (X) : (Y))
00066 #define MAX(X,Y) (((X)>(Y))? (X) : (Y))
00067 
00070 static const float MAX_ENERGY = 150.f; 
00071 
00072 
00073 
00075 
00076 VolMapCreate::VolMapCreate(VMDApp *the_app, AtomSel *the_sel, float resolution) {
00077 volmap = NULL;
00078 app = the_app;
00079 sel = the_sel;
00080 delta = resolution;
00081 computed_frames = 0;
00082 checkpoint_freq = 0;
00083 checkpoint_name = NULL;
00084 
00085 char dataname[1];
00086 strcpy(dataname, ""); // null-terminated empty string
00087 float zerovec[3];
00088 memset(zerovec, 0, 3L*sizeof(float));
00089 volmap = new VolumetricData(dataname, zerovec,
00090 zerovec, zerovec, zerovec, 0, 0, 0, NULL);
00091 user_minmax = false;
00092 }
00093 
00094 
00095 VolMapCreate::~VolMapCreate() {
00096 if (volmap) delete volmap;
00097 if (checkpoint_name) delete[] checkpoint_name;
00098 }
00099 
00100 
00101 void VolMapCreate::set_minmax (float minx, float miny, float minz, float maxx, float maxy, float maxz) {
00102 user_minmax = true;
00103 min_coord[0] = minx;
00104 min_coord[1] = miny;
00105 min_coord[2] = minz;
00106 max_coord[0] = maxx;
00107 max_coord[1] = maxy;
00108 max_coord[2] = maxz;
00109 }
00110 
00111 
00112 void VolMapCreate::set_checkpoint (int checkpointfreq_tmp, char *checkpointname_tmp) {
00113 if (checkpointfreq_tmp > -1) checkpoint_freq = checkpointfreq_tmp;
00114 if (!checkpointname_tmp) return;
00115 
00116 if (checkpoint_name) delete[] checkpoint_name;
00117 checkpoint_name = new char[strlen(checkpointname_tmp)+1];
00118 strcpy(checkpoint_name, checkpointname_tmp);
00119 }
00120 
00121 
00124 int VolMapCreate::calculate_minmax (float *my_min_coord, float *my_max_coord) {
00125 DrawMolecule *mol = app->moleculeList->mol_from_id(sel->molid());
00126 int numframes = app->molecule_numframes(sel->molid()); // XXX need a frame selection object
00127 
00128 float frame_min_coord[3], frame_max_coord[3], *coords;
00129 
00130 msgInfo << "volmap: Computing bounding box coordinates" << sendmsg;
00131 
00132 int save_frame = sel->which_frame;
00133 int frame;
00134 for (frame=0; frame<numframes; frame++) {
00135 sel->which_frame=frame;
00136 sel->change(NULL,mol);
00137 coords = sel->coordinates(app->moleculeList);
00138 if (!coords) continue;
00139 
00140 int err = measure_minmax(sel->num_atoms, sel->on, coords, NULL, 
00141 frame_min_coord, frame_max_coord);
00142 if (err != MEASURE_NOERR) {
00143 sel->which_frame = save_frame;
00144 return err;
00145 }
00146 
00147 int i;
00148 for (i=0; i<3; i++) {
00149 if (!frame || frame_min_coord[i] < my_min_coord[i]) my_min_coord[i] = frame_min_coord[i];
00150 if (!frame || frame_max_coord[i] > my_max_coord[i]) my_max_coord[i] = frame_max_coord[i];
00151 }
00152 }
00153 sel->which_frame = save_frame;
00154 
00155 return 0;
00156 }
00157 
00158 
00160 int VolMapCreate::calculate_max_radius (float &max_rad) {
00161 DrawMolecule *mol = app->moleculeList->mol_from_id(sel->molid());
00162 if (!mol) return -1;
00163 const float *radius = mol->extraflt.data("radius");
00164 if (!radius) return MEASURE_ERR_NORADII;
00165 
00166 max_rad = 0.f;
00167 int i;
00168 for (i=sel->firstsel; i<=sel->lastsel; i++) 
00169 if (sel->on[i] && radius[i] > max_rad) max_rad = radius[i];
00170 
00171 return 0;
00172 }
00173 
00174 // Combo routines are used to combine the different frame maps together into a
00175 // final entity
00176 
00177 // Initialize the frame combination buffer
00178 void VolMapCreate::combo_begin(CombineType method, void **customptr, void *params) {
00179 int gridsize = volmap->xsize*volmap->ysize*volmap->zsize;
00180 
00181 *customptr = NULL;
00182 memset(volmap->data, 0, sizeof(float)*gridsize);
00183 computed_frames = 0;
00184 
00185 // these combine types need additional storage
00186 if (method == COMBINE_STDEV) {
00187 float *voldata2 = new float[gridsize];
00188 memset(voldata2, 0, gridsize*sizeof(float));
00189 *customptr = (void*) voldata2;
00190 }
00191 }
00192 
00193 // Add a frame to the combination buffer
00194 void VolMapCreate::combo_addframe(CombineType method, float *voldata, void *customptr, float *frame_voldata) {
00195 float *voldata2 = (float*) customptr;
00196 int gridsize = volmap->xsize*volmap->ysize*volmap->zsize;
00197 int n;
00198 
00199 computed_frames++;
00200 
00201 if (computed_frames == 1) { // FIRST FRAME
00202 switch (method) {
00203 case COMBINE_AVG:
00204 case COMBINE_MAX: 
00205 case COMBINE_MIN: 
00206 memcpy(voldata, frame_voldata, gridsize*sizeof(float));
00207 break;
00208 case COMBINE_PMF:
00209 memcpy(voldata, frame_voldata, gridsize*sizeof(float));
00210 break;
00211 case COMBINE_STDEV: 
00212 memcpy(voldata, frame_voldata, gridsize*sizeof(float));
00213 for (n=0; n<gridsize; n++) voldata2[n] = frame_voldata[n]*frame_voldata[n];
00214 break;
00215 }
00216 
00217 return;
00218 }
00219 
00220 // THE FOLLOWING ONLY APPLIES TO OTHER FRAMES THAN FIRST
00221 switch (method) {
00222 case COMBINE_AVG:
00223 for (n=0; n<gridsize; n++) voldata[n] += frame_voldata[n];
00224 break;
00225 case COMBINE_PMF:
00226 for (n=0; n<gridsize; n++) voldata[n] = (float) -log(exp(-voldata[n]) + exp(-frame_voldata[n]));
00227 break;
00228 case COMBINE_MAX: 
00229 for (n=0; n<gridsize; n++) voldata[n] = MAX(voldata[n], frame_voldata[n]);
00230 break;
00231 case COMBINE_MIN: 
00232 for (n=0; n<gridsize; n++) voldata[n] = MIN(voldata[n], frame_voldata[n]);
00233 break;
00234 case COMBINE_STDEV: 
00235 for (n=0; n<gridsize; n++) voldata[n] += frame_voldata[n];
00236 for (n=0; n<gridsize; n++) voldata2[n] += frame_voldata[n]*frame_voldata[n];
00237 break;
00238 }
00239 }
00240 
00241 
00246 void VolMapCreate::combo_export(CombineType method, float *voldata, void *customptr) {
00247 float *voldata2 = (float*) customptr;
00248 int gridsize = volmap->xsize*volmap->ysize*volmap->zsize;
00249 int n;
00250 
00251 switch (method) {
00252 case COMBINE_AVG:
00253 for (n=0; n<gridsize; n++)
00254 volmap->data[n] = voldata[n]/computed_frames;
00255 break;
00256 case COMBINE_PMF:
00257 for (n=0; n<gridsize; n++) {
00258 float val = voldata[n] + logf((float) computed_frames);
00259 if (val > MAX_ENERGY) val = MAX_ENERGY; // weed out outlying data
00260 volmap->data[n] = val;
00261 }
00262 break;
00263 case COMBINE_MAX:
00264 case COMBINE_MIN: 
00265 memcpy(volmap->data, voldata, gridsize*sizeof(float)); 
00266 break;
00267 case COMBINE_STDEV: 
00268 for (n=0; n<gridsize; n++) {
00269 volmap->data[n] = voldata[n]/computed_frames;
00270 volmap->data[n] = sqrtf(voldata2[n]/computed_frames - volmap->data[n]*volmap->data[n]); 
00271 }
00272 break;
00273 }
00274 }
00275 
00276 
00278 void VolMapCreate::combo_end(CombineType method, void *customptr) {
00279 if (method == COMBINE_STDEV) {
00280 float *voldata2 = (float*) customptr;
00281 delete[] voldata2;
00282 }
00283 }
00284 
00285 
00290 int VolMapCreate::compute_all (bool allframes, CombineType method, void *params) {
00291 int err = this->compute_init();
00292 if (err) return err;
00293 
00294 int gridsize = volmap->xsize*volmap->ysize*volmap->zsize;
00295 
00296 // Special case: if only have one frame do it here the fast way
00297 if (!allframes) {
00298 if (volmap->data) delete[] volmap->data;
00299 volmap->data = new float[gridsize]; // final exported voldata
00300 
00301 msgInfo << "volmap: grid size = " << volmap->xsize
00302 <<"x"<< volmap->ysize <<"x"<< volmap->zsize;
00303 char tmp[64];
00304 sprintf(tmp, " (%.1f MB)", sizeof(float)*gridsize/(1024.*1024.));
00305 msgInfo << tmp << sendmsg;
00306 
00307 // only compute the current frame of the given selection
00308 this->compute_frame(sel->which_frame, volmap->data); 
00309 
00310 //err = this->compute_end();
00311 //if (err) return err;
00312 return 0; // no error
00313 }
00314 
00315 
00316 msgInfo << "volmap: grid size = " << volmap->xsize
00317 <<"x"<< volmap->ysize <<"x"<< volmap->zsize;
00318 char tmp[64];
00319 sprintf(tmp, " (%.1f MB)", sizeof(float)*gridsize/(1024.*1024.));
00320 msgInfo << tmp << sendmsg;
00321 
00322 int numframes = app->molecule_numframes(sel->molid());
00323 msgInfo << "volmap: Computing " << numframes << " frames in total..." << sendmsg;
00324 
00325 if (volmap->data) delete[] volmap->data;
00326 volmap->data = new float[gridsize]; // final exported voldata
00327 float *frame_voldata = new float[gridsize]; // individual frame voldata
00328 float *voldata = new float[gridsize]; // combo cache voldata
00329 
00330 void *customptr = NULL;
00331 combo_begin(method, &customptr, params);
00332 wkf_timerhandle timer = wkf_timer_create();
00333 
00334 // Combine frame_voldata into voldata, one frame at a time, starting with 1st frame
00335 int frame;
00336 for (frame=0; frame<numframes; frame++) { 
00337 // XXX to-do, only take frames from a frame selection
00338 msgInfo << "volmap: frame " << frame << "/" << numframes;
00339 #ifdef TIMING
00340 msgInfo << sendmsg;
00341 #else
00342 msgInfo << " ";
00343 #endif
00344 
00345 wkf_timer_start(timer);
00346 
00347 this->compute_frame(frame, frame_voldata);
00348 wkf_timer_stop(timer);
00349 msgInfo << "Total time = " << wkf_timer_time(timer) << " s" << sendmsg;
00350 
00351 combo_addframe(method, voldata, customptr, frame_voldata);
00352 if (checkpoint_freq && computed_frames && !(computed_frames%checkpoint_freq)) {
00353 combo_export(method, voldata, customptr);
00354 const char *filename;
00355 if (checkpoint_name) filename=checkpoint_name;
00356 else filename = "checkpoint.dx";
00357 write_map(filename);
00358 }
00359 }
00360 
00361 wkf_timer_destroy(timer);
00362 
00363 delete[] frame_voldata;
00364 
00365 // All frames have been combined, perform finishing steps here
00366 combo_export(method, voldata, customptr);
00367 combo_end (method, customptr);
00368 delete[] voldata;
00369 
00370 return 0; // no error
00371 }
00372 
00373 
00374 
00375 // compute_init() sets up the grid coordinate system and dimensions
00376 // If the user did not specify the grid's minmax boundary, it is
00377 // defaulted to the trajectory's minmax coordinates, to which "padding"
00378 // is added in all dimensions. 
00379 int VolMapCreate::compute_init (float padding) {
00380 if (!sel) return MEASURE_ERR_NOSEL;
00381 if (sel->num_atoms == 0) return MEASURE_ERR_NOATOMS;
00382 
00383 int err, i;
00384 
00385 if (!volmap) return -1;
00386 
00387 if (user_minmax)
00388 padding = 0.; // don't want to pad user's defaults
00389 else {
00390 // The minmax coordinates over all frames
00391 err = calculate_minmax(min_coord, max_coord);
00392 if (err) return err;
00393 }
00394 
00395 // Depending on the volmap type we are computing we add some padding
00396 // to the dimensions of the grid. This function is called by the 
00397 // compute_init() methods of the VolMapCreate<type> subclasses. 
00398 // The caller provides a reasonable map type specific padding value.
00399 // The padding can for example be the radius of the largest atom or
00400 // an interaction distance cutoff.
00401 // After padding we align the grid with integer coordinates.
00402 for (i=0; i<3; i++) {
00403 //adjust padding and ensure that different maps are properly aligned
00404 min_coord[i] = (float) floor((min_coord[i] - padding)/delta)*delta; 
00405 max_coord[i] = (float) ceil((max_coord[i] + padding)/delta)*delta;
00406 }
00407 
00408 volmap->xsize = MAX((int)((max_coord[0] - min_coord[0])/delta), 0);
00409 volmap->ysize = MAX((int)((max_coord[1] - min_coord[1])/delta), 0);
00410 volmap->zsize = MAX((int)((max_coord[2] - min_coord[2])/delta), 0);
00411 
00412 char tmpstr[256] = { 0 };
00413 sprintf(tmpstr, "{%g %g %g} {%g %g %g}\n", min_coord[0], min_coord[1], min_coord[2], max_coord[0], max_coord[1], max_coord[2]);
00414 msgInfo << "volmap: grid minmax = " << tmpstr << sendmsg; 
00415 
00416 float cellx[3], celly[3], cellz[3];
00417 cellx[0] = delta;
00418 cellx[1] = 0.f;
00419 cellx[2] = 0.f;
00420 celly[0] = 0.f;
00421 celly[1] = delta;
00422 celly[2] = 0.f;
00423 cellz[0] = 0.f;
00424 cellz[1] = 0.f;
00425 cellz[2] = delta;
00426 
00427 // define origin by shifting to middle of each cell,
00428 // compute_frame() needs to take this into account
00429 for (i=0; i<3; i++) 
00430 volmap->origin[i] = min_coord[i] + \
00431 0.5f*(cellx[i] + celly[i] + cellz[i]);
00432 int d;
00433 for (d=0; d<3; d++) {
00434 volmap->xaxis[d] = cellx[d]*(volmap->xsize-1);
00435 volmap->yaxis[d] = celly[d]*(volmap->ysize-1);
00436 volmap->zaxis[d] = cellz[d]*(volmap->zsize-1);
00437 }
00438 
00439 if ((volmap->xsize*volmap->ysize*volmap->zsize) == 0)
00440 return MEASURE_ERR_ZEROGRIDSIZE;
00441 
00442 return 0; // no error
00443 }
00444 
00445 
00446 
00448 
00453 
00454 int VolMapCreateMask::compute_init() {
00455 char tmpstr[255] = { 0 };
00456 sprintf(tmpstr, "mask (%s.200)", sel->cmdStr);
00457 volmap->set_name(tmpstr);
00458 
00459 return VolMapCreate::compute_init(atomradius+0.5f);
00460 }
00461 
00462 
00463 int VolMapCreateMask::compute_frame (int frame, float *voldata) {
00464 DrawMolecule *mol = app->moleculeList->mol_from_id(sel->molid());
00465 if (!mol) return -1;
00466 
00467 int i;
00468 int GRIDSIZEX = volmap->xsize;
00469 int GRIDSIZEY = volmap->ysize;
00470 int GRIDSIZEZ = volmap->zsize;
00471 int gridsize = volmap->xsize*volmap->ysize*volmap->zsize;
00472 
00473 //create volumetric mask grid
00474 memset(voldata, 0, gridsize*sizeof(float));
00475 int save_frame = sel->which_frame;
00476 sel->which_frame=frame;
00477 sel->change(NULL,mol);
00478 
00479 const float *coords = sel->coordinates(app->moleculeList);
00480 if (!coords) {
00481 sel->which_frame = save_frame;
00482 return -1;
00483 }
00484 
00485 float cellx[3], celly[3], cellz[3];
00486 volmap->cell_axes(cellx, celly, cellz);
00487 
00488 float min_coords[3];
00489 for (i=0; i<3; i++)
00490 min_coords[i] = float(volmap->origin[i] - 0.5f*(cellx[i] + celly[i] + cellz[i]));
00491 
00492 // paint atomic spheres on map
00493 int gx, gy, gz;
00494 for (i=sel->firstsel; i<=sel->lastsel; i++) { 
00495 if (!sel->on[i]) continue; //atom is not selected
00496 
00497 gx = (int) ((coords[3L*i ] - min_coords[0])/delta);
00498 gy = (int) ((coords[3L*i+1] - min_coords[1])/delta);
00499 gz = (int) ((coords[3L*i+2] - min_coords[2])/delta);
00500 
00501 int steps = (int)(atomradius/delta)+1;
00502 int iz, iy, ix;
00503 for (iz=MAX(gz-steps,0); iz<=MIN(gz+steps,GRIDSIZEZ-1); iz++)
00504 for (iy=MAX(gy-steps,0); iy<=MIN(gy+steps,GRIDSIZEY-1); iy++)
00505 for (ix=MAX(gx-steps,0); ix<=MIN(gx+steps,GRIDSIZEX-1); ix++) {
00506 int n = ix + iy*GRIDSIZEX + iz*GRIDSIZEY*GRIDSIZEX;
00507 float dx = float(coords[3L*i ] - volmap->origin[0] - ix*delta);
00508 float dy = float(coords[3L*i+1] - volmap->origin[1] - iy*delta);
00509 float dz = float(coords[3L*i+2] - volmap->origin[2] - iz*delta);
00510 float dist2 = dx*dx+dy*dy+dz*dz;
00511 if (dist2 <= atomradius*atomradius) voldata[n] = 1.f;
00512 }
00513 }
00514 
00515 sel->which_frame = save_frame;
00516 
00517 return 0;
00518 } 
00519 
00520 
00521 
00523 
00531 
00532 int VolMapCreateDensity::compute_init () {
00533 char tmpstr[255] = { 0 };
00534 sprintf(tmpstr, "density (%.200s) [A^-3]", sel->cmdStr);
00535 volmap->set_name(tmpstr);
00536 
00537 float max_rad=0.f;
00538 calculate_max_radius(max_rad);
00539 
00540 return VolMapCreate::compute_init(MAX(3.f*radius_scale*max_rad,10.f));
00541 }
00542 
00543 
00544 int VolMapCreateDensity::compute_frame (int frame, float *voldata) {
00545 if (!weight) return MEASURE_ERR_NOWEIGHT;
00546 
00547 DrawMolecule *mol = app->moleculeList->mol_from_id(sel->molid());
00548 if (!mol) return -1;
00549 int i;
00550 
00551 const float *radius = mol->extraflt.data("radius");
00552 if (!radius) return MEASURE_ERR_NORADII;
00553 
00554 int GRIDSIZEX = volmap->xsize;
00555 int GRIDSIZEY = volmap->ysize;
00556 int GRIDSIZEZ = volmap->zsize;
00557 int gridsize = volmap->xsize*volmap->ysize*volmap->zsize;
00558 
00559 //create volumetric density grid
00560 memset(voldata, 0, gridsize*sizeof(float));
00561 int save_frame = sel->which_frame;
00562 sel->which_frame = frame;
00563 sel->change(NULL,mol);
00564 const float *coords = sel->coordinates(app->moleculeList);
00565 if (!coords) {
00566 sel->which_frame = save_frame;
00567 return -1;
00568 }
00569 
00570 if (weight_mutable) {
00571 if (weight_string) {
00572 get_weights_from_attribute(app, sel, weight_string, weight);
00573 } else {
00574 atomsel_default_weights(sel, weight);
00575 }
00576 }
00577 
00578 float cellx[3], celly[3], cellz[3];
00579 volmap->cell_axes(cellx, celly, cellz);
00580 
00581 float min_coords[3];
00582 for (i=0; i<3; i++) 
00583 min_coords[i] = float(volmap->origin[i] - 0.5f*(cellx[i] + celly[i] + cellz[i]));
00584 
00585 int gx, gy, gz; // grid coord indices
00586 for (i=sel->firstsel; i<=sel->lastsel; i++) { 
00587 if (!sel->on[i]) continue; //atom is not selected
00588 
00589 gx = (int) ((coords[3L*i ] - min_coords[0])/delta);
00590 gy = (int) ((coords[3L*i+1] - min_coords[1])/delta);
00591 gz = (int) ((coords[3L*i+2] - min_coords[2])/delta);
00592 
00593 float scaled_radius = 0.5f*radius_scale*radius[i];
00594 float exp_factor = 1.0f/(2.0f*scaled_radius*scaled_radius);
00595 float norm = weight[i]/(sqrtf((float) (8.0f*VMD_PI*VMD_PI*VMD_PI))*scaled_radius*scaled_radius*scaled_radius);
00596 
00597 int steps = (int)(4.1f*scaled_radius/delta);
00598 int iz, iy, ix;
00599 for (iz=MAX(gz-steps,0); iz<=MIN(gz+steps,GRIDSIZEZ-1); iz++)
00600 for (iy=MAX(gy-steps,0); iy<=MIN(gy+steps,GRIDSIZEY-1); iy++)
00601 for (ix=MAX(gx-steps,0); ix<=MIN(gx+steps,GRIDSIZEX-1); ix++) {
00602 int n = ix + iy*GRIDSIZEX + iz*GRIDSIZEY*GRIDSIZEX;
00603 float dx = float(coords[3L*i ] - volmap->origin[0] - ix*delta);
00604 float dy = float(coords[3L*i+1] - volmap->origin[1] - iy*delta);
00605 float dz = float(coords[3L*i+2] - volmap->origin[2] - iz*delta);
00606 float dist2 = dx*dx+dy*dy+dz*dz;
00607 voldata[n] += norm * expf(-dist2*exp_factor);
00608 // Uncomment the following line for a much faster implementation
00609 // This is useful is all you care about is the smooth visual appearance
00610 // voldata[n] += exp_factor/(dist2+10.f);
00611 }
00612 }
00613 
00614 sel->which_frame = save_frame;
00615 
00616 return 0;
00617 } 
00618 
00619 
00620 
00621 
00623 
00627 
00628 int VolMapCreateInterp::compute_init () {
00629 char tmpstr[255] = { 0 };
00630 sprintf(tmpstr, "interp (%.200s) [A^-3]", sel->cmdStr);
00631 volmap->set_name(tmpstr);
00632 
00633 return VolMapCreate::compute_init(delta+0.5f);
00634 }
00635 
00636 
00637 int VolMapCreateInterp::compute_frame (int frame, float *voldata) {
00638 if (!weight) return MEASURE_ERR_NOWEIGHT;
00639 
00640 DrawMolecule *mol = app->moleculeList->mol_from_id(sel->molid());
00641 if (!mol) return -1;
00642 int i;
00643 
00644 int GRIDSIZEX = volmap->xsize;
00645 int GRIDSIZEY = volmap->ysize;
00646 int GRIDSIZEXY = GRIDSIZEX * GRIDSIZEY;
00647 int gridsize = volmap->xsize*volmap->ysize*volmap->zsize;
00648 
00649 // create volumetric density grid
00650 memset(voldata, 0, gridsize*sizeof(float));
00651 int save_frame = sel->which_frame;
00652 sel->which_frame = frame;
00653 sel->change(NULL,mol);
00654 const float *coords = sel->coordinates(app->moleculeList);
00655 if (!coords) {
00656 sel->which_frame = save_frame;
00657 return -1;
00658 }
00659 
00660 if (weight_mutable) {
00661 if (weight_string) {
00662 get_weights_from_attribute(app, sel, weight_string, weight);
00663 } else {
00664 atomsel_default_weights(sel, weight);
00665 }
00666 }
00667 
00668 int gx, gy, gz; // grid coord indices
00669 float fgx, fgy, fgz; // fractional grid coord indices
00670 float dx, dy, dz; // to measure distances
00671 for (i=sel->firstsel; i<=sel->lastsel; i++) { 
00672 if (!sel->on[i]) continue; //atom is not selected
00673 
00674 // Find position of the atom within the map ("fractional indices")
00675 fgx = float(coords[3L*i ] - volmap->origin[0])/delta;
00676 fgy = float(coords[3L*i+1] - volmap->origin[1])/delta;
00677 fgz = float(coords[3L*i+2] - volmap->origin[2])/delta;
00678 
00679 // Find nearest voxel with lowest indices
00680 gx = (int) fgx;
00681 gy = (int) fgy;
00682 gz = (int) fgz;
00683 
00684 // Calculate distance between atom and each voxel
00685 dx = fgx - gx;
00686 dy = fgy - gy;
00687 dz = fgz - gz;
00688 
00689 // Perform trilinear interpolation
00690 
00691 voldata[ gx + gy*GRIDSIZEX + gz*GRIDSIZEXY ] \
00692 += (1.0f - dx) * (1.0f - dy) * (1.0f - dz) * weight[i];
00693 
00694 voldata[ (gx+1) + (gy+1)*GRIDSIZEX + (gz+1)*GRIDSIZEXY ] \
00695 += dx * dy * dz * weight[i];
00696 
00697 voldata[ (gx+1) + (gy+1)*GRIDSIZEX + gz*GRIDSIZEXY ] \
00698 += dx * dy * (1.0f - dz) * weight[i];
00699 
00700 voldata[ gx + gy*GRIDSIZEX + (gz+1)*GRIDSIZEXY ] \
00701 += (1.0f - dx) * (1.0f - dy) * dz * weight[i];
00702 
00703 voldata[ (gx+1) + gy*GRIDSIZEX + gz*GRIDSIZEXY ] \
00704 += dx * (1.0f - dy) * (1.0f - dz) * weight[i];
00705 
00706 voldata[ gx + (gy+1)*GRIDSIZEX + (gz+1)*GRIDSIZEXY ] \
00707 += (1.0f - dx) * dy * dz * weight[i];
00708 
00709 voldata[ gx + (gy+1)*GRIDSIZEX + gz*GRIDSIZEXY ] \
00710 += (1.0f - dx) * dy * (1.0f - dz) * weight[i];
00711 
00712 voldata[ (gx+1) + gy*GRIDSIZEX + (gz+1)*GRIDSIZEXY ] \
00713 += dx * (1.0f - dy) * dz * weight[i];
00714 }
00715 
00716 sel->which_frame = save_frame;
00717 
00718 return 0;
00719 } 
00720 
00721 
00722 
00723 
00724 
00726 
00732 
00733 int VolMapCreateOccupancy::compute_init () {
00734 char tmpstr[255] = { 0 };
00735 sprintf(tmpstr, "occupancy (%.200s)", sel->cmdStr);
00736 volmap->set_name(tmpstr);
00737 
00738 float max_rad=0.f; 
00739 if (use_points)
00740 max_rad = 1.f;
00741 else
00742 calculate_max_radius(max_rad);
00743 
00744 return VolMapCreate::compute_init(max_rad);
00745 }
00746 
00747 
00748 int VolMapCreateOccupancy::compute_frame(int frame, float *voldata) { 
00749 DrawMolecule *mol = app->moleculeList->mol_from_id(sel->molid());
00750 if (!mol) return -1;
00751 
00752 int GRIDSIZEX = volmap->xsize;
00753 int GRIDSIZEY = volmap->ysize;
00754 int GRIDSIZEZ = volmap->zsize;
00755 int gridsize = volmap->xsize*volmap->ysize*volmap->zsize;
00756 int i;
00757 
00758 //create volumetric density grid
00759 memset(voldata, 0, gridsize*sizeof(float));
00760 int save_frame = sel->which_frame;
00761 sel->which_frame=frame;
00762 sel->change(NULL,mol);
00763 const float *coords = sel->coordinates(app->moleculeList);
00764 
00765 if (!coords) {
00766 sel->which_frame = save_frame;
00767 return -1;
00768 }
00769 
00770 float cellx[3], celly[3], cellz[3];
00771 volmap->cell_axes(cellx, celly, cellz);
00772 
00773 float min_coords[3];
00774 for (i=0; i<3; i++) 
00775 min_coords[i] = float(volmap->origin[i] - 0.5f*(cellx[i] + celly[i] + cellz[i]));
00776 
00777 int gx, gy, gz;
00778 
00779 if (use_points) { // draw single points
00780 for (i=sel->firstsel; i<=sel->lastsel; i++) { 
00781 if (!sel->on[i]) continue; //atom is not selected
00782 
00783 gx = (int) ((coords[3L*i ] - min_coords[0])/delta);
00784 if (gx<0 || gx>=GRIDSIZEX) continue;
00785 gy = (int) ((coords[3L*i+1] - min_coords[1])/delta);
00786 if (gy<0 || gy>=GRIDSIZEY) continue;
00787 gz = (int) ((coords[3L*i+2] - min_coords[2])/delta);
00788 if (gz<0 || gz>=GRIDSIZEZ) continue;
00789 
00790 voldata[gx+GRIDSIZEX*gy+GRIDSIZEX*GRIDSIZEY*gz] = 1.f; 
00791 }
00792 }
00793 else { // paint atomic spheres on map
00794 const float *radius = mol->extraflt.data("radius");
00795 if (!radius) {
00796 sel->which_frame = save_frame;
00797 return MEASURE_ERR_NORADII;
00798 }
00799 
00800 for (i=sel->firstsel; i<=sel->lastsel; i++) { 
00801 if (!sel->on[i]) continue; //atom is not selected
00802 
00803 gx = (int) ((coords[3L*i ] - min_coords[0])/delta);
00804 gy = (int) ((coords[3L*i+1] - min_coords[1])/delta);
00805 gz = (int) ((coords[3L*i+2] - min_coords[2])/delta);
00806 
00807 int steps = (int)(radius[i]/delta)+1;
00808 int iz, iy, ix;
00809 for (iz=MAX(gz-steps,0); iz<=MIN(gz+steps,GRIDSIZEZ-1); iz++)
00810 for (iy=MAX(gy-steps,0); iy<=MIN(gy+steps,GRIDSIZEY-1); iy++)
00811 for (ix=MAX(gx-steps,0); ix<=MIN(gx+steps,GRIDSIZEX-1); ix++) {
00812 int n = ix + iy*GRIDSIZEX + iz*GRIDSIZEY*GRIDSIZEX;
00813 float dx = float(coords[3L*i ] - volmap->origin[0] - ix*delta);
00814 float dy = float(coords[3L*i+1] - volmap->origin[1] - iy*delta);
00815 float dz = float(coords[3L*i+2] - volmap->origin[2] - iz*delta);
00816 float dist2 = dx*dx+dy*dy+dz*dz;
00817 if (dist2 <= radius[i]*radius[i]) voldata[n] = 1.f;
00818 }
00819 }
00820 }
00821 
00822 sel->which_frame = save_frame;
00823 
00824 return 0;
00825 }
00826 
00827 
00828 
00830 
00836 
00837 int VolMapCreateDistance::compute_init () {
00838 char tmpstr[255] = { 0 };
00839 sprintf(tmpstr, "distance (%.200s) [A]", sel->cmdStr);
00840 volmap->set_name(tmpstr);
00841 
00842 float max_rad=0.f;
00843 calculate_max_radius(max_rad);
00844 
00845 return VolMapCreate::compute_init(max_rad+max_dist);
00846 }
00847 
00848 
00851 int VolMapCreateDistance::compute_frame(int frame, float *voldata) { 
00852 int i, n; 
00853 DrawMolecule *mol = app->moleculeList->mol_from_id(sel->molid());
00854 if (!mol) return -1;
00855 const float *radius = mol->extraflt.data("radius");
00856 if (!radius) return MEASURE_ERR_NORADII;
00857 
00858 int GRIDSIZEX = volmap->xsize;
00859 int GRIDSIZEY = volmap->ysize;
00860 int gridsize = volmap->xsize*volmap->ysize*volmap->zsize;
00861 
00862 float dx, dy, dz;
00863 float dist, mindist, r;
00864 
00865 float max_rad=0.f;
00866 calculate_max_radius(max_rad);
00867 
00868 // 1. Create a fake "molecule" containing all of the grid points
00869 // this is quite memory intensive but _MUCH_ faster doing it point-by point!
00870 
00871 float *gridpos = new float[3L*gridsize]; 
00872 int *gridon = new int[gridsize]; 
00873 for (n=0; n<gridsize; n++) {
00874 gridpos[3L*n ] = float((n%GRIDSIZEX)*delta + volmap->origin[0]); //position of grid cell's center
00875 gridpos[3L*n+1] = float(((n/GRIDSIZEX)%GRIDSIZEY)*delta + volmap->origin[1]);
00876 gridpos[3L*n+2] = float((n/(GRIDSIZEX*GRIDSIZEY))*delta + volmap->origin[2]); 
00877 gridon[n] = 1;
00878 }
00879 
00880 GridSearchPair *pairlist, *p;
00881 
00882 int save_frame = sel->which_frame;
00883 sel->which_frame = frame;
00884 sel->change(NULL,mol);
00885 const float *coords = sel->coordinates(app->moleculeList);
00886 if (!coords) {
00887 sel->which_frame = save_frame;
00888 return -1;
00889 }
00890 
00891 // initialize all grid points to be the maximal allowed distance = cutoff
00892 for (n=0; n<gridsize; n++) voldata[n] = max_dist;
00893 
00894 // 2. Create a list of all bonds between the grid and the real molecule
00895 // which are within the user-set cutoff distance 
00896 // (the use of a cutoff is purely to speed this up tremendously)
00897 
00898 pairlist = vmd_gridsearch3(gridpos, gridsize, gridon, coords,
00899 sel->num_atoms, sel->on, max_dist+max_rad, true, -1);
00900 for (p=pairlist; p; p=p->next) {
00901 n = p->ind1;
00902 // if a grid point is already known to be inside an atom, skip it and save some time
00903 if ((mindist = voldata[n]) == 0.f) continue;
00904 i = p->ind2;
00905 r = radius[i];
00906 dx = gridpos[3L*n ] - coords[3L*i];
00907 dy = gridpos[3L*n+1] - coords[3L*i+1];
00908 dz = gridpos[3L*n+2] - coords[3L*i+2];
00909 
00910 // 3. At each grid point, store the _smallest_ recorded distance
00911 // to a nearby atomic surface
00912 
00913 dist = sqrtf(dx*dx+dy*dy+dz*dz) - r;
00914 if (dist < 0) dist = 0.f;
00915 if (dist < mindist) voldata[n] = dist;
00916 }
00917 
00918 // delete pairlist
00919 for (p=pairlist; p;) {
00920 GridSearchPair *tmp = p;
00921 p = p->next;
00922 free(tmp);
00923 } 
00924 
00925 delete [] gridpos; 
00926 delete [] gridon; 
00927 
00928 sel->which_frame = save_frame;
00929 
00930 return MEASURE_NOERR; 
00931 }
00932 
00933 
00934 
00935 
00937 
00938 int VolMapCreateCoulombPotential::compute_init () {
00939 char tmpstr[255] = { 0 };
00940 sprintf(tmpstr, "Potential (kT/e at 298.15K) (%.200s)", sel->cmdStr);
00941 volmap->set_name(tmpstr);
00942 
00943 float max_rad;
00944 calculate_max_radius(max_rad);
00945 
00946 // init object, no extra padding by default 
00947 return VolMapCreate::compute_init(0.f);
00948 }
00949 
00950 
00951 int VolMapCreateCoulombPotential::compute_frame(int frame, float *voldata) {
00952 DrawMolecule *mol = app->moleculeList->mol_from_id(sel->molid());
00953 if (!mol) return -1;
00954 int i;
00955 
00956 const float *charge = mol->extraflt.data("charge");
00957 if (!charge) return MEASURE_ERR_NORADII; // XXX fix this later
00958 
00959 int gridsize=volmap->xsize*volmap->ysize*volmap->zsize;
00960 
00961 // create volumetric density grid
00962 memset(voldata, 0, gridsize*sizeof(float));
00963 int save_frame = sel->which_frame;
00964 sel->which_frame=frame;
00965 sel->change(NULL,mol);
00966 const float *coords = sel->coordinates(app->moleculeList);
00967 if (!coords) {
00968 sel->which_frame = save_frame;
00969 return -1;
00970 }
00971 
00972 float cellx[3], celly[3], cellz[3];
00973 volmap->cell_axes(cellx, celly, cellz);
00974 
00975 float min_coords[3];
00976 for (i=0; i<3; i++) 
00977 min_coords[i] = float(volmap->origin[i] - 0.5f*(cellx[i] + celly[i] + cellz[i]));
00978 
00979 // copy selected atom coordinates and charges to a contiguous memory
00980 // buffer and translate them to the starting corner of the map.
00981 float *xyzq = (float *) malloc(sel->selected * 4L * sizeof(float));
00982 float *curatom = xyzq;
00983 for (i=sel->firstsel; i<=sel->lastsel; i++) { 
00984 if (sel->on[i]) {
00985 curatom[0] = coords[3L*i ] - min_coords[0];
00986 curatom[1] = coords[3L*i+1] - min_coords[1];
00987 curatom[2] = coords[3L*i+2] - min_coords[2];
00988 curatom[3] = charge[i] * float(POT_CONV);
00989 curatom += 4;
00990 }
00991 }
00992 
00993 vol_cpotential(sel->selected, xyzq, voldata, 
00994 volmap->zsize, volmap->ysize, volmap->xsize, delta);
00995 
00996 free(xyzq);
00997 
00998 sel->which_frame = save_frame;
00999 
01000 return 0;
01001 } 
01002 
01004 
01005 #if defined(VMDUSEMSMPOT)
01006 int VolMapCreateCoulombPotentialMSM::compute_init () {
01007 char tmpstr[255] = { 0 };
01008 sprintf(tmpstr, "Potential (kT/e at 298.15K) (%.200s)", sel->cmdStr);
01009 volmap->set_name(tmpstr);
01010 
01011 float max_rad;
01012 calculate_max_radius(max_rad);
01013 
01014 // init object, no extra padding by default
01015 // Note: padding would create serious problems for the periodic case 
01016 return VolMapCreate::compute_init(0.f);
01017 }
01018 
01019 
01020 int VolMapCreateCoulombPotentialMSM::compute_frame(int frame, float *voldata) {
01021 DrawMolecule *mol = app->moleculeList->mol_from_id(sel->molid());
01022 if (!mol) return -1;
01023 int i;
01024 
01025 int usepbc = 0;
01026 
01027 const float *charge = mol->extraflt.data("charge");
01028 if (!charge) return MEASURE_ERR_NORADII; // XXX fix this later
01029 
01030 int gridsize=volmap->xsize*volmap->ysize*volmap->zsize;
01031 
01032 // create volumetric density grid
01033 memset(voldata, 0, gridsize*sizeof(float));
01034 int save_frame = sel->which_frame;
01035 sel->which_frame=frame;
01036 sel->change(NULL,mol);
01037 const float *coords = sel->coordinates(app->moleculeList);
01038 const Timestep *ts = sel->timestep(app->moleculeList); 
01039 if (!coords) {
01040 sel->which_frame = save_frame;
01041 return -1;
01042 }
01043 if (!ts) {
01044 return -1;
01045 }
01046 
01047 float cellx[3], celly[3], cellz[3];
01048 volmap->cell_axes(cellx, celly, cellz);
01049 
01050 float min_coords[3];
01051 for (i=0; i<3; i++) 
01052 min_coords[i] = float(volmap->origin[i] - 0.5f*(cellx[i] + celly[i] + cellz[i]));
01053 
01054 // copy selected atom coordinates and charges to a contiguous memory
01055 // buffer and translate them to the starting corner of the map.
01056 float *xyzq = (float *) malloc(sel->selected * 4L * sizeof(float));
01057 float *curatom = xyzq;
01058 for (i=sel->firstsel; i<=sel->lastsel; i++) { 
01059 if (sel->on[i]) {
01060 curatom[0] = coords[3L*i ] - min_coords[0];
01061 curatom[1] = coords[3L*i+1] - min_coords[1];
01062 curatom[2] = coords[3L*i+2] - min_coords[2];
01063 curatom[3] = charge[i] * float(POT_CONV);
01064 curatom += 4;
01065 }
01066 }
01067 
01068 Msmpot *msm = Msmpot_create(); // create a multilevel summation object
01069 #if 0
01070 int msmrc;
01071 int mx = volmap->xsize; /* map lattice dimensions */
01072 int my = volmap->ysize;
01073 int mz = volmap->zsize;
01074 float lx = delta*mx; /* map lattice lengths */
01075 float ly = delta*my;
01076 float lz = delta*mz;
01077 float x0=0, y0=0, z0=0; /* map origin */
01078 float vx=0, vy=0, vz=0; /* periodic domain lengths (0 for nonperiodic) */
01079 
01080 if (getenv("MSMPOT_NOCUDA")) {
01081 /* turn off use of CUDA (with 0 in last parameter) */
01082 Msmpot_configure(msm, 0, 0, 0, 0, 0, 0, 0, 0, 0);
01083 }
01084 
01085 if (getenv("MSMPOT_PBCON")) {
01086 vx = lx, vy = ly, vz = lz; /* use periodic boundary conditions */
01087 }
01088 
01089 if (getenv("MSMPOT_EXACT")) { 
01090 msmrc = Msmpot_compute_exact(msm, voldata, 
01091 mx, my, mz, lx, ly, lz, x0, y0, z0, vx, vy, vz,
01092 xyzq, sel->selected);
01093 }
01094 else {
01095 msmrc = Msmpot_compute(msm, voldata, 
01096 mx, my, mz, lx, ly, lz, x0, y0, z0, vx, vy, vz,
01097 xyzq, sel->selected);
01098 }
01099 if (msmrc != MSMPOT_SUCCESS) {
01100 printf("MSM return code: %d\n", msmrc);
01101 printf("MSM error string: '%s'\n", Msmpot_error_string(msmrc));
01102 } 
01103 #endif
01104 #if 1
01105 // New MSM API: both non-periodic and periodic MSM calcs
01106 int msmrc;
01107 
01108 // XXX hack for ease of initial testing
01109 if (getenv("VMDMSMUSEPBC"))
01110 usepbc = 1;
01111 
01112 if (usepbc) {
01113 // get periodic cell information for current frame
01114 float a, b, c, alpha, beta, gamma;
01115 a = ts->a_length;
01116 b = ts->b_length;
01117 c = ts->c_length;
01118 alpha = ts->alpha;
01119 beta = ts->beta;
01120 gamma = ts->gamma;
01121 
01122 // check validity of PBC cell side lengths
01123 if (fabsf(a*b*c) < 0.0001) {
01124 msgErr << "volmap coulombmsm: unit cell volume is zero." << sendmsg;
01125 return -1;
01126 }
01127 
01128 // check PBC unit cell shape to select proper low level algorithm.
01129 if ((alpha != 90.0) || (beta != 90.0) || (gamma != 90.0)) {
01130 msgErr << "volmap coulombmsm: unit cell is non-orthogonal." << sendmsg;
01131 return -1;
01132 }
01133 
01134 #ifdef MSMPOT_COMPUTE_EXACT
01135 if (getenv("MSMPOT_EXACT")) {
01136 // XXX the current PBC code will currently use the initially specified 
01137 // map dimensions and coordinates for all frames in a time-averaged
01138 // calculation. In the case that one would prefer the map to cover
01139 // a fixed region of the unit cell in reciprocal space, we will need
01140 // to change this code to update the effective map geometry on-the-fly.
01141 msgInfo << "Running EXACT periodic MSM calculation..." << sendmsg;
01142 msmrc = Msmpot_compute_exact(msm, voldata, 
01143 volmap->xsize, volmap->ysize, volmap->zsize,
01144 volmap->xsize * delta, 
01145 volmap->ysize * delta, 
01146 volmap->zsize * delta, 
01147 0, 0, 0, // origin, already translated to min 
01148 a, b, c, // pbc cell length 0 == nonperiodic calc
01149 xyzq, sel->selected);
01150 } else {
01151 #endif
01152 // XXX the current PBC code will currently use the initially specified 
01153 // map dimensions and coordinates for all frames in a time-averaged
01154 // calculation. In the case that one would prefer the map to cover
01155 // a fixed region of the unit cell in reciprocal space, we will need
01156 // to change this code to update the effective map geometry on-the-fly.
01157 msgInfo << "Running periodic MSM calculation..." << sendmsg;
01158 msmrc = Msmpot_compute(msm, voldata, 
01159 volmap->xsize, volmap->ysize, volmap->zsize,
01160 volmap->xsize * delta, 
01161 volmap->ysize * delta, 
01162 volmap->zsize * delta, 
01163 0, 0, 0, // origin, already translated to min 
01164 a, b, c, // pbc cell length 0 == nonperiodic calc
01165 xyzq, sel->selected);
01166 #ifdef MSMPOT_COMPUTE_EXACT
01167 }
01168 #endif
01169 
01170 } else {
01171 
01172 #ifdef MSMPOT_COMPUTE_EXACT
01173 if (getenv("MSMPOT_EXACT")) {
01174 msgInfo << "Running EXACT non-periodic MSM calculation..." << sendmsg;
01175 msmrc = Msmpot_compute_exact(msm, voldata, 
01176 volmap->xsize, volmap->ysize, volmap->zsize,
01177 volmap->xsize * delta, 
01178 volmap->ysize * delta, 
01179 volmap->zsize * delta, 
01180 0, 0, 0, // origin, already translated to min 
01181 0, 0, 0, // pbc cell length 0 == nonperiodic calc
01182 xyzq, sel->selected);
01183 } else {
01184 #endif
01185 msgInfo << "Running non-periodic MSM calculation..." << sendmsg;
01186 msmrc = Msmpot_compute(msm, voldata, 
01187 volmap->xsize, volmap->ysize, volmap->zsize,
01188 volmap->xsize * delta, 
01189 volmap->ysize * delta, 
01190 volmap->zsize * delta, 
01191 0, 0, 0, // origin, already translated to min 
01192 0, 0, 0, // pbc cell length 0 == nonperiodic calc
01193 xyzq, sel->selected);
01194 #ifdef MSMPOT_COMPUTE_EXACT
01195 }
01196 #endif
01197 
01198 }
01199 
01200 if (msmrc != MSMPOT_SUCCESS) {
01201 printf("MSM return code: %d\n", msmrc);
01202 printf("MSM error string: '%s'\n", Msmpot_error_string(msmrc));
01203 } 
01204 #else
01205 // old MSM API: non-periodic MSM calcs only
01206 int msmrc = Msmpot_compute(msm, voldata, 
01207 volmap->xsize, volmap->ysize, volmap->zsize,
01208 delta, delta, delta, 
01209 0, 0, 0,
01210 0, 0, 0, // origin, already translated to min 
01211 xyzq, sel->selected);
01212 if (msmrc != MSMPOT_ERROR_NONE) {
01213 printf("MSM return code: %d\n", msmrc);
01214 printf("MSM error string: '%s'\n", Msmpot_error_string(msmrc));
01215 } 
01216 #endif
01217 Msmpot_destroy(msm);
01218 
01219 free(xyzq);
01220 
01221 sel->which_frame = save_frame;
01222 
01223 return 0;
01224 } 
01225 
01226 #endif
01227 
01228 
01229 
01230 
01231 
01232 // Write the map as a DX file.
01233 // This is the default base class function which can be
01234 // overridden by the derived classes. 
01235 // E.g. VolMapCreateFastEnergy defines its own write_map().
01236 void VolMapCreate::write_map(const char *filename) {
01237 volmap_write_dx_file(volmap, filename);
01238 }
01239 
01240 
01241 int volmap_write_dx_file (VolumetricData *volmap, const char *filename) {
01242 if (!volmap->data) return -1; // XXX is this a good random error code?
01243 int i;
01244 int xsize = volmap->xsize;
01245 int ysize = volmap->ysize;
01246 int zsize = volmap->zsize;
01247 int gridsize = xsize*ysize*zsize;
01248 
01249 float cellx[3], celly[3], cellz[3];
01250 volmap->cell_axes(cellx, celly, cellz);
01251 
01252 
01253 msgInfo << "volmap: writing file \"" << filename << "\"." << sendmsg;
01254 
01255 FILE *fout = fopen(filename, "w");
01256 if (!fout) {
01257 msgErr << "volmap: Cannot open file \"" << filename
01258 << "\" for writing." << sendmsg;
01259 return errno;
01260 };
01261 
01262 fprintf(fout, "# Data calculated by the VMD volmap function\n");
01263 
01264 // Since the data origin and the grid origin are aligned we have
01265 // grid centered data, even though we were thinking in terms of
01266 // voxels centered in datapoints. VMD treats all dx file maps as
01267 // grid centered data, so this is right.
01268 fprintf(fout, "object 1 class gridpositions counts %d %d %d\n", xsize, ysize, zsize);
01269 fprintf(fout, "origin %g %g %g\n", volmap->origin[0], volmap->origin[1], volmap->origin[2]);
01270 fprintf(fout, "delta %g %g %g\n", cellx[0], cellx[1], cellx[2]);
01271 fprintf(fout, "delta %g %g %g\n", celly[0], celly[1], celly[2]);
01272 fprintf(fout, "delta %g %g %g\n", cellz[0], cellz[1], cellz[2]);
01273 fprintf(fout, "object 2 class gridconnections counts %d %d %d\n", xsize, ysize, zsize);
01274 fprintf(fout, "object 3 class array type double rank 0 items %d data follows\n", gridsize);
01275 
01276 // This reverses the ordering from x fastest to z fastest changing variable
01277 float val1,val2,val3;
01278 int gx=0, gy=0, gz=-1;
01279 for (i=0; i < (gridsize/3)*3; i+=3) {
01280 if (++gz >= zsize) {
01281 gz=0;
01282 if (++gy >= ysize) {gy=0; gx++;}
01283 }
01284 val1 = volmap->voxel_value(gx,gy,gz);
01285 if (++gz >= zsize) {
01286 gz=0;
01287 if (++gy >= ysize) {gy=0; gx++;}
01288 }
01289 val2 = volmap->voxel_value(gx,gy,gz);
01290 if (++gz >= zsize) {
01291 gz=0;
01292 if (++gy >= ysize) {gy=0; gx++;}
01293 }
01294 val3 = volmap->voxel_value(gx,gy,gz); 
01295 fprintf(fout, "%g %g %g\n", val1, val2, val3);
01296 }
01297 for (i=(gridsize/3)*3; i < gridsize; i++) {
01298 if (++gz >= zsize) {
01299 gz=0;
01300 if (++gy >= ysize) {gy=0; gx++;}
01301 }
01302 fprintf(fout, "%g ", volmap->voxel_value(gx,gy,gz));
01303 }
01304 if (gridsize%3) fprintf(fout, "\n");
01305 fprintf(fout, "\n");
01306 
01307 // Replace any double quotes (") by single quotes (') in the 
01308 // dataname string to make sure that we don't prematurely
01309 // terminate the string in the dx file.
01310 char *squotes = new char[strlen(volmap->name)+1];
01311 strcpy(squotes, volmap->name);
01312 char *s = squotes;
01313 while((s=strchr(s, '"'))) *s = '\'';
01314 
01315 if (volmap->name) {
01316 fprintf(fout, "object \"%s\" class field\n", squotes);
01317 } else {
01318 char dataname[10];
01319 strcpy(dataname, "(no name)");
01320 fprintf(fout, "object \"%s\" class field\n", dataname);
01321 }
01322 
01323 delete [] squotes;
01324 
01325 fclose(fout);
01326 return 0;
01327 }