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

VolMapCreateILS.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  * RCS INFORMATION:
00010  *
00011  *      $RCSfile: VolMapCreateILS.C,v $
00012  *      $Author: johns $        $Locker:  $             $State: Exp $
00013  *      $Revision: 1.172 $      $Date: 2020/07/08 04:40:23 $
00014  *
00015  ***************************************************************************/
00016 
00017 #include <math.h>
00018 #include <stdio.h>
00019 #include "VolMapCreate.h"
00020 #include "MoleculeList.h"
00021 #include "VolumetricData.h"
00022 #include "utilities.h"
00023 #include "WKFUtils.h"
00024 
00025 /* avoid parameter name collisions with AIX5 "hz" macro */
00026 #undef hz
00027 
00028 //#define DEBUG 1
00029 #include "VMDApp.h"
00030 #if defined(DEBUG)
00031 #include "MoleculeGraphics.h"  // needed only for debugging
00032 #endif
00033 
00034 #include "Measure.h"
00035 #include "Inform.h"
00036 
00037 #if defined(VMDCUDA)
00038 #include "CUDAKernels.h"
00039 #endif
00040 
00041 #define TIMING
00042 
00043 typedef union flint_t {
00044   float f;
00045   int i;
00046   char c[4];
00047 } flint;
00048 
00049 #define BIN_DEPTH     8  // number of slots per bin
00050 #define BIN_SLOTSIZE  4  // slot permits x, y, z, vdwtype  (4 elements)
00051 #define BIN_SIZE      (BIN_DEPTH * BIN_SLOTSIZE)  // size given in "flints"
00052 
00053 
00054 #define BIN_DEPTH  8
00055 
00056 
00057 struct AtomPosType {
00058   float x, y, z;                   // position coordinates of atom
00059   int vdwtype;                     // type index for van der Waals params
00060 };
00061 
00062 
00063 struct BinOfAtoms {
00064   AtomPosType atom[BIN_DEPTH];     // use fixed bin depth for CUDA
00065 };
00066 
00067 
00068 typedef struct ComputeOccupancyMap_t {
00069 
00070   // these are initialized by caller (pointers to existing memory allocations)
00071 
00072   float *map;                    // buffer space for occupancy map
00073   int mx, my, mz;                // map dimensions
00074   float lx, ly, lz;              // map lengths
00075   float x0, y0, z0;              // map origin
00076   float ax[3], ay[3], az[3];     // map basis vectors, aligned
00077   float alignmat[16];            // 4x4 matrix used for alignment
00078   int num_coords;                // number of atoms
00079   const float *coords;           // atom coords x/y/z, length 3*num_coords
00080 
00081   const int *vdw_type;           // type index for each atom, length num_coords
00082 
00083   const float *vdw_params;       // vdw parameters for system atoms, listed as
00084                                  // (epsilon, rmin) pairs for each type number
00085 
00086   const float *probe_vdw_params; // vdw parameters for probe atoms, listed as
00087                                  // (epsilon, rmin) pairs for each probe atom,
00088                                  // length 2*num_probes
00089 
00090   const float *conformers;       // length 3*num_probes*num_conformers,
00091                                  // with probe atoms listed consecutively
00092                                  // for each conformer
00093 
00094   int num_probes;                // number of probe atoms
00095   int num_conformers;            // number of conformers
00096 
00097   float cutoff;                  // cutoff distance
00098   float extcutoff;               // extended cutoff, includes probe radius
00099 
00100   float excl_dist;               // exclusion distance threshold
00101   float excl_energy;             // exclusion energy threshold
00102 
00103   int kstart, kstop;             // z-direction indexing for multi-threaded
00104                                  // slab decomposition of maps
00105                                  // kstart = starting index of slab
00106                                  // kstop = last index + 1
00107 
00108   // internal data (memory is allocated)
00109 
00110   float hx, hy, hz;              // derived map spacings
00111   float bx, by, bz;              // atom bin lengths
00112   float bx_1, by_1, bz_1;        // reciprocal bin lengths
00113   int mpblx, mpbly, mpblz;       // map points per bin side in x, y, z
00114   int cpu_only;                  // flag indicates that CUDA cannot be used
00115 
00116   BinOfAtoms *bin;               // padded bin of atoms
00117   BinOfAtoms *bin_zero;          // bin pointer shifted to (0,0,0)-index
00118   char *bincnt;                  // count occupied slots in each bin
00119   char *bincnt_zero;             // bincnt pointer shifted to (0,0,0)-index
00120 
00121   int nbx, nby, nbz;             // number of bins in x, y, z directions
00122   int padx, pady, padz;          // bin padding
00123 
00124   char *bin_offsets;             // bin neighborhood index offset
00125   int num_bin_offsets;           // length of tight 
00126 
00127   AtomPosType *extra;            // extra atoms that over fill bins
00128   int num_extras;                // number of extra atoms
00129 
00130   char *exclusions;              // same dimensions as map
00131 
00132   // data for CUDA goes here
00133 
00134 } ComputeOccupancyMap;
00135 
00136 
00137 #define DEFAULT_EXCL_DIST      1.f
00138 #define DEFAULT_EXCL_ENERGY   87.f
00139 
00140 #define DEFAULT_BIN_LENGTH     3.f
00141 #define MAX_BIN_VOLUME        27.f
00142 #define MIN_BIN_VOLUME         8.f
00143 
00144 
00145 // Must already have ComputeOccupancyMap initialized.
00146 // Performs geometric hashing of atoms into bins, along with all related
00147 // memory management.  Also allocates memory for map exclusions.
00148 static int ComputeOccupancyMap_setup(ComputeOccupancyMap *);
00149 
00150 // Calculate occupancy for slab of map indexed from kstart through kstop.
00151 // Starts by finding the exclusions, then continues by calculating the
00152 // occupancy for the given probe and conformers.
00153 static int ComputeOccupancyMap_calculate_slab(ComputeOccupancyMap *);
00154 
00155 // Cleanup memory allocations.
00156 static void ComputeOccupancyMap_cleanup(ComputeOccupancyMap *);
00157 
00158 
00159 // Write bin histogram into a dx map
00160 static void write_bin_histogram_map(
00161     const ComputeOccupancyMap *,
00162     const char *filename
00163     );
00164 
00165 static void atom_bin_stats(const ComputeOccupancyMap *);
00166 
00167 
00168 // XXX slow quadratic complexity algorithm for checking correctness
00169 //     for every map point, for every probe atom in all conformers,
00170 //     iterate over all atoms
00171 static void compute_allatoms(
00172     float *map,                    // return calculated occupancy map
00173     int mx, int my, int mz,        // dimensions of map
00174     float lx, float ly, float lz,  // lengths of map
00175     const float origin[3],         // origin of map
00176     const float axes[9],           // basis vectors of aligned map
00177     const float alignmat[16],      // 4x4 alignment matrix
00178     int num_coords,                // number of atoms
00179     const float *coords,           // atom coordinates, length 3*num_coords
00180     const int *vdw_type,           // vdw type numbers, length num_coords
00181     const float *vdw_params,       // scaled vdw parameters for atoms
00182     float cutoff,                  // cutoff distance
00183     const float *probe_vdw_params, // scaled vdw parameters for probe atoms
00184     int num_probe_atoms,           // number of atoms in probe
00185     int num_conformers,            // number of conformers
00186     const float *conformers,       // length 3*num_probe_atoms*num_conformers
00187     float excl_energy              // exclusion energy threshold
00188     );
00189 
00190 
00192 
00193 // Computes free energy maps from implicit ligand sampling.
00194 // I.e. it creates a map of the estimated potential of mean force
00195 // (in units of k_B*T at the specified temperature) of placing a 
00196 // weakly-interacting gas monoatomic or diatomic ligand in every 
00197 // voxel. Each voxel can be divided into a subgrid in order to obtain
00198 // a potential value for that voxel that is averaged over its subgrid 
00199 // positions. For diatomic ligands for each (sub-)gridpoint one can
00200 // (and should) also average over different random rotamers.
00201 
00202 // For the computation the trajectory frames must be aligned. One can
00203 // either manually align them beforehand or provide a selection
00204 // (AtomSel *alignsel) on which automatic alignment should be based.
00205 // If such a selection was provided then it is used to align all
00206 // trajectory frames to the first frame.
00207 // Note that the results are slightly different from what you get when
00208 // you align all frames manually prior to the computation.
00209 // The differences are quite noticable when comparing the files as
00210 // ascii texts but actually they are of numerical nature and when you
00211 // display the maps they look the same.
00212 // Suppose you want to align your trajectory to a reference frame from a
00213 // different molecule then you can specify the alignment matrix (as
00214 // returned by "measure fit") that would align the first frame to the
00215 // reference. The corresponding member variable is Matrix4 transform.
00216 
00217 VolMapCreateILS::VolMapCreateILS(VMDApp *_app,
00218                                  int mol, int firstframe, int lastframe,
00219                                  float T, float res, int nsub,
00220                                  float cut, int mask) : 
00221   app(_app), molid(mol), cutoff(cut),
00222   temperature(T), first(firstframe), last(lastframe),
00223   maskonly(mask)
00224 {
00225   compute_elec = false;
00226 
00227   num_probe_atoms = 0;
00228   num_conformers = 0;
00229   conformer_freq = 1;
00230   conformers = NULL;
00231   probe_coords = NULL;
00232   probe_vdw    = NULL;
00233   probe_charge = NULL;
00234   probe_axisorder1 = 1;
00235   probe_axisorder2 = 1;
00236   probe_tetrahedralsymm = 0;
00237   probe_effsize = 0.f;
00238   extcutoff = cutoff;
00239 
00240   max_energy = DEFAULT_EXCL_ENERGY;
00241   min_occup = expf(-max_energy);
00242 
00243   nsubsamp = 1;
00244   if (nsub>1) nsubsamp = nsub;
00245   delta = res/float(nsubsamp);
00246 
00247   alignsel = NULL;
00248   pbc = false;
00249   pbcbox = false;
00250   pbccenter[0] = pbccenter[1] = pbccenter[2] = 0.f;
00251   
00252   num_unique_types = 0;
00253   vdwparams  = NULL;
00254   atomtypes  = NULL;
00255 
00256   minmax[0] = 0.f;
00257   minmax[1] = 0.f;
00258   minmax[2] = 0.f;
00259   minmax[3] = 0.f;
00260   minmax[4] = 0.f;
00261   minmax[5] = 0.f;
00262   
00263   volmap  = NULL;
00264   volmask = NULL;
00265 }
00266 
00267 
00268 VolMapCreateILS::~VolMapCreateILS() {
00269   if (probe_coords) delete [] probe_coords;
00270   if (probe_vdw)    delete [] probe_vdw;
00271   if (probe_charge) delete [] probe_charge;
00272 
00273   if (conformers) delete [] conformers;
00274   if (vdwparams)  delete [] vdwparams;
00275   if (atomtypes)  delete [] atomtypes;
00276 
00277   if (volmap) delete volmap;
00278 
00279   //  if (volmask) delete volmask;
00280 }
00281 
00282 
00285 int VolMapCreateILS::write_map(const char *filename) {
00286   // Override the name string:
00287   // We include the number of frames used for sampling and the temperature
00288   // of the underlying MD simulation in the dataset name string.
00289   // This way we can use the VolumetricData class without adding new
00290   // non-general fields to it.
00291   char tmpstr[256] = { 0 };
00292   if (maskonly) {
00293     sprintf(tmpstr, "ligand pmf mask, %i frames, cutoff = %.2fA",
00294             computed_frames, cutoff);
00295   } else {
00296     sprintf(tmpstr, "ligand pmf [kT], %i frames, T = %.2f Kelvin, maxenergy = %g, cutoff = %.2fA",
00297             computed_frames, temperature, max_energy, cutoff);
00298   }
00299   volmap->set_name(tmpstr);
00300 
00301 
00302   // Add volmap to a molecule so that we can use
00303   // the plugin interface to write the dx map.
00304   if (!add_map_to_molecule()) return 0;
00305 
00306   Molecule *mol = app->moleculeList->mol_from_id(molid);
00307   FileSpec spec;
00308   spec.nvolsets = 1;        // one volumeset to write
00309   spec.setids = new int[1]; // ID of the volumeset
00310   spec.setids[0] = mol->num_volume_data()-1;
00311 
00312   if (!app->molecule_savetrajectory(molid, filename, "dx", &spec)) {
00313     msgErr << "Couldn't write dx file!" << sendmsg;
00314   }
00315 
00316   return 1;
00317 }
00318 
00319 
00320 // Add volumetric data to the molecule
00321 int VolMapCreateILS::add_map_to_molecule() {
00322   Molecule *mol = app->moleculeList->mol_from_id(molid);
00323   float origin[3], xaxis[3], yaxis[3], zaxis[3];
00324   int i;
00325   for (i=0; i<3; i++) {
00326     origin[i] = (float) volmap->origin[i];
00327     xaxis[i] = (float) volmap->xaxis[i];
00328     yaxis[i] = (float) volmap->yaxis[i];
00329     zaxis[i] = (float) volmap->zaxis[i];
00330   }
00331   float *voldata = volmap->data;
00332 
00333   int err = app->molecule_add_volumetric(molid, 
00334               volmap->name, origin, xaxis, yaxis, zaxis,
00335               volmap->xsize, volmap->ysize, volmap->zsize,
00336               voldata);
00337   if (err != 1) {
00338     msgErr << "ERROR: Adding volmap " << mol->num_volume_data()-1
00339            << " to molecule " << molid << " was unsuccessful!"
00340            << sendmsg;
00341     return 0;
00342   }
00343 
00344   // Avoid data being deleted by volmap's destructor
00345   // since it is now owned by the molecule and will be
00346   // freed by its destructor .
00347   volmap->data = NULL;
00348 
00349   msgInfo << "Added volmap " << mol->num_volume_data()-1
00350           << " to molecule " << molid << "." << sendmsg;
00351 
00352   return 1;
00353 }
00354 
00355 // Set probe coordinates,charges and VDW parameters.
00356 // Currently we only support up to 2 probe atoms but
00357 // we should change that later.
00358 void VolMapCreateILS::set_probe(int numatoms, int num_conf,
00359                                 const float *pcoords,
00360                                 const float *vdwrmin,
00361                                 const float *vdweps,
00362                                 const float *charge) {
00363   if (numatoms<=0) return;
00364   if (numatoms==1 && num_conf>0) num_conf = 0;
00365 
00366   conformer_freq  = num_conf;
00367   num_probe_atoms = numatoms;
00368   probe_coords = new float[3*numatoms];
00369   probe_vdw    = new float[2*numatoms];
00370   probe_charge = new float[numatoms];
00371 
00372   // Thermopdynamic beta = 1/kT factor (k = Boltzmann const.)
00373   // 1/k [kcal/mol] = 1.0/(1.38066*6.022/4184) = 503.2206
00374   const float beta = 503.2206f/temperature;
00375 
00376   // The combination rules for VDW parameters of
00377   // two atoms are:
00378   //
00379   // eps(i,j)  = sqrt(eps(i) * eps(j))
00380   // rmin(i,j) = rmin(i)/2 + rmin(j)/2
00381   //
00382   // where rmin(i,j) is the distance of the two atoms at the
00383   // energy minimum. The parameters are provided from the
00384   // user interface in form the contribution of each atom
00385   // rmin(i)/2 (which can be interpreted as the VDW radius).
00386 
00387   // We take the sqrt(eps) here already so that during
00388   // the computation we merely have to multiply the two
00389   // parameters.
00390   // We are computing the occupancy
00391   //   rho = sum_i exp(-1/kT * U[i]) 
00392   //     U = eps*[(rmin/r)^12 - 2*(rmin/r)^6]
00393   // We include beta = 1/kT into the probe eps parameter
00394   // for speed.
00395   int i;
00396   for (i=0; i<num_probe_atoms; i++) {
00397     probe_vdw[2*i  ] = beta*sqrtf(-vdweps[i]);
00398     probe_vdw[2*i+1] = vdwrmin[i];
00399   }
00400 
00401   if (charge) {
00402     memcpy(probe_charge, charge, numatoms*sizeof(float));
00403   } else {
00404     memset(probe_charge, 0, numatoms*sizeof(float));
00405   }
00406 
00407   // Get geometric center:
00408   float cgeom[3];
00409   vec_zero(cgeom);
00410 
00411   for (i=0; i<num_probe_atoms; i++) {
00412     vec_add(cgeom, cgeom, &pcoords[3*i]);
00413   }
00414   vec_scale(cgeom, 1.f/(float)num_probe_atoms, cgeom);
00415 
00416   // Shift geometric center to origin:
00417   for (i=0; i<num_probe_atoms; i++) {
00418     vec_sub(&probe_coords[3*i], &pcoords[3*i], cgeom);
00419   }
00420 }
00421 
00422 
00423 // Set the two highest symmetry axes for the probe and a flag
00424 // telling if we have a tetrahedral symmetry.
00425 void VolMapCreateILS::set_probe_symmetry(
00426         int order1, const float *axis1,
00427         int order2, const float *axis2,
00428         int tetrahedral) {
00429   probe_tetrahedralsymm = tetrahedral;
00430 
00431   if (axis1 && order1>1) {
00432     probe_axisorder1 = order1;
00433     vec_copy(probe_symmaxis1, axis1);
00434     vec_normalize(probe_symmaxis1);
00435   }
00436   if (axis2 && order2>1) {
00437     probe_axisorder2 = order2;
00438     vec_copy(probe_symmaxis2, axis2);
00439     vec_normalize(probe_symmaxis2);
00440   }
00441   if (!tetrahedral && probe_axisorder1>1 && probe_axisorder2>1 &&
00442       dot_prod(probe_symmaxis1, probe_symmaxis2) > 0.05) {
00443     // Axes not orthogonal, drop lower order axis
00444     if (probe_axisorder1<probe_axisorder2) {
00445       probe_axisorder1 = probe_axisorder2;
00446     }
00447     probe_axisorder2 = 1;
00448   }
00449 }
00450 
00451 
00452 // Set minmax coordinates of rectangular grid bounding box
00453 void VolMapCreateILS::set_minmax (float minx, float miny, float minz,
00454                                   float maxx, float maxy, float maxz) {
00455   minmax[0] = minx;
00456   minmax[1] = miny;
00457   minmax[2] = minz;
00458   minmax[3] = maxx;
00459   minmax[4] = maxy;
00460   minmax[5] = maxz;
00461 }
00462 
00463 
00464 // Request PBC aware computation.
00465 // Optionally one can set the PBC cell center
00466 // (default is {0 0 0}).
00467 // If bbox is TRUE then the bounding box of the grid will
00468 // be chosen as the bounding box of the (possibly rhombic)
00469 // PBC cell.
00470 void VolMapCreateILS::set_pbc(float center[3], int bbox) {
00471   pbc    = 1;
00472   pbcbox = bbox;
00473   if (center) vec_copy(pbccenter, center);
00474 }
00475 
00476 
00477 // Set maximum energy considered in the calculation.
00478 void VolMapCreateILS::set_maxenergy(float maxenergy) {
00479   max_energy = maxenergy;
00480   if (max_energy>DEFAULT_EXCL_ENERGY) {
00481     max_energy = DEFAULT_EXCL_ENERGY;
00482   }
00483 }
00484 
00485 
00486 // Set selection to be used for alignment.
00487 void VolMapCreateILS::set_alignsel(AtomSel *asel) {
00488   if (asel) alignsel = asel;
00489 }
00490 
00491 
00492 // Set transformation matrix that was used for the
00493 // alignment of the first frame.
00494 void VolMapCreateILS::set_transform(const Matrix4 *mat) {
00495   transform = *mat;
00496 }
00497 
00498 
00499 // Set grid dimensions to the minmax coordinates and
00500 // pad the positive side of each dimension.
00501 //
00502 //        lattice
00503 //   +---+---+---+---+
00504 //  _|___|___|___|_  |  _
00505 // | |   |   |   | | |  |
00506 // | o---o---o---o---+  |
00507 // | |   |   |   | | |  |
00508 // | |   |   |   | | |  |
00509 // | o---o---o---o---+  |minmax
00510 // |_|_  |   |   | | |  |
00511 // | | | |   |   | | |  |
00512 // | o---o---o---o---+  |
00513 // |___|___________|    L
00514 //  
00515 //   |---|       
00516 //   delta
00517 //               
00518 //   |-----------|
00519 //   xaxis
00520 //
00521 int VolMapCreateILS::set_grid() {
00522   // Number of samples in the final downsampled map.
00523   int nx = (int) ceilf((minmax[3]-minmax[0])/(delta*nsubsamp));
00524   int ny = (int) ceilf((minmax[4]-minmax[1])/(delta*nsubsamp));
00525   int nz = (int) ceilf((minmax[5]-minmax[2])/(delta*nsubsamp));
00526 
00527   // Number of samples used during computation.
00528   nsampx = nx*nsubsamp;
00529   nsampy = ny*nsubsamp;
00530   nsampz = nz*nsubsamp;
00531 
00532   // For volumetric maps the origin is the center of the
00533   // first cell (i.e. the location of the first sample)
00534   // rather than the lower corner of the minmax box.
00535 
00536   // Origin for final downsampled map
00537   float origin[3];
00538   origin[0] = minmax[0] + 0.5f*delta*nsubsamp;
00539   origin[1] = minmax[1] + 0.5f*delta*nsubsamp;
00540   origin[2] = minmax[2] + 0.5f*delta*nsubsamp;
00541 
00542   // Origin for the highres map used during computation
00543   gridorigin[0] = minmax[0] + 0.5f*delta;
00544   gridorigin[1] = minmax[1] + 0.5f*delta;
00545   gridorigin[2] = minmax[2] + 0.5f*delta;
00546   
00547 
00548   // Cell spanning vectors,
00549   // (delta is the distance between two samples)
00550   float cellx[3], celly[3], cellz[3];
00551   cellx[0] = delta*nsubsamp;
00552   cellx[1] = 0.f;
00553   cellx[2] = 0.f;
00554   celly[0] = 0.f;
00555   celly[1] = delta*nsubsamp;
00556   celly[2] = 0.f;
00557   cellz[0] = 0.f;
00558   cellz[1] = 0.f;
00559   cellz[2] = delta*nsubsamp;
00560 
00561   // The axes that span the whole lattice i.e. the vector
00562   // between the first and the last sample point in each
00563   // dimension.
00564   float xaxis[3], yaxis[3], zaxis[3];
00565   int i;
00566   for (i=0; i<3; i++) {
00567     xaxis[i] = cellx[i]*(nx-1);
00568     yaxis[i] = celly[i]*(ny-1);
00569     zaxis[i] = cellz[i]*(nz-1);
00570   }
00571 
00572   // Initialize the final downsampled map:
00573   float *data = new float[nx*ny*nz];
00574   if (maskonly) {
00575     // Fill mask map with ones
00576     for (i=0; i<nx*ny*nz; i++) {
00577       data[i] = 1.f;
00578     }    
00579   } else {
00580     memset(data, 0, nx*ny*nz*sizeof(float));
00581   }
00582 
00583   volmap = new VolumetricData("\0", origin,
00584                   xaxis, yaxis, zaxis, nx, ny, nz, data);
00585 
00586   if (!volmap->gridsize())
00587     return MEASURE_ERR_ZEROGRIDSIZE;
00588 
00589   return 0;
00590 }
00591 
00592 
00593 // Initialize the ILS calculation
00594 int VolMapCreateILS::initialize() {
00595   if (!app->molecule_numframes(molid))
00596     return MEASURE_ERR_NOFRAMES;
00597 
00598   msgInfo << "\n-- Implicit Ligand Sampling --\n\n"
00599           << "This computes the potential of mean force (free energy) over\n"
00600           << "a 3-D grid, for a small ligand.\n\n" << sendmsg;
00601     //    << "If you use this method in your work, please cite:\n\n"
00602     //    << "  J COHEN, A ARKHIPOV, R BRAUN and K SCHULTEN, \"Imaging the\n"
00603     //    << "  migration pathways for O2, CO, NO, and Xe inside myoglobin\".\n"
00604     //    << "  Biophysical Journal 91:1844-1857, 2006.\n\n" << sendmsg;
00605   char tmpstr[256] = { 0 };
00606   sprintf(tmpstr, "Temperature:    %g K", temperature);
00607   msgInfo << tmpstr << sendmsg;
00608   sprintf(tmpstr, "Energy cutoff:  %g kT", max_energy);
00609   msgInfo << tmpstr << sendmsg;
00610 
00611   // Compute electrostatics only if probe charges are nonzero
00612   if (compute_elec && probe_charge[0]==0.0 && probe_charge[1]==0.0) {
00613     compute_elec = 0;
00614   }
00615   msgInfo << "Electrostatics: ";
00616   if (compute_elec) msgInfo << "on"  << sendmsg;
00617   else              msgInfo << "off" << sendmsg;
00618 
00619   sprintf(tmpstr, "Map resolution: %g Angstrom", delta);  
00620   msgInfo << tmpstr << sendmsg;
00621   sprintf(tmpstr, "Subgrid res:    %d points", nsubsamp);
00622   msgInfo << tmpstr << sendmsg;
00623     
00624   // Initialize parameters
00625   if (num_probe_atoms==1) {
00626     sprintf(tmpstr, "VDW cutoff:     %6.3f Angstrom", cutoff);
00627     msgInfo << tmpstr << sendmsg;
00628   }
00629   else {
00630     // Get the max probe radius, i.e. the largest distance
00631     // of an atom to the center of mass
00632 
00633     float max_proberad = 0.f;
00634     int i;
00635     for (i=0; i<num_probe_atoms; i++) {
00636       float dist = norm(&probe_coords[3*i]);
00637       if (dist>max_proberad) max_proberad = dist;
00638     }
00639 
00640     sprintf(tmpstr, "VDW cutoff:     %6.3f Angstrom (%6.3f + probe radius)",
00641             cutoff+max_proberad, cutoff);
00642     msgInfo << tmpstr << sendmsg;
00643 
00644     extcutoff = cutoff + max_proberad;
00645   }
00646 
00647   if (alignsel) {
00648     // Get coordinates of the alignment reference frame (the current frame of alignsel)
00649     alignrefpos = alignsel->coordinates(app->moleculeList);
00650 
00651   } else if (!alignsel && last-first>1) {
00652     msgWarn << sendmsg;
00653     msgWarn << "Use of periodic boundaries requested (-pbc) but" << sendmsg;
00654     msgWarn << "no alignment matrix specified (you didn't use -alignsel)." << sendmsg;
00655     msgWarn << "Have you aligned your structure prior to this calculation?" << sendmsg;
00656     msgWarn << "Hopefully not, since it will have totally messed" << sendmsg;
00657     msgWarn << "up the definition of your PBC cells. Instead you" << sendmsg;
00658     msgWarn << "should use the -alignsel option and let volmap handle" << sendmsg;
00659     msgWarn << "the alignment." << sendmsg;
00660   }
00661 
00662   if (pbc && pbcbox) {
00663     // Compute minmax based on the PBC cell of the first frame:
00664     // Get the smallest rectangular box that encloses the
00665     // entire (possibly nonorthogonal) PBC cell
00666     int err = compute_pbcminmax(app->moleculeList, molid, first,
00667                    pbccenter, &transform, &minmax[0], &minmax[3]);
00668     if (err) return err;
00669   }
00670 
00671   sprintf(tmpstr, "{%g %g %g} {%g %g %g}\n", minmax[0], minmax[1], minmax[2], minmax[3], minmax[4], minmax[5]);
00672   msgInfo << "Grid minmax = " << tmpstr << sendmsg; 
00673 
00674   // Automatically add the force field cutoff to the system
00675   float ffcutoff[3];
00676   ffcutoff[0] = ffcutoff[1] = ffcutoff[2] = extcutoff;
00677   float cutminmax[6];
00678   vec_sub(cutminmax,   minmax,   ffcutoff);
00679   vec_add(cutminmax+3, minmax+3, ffcutoff);
00680   
00681   if (!box_inside_pbccell(first, minmax)) {
00682     if (!pbc) {
00683       msgWarn << sendmsg;
00684       msgWarn << "Selected grid exceeds periodic cell boundaries." << sendmsg;
00685       msgWarn << "Parts of the map will be undefined!" << sendmsg;
00686       msgWarn << "(Consider using -pbc flag to ensure map is valid over entire grid.)"
00687               << sendmsg << sendmsg;
00688     } else {
00689       msgInfo << "Selected grid exceeds periodic cell boundaries." << sendmsg;
00690       msgInfo << "Using PBC image atoms to ensure map is valid over entire grid."
00691               << sendmsg;
00692     }
00693   } else if (!box_inside_pbccell(first, cutminmax)) {
00694     if (!pbc) {
00695       msgWarn << sendmsg;
00696       msgWarn << "Nonbonded interaction cutoff region needed around" << sendmsg;
00697       msgWarn << "selected grid exceeds periodic cell boundaries." << sendmsg;
00698       msgWarn << "Parts of the map will be ill-defined!" << sendmsg;
00699       msgWarn << "(Consider using -pbc flag to ensure map is valid over entire grid.)"
00700               << sendmsg << sendmsg;;
00701     } else {
00702       msgInfo << "Nonbonded interaction cutoff region needed around" << sendmsg;
00703       msgInfo << "selected grid exceeds periodic cell boundaries." << sendmsg;
00704       msgInfo << "Using PBC image atoms to ensure map is valid over entire grid."
00705               << sendmsg;
00706     }
00707   }
00708 
00709   // Compute and set the grid dimensions
00710   set_grid();
00711   msgInfo << "Grid origin = {"
00712           << volmap->origin[0] << " "
00713           << volmap->origin[1] << " "
00714           << volmap->origin[2] << "}" << sendmsg;
00715 
00716   char tmp[64];
00717   sprintf(tmp, "Grid size   = %dx%dx%d (%.1f MB)",
00718           nsampx, nsampy, nsampz,
00719           sizeof(float)*nsampx*nsampy*nsampz/(1024.*1024.));
00720   msgInfo << tmp << sendmsg;
00721 
00722   if (nsubsamp>1) {
00723     sprintf(tmp, "Downsampling final map to %dx%dx%d  (%.1f MB)",
00724             volmap->xsize, volmap->ysize, volmap->zsize,
00725             sizeof(float)*volmap->gridsize()/(1024.*1024.));
00726     msgInfo << tmp << sendmsg;
00727   }
00728 
00729   Molecule *mol = app->moleculeList->mol_from_id(molid);
00730   num_atoms = mol->nAtoms;
00731 
00732   msgInfo << "Global transformation for all frames:" << sendmsg;
00733   print_Matrix4(&transform);
00734 
00735   if (alignsel) msgInfo << "Aligning all frames to the first one." << sendmsg;
00736   else          msgInfo << "Assuming all frames are aligned." << sendmsg;
00737 
00738   if (maskonly) {
00739     msgInfo << sendmsg << "Masking mode:" << sendmsg
00740             << "Generating a mask map containing 1 for grid points that" << sendmsg
00741             << "have a valid contribution from each frame and 0 otherwise."
00742             << sendmsg << sendmsg;
00743     return MEASURE_NOERR;
00744   }
00745 
00746 
00747   // Create list of unique VDW parameters which can be accessed
00748   // through an index list.
00749   create_unique_paramlist();
00750 
00751   // Find smallest VDW rmin parameter for all system atoms 
00752   float min_sysrmin = vdwparams[1];
00753   int i;
00754   for (i=1; i<num_unique_types; i++) {
00755     if (vdwparams[2*i+1]<min_sysrmin) min_sysrmin = vdwparams[2*i+1];
00756   }
00757 
00758   // Find largest VDW rmin parameter for all probe atoms
00759   float min_probermin = probe_vdw[1];
00760   for (i=1; i<num_probe_atoms; i++) {
00761     if (probe_vdw[2*i+1]<min_probermin) min_probermin = probe_vdw[2*i+1];
00762   }
00763   
00764 
00765   const float invbeta = temperature/503.2206f;
00766   msgInfo << "Probe with "<<num_probe_atoms<<" atoms:" << sendmsg;
00767   for (i=0; i<num_probe_atoms; i++) {
00768     sprintf(tmpstr, "  atom %d: epsilon = %g, rmin/2 = %g, charge = % .3f",
00769             i, -pow(invbeta*probe_vdw[2*i],2),
00770             probe_vdw[2*i+1], probe_charge[i]);
00771     msgInfo << tmpstr << sendmsg;
00772   }
00773 
00774   // Create conformers for multiatom probes
00775   if (conformer_freq && num_probe_atoms>1) {
00776     initialize_probe();
00777   } else {
00778     msgInfo << "Ignoring orientations for monoatomic probe." << sendmsg;
00779   }
00780 
00781   get_eff_proberadius();
00782 
00783 
00784   // Define a cutoff for identifying obvious atom clashes:
00785   sprintf(tmpstr, "Clash exclusion distance:  %.3f Angstrom", excl_dist);
00786   msgInfo << tmpstr << sendmsg;
00787 
00788 //   volmask = new VolumetricData(volmap->name, volmap->origin,
00789 //                     volmap->xaxis, volmap->yaxis, volmap->zaxis,
00790 //                     volmap->xsize, volmap->ysize, volmap->zsize, NULL);
00791   //volmask = new VolumetricData(*volmap);
00792   //volmask->data = new float[gridsize*sizeof(float)];
00793   //memset(volmask->data, 0, sizeof(float)*gridsize);
00794   //char name[8];
00795   //strcpy(name, "mask");
00796   //volmask->set_name(name);
00797 
00798   return MEASURE_NOERR;
00799 }
00800 
00801 
00802 void VolMapCreateILS::get_eff_proberadius() {
00803   int numconf, numorient, numrot;
00804   float *conf;
00805   if (probe_tetrahedralsymm) {
00806     numconf = gen_conf_tetrahedral(conf, 6, numorient, numrot);
00807   } else {
00808     numconf = gen_conf(conf, 8, numorient, numrot);
00809   }
00810 
00811   int t, i, j, k;
00812   float max_proberad = 0.f;
00813   for (k=0; k<num_probe_atoms; k++) {
00814     float dist = norm(&probe_coords[3*k]);
00815     float rmin = probe_vdw[2*k+1];
00816     if (dist+rmin>max_proberad) max_proberad = dist+rmin;
00817   }
00818 
00819   float stepsize = 0.01f;
00820   float *effrad = new float[num_unique_types];
00821   excl_dist = 999999.f;
00822 
00823   for (t=0; t<num_unique_types; t++) {
00824     float vdweps  = vdwparams[2*t  ];
00825     float vdwrmin = vdwparams[2*t+1];
00826     //printf("vdweps=%f, vdwrmin=%f\n", vdweps, vdwrmin);
00827     float begin = max_proberad + vdwrmin;
00828     int maxnumstep = int(0.5+begin/stepsize);
00829     //printf("maxproberad=%.2f\n", max_proberad);
00830     float Wmin = 0.f;
00831 //    float Ropt = 0.f;
00832 
00833     for (i=0; i<maxnumstep; i++) {
00834       float dist = begin - float(i)*stepsize;
00835       if (dist<=0.0f) break;
00836 
00837       float avgocc = 0.f;
00838       for (j=0; j<numconf; j++) {
00839         float *coor = &conf[3*num_probe_atoms*j];
00840         float u = 0.f;
00841         for (k=0; k<num_probe_atoms; k++) {
00842           float dx = dist-coor[3*k];
00843           float dy = coor[3*k+1];
00844           float dz = coor[3*k+2];
00845           float r2 = dx*dx + dy*dy + dz*dz;
00846           float epsilon = vdweps  * probe_vdw[2*k];
00847           float rmin    = vdwrmin + probe_vdw[2*k+1];
00848           float rm6 = rmin*rmin / r2;
00849           rm6 = rm6 * rm6 * rm6;
00850           u += epsilon * rm6 * (rm6 - 2.f);  // sum vdw contribution
00851           //printf("     u[%d] = %f, r2=%f\n", k, epsilon * rm6 * (rm6 - 2.f), r2);
00852         }
00853         //printf("  u[%d] = %f\n", j, u);
00854 
00855         float occ = expf(-float(u));
00856 
00857         avgocc += occ;
00858       }
00859       avgocc /= float(numconf);
00860       float W = -logf(avgocc);
00861       //printf("dist=%.2f occ=%g; dG=%g\n", dist, avgocc, -logf(avgocc));
00862       if (W<Wmin) { 
00863         Wmin = W;
00864 //        Ropt = dist; 
00865       }
00866       if (W>max_energy+5.f) {
00867         effrad[t] = dist;
00868         break;
00869       }
00870     }
00871 
00872     //printf("effrad[%d]=%.3f Ropt=%.3f, Wmin=%f\n", t, effrad[t], Ropt, Wmin);
00873     if (effrad[t]<excl_dist) excl_dist = effrad[t];
00874   }
00875 
00876   delete [] effrad;
00877   delete [] conf;
00878 }
00879 
00880 
00881 // Check if two vectors are collinear within a given tolerance.
00882 // Assumes that the vectors are normalized.
00883 static bool collinear(const float *axis1, const float *axis2, float tol) {
00884   if (fabs(dot_prod(axis1, axis2)) > 1.f-DEGTORAD(tol)) return 1;
00885   return 0;
00886 }
00887 
00888 
00889 // Check if probe is a linear molecule and returns
00890 // the Cinf axis.
00891 int VolMapCreateILS::is_probe_linear(float *axis) {
00892   if (num_probe_atoms==1) return 0;
00893 
00894   float vec0[3], vec1[3];
00895   vec_sub(vec0, &probe_coords[3], &probe_coords[0]);
00896   vec_copy(axis, vec0);
00897   if (num_probe_atoms==2) return 1;
00898 
00899   float norm0 = norm(vec0);
00900   int i;
00901   for (i=2; i<num_probe_atoms; i++) {
00902     vec_sub(vec1, &probe_coords[3*i], &probe_coords[0]);
00903     float dot = dot_prod(vec0, vec1)/(norm0*norm(vec1));
00904     if (fabs(dot) < 0.95f) return 0;
00905 
00906     if (dot>0.f) vec_add(axis, axis, vec1);
00907     else         vec_sub(axis, axis, vec1);
00908   }
00909   vec_normalize(axis);
00910   return 1;
00911 }
00912 
00913 
00914 // Subdivide a triangle specified by the points pole, eq1 and eq2.
00915 // A freqency of 1 means no partitioning, 2 signifies that each edge
00916 // is subdivided into 2 segments and so on.
00917 // The resulting vertices are returned in v.
00918 static int triangulate(const float *pole, const float *eq1,
00919                         const float *eq2, int freq, float *v) {
00920   if (freq==0) {
00921     vec_copy(v, pole);
00922     return 1;
00923   }
00924 
00925   float meridian[3], parallel[3];
00926   vec_sub(meridian, eq1, pole);
00927   vec_sub(parallel, eq2, eq1);
00928   float mlen = norm(meridian);
00929   float plen = norm(parallel);
00930   vec_normalize(meridian);
00931   vec_normalize(parallel);
00932   int i, k = 0;
00933   for (i=0; i<=freq; i++) {
00934     float latitude = float(i)/float(freq)*mlen;
00935     float p[3], p0[3];
00936     vec_copy(p0, pole);
00937     vec_scaled_add(p0, latitude, meridian);
00938     int j;
00939     for (j=0; j<=i; j++) {
00940       float longitude = float(j)/float(freq)*plen;
00941       vec_copy(p, p0);
00942       vec_scaled_add(p, longitude, parallel);
00943       vec_copy(&v[3*k], p);
00944       k++;
00945     }
00946   }
00947   return k;
00948 }
00949 
00950 
00951 // Generate the 6 vertices of an octahedron
00952 // (or only 3 if the symmetry flag was set)
00953 static void octahedron(float *vertices, int C2symm) {
00954   const float v[] = {1.f, 0.f, 0.f,
00955                      0.f, 1.f, 0.f,
00956                      0.f, 0.f, 1.f};
00957   memcpy(vertices, v, 9*sizeof(float));
00958   if (!C2symm) {
00959     int i;
00960     for (i=0; i<3; i++) {
00961       vec_negate(&vertices[9+3*i], &vertices[3*i]);
00962     }
00963   }
00964 }
00965 
00966 // Generate the 8 vertices of a hexahedron
00967 // (or only 4 if the symmetry flag was set)
00968 static void hexahedron(float *vertices, int C2symm) {
00969   const float v[] = {1.f,  1.f,  1.f,
00970                      1.f, -1.f,  1.f,
00971                      1.f,  1.f, -1.f,
00972                      1.f, -1.f, -1.f};
00973   memcpy(vertices, v, 12*sizeof(float));
00974 
00975   if (!C2symm) {
00976     int i;
00977     for (i=0; i<4; i++) {
00978       vec_negate(&vertices[12+3*i], &vertices[3*i]);
00979     }
00980   }
00981 }
00982 
00983 // Generate normal vectors for the 12 dodecahedral faces
00984 // (or only 6 if the symmetry flag was set) and the 20
00985 // vertices. The vertices are at the same time the faces
00986 // of the icosahedron, the dual of the dodecahedron.
00987 // XXX I know that this takes more space than just listing
00988 //     the hardcoded points...
00989 static void dodecahedron(float *faces, float *vertices, int C2symm) {
00990   // Angle between two faces
00991   const float dihedral = float(180.f-RADTODEG(acos(-1.f/sqrtf(5.f))));
00992 
00993   // Faces:
00994   float x[3];
00995   vec_zero(x); x[0] = 1.f;
00996   vec_copy(&faces[0], x);
00997   // Contruct first point in ring
00998   Matrix4 rot;
00999   rot.rot(dihedral, 'z');
01000   rot.multpoint3d(&faces[0], &faces[3]);
01001   // Get other points by rotation
01002   int i;
01003   for (i=1; i<5; i++) {
01004     rot.identity();
01005     rot.rot(float(i)*72.f, 'x');
01006     rot.multpoint3d(&faces[3], &faces[3+3*i]);
01007   }
01008   if (!C2symm) {
01009     for (i=0; i<6; i++) {
01010       vec_negate(&faces[18+3*i], &faces[3*i]);
01011     }
01012   }
01013 
01014   // Vertices
01015   // center and first ring
01016   for (i=0; i<5; i++) {
01017     vec_copy(&vertices[3*i], &faces[0]);
01018     vec_add(&vertices[3*i], &vertices[3*i], &faces[3*(i+1)]);
01019     vec_add(&vertices[3*i], &vertices[3*i], &faces[3*((i+1)%5+1)]);
01020     vec_normalize(&vertices[3*i]);
01021   }
01022   // second ring
01023   vec_copy(&vertices[3*5], &faces[3*1]);
01024   vec_add(&vertices[3*5], &vertices[3*5], &faces[3*2]);
01025   vec_normalize(&vertices[3*5]);
01026   float cross[3];
01027   cross_prod(cross, &vertices[3*5], &faces[0]);
01028   vec_normalize(cross);
01029   float phi = angle(&vertices[3*5], &vertices[0]);
01030   rot.identity();
01031   rot.rotate_axis(cross, float(DEGTORAD(-phi)));
01032   rot.multpoint3d(&vertices[3*5], &vertices[3*5]);
01033   for (i=1; i<5; i++) {
01034     rot.identity();
01035     rot.rot(float(i)*72.f, 'x');
01036     rot.multpoint3d(&vertices[3*5], &vertices[3*5+3*i]);
01037   }
01038 
01039   // opposite orientations
01040   if (!C2symm) {
01041     for (i=0; i<10; i++) {
01042       vec_negate(&vertices[30+3*i], &vertices[3*i]);
01043     }
01044   }
01045 }
01046 
01047 
01048 // Geodesic triangulation of an icosahedron.
01049 // Parameter freq is the geodesic frequency, the 
01050 // flag C2symm signifies if the molecule is symmetric
01051 // and we can omit one hemisphere.
01052 // Allocates memory for the orientation vectors
01053 // and returns the number of generated orientations.
01054 static int icosahedron_geodesic(float *(&orientations),
01055                                 int C2symm, int freq) {
01056   int i;
01057   float faces[3*12];
01058   float junk[3*20];
01059   dodecahedron(faces, junk, 0);
01060   float meridian[3], parallel[3];
01061 
01062   int numvertex = 0; // number of triangle vertices
01063   for (i=1; i<=freq+1; i++) numvertex += i;
01064   int symmfac = C2symm ? 2 : 1;
01065   int numorient = (10*freq*freq + 2)/symmfac;
01066   orientations = new float[3*numorient];
01067 
01068   // Add pole to array of orientations
01069   vec_copy(&orientations[0], &faces[0]);
01070   int k = 1;
01071 
01072   for (i=0; i<5; i++) {
01073     // First ring
01074     float p0[3], p1[3], p2[3];
01075     vec_sub(parallel, &faces[3+3*((i+1)%5)], &faces[3*(i+1)]);
01076     float edgelen = norm(parallel);
01077     vec_normalize(parallel);
01078     vec_sub(meridian, &faces[3*(i+1)], &faces[0]);
01079     vec_normalize(meridian);
01080     vec_copy(p0, &faces[0]);
01081     vec_copy(p1, &faces[3*(i+1)]);
01082     vec_copy(p2, &faces[3+3*((i+1)%5)]);
01083     vec_scaled_add(p0,  1.f/float(freq)*edgelen, meridian);
01084     vec_scaled_add(p2, -1.f/float(freq)*edgelen, parallel);
01085     triangulate(p0, p1, p2, freq-1, &orientations[3*k]);
01086     k += numvertex-(freq+1);
01087 
01088     // Second ring
01089     vec_sub(meridian, &faces[3*(i+1)], &faces[21+3*((i+3)%5)]);
01090     vec_normalize(meridian);
01091     vec_copy(p0, &faces[21+3*((i+3)%5)]);
01092     vec_copy(p1, &faces[3*(i+1)]);
01093     vec_scaled_add(p0,  1.f/float(freq)*edgelen, meridian);
01094     vec_scaled_add(p1, -1.f/float(freq)*edgelen, meridian);
01095     vec_copy(p2, p1);
01096     vec_scaled_add(p2, float(freq-2)/float(freq)*edgelen, parallel);
01097     triangulate(p0, p1, p2, freq-2, &orientations[3*k]);
01098     k += numvertex-(freq+1)-freq;
01099   } 
01100 
01101   if (!C2symm) {
01102     for (i=0; i<numorient/2; i++) {
01103       vec_negate(&orientations[3*numorient/2+3*i], &orientations[3*i]);
01104     }
01105   }
01106 
01107   return numorient;
01108 }
01109 
01110 float VolMapCreateILS::dimple_depth(float phi) {
01111   int i;
01112   phi = 0.5f*phi;
01113   // Find smallest system atom
01114   float min_syseps  = vdwparams[0];
01115   float min_sysrmin = vdwparams[1];
01116   for (i=1; i<num_unique_types; i++) {
01117     if (vdwparams[2*i+1]<min_sysrmin) {
01118       min_syseps  = vdwparams[2*i  ];
01119       min_sysrmin = vdwparams[2*i+1];
01120     }
01121   }
01122   float maxdepth = 0.f;
01123 
01124   // Check all probe atoms for dimple depth
01125   for (i=0; i<num_probe_atoms; i++) {
01126     float d = norm(&probe_coords[3*i]);
01127     float sphi, cphi;
01128     sincosf(float(DEGTORAD(phi)), &sphi, &cphi);
01129     float a = d*sphi;
01130     float m = a/cphi;
01131     if (phi == 90.f) m = a;
01132     float c = probe_vdw[2*i+1] + min_sysrmin;
01133     if (m>c) {
01134       maxdepth = d;
01135       break;
01136     }
01137     float b = sqrtf(c*c-m*m);
01138     float depth = d + c - d*cosf(float(DEGTORAD(phi))) - b;
01139     //printf("d=%f, rp=%.3f, rs=%.3f, g=%f, b=%f, a=%f, c=%f\n",
01140     //       d, probe_vdw[2*i+1], min_sysrmin, d*cosf(DEGTORAD(phi)), b, m, c);
01141 
01142     float epsilon = min_syseps * probe_vdw[2*i];
01143     float rmin = min_sysrmin + probe_vdw[2*i+1];
01144     // Get energy in dimple (atoms are touching)
01145     float r2 = c*c;
01146     float rm6 = rmin*rmin / r2;
01147     rm6 = rm6 * rm6 * rm6;
01148     float u0 = epsilon * rm6 * (rm6 - 2.f);
01149     // Get energy for outer radius
01150     r2 = (c+depth)*(c+depth);
01151     rm6 = rmin*rmin / r2;
01152     rm6 = rm6 * rm6 * rm6;
01153     float u1 = epsilon * rm6 * (rm6 - 2.f);
01154     // Get energy for outer radius
01155     r2 = (c-depth)*(c-depth);
01156     rm6 = rmin*rmin / r2;
01157     rm6 = rm6 * rm6 * rm6;
01158     float u2 = epsilon * rm6 * (rm6 - 2.f);
01159     float du1 = u1-u0;
01160     float du2 = u2-u0;
01161     printf("phi = %.2f: %d dimple depth = %f = %5.2f%%, dU1 = %fkT = %5.2f%%; dU1 = %fkT = %5.2f%%\n",
01162            phi, i, depth, 100.f*depth/(d+probe_vdw[2*i]), du1, fabs(100.f*du1/u0), du2, fabs(100.f*du2/u0));
01163 
01164     if (depth>maxdepth) maxdepth = depth;
01165   }
01166   return maxdepth;
01167 }
01168 
01169 
01170 // Generate conformers for tetrahedral symmetries.
01171 // Allocates memory for *conform array and returns the number
01172 // of generated conformers.
01173 // numorient:  # symmetry unique orientations
01174 // numrot:     # rotamers per orientation
01175 //
01176 // Rotate around the 4 axes defined by the corners of
01177 // the tetrahedron and its dual (also a tetrahedron).
01178 // XXX:
01179 // This approach exploits the probe symmetry very well
01180 // but for higher frequencies you will start seeing 
01181 // "holes" in the pattern. These holes are in the middle
01182 // of the 12 corners of the cube spanned by the vertices
01183 // of the two dual tetrahedra
01184 // One idea how to fix this would be to generate two
01185 // extra sets of conformations where the basic tetrahedra
01186 // are rotated such that their vertices are in the holes.
01187 int VolMapCreateILS::gen_conf_tetrahedral(float *(&conform),
01188                      int freq, int &numorient, int &numrot) {
01189   // Generate the 4 corners of the tetrahedron
01190   float tetra0[3], tetra1[3], tetra2[3], tetra3[3];
01191   vec_zero(tetra0);
01192   tetra0[0] = 1.f;
01193   Matrix4 rot;
01194   rot.rot(109.47122f, 'z');
01195   rot.multpoint3d(tetra0, tetra1);
01196   rot.identity();
01197   rot.rot(120.f, 'x');
01198   rot.multpoint3d(tetra1, tetra2);
01199   rot.multpoint3d(tetra2, tetra3);
01200 
01201 #if defined(DEBUG)
01202   Molecule *mol = app->moleculeList->mol_from_id(molid);
01203   MoleculeGraphics *gmol = mol->moleculeGraphics();
01204   gmol->use_color(8);
01205   gmol->add_line(tetra0, tetra1, 0, 1);
01206   gmol->add_line(tetra0, tetra2, 0, 1);
01207   gmol->add_line(tetra0, tetra3, 0, 1);
01208   gmol->add_line(tetra1, tetra2, 0, 1);
01209   gmol->add_line(tetra1, tetra3, 0, 1);
01210   gmol->add_line(tetra2, tetra3, 0, 1);
01211 #endif
01212 
01213   // array of probe orientation vectors
01214   float *orientations;
01215   orientations = new float[3*8];
01216 
01217   vec_copy(&orientations[3*0], tetra0);
01218   vec_copy(&orientations[3*1], tetra1);
01219   vec_copy(&orientations[3*2], tetra2);
01220   vec_copy(&orientations[3*3], tetra3);
01221   float face[3];
01222   vec_copy(face, tetra0);
01223   vec_add(face, face, tetra1);
01224   vec_add(face, face, tetra2);
01225   vec_copy(&orientations[3*4], face);
01226   vec_copy(face, tetra0);
01227   vec_add(face, face, tetra1);
01228   vec_add(face, face, tetra3);
01229   vec_copy(&orientations[3*5], face);
01230   vec_copy(face, tetra0);
01231   vec_add(face, face, tetra2);
01232   vec_add(face, face, tetra3);
01233   vec_copy(&orientations[3*6], face);
01234   vec_copy(face, tetra1);
01235   vec_add(face, face, tetra2);
01236   vec_add(face, face, tetra3);
01237   vec_copy(&orientations[3*7], face);
01238 
01239   numrot = freq;  // # rotamers per orientation
01240   numorient = 8;  // the tetrahedron and its dual
01241   int numconf = numorient*numrot-6;
01242   conform = new float[3*num_probe_atoms*numconf];
01243   memset(conform, 0, 3*num_probe_atoms*numconf*sizeof(float));
01244 
01245   int conf = 0;
01246   int i;
01247   for (i=0; i<numorient; i++) {
01248     float *dir = &orientations[3*i];
01249     vec_normalize(dir);
01250 
01251     // Apply rotations around dir
01252     int j;
01253     for (j=0; j<numrot; j++) {
01254       // The non-rotated orientations are equivalent
01255       // for each tetrahedral dual so we skip them.
01256       if (i%4 && j==0) continue;
01257 
01258       float psi = float(j)/float(numrot)*360.f/float(probe_axisorder1);
01259       Matrix4 rot;
01260       if (i>=numorient/2) {
01261         // Flip orientation by 180 deg in order to get the
01262         // dual of the tetrahedron which is also a tetrahedron
01263         // with corners and faces exchanged.
01264         float z[3];
01265         vec_zero(z); z[2]=1.f;
01266         rot.rotate_axis(z, float(DEGTORAD(180.f)));
01267       }
01268 
01269       // Rotate around dir
01270       rot.rotate_axis(dir, float(DEGTORAD(psi)));
01271     
01272       // Apply rotation to all probe atoms
01273       int k;
01274       for (k=0; k<num_probe_atoms; k++) {
01275         rot.multpoint3d(&probe_coords[3*k],
01276                         &conform[3*num_probe_atoms*conf + 3*k]);
01277       }
01278 
01279       conf++;
01280     }
01281   }
01282 
01283   delete [] orientations;
01284 
01285   // Return the number of generated conformers
01286   return numconf;
01287 }
01288 
01289 
01290 // Compute angle that rotates cross(axis, v1) into cross(axis, v2)
01291 // with respect to rotation around axis.
01292 static float signed_angle(const float *axis,
01293                           const float *v1, const float *v2) {
01294   float normaxis[3], cross1[3], cross2[3], cross3[3];
01295   cross_prod(cross1, axis, v1);
01296   cross_prod(cross2, axis, v2);
01297   cross_prod(cross3, v1, v2);
01298   vec_normalize(cross3);
01299   vec_copy(normaxis, axis);
01300   vec_normalize(normaxis);
01301   float phi = angle(cross1, cross2);
01302   if (dot_prod(axis, cross3)<0) {
01303     phi = -phi;
01304   }
01305   return phi;
01306 }
01307 
01308 
01309 // Generate conformers for all non-tetrahedral symmetries.
01310 // Allocates memory for *conform array and returns the number
01311 // of generated conformers.
01312 // numorient:  # symmetry unique orientations
01313 // numrot:     # rotamers per orientation
01314 int VolMapCreateILS::gen_conf(float *(&conform), int freq,
01315                               int &numorient, int &numrot) {
01316   int i;
01317   float *orientations = NULL; // array of probe orientation vectors
01318   int C2symm = (probe_axisorder2==2 ? 1 : 0);
01319   int symmfac = C2symm ? 2 : 1;
01320   float anglespacing = 360.f;
01321 
01322   switch (freq) {
01323   case 1:
01324     numorient = 1;
01325     orientations = new float[3*numorient];
01326     break;
01327   case 2:
01328     // 6 octahedral vertices
01329     numorient = 6/symmfac;
01330     orientations = new float[3*numorient];
01331     octahedron(orientations, C2symm);
01332     anglespacing = 90.f;
01333     break;
01334   case 3:
01335     // 8 hexahedral vertices
01336     numorient = 8/symmfac;
01337     orientations = new float[3*numorient];
01338     hexahedron(orientations, C2symm);
01339     anglespacing = 109.47122f;
01340     break;
01341   case 4:
01342     // 12 dodecahedral faces 
01343     numorient = 12/symmfac;
01344     orientations = new float[3*numorient];
01345     float vertices[3*20]; // junk
01346     dodecahedron(orientations, vertices, C2symm);
01347     anglespacing = float(180.f-RADTODEG(acos(-1.f/sqrtf(5.f))));
01348       break;
01349   case 5:
01350     // 20 dodecahedral vertices
01351     numorient = 20/symmfac;
01352     orientations = new float[3*numorient];
01353     float faces[3*12]; // junk
01354     dodecahedron(faces, orientations, C2symm);
01355     anglespacing = 180.f-138.189685f;
01356     break;
01357   case 6:
01358     // 12 faces and 20 vertices of a dodecahedron
01359     numorient = 32/symmfac;
01360     orientations = new float[3*numorient];
01361     dodecahedron(&orientations[0], &orientations[3*12/symmfac], C2symm);
01362     anglespacing = 37.377380f;
01363     break;
01364   default:
01365     // Triangulate icosahedral faces
01366     freq -= 5;
01367  
01368     numorient = icosahedron_geodesic(orientations, C2symm, freq);
01369 
01370     anglespacing = (180.f-138.189685f)/freq;
01371     break;
01372   }
01373 
01374   // Number of rotamers per orientation
01375   // Chosen such that the rotation stepping angle is as close
01376   // as possible to the angle between two adjacent orientations.
01377   numrot = 1;
01378   if (probe_axisorder1>=0) {
01379     numrot = int(floorf(360.f/probe_axisorder1/anglespacing + 0.5f));
01380   }
01381 
01382   int numconf = numorient*numrot;
01383   //printf("numorient=%d, numrot=%d, numconf=%d, num_probe_atoms=%d\n",numorient,numrot,numconf,num_probe_atoms);
01384   conform = new float[3*num_probe_atoms*numconf];
01385   memset(conform, 0, 3*num_probe_atoms*numconf*sizeof(float));
01386 
01387 #if defined(DEBUG)
01388   Molecule *mol = app->moleculeList->mol_from_id(molid);
01389   MoleculeGraphics *gmol = mol->moleculeGraphics();
01390   for (i=0; i<numorient; i++) {
01391     float dir[3];
01392     vec_scale(dir, 0.8, &orientations[3*i]);
01393     gmol->use_color(7);
01394     if (i==0) gmol->use_color(4);
01395     if (i==1) gmol->use_color(8);
01396     if (i==2) gmol->use_color(9);
01397     gmol->add_sphere(dir, 0.1, 8);
01398   }
01399 #endif
01400 
01401   // Generate conformer coordinates
01402   int conf = 0;
01403   for (i=0; i<numorient; i++) {
01404     float *dir = &orientations[3*i];
01405     vec_normalize(dir);
01406 
01407     float cross[3], x[3], y[3], z[3];
01408     vec_zero(x); x[0]=1.f;
01409     vec_zero(y); y[1]=1.f;
01410     vec_zero(z); z[2]=1.f;
01411     float phi = 0.f;
01412     float theta = 0.f;
01413 
01414     if (!collinear(x, dir, 2.f)) {
01415       // Get rotation axis and angle phi to rotate x-axis into dir
01416       cross_prod(cross, x, dir);
01417       phi = angle(dir, x);
01418       // Get rotation around X so that Y would be in the plane
01419       // spanned by X and dir. (If we have a second symmetry axis
01420       // then this rotates axis2 into that plane because we have
01421       // previously aligned axis2 with Y.)
01422       float cross2[3];
01423       cross_prod(cross2, x, y);
01424       theta = signed_angle(x, cross2, cross);
01425     } else if (dot_prod(x, dir)<0.f) {
01426       // dir and x are antiparallel
01427       phi = 180.f;
01428     }
01429     //printf("dir[%d] = {%.2f %.2f %.2f},  phi=%.2f, theta=%.2f\n", i, dir[0], dir[1], dir[2], phi, theta);
01430 
01431     
01432     // Apply rotations around dir
01433     int j;
01434     for (j=0; j<numrot; j++) {
01435       Matrix4 m;
01436       float psi = float(j)/float(numrot)*360.f/float(probe_axisorder1);
01437 
01438       // Rotate around dir
01439       m.rotate_axis(dir, float(DEGTORAD(psi)));
01440 
01441       // Rotate around X
01442       m.rot(theta, 'x');
01443       
01444       // Tilt X into dir
01445       m.rotate_axis(z, float(DEGTORAD(phi)));
01446 
01447       // Apply rotation to all probe atoms
01448       int k;
01449       for (k=0; k<num_probe_atoms; k++) {
01450         m.multpoint3d(&probe_coords[3*k],
01451                       &conform[3*num_probe_atoms*conf + 3*k]);
01452       }
01453 
01454       conf++;
01455     }
01456   }
01457 
01458   delete [] orientations;
01459 
01460   // Return the number of generated conformers
01461   return numconf;
01462 }
01463 
01464 
01465 // Perform simple symmetry check on the probe molecule.
01466 // Checks if the molecule is linear molecules and if it has an
01467 // additional horizontal symmetry plane (Cinfv vs. Dinfh pointgroup)
01468 // Also recognizes sp3 centers with 4 identical ligands like
01469 // methane (Td pointgroup).
01470 void VolMapCreateILS::check_probe_symmetry() {
01471   float principal[3];
01472   if (is_probe_linear(principal)) {
01473     probe_axisorder1 = -1;
01474     vec_copy(probe_symmaxis1, principal);
01475 
01476     // Check if there is an additional C2 axis, i.e.
01477     // is there an identical image for each atom.
01478     // This means the pointgroup would be Dinfv, 
01479     // otherwise we just have Cinfv.
01480     int Dinfv = 1;
01481     int i, j;
01482     for (i=0; i<num_probe_atoms; i++) {
01483       float image[3];
01484       vec_negate(image, &probe_coords[3*i]);
01485       int match = 1;
01486       for (j=i+1; j<num_probe_atoms; j++) {
01487         if (distance(&probe_coords[3*j], image)>0.05) continue;
01488         if (probe_vdw[2*i  ]!=probe_vdw[2*j  ] ||
01489             probe_vdw[2*i+1]!=probe_vdw[2*j+1] ||
01490             probe_charge[i]!=probe_charge[j]) {
01491           match = 0;
01492           break;
01493         }
01494       }
01495       if (!match) {
01496         Dinfv = 0;
01497         break;
01498       }
01499     }
01500 
01501     if (Dinfv) {
01502       // Construct perpendicular C2 symmetry axis
01503       float v[3]; // helper vector used to construct orthogonal
01504       vec_zero(v); v[0] = 1.f;
01505       if (fabs(dot_prod(probe_symmaxis1, v))>0.95) {
01506         // Almost parallel, choose a different vector
01507         v[0] = 0.f; v[1] = 1.f;
01508       }
01509       float cross[3];
01510       cross_prod(cross, probe_symmaxis1, v);
01511       cross_prod(probe_symmaxis2, cross, probe_symmaxis1);
01512       probe_axisorder2 = 2;
01513     }
01514   }
01515   else if (num_probe_atoms==5) {
01516     // Try a very simple check for tetrahedral symmetry:
01517     // It will recognize molecules with a central atom and
01518     // 4 equivalent atoms in the corners of a tetrahedron
01519     // such as methane.
01520 
01521     // Look for central atom
01522     int i, icenter = -1;
01523     float zero[3];
01524     vec_zero(zero);
01525     for (i=0; i<num_probe_atoms; i++) {
01526       if (distance(&probe_coords[3*i], zero)<0.05) {
01527         icenter = i;
01528         break;
01529       }
01530     }
01531 
01532     if (icenter>=0) {
01533       float corner[12];
01534       float vdweps=0.f, vdwrmin=0.f, charge=0.f, dist=0.f;
01535       int match = 1;
01536       int j = 0;
01537 
01538       // Check if all ligand atom have the same type and
01539       // build a coordinate list.
01540       for (i=0; i<num_probe_atoms; i++) {
01541         if (i==icenter) continue;
01542         if (j==0) {
01543           vdweps  = probe_vdw[2*i  ];
01544           vdwrmin = probe_vdw[2*i+1];
01545           charge  = probe_charge[i];
01546           dist = norm(&probe_coords[3*i]);
01547         }
01548         else if (probe_vdw[2*i  ] != vdweps  ||
01549                  probe_vdw[2*i+1] != vdwrmin ||
01550                  probe_charge[i]  != charge   ||
01551                  norm(&probe_coords[3*i])-dist > 0.05) {
01552           match = 0;
01553           break;
01554         }
01555 
01556         vec_copy(&corner[3*j], &probe_coords[3*i]);
01557         j++;
01558       }
01559 
01560       // Check the tetrahedral angles
01561       if (match &&
01562           angle(&corner[0], &corner[3])-109.47122f < 5.f &&
01563           angle(&corner[0], &corner[6])-109.47122f < 5.f &&
01564           angle(&corner[0], &corner[9])-109.47122f < 5.f &&
01565           angle(&corner[3], &corner[6])-109.47122f < 5.f &&
01566           angle(&corner[3], &corner[9])-109.47122f < 5.f &&
01567           angle(&corner[6], &corner[9])-109.47122f < 5.f) {
01568         probe_tetrahedralsymm = 1;
01569         probe_axisorder1 = 3;
01570         probe_axisorder2 = 3;
01571         vec_copy(probe_symmaxis1, &corner[0]);
01572         vec_copy(probe_symmaxis2, &corner[3]);
01573       }
01574     }
01575   }
01576 }
01577 
01578 
01579 // Generate probe conformations with uniformly distributed
01580 // orientations and rotations around the orientation vector
01581 // and store the resulting probe coordinates in *conformers.
01582 void VolMapCreateILS::initialize_probe() {
01583   // Perform simple check on probe symmetry and determine
01584   // symmetry axes and order.
01585   check_probe_symmetry();
01586 
01587   // We can make use of up to two symmetry axes in two
01588   // independent operations: The orientation of the probe's
01589   // principal axis and the rotation around it.
01590   // If only one axis was found or specified we will use it to
01591   // exploit symmetry during rotation of the probe around the
01592   // orientation vectors.
01593   // In case we have an additional (orthogonal) 2-fold rotary
01594   // axis we can omit half of the orientations. The idea is that
01595   // if a 180 deg rotation turns the molecule into an identical
01596   // image then we don't have to generate the conformer corresponding
01597   // to the orientation vector pointing in the opposite direction.
01598   // Actually this applies only to linear symmetric molecules like
01599   // oxygen, but since for each orientation the molecule will also be
01600   // rotated around the orientation vector we can extend the concept
01601   // to cases where such an additional rotation turns the flipped
01602   // molecule into the identical image. Strictly, this rotation
01603   // should correspond to one of the generated rotamers but because
01604   // the phase for generating the rotamers was chosen arbitrarily
01605   // anyway we don't need this correspondence.
01606   //
01607   // probe_axisorder1: for probe rotation
01608   // probe_axisorder2: for probe orientation
01609   if (probe_axisorder1==1 || 
01610       (probe_axisorder2!=1 && probe_axisorder2!=2)) {
01611     // Swap the axes
01612     int tmpord = probe_axisorder1;
01613     probe_axisorder1 = probe_axisorder2;
01614     probe_axisorder2 = tmpord;
01615 
01616     float tmp[3];
01617     vec_copy(tmp, probe_symmaxis1);
01618     vec_copy(probe_symmaxis1, probe_symmaxis2);
01619     vec_copy(probe_symmaxis2, tmp);
01620   }
01621 
01622 
01623   // Rotate the probe so that symmetry axis 1 is along X
01624   Matrix4 rot;
01625   rot.transvecinv(probe_symmaxis1[0], probe_symmaxis1[1], probe_symmaxis1[2]);
01626   int i;
01627   for (i=0; i<num_probe_atoms; i++) {
01628     rot.multpoint3d(&probe_coords[3*i], &probe_coords[3*i]);
01629   }
01630   rot.multpoint3d(probe_symmaxis1, probe_symmaxis1);
01631   rot.multpoint3d(probe_symmaxis2, probe_symmaxis2);
01632 
01633   // Rotate the probe so that symmetry axis 2 is along Y
01634   if (probe_axisorder2>1) {
01635     float cross[3], z[3];
01636     vec_zero(z); z[2] = 1.f;
01637     cross_prod(cross, probe_symmaxis1, probe_symmaxis2);
01638 
01639     float phi = angle(cross, z);
01640     rot.identity();
01641     rot.rotate_axis(probe_symmaxis1, float(DEGTORAD(phi)));
01642 
01643     for (i=0; i<num_probe_atoms; i++) {
01644       rot.multpoint3d(&probe_coords[3*i], &probe_coords[3*i]);
01645     }
01646     rot.multpoint3d(probe_symmaxis1, probe_symmaxis1);
01647     rot.multpoint3d(probe_symmaxis2, probe_symmaxis2);
01648   }
01649 
01650 
01651   if (getenv("VMDILSNOSYMM")) {
01652     // Omit any symmetry in the probe conformer generation
01653     // (useful for benchmarking).
01654     probe_axisorder1 = 1;
01655     probe_axisorder2 = 1;
01656     probe_tetrahedralsymm = 0;
01657     msgWarn << "env(VMDILSNOSYMM) is set: Ignoring probe symmetry!" << sendmsg;
01658   }
01659 
01660   num_orientations = 0;  // # symmetry unique orientations
01661   num_rotations    = 1;  // # rotamers per orientation
01662 
01663   if (probe_tetrahedralsymm) {
01664     msgInfo << "Probe symmetry: tetrahedral" << sendmsg;
01665 
01666     num_conformers = gen_conf_tetrahedral(conformers,
01667                          conformer_freq, num_orientations,
01668                          num_rotations);
01669   }
01670 
01671   else {
01672     // Probe is not tetrahedral, generate geodesic orientations
01673     // based on platonic solids.
01674 
01675     if (probe_axisorder1<=1 && probe_axisorder1<=1) {
01676       msgInfo << "Probe symmetry: none" << sendmsg;
01677     }
01678 
01679     else if (probe_axisorder1==-1) {
01680       if (probe_axisorder2==2) {
01681         msgInfo << "Probe symmetry: Dinfh (linear, C2)" << sendmsg;
01682       } else {
01683         msgInfo << "Probe symmetry: Cinfv (linear)" << sendmsg;
01684       }
01685       msgInfo << "  Probe is linear, generating probe orientations only." << sendmsg;
01686     }
01687     else if (probe_axisorder1>1) {
01688       if (probe_axisorder2==2) {
01689         msgInfo << "Probe symmetry: C" << probe_axisorder1
01690                 << ", C2" << sendmsg;
01691       } else {
01692         msgInfo << "Probe symmetry: C" << probe_axisorder1 << sendmsg;
01693       }
01694       msgInfo << "  Exploiting C" << probe_axisorder1
01695               << " rotary symmetry for rotation of the oriented probe." << sendmsg;
01696     }
01697 
01698     if (probe_axisorder2==2) {
01699       msgInfo << "  Exploiting C2 rotary symmetry for probe orientations" 
01700               << sendmsg;
01701     }
01702 
01703     // dimple_depth(180.f); // single orientation
01704     // dimple_depth(90.f); // hexahedron
01705     // dimple_depth(180.f-109.47122f); // octahedron
01706     // dimple_depth(180.f-116.56557f); // dodecahedron
01707     // dimple_depth(180.f-138.18966f); // icosahedron
01708     // dimple_depth(25.f); // icosahedron faces+vertices
01709 
01710     num_conformers = gen_conf(conformers, conformer_freq,
01711                               num_orientations, num_rotations);
01712   }
01713 
01714   msgInfo << "Probe orientations:       " << num_orientations
01715           << sendmsg;
01716   msgInfo << "Rotamers per orientation: " << num_rotations
01717           << sendmsg;
01718   msgInfo << "Conformers generated:     " << num_conformers
01719           << sendmsg << sendmsg;
01720 }
01721 
01722 
01723 // Check if the box given by the minmax coordinates is located
01724 // entirely inside the PBC unit cell of the given frame and in
01725 // this case return 1, otherwise return 0. 
01726 int VolMapCreateILS::box_inside_pbccell(int frame, float *minmaxcoor) {
01727   Matrix4 mat;
01728 
01729   // Get the PBC --> orthonormal unitcell transformation
01730   measure_pbc2onc(app->moleculeList, molid,
01731                   frame, pbccenter, mat);
01732 
01733   // Get vectors describing the box edges
01734   float box[9];
01735   memset(box, 0, 9*sizeof(float));
01736   box[0] = minmaxcoor[3]-minmaxcoor[0];
01737   box[4] = minmaxcoor[4]-minmaxcoor[1];
01738   box[8] = minmaxcoor[5]-minmaxcoor[2];
01739   // printf("box = {%g %g %g}\n", box[0], box[1], box[2]);
01740   // printf("box = {%g %g %g}\n", box[3], box[4], box[5]);
01741   // printf("box = {%g %g %g}\n", box[6], box[7], box[8]);
01742 
01743   // Create coordinates for the 8 corners of the box
01744   // and transform them into the system of the orthonormal
01745   // PBC unit cell.
01746   float node[8*3];
01747   memset(node, 0, 8*3*sizeof(float));
01748   int n=0;
01749   int i, j, k;
01750   for (i=0; i<=1; i++) {
01751     for (j=0; j<=1; j++) {
01752       for (k=0; k<=1; k++) {
01753         vec_copy(node+3*n, &minmaxcoor[0]);
01754         vec_scaled_add(node+3*n, float(i), &box[0]);
01755         vec_scaled_add(node+3*n, float(j), &box[3]);
01756         vec_scaled_add(node+3*n, float(k), &box[6]);
01757         // Apply the PBC --> orthonormal unitcell transformation
01758         // to the current test point.
01759         mat.multpoint3d(node+3*n, node+3*n);
01760         n++;
01761       }
01762     }
01763   }
01764 
01765   // Check if corners lie inside the orthonormal unitcell
01766   for (n=0; n<8; n++) {
01767     //printf("node[%i] = {%g %g %g}\n", n, node[3*n], node[3*n+1], node[3*n+2]);
01768     if (node[3*n  ]<0.f) return 0;
01769     if (node[3*n+1]<0.f) return 0;
01770     if (node[3*n+2]<0.f) return 0;
01771     if (node[3*n  ]>1.f) return 0;
01772     if (node[3*n+1]>1.f) return 0;
01773     if (node[3*n+2]>1.f) return 0;
01774   }
01775   return 1;
01776 }
01777 
01778 
01779 // Check if the entire volmap grid is located entirely inside
01780 // the PBC unit cell of the given frame (taking the alignment
01781 // into account) and in this case return 1, otherwise return 0.
01782 // Also sets the the gridpoints of volmap that are outside the
01783 // PBC cell to zero (used in the maskonly mode).
01784 int VolMapCreateILS::grid_inside_pbccell(int frame, float *maskvoldata, 
01785                                          const Matrix4 &alignment) {
01786   Matrix4 AA, BB, CC;
01787 
01788   Molecule *mol = app->moleculeList->mol_from_id(molid);
01789   mol->get_frame(frame)->get_transforms(AA, BB, CC);
01790 
01791   // Construct the cell spanning vectors
01792   float cella[3], cellb[3], cellc[3];
01793   cella[0] = AA.mat[12];
01794   cella[1] = AA.mat[13];
01795   cella[2] = AA.mat[14];
01796   cellb[0] = BB.mat[12];
01797   cellb[1] = BB.mat[13];
01798   cellb[2] = BB.mat[14];
01799   cellc[0] = CC.mat[12];
01800   cellc[1] = CC.mat[13];
01801   cellc[2] = CC.mat[14];
01802   // Construct the normals of the 6 cell boundary planes
01803   float normals[3*6];
01804   cross_prod(&normals[0], cella, cellb);
01805   cross_prod(&normals[3], cella, cellc);
01806   cross_prod(&normals[6], cellb, cellc);
01807   vec_normalize(&normals[0]);
01808   vec_normalize(&normals[3]);
01809   vec_normalize(&normals[6]);
01810   vec_scale(&normals[0], cutoff, &normals[0]);
01811   vec_scale(&normals[3], cutoff, &normals[3]);
01812   vec_scale(&normals[6], cutoff, &normals[6]);
01813   vec_negate(&normals[9],  &normals[0]);
01814   vec_negate(&normals[12], &normals[3]);
01815   vec_negate(&normals[15], &normals[6]);
01816 
01817   Matrix4 pbc2onc;
01818   int allinside = 1;
01819 
01820   // Get the PBC --> orthonormal unitcell transformation
01821   measure_pbc2onc(app->moleculeList, molid, frame, pbccenter, pbc2onc);
01822 
01823   // In order to transform a point P into the orthonormal cell (P') it 
01824   // first has to be unaligned (the inverse of the alignment):
01825   // P' = M_norm * (alignment^-1) * P
01826   Matrix4 alignmentinv(alignment);
01827   alignmentinv.inverse();
01828 
01829   Matrix4 coretransform(pbc2onc);
01830   coretransform.multmatrix(alignmentinv);
01831 
01832   float testpos[3], gridpos[3], extgridpos[3];
01833 
01834   int n;
01835   for (n=0; n<nsampx*nsampy*nsampz; n++) {
01836     // Position of grid cell's center
01837     gridpos[0] = float((n%nsampx)         *delta + gridorigin[0]);
01838     gridpos[1] = float(((n/nsampx)%nsampy)*delta + gridorigin[1]);
01839     gridpos[2] = float((n/(nsampx*nsampy))*delta + gridorigin[2]); 
01840 
01841     // Construct 6 test points that are at cutoff distance along
01842     // the 6 cell normal vectors. The closest system boundary
01843     // must lie along one of these normals. If all 6 points are
01844     // within the PBC cell then all possible interaction partner
01845     // will be within the cell, too.
01846     int i;
01847     for (i=0; i<6; i++) {
01848       vec_add(extgridpos, gridpos, &normals[3*i]);
01849       // Project into an orthonormal system that is convenient
01850       // for testing if a point is outside the cell:
01851       coretransform.multpoint3d(extgridpos, testpos);
01852       if (testpos[0]<0.f || testpos[0]>1.f ||
01853           testpos[1]<0.f || testpos[1]>1.f ||
01854           testpos[2]<0.f || testpos[2]>1.f) {
01855         // The test point is outside the PBC cell
01856         maskvoldata[n] = 0.f;
01857         allinside = 0;
01858         i = 6;
01859       }
01860     }
01861   }
01862 
01863   return allinside;
01864 }
01865 
01866 
01875 int VolMapCreateILS::create_unique_paramlist() {
01876   Molecule *mol = app->moleculeList->mol_from_id(molid);
01877 
01878   // typecast pointers to "flint" so that we can do int compare below
01879   const flint *radius = (flint *) mol->extraflt.data("radius");
01880   if (!radius) return MEASURE_ERR_NORADII;
01881   const flint *occupancy = (flint *) mol->extraflt.data("occupancy");
01882   if (!occupancy) return MEASURE_ERR_NORADII;
01883 
01884   int i, j;
01885 
01886 #define MAX_UNIQUE_TYPES  200
01887   // Any sane data set should have no more than about 25 unique types.
01888   // We guard against malformed data by setting an upper bound on the
01889   // number of types, preventing our O(N) search to devolve into O(N^2).
01890 
01891   atomtypes = new int[mol->nAtoms];  // each atom stores its type index
01892   atomtypes[0] = 0;  // first atom is automatically assigned first type
01893   flint *unique_occ = new flint[MAX_UNIQUE_TYPES];
01894   flint *unique_rad = new flint[MAX_UNIQUE_TYPES];
01895   unique_occ[0].f = occupancy[0].f;
01896   unique_rad[0].f = radius[0].f;
01897   num_unique_types = 1;
01898 
01899   for (i=1; i<mol->nAtoms; i++) {
01900     int found = 0;
01901     // Compare VDW params of current atom against all
01902     // existing unique VDW types.
01903     for (j=0; j<num_unique_types; j++) {
01904       // perform == test on ints because it's safer
01905       if (occupancy[i].i==unique_occ[j].i && radius[i].i==unique_rad[j].i) {
01906         found = 1;
01907         break;
01908       }
01909     }
01910     if (!found) {
01911       if (MAX_UNIQUE_TYPES==num_unique_types) {
01912         msgErr << "Exceeded maximum number " << MAX_UNIQUE_TYPES
01913                << " of unique atom parameter types" << sendmsg;
01914         return -1;
01915       }
01916       // No matching VDW type found, create a new one
01917       unique_occ[j].f = occupancy[i].f;
01918       unique_rad[j].f = radius[i].f;
01919       num_unique_types++;
01920     }
01921     atomtypes[i] = j;
01922   }
01923 
01924   vdwparams  = new float[2*num_unique_types];
01925   for (j=0; j<num_unique_types; j++) {
01926     // check validity of VDW parameters
01927     if ( !(unique_occ[j].f <= 0.f && unique_rad[j].f > 0.f) ) {
01928       msgErr << "Found invalid VDW parameters " << j
01929         << ": occupancy=" << unique_occ[j].f
01930         << ", radius=" << unique_rad[j].f
01931         << sendmsg;
01932       return -1;
01933     }
01934     // The combination rule for VDW eps of 2 atoms is:
01935     // eps(i,j)  = sqrt(eps(i)) * sqrt(eps(j))
01936     // so we take the sqrt here already.
01937     vdwparams[2*j  ] = sqrtf(-unique_occ[j].f);
01938     vdwparams[2*j+1] = unique_rad[j].f;
01939   }
01940   delete [] unique_occ;
01941   delete [] unique_rad;
01942 
01943   msgInfo << "Number of atom types: " << num_unique_types << sendmsg;
01944 
01945 #if 0
01946   float *epsilon = new float[mol->nAtoms];
01947   for (i=0; i<mol->nAtoms; i++) {
01948     // The combination rule for VDW eps of 2 atoms is:
01949     // eps(i,j)  = sqrt(eps(i)) * sqrt(eps(j))
01950     // so we take the sqrt here already.
01951     epsilon[i] = sqrtf(-occupancy[i]);
01952   }  
01953 
01954   atomtypes = new int[mol->nAtoms];
01955   atomtypes[0] = 0;
01956   float *unique_eps = new float[mol->nAtoms];
01957   float *unique_rad = new float[mol->nAtoms];
01958   unique_eps[0] = epsilon[0];
01959   unique_rad[0] = radius[0];
01960   num_unique_types = 1;
01961 
01962   for (i=1; i<mol->nAtoms; i++) {
01963     int found = 0;
01964     // Compare VDW params of current atom against all
01965     // existing unique VDW types.
01966     for (j=0; j<num_unique_types; j++) {
01967       if (epsilon[i]==unique_eps[j] && radius[i]==unique_rad[j]) {
01968         found = 1;
01969         break;
01970       }
01971     }
01972     if (!found) {
01973       // No matching VDW type found, create a new one
01974       unique_eps[j] = epsilon[i];
01975       unique_rad[j] = radius[i];
01976       num_unique_types++;
01977     }
01978     atomtypes[i] = j;
01979   }
01980 
01981   vdwparams  = new float[2*num_unique_types];
01982   for (j=0; j<num_unique_types; j++) {
01983     vdwparams[2*j  ] = unique_eps[j];
01984     vdwparams[2*j+1] = unique_rad[j];
01985     //printf("eps=%f, rmin=%f\n", unique_eps[j], unique_rad[j]);
01986   }
01987   delete [] epsilon;
01988   delete [] unique_eps;
01989   delete [] unique_rad;
01990 
01991   msgInfo << "Number of atom types: " << num_unique_types << sendmsg;
01992 #endif
01993 
01994   return 0;
01995 }
01996 
01997 
01998 // Perform ILS calculation for all specified frames.
01999 int VolMapCreateILS::compute() {
02000   int numframes = app->molecule_numframes(molid);
02001   if (first<0 || last>=numframes) return -1;
02002 
02003   int err = initialize();
02004   if (err) return err;
02005   
02006   int n, frame;
02007   int gridsize = nsampx*nsampy*nsampz;
02008   float *frame_voldata = new float[gridsize]; // individual frame voldata
02009   float *combo_voldata = new float[gridsize]; // combo cache voldata
02010   if (maskonly) {
02011     // Fill mask map with ones
02012     for (n=0; n<gridsize; n++) {
02013       combo_voldata[n] =  1.f;
02014       frame_voldata[n] =  1.f;
02015     }    
02016   } else {
02017     memset(combo_voldata, 0, gridsize*sizeof(float));
02018   }
02019 
02020   msgInfo << "Processing frames " << first << "-" << last
02021           << ", " << last-first+1 << " frames in total..." << sendmsg;
02022 
02023   computed_frames = 0;
02024 
02025   wkf_timerhandle timer = wkf_timer_create();
02026   wkf_timerhandle alltimer = wkf_timer_create();
02027   wkf_timer_start(alltimer);
02028 
02029   // Combine frame_voldata into combo_voldata, one frame at a time
02030   for (frame=first; frame<=last; frame++) { 
02031     msgInfo << "ILS frame " << frame-first+1 << "/" << last-first+1;
02032 
02033 #ifdef TIMING
02034     msgInfo << sendmsg;
02035 #else
02036     msgInfo << "   ";
02037 #endif
02038 
02039     wkf_timer_start(timer);
02040 
02041     // Perform the actual ILS calculation for this frame
02042     compute_frame(frame, frame_voldata);
02043 
02044     msgInfo << "Total frame time = " << wkf_timer_timenow(timer) << " s" << sendmsg;
02045 
02046     if (maskonly) {
02047       for (n=0; n<gridsize; n++) {
02048         combo_voldata[n] *= frame_voldata[n];
02049       }
02050     } else {
02051       // For each cell combine occupancies of the new frame with the
02052       // sum of the existing ones (we will divide by the number of
02053       // frames later to get the average).
02054       int numexcl = 0;
02055       for (n=0; n<gridsize; n++) {
02056         combo_voldata[n] += frame_voldata[n];
02057         if (frame_voldata[n]<=min_occup) numexcl++;
02058       }
02059       //printf("numexcl = %d/%d\n", numexcl, gridsize);
02060     }
02061 
02062     computed_frames++;
02063   }
02064 
02065   double allframetime = wkf_timer_timenow(alltimer);
02066 
02067   // Downsampling of the final map
02068   if (nsubsamp>1||1) {
02069     int ndownsampx = volmap->xsize;
02070     int ndownsampy = volmap->ysize;
02071     int ix, iy, iz;
02072 
02073     if (maskonly) {
02074       for (iz=0; iz<nsampz; iz++) {
02075         int isubz = iz/nsubsamp*ndownsampy*ndownsampx;
02076         for (iy=0; iy<nsampy; iy++) {
02077           int isuby = iy/nsubsamp*ndownsampx;
02078           for (ix=0; ix<nsampx; ix++) {
02079             n = iz*nsampy*nsampx + iy*nsampx + ix;
02080             float val = combo_voldata[n];
02081             int m = isubz + isuby + ix/nsubsamp; 
02082             // If any of the subsamples where zero,
02083             // the downsampled voxel will be zero:
02084             volmap->data[m] *= val; 
02085           }
02086         }
02087       }
02088 
02089     } else {
02090       for (iz=0; iz<nsampz; iz++) {
02091         int isubz = iz/nsubsamp*ndownsampy*ndownsampx;
02092         for (iy=0; iy<nsampy; iy++) {
02093           int isuby = iy/nsubsamp*ndownsampx;
02094           for (ix=0; ix<nsampx; ix++) {
02095             n = iz*nsampy*nsampx + iy*nsampx + ix;
02096             float val = combo_voldata[n];
02097             int m = isubz + isuby + ix/nsubsamp; 
02098             //printf("%d: val[%2d,%2d,%2d]=%g -->%d\n", n, ix, iy, iz, val, m);
02099             volmap->data[m] += val;
02100           }
02101         }
02102       }
02103 
02104       // Finally, we have to divide by the number of frames
02105       // and by the number of subsamples.
02106       float nsamppercell = float(nsubsamp*nsubsamp*nsubsamp*computed_frames);
02107       for (n=0; n<volmap->gridsize(); n++) {
02108         volmap->data[n] = volmap->data[n]/nsamppercell;
02109       }
02110     }
02111 
02112     if (!maskonly) {
02113       // Our final maps contain the PMF W which is related to the
02114       // occupancy rho (the probability to find a particle at
02115       // that point) by
02116       //
02117       // W = -ln(rho);  [in units of kT]
02118       //
02119       // Additionally we clamp large energies to the user-provided
02120       // max_energy value.
02121       //
02122       for (n=0; n<volmap->gridsize(); n++) {
02123         float val = volmap->data[n];
02124         if (val<=min_occup) {
02125           volmap->data[n] = max_energy;
02126         } else {
02127           val = -logf(val);
02128           if (val > max_energy) val = max_energy;
02129           volmap->data[n] = val;
02130         }
02131       }
02132     }
02133   }
02134 
02135   delete[] frame_voldata;
02136   delete[] combo_voldata;
02137 
02138   msgInfo << "#################################################"
02139           << sendmsg << sendmsg;
02140   msgInfo << "Total time for all frames = "
02141           << allframetime << " s" << sendmsg;
02142   msgInfo << "Avg time per frame        = " 
02143           << allframetime/(last-first+1) << " s" << sendmsg;
02144   msgInfo << "Downsampling              = "
02145           << wkf_timer_timenow(alltimer)-allframetime << " s"
02146           << sendmsg << sendmsg;
02147   msgInfo << "#################################################"
02148           << sendmsg << sendmsg;
02149 
02150   wkf_timer_destroy(timer);
02151   wkf_timer_destroy(alltimer);
02152 
02153   return 0;
02154 }
02155 
02156 
02157 
02158 // Align current frame to the reference
02159 void VolMapCreateILS::align_frame(Molecule *mol, int frame, float *coords,
02160                                   Matrix4 &alignment) {
02161   // In case alignsel is not NULL we align the current frame to the
02162   // first frame according to the provided selection.
02163   if (alignsel) {
02164     int i;
02165     int save_frame_alignsel = alignsel->which_frame;
02166     
02167     alignsel->which_frame = frame;
02168     alignsel->change(NULL, mol);
02169 
02170     float *weight = new float[alignsel->selected]; 
02171     for (i=0; i<alignsel->selected; i++) weight[i] = 1.0;
02172     
02173     measure_fit(alignsel, alignsel, coords, alignrefpos, weight, NULL, &alignment);
02174     delete[] weight;
02175     
02176     if (!getenv("VMDILSALIGNMAPS")) {
02177       // Do the alignment
02178       // (For the neighboring pbc coordinates the alignment is 
02179       // implicitely done below).
02180       for (i=0; i < mol->nAtoms; i++) { 
02181         alignment.multpoint3d(coords, coords);
02182         coords += 3;
02183       }
02184     }
02185     
02186     alignsel->which_frame = save_frame_alignsel;
02187   }
02188 
02189   // Combine frame alignment with general transform (global alignment)
02190   alignment.multmatrix(transform);
02191 }
02192 
02193 
02194 // Get array of coordinates of selected atoms.
02195 // If the pbc flag was set we also generate coordinates
02196 // for atoms within a cutoff in adjacent pbc image cells.
02197 // The image coordinates can be related to the according atoms
02198 // in the main cell through the indexmap.
02199 int VolMapCreateILS::get_atom_coordinates(int frame, Matrix4 &alignment,
02200                                           int *(&vdwtypes),
02201                                           float *(&coords)) {
02202   wkf_timerhandle timer = wkf_timer_create();
02203   wkf_timer_start(timer);
02204 
02205   // Select all atoms within the extended cutoff of the
02206   // user specified grid minmax box.
02207   int *selon = new int[num_atoms];
02208   memset(selon, 0, num_atoms*sizeof(int));
02209 
02210   float minx = minmax[0]-extcutoff;
02211   float miny = minmax[1]-extcutoff;
02212   float minz = minmax[2]-extcutoff;
02213   float maxx = minmax[3]+extcutoff;
02214   float maxy = minmax[4]+extcutoff;
02215   float maxz = minmax[5]+extcutoff;
02216 
02217   int numselected = 0;
02218   int i;
02219   for (i=0; i<num_atoms; i++) {
02220     float x = coords[3*i  ];
02221     float y = coords[3*i+1];
02222     float z = coords[3*i+2];
02223     if (x>=minx && x<=maxx &&
02224         y>=miny && y<=maxy &&
02225         z>=minz && z<=maxz) {
02226       selon[i] = 1;
02227       numselected++;
02228     }
02229   }
02230 
02231   int numcoords = numselected;
02232 
02233   float *selcoords = NULL;
02234 
02235   // If pbc is set the user requests a PBC aware computation and we
02236   // must generate the extended coordinates, i.e. the positions of
02237   // the atoms in the neighboring pbc cells.
02238   if (pbc) {
02239     // Automatically add the force field cutoff to the system.
02240     float ffcutoff[3];
02241     ffcutoff[0] = ffcutoff[1] = ffcutoff[2] = cutoff;
02242 
02243     // Positions of the atoms in the neighboring pbc cells.
02244     ResizeArray<float> extcoord_array;
02245 
02246     // The indexmap_array contains the index of the according
02247     // unitcell atom for each extended atom.
02248     ResizeArray<int>   indexmap_array;
02249 
02250     // Generate coordinates for atoms in the neighboring cells.
02251     // The indexmap_array tells to which atom the extended coordinate
02252     // corresponds.
02253     // We have to use NULL instead of sel (2nd parameter) in order
02254     // to return all PBC neighbors and not only the ones within
02255     // ffcutoff of sel. The reason is that we have no atomselection
02256     // and are computing the interaction between the all system atoms
02257     // and the probe located at the gridpoints.
02258     measure_pbc_neighbors(app->moleculeList, NULL, molid, frame,
02259                           &alignment, pbccenter, ffcutoff, minmax, 
02260                           &extcoord_array, &indexmap_array);
02261 
02262     numcoords = numselected+indexmap_array.num();
02263 
02264     selcoords = new float[3*numcoords];
02265     vdwtypes  = new int[numcoords];
02266 
02267     int j = numselected;
02268     for (i=0; i<indexmap_array.num(); i++) {
02269       selcoords[3*j  ] = extcoord_array[3*i  ];
02270       selcoords[3*j+1] = extcoord_array[3*i+1];
02271       selcoords[3*j+2] = extcoord_array[3*i+2];
02272       vdwtypes[j] = atomtypes[indexmap_array[i]];
02273       j++;
02274     }
02275 
02276     //printf("volmap: considering %d PBC neighbor atoms.\n", indexmap_array.num());
02277   } else {
02278     selcoords = new float[3*numcoords];
02279     vdwtypes  = new int[numcoords];
02280   }
02281 
02282   // Get the core coordinates (selected atoms in the main PBC cell)
02283   int j=0;
02284   for (i=0; i<num_atoms; i++) { 
02285     if (!selon[i]) continue; //atom is not selected
02286     selcoords[3*j  ] = coords[3*i  ];
02287     selcoords[3*j+1] = coords[3*i+1];
02288     selcoords[3*j+2] = coords[3*i+2];
02289     vdwtypes[j] = atomtypes[i];
02290     j++;
02291   }
02292   //printf("volmap: considering %d core atoms.\n", j);
02293 
02294   coords = selcoords;
02295 
02296   delete [] selon;
02297 
02298 #ifdef TIMING
02299   msgInfo << "Coord setup: " << wkf_timer_timenow(timer) << " s" << sendmsg;
02300 #endif
02301   wkf_timer_destroy(timer);
02302 
02303   return numcoords;
02304 }
02305 
02306 
02307 
02309 //  This is the function driving ILS for each frame    //
02311 
02312 // Computes, for each gridpoint, the VdW energy to the nearest atoms
02313 int VolMapCreateILS::compute_frame(int frame, float *voldata) { 
02314   Matrix4 alignment;
02315   float *coords;
02316 
02317   Molecule *mol = app->moleculeList->mol_from_id(molid);
02318   if (!mol) return -1;
02319 
02320 #ifdef TIMING
02321   char report[128] = { 0 };
02322   wkf_timerhandle timer = wkf_timer_create();
02323 #endif
02324 
02325   // Advance to next frame
02326   coords = mol->get_frame(frame)->pos;
02327 
02328   // In case alignsel is not NULL, align the current frame to the
02329   // first frame according to the provided selection and get the
02330   // alignment transformation matrix.
02331   align_frame(mol, frame, coords, alignment);
02332 
02333   if (maskonly) {
02334 #ifdef TIMING
02335     wkf_timer_start(timer);
02336 #endif
02337 
02338     // We are only creating a mask map that defines the
02339     // gridpoints that (after alignment) still overlap in 
02340     // each frame with the reference grid and thus are not
02341     // undersampled.
02342     grid_inside_pbccell(frame, voldata, alignment);
02343 
02344 #ifdef TIMING
02345     msgInfo << "Masking:     " << wkf_timer_timenow(timer) << " s" << sendmsg;
02346 #endif
02347 
02348     return MEASURE_NOERR;
02349   }
02350 
02351 
02352   // Get array of coordinates of selected atoms and their
02353   // neighbors (within a cutoff) in the PBC images.
02354   // Memory for *vdwtypes and *coords will be allocated.
02355   int *vdwtypes = NULL;
02356   int numcoords;
02357   float originalign[3];
02358   float axesalign[9];
02359   float gridaxes[9];
02360   memset(gridaxes, 0, 9*sizeof(float));
02361   gridaxes[0] = gridaxes[4] = gridaxes[8] = 1.f;
02362 
02363   if (getenv("VMDILSALIGNMAPS")) {
02364     // We use all atoms unaligned, but be align the map instead
02365     numcoords = num_atoms;
02366     vdwtypes = atomtypes;
02367     alignment.multpoint3d(gridorigin, originalign);
02368     alignment.multpoint3d(&gridaxes[0], &axesalign[0]);
02369     alignment.multpoint3d(&gridaxes[3], &axesalign[3]);
02370     alignment.multpoint3d(&gridaxes[6], &axesalign[6]);
02371     msgInfo << "Aligning maps." << sendmsg;
02372   } else {
02373     // Get extended list of aligned atom coordinates
02374     numcoords = get_atom_coordinates(frame, alignment,
02375                                      vdwtypes, coords);
02376     memcpy(originalign, gridorigin, 3*sizeof(float));
02377     memcpy(axesalign,   gridaxes,   9*sizeof(float));
02378     msgInfo << "Aligning frames." << sendmsg;
02379   }
02380   
02381   if (getenv("VMDALLATOMILS")) {
02382 #ifdef TIMING
02383     wkf_timer_start(timer);
02384 #endif
02385     
02386     // Assuming the grid is aligned with the coordinate axes:
02387     float lenx = float(nsampx*delta);
02388     float leny = float(nsampy*delta);
02389     float lenz = float(nsampz*delta);
02390     
02391     compute_allatoms(voldata, nsampx, nsampy, nsampz,
02392                      lenx, leny, lenz, originalign, axesalign,
02393                      alignment.mat, numcoords, coords,
02394                      vdwtypes, vdwparams, cutoff, probe_vdw, num_probe_atoms,
02395                      num_conformers, conformers, max_energy); 
02396     
02397 #ifdef TIMING
02398     sprintf(report, "compute_allatoms()                                     "
02399         "%f s\n", wkf_timer_timenow(timer));
02400     msgInfo << report << sendmsg;
02401 #endif
02402 
02403   } else {
02404 
02405 #ifdef TIMING
02406     wkf_timer_start(timer);
02407 #endif
02408 
02409     // Assuming the grid is aligned with the coordinate axes:
02410     float lenx = float(nsampx*delta);
02411     float leny = float(nsampy*delta);
02412     float lenz = float(nsampz*delta);
02413 
02414     int retval;
02415 
02416     ComputeOccupancyMap om;
02417 
02418     // must be set by caller
02419     om.map = voldata;
02420     om.mx = nsampx;
02421     om.my = nsampy;
02422     om.mz = nsampz;
02423     om.lx = lenx;
02424     om.ly = leny;
02425     om.lz = lenz;
02426     om.x0 = originalign[0];
02427     om.y0 = originalign[1];
02428     om.z0 = originalign[2];
02429     memcpy(om.ax, &axesalign[0], 3*sizeof(float));
02430     memcpy(om.ay, &axesalign[3], 3*sizeof(float));
02431     memcpy(om.az, &axesalign[6], 3*sizeof(float));
02432     memcpy(om.alignmat, alignment.mat, 16*sizeof(float));
02433     om.num_coords = numcoords;
02434     om.coords = coords;
02435     om.vdw_type = vdwtypes;
02436     om.vdw_params = vdwparams;
02437     om.probe_vdw_params = probe_vdw;
02438     om.conformers = conformers;
02439     om.num_probes = num_probe_atoms;
02440     om.num_conformers = num_conformers;
02441     om.cutoff = cutoff;
02442     om.extcutoff = extcutoff;
02443     om.excl_dist = excl_dist; 
02444     om.excl_energy = max_energy;
02445 
02446     // single threaded version calculates one largest slab
02447     om.kstart = 0;
02448     om.kstop = om.mz;
02449 
02450     retval = ComputeOccupancyMap_setup(&om);
02451 
02452 #ifdef TIMING
02453     sprintf(report, "ComputeOccupancyMap_setup()                            "
02454         "%f s\n", wkf_timer_timenow(timer));
02455     msgInfo << report << sendmsg;
02456 #endif
02457 
02458     if (getenv("VMDILSVERBOSE")) { // XXX debugging
02459       atom_bin_stats(&om
02460           /*
02461           gridorigin, lenx, leny, lenz, extcutoff, cutoff,
02462           om.bincnt, om.nbx, om.nby, om.nbz, om.padx,
02463           om.num_bin_offsets, om.num_extras */);
02464     }
02465 
02466     if (retval != 0) {
02467       if (getenv("VMDILSVERBOSE")) { // XXX debugging
02468         int i, j, k;
02469         int total_extra_atoms = 0;
02470         for (k = 0;  k < om.nbz;  k++) {
02471           for (j = 0;  j < om.nby;  j++) {
02472             for (i = 0;  i < om.nbx;  i++) {
02473               int index = (k*om.nby + j)*om.nbx + i;
02474               if (om.bincnt[index] > BIN_DEPTH) {
02475                 printf("*** bin[%d,%d,%d] tried to fill with %d atoms\n",
02476                     i, j, k, om.bincnt[index]);
02477                 total_extra_atoms += om.bincnt[index] - BIN_DEPTH;
02478               }
02479             }
02480           }
02481         }
02482         // XXX should have total_extra_atoms > num_extra_atoms
02483         printf("*** can't handle total of %d extra atoms\n", total_extra_atoms);
02484       }
02485       ComputeOccupancyMap_cleanup(&om);
02486       return -1;
02487     }
02488 
02489 #if defined(VMDCUDA)
02490     if (getenv("VMDILSVERBOSE")) {
02491       printf("*** cpu only = %d\n", om.cpu_only);
02492     }
02493     if (!getenv("VMDNOCUDA") && !(om.cpu_only) &&
02494         (retval=
02495          vmd_cuda_evaluate_occupancy_map(om.mx, om.my, om.mz, om.map,
02496            om.excl_energy, om.cutoff, om.hx, om.hy, om.hz,
02497            om.x0, om.y0, om.z0, om.bx_1, om.by_1, om.bz_1, 
02498            om.nbx, om.nby, om.nbz, (float *) om.bin, (float *) om.bin_zero,
02499            om.num_bin_offsets, om.bin_offsets,
02500            om.num_extras, (float *) om.extra,
02501            num_unique_types, om.vdw_params,
02502            om.num_probes, om.probe_vdw_params,
02503            om.num_conformers, om.conformers)) == 0) {
02504       // successfully ran ILS with CUDA, otherwise fall back on CPU
02505     } else {
02506       if (retval != 0) {
02507         msgInfo << "vmd_cuda_evaluate_occupancy_map() FAILED,"
02508           " using CPU for calculation\n" << sendmsg;
02509       }
02510 #endif /* CUDA... */
02511 
02512 #ifdef TIMING
02513       wkf_timer_start(timer);
02514 #endif
02515 
02516       retval = ComputeOccupancyMap_calculate_slab(&om);
02517 
02518 #ifdef TIMING
02519       sprintf(report, "ComputeOccupancyMap_calculate_slab()                   "
02520               "%f s\n", wkf_timer_timenow(timer));
02521       msgInfo << report << sendmsg;
02522 #endif
02523 
02524       if (retval != 0) {
02525         if (getenv("VMDILSVERBOSE")) { // XXX debugging
02526           printf("*** ComputeOccupancyMap_calculate_slab() failed\n");
02527         }
02528         ComputeOccupancyMap_cleanup(&om);
02529         return -1;
02530       }
02531 
02532 #if defined(VMDCUDA)
02533     } // end else not VMDCUDAILS
02534 #endif
02535 
02536     ComputeOccupancyMap_cleanup(&om);
02537   }
02538 
02539   if (!getenv("VMDILSALIGNMAPS")) {
02540     delete[] coords;
02541     delete[] vdwtypes;
02542   }
02543 
02544 #ifdef TIMING
02545   wkf_timer_destroy(timer);
02546 #endif
02547       
02548   return MEASURE_NOERR; 
02549 }
02550 
02551 
02553 // Here follows the new implementation.                //
02555 
02556 static int fill_atom_bins(ComputeOccupancyMap *p);
02557 static void tighten_bin_neighborhood(ComputeOccupancyMap *p);
02558 static void find_distance_exclusions(ComputeOccupancyMap *p);
02559 static void find_energy_exclusions(ComputeOccupancyMap *p);
02560 static void compute_occupancy_monoatom(ComputeOccupancyMap *p);
02561 static void compute_occupancy_multiatom(ComputeOccupancyMap *p);
02562 
02563 
02564 int ComputeOccupancyMap_setup(ComputeOccupancyMap *p) {
02565 
02566   // initialize pointer fields to zero
02567   p->bin = NULL;
02568   p->bin_zero = NULL;
02569   p->bincnt = NULL;
02570   p->bincnt_zero = NULL;
02571   p->bin_offsets = NULL;
02572   p->extra = NULL;
02573   p->exclusions = NULL;
02574 
02575   // initialize occupancy map, allocate and initialize exclusion map
02576   int mtotal = p->mx * p->my * p->mz;
02577   memset(p->map, 0, mtotal * sizeof(float));  // zero occupancy by default
02578   p->exclusions = new char[mtotal];
02579   memset(p->exclusions, 0, mtotal * sizeof(char));  // no exclusions yet
02580 
02581   // derive map spacing based on length and number of points
02582   p->hx = p->lx / p->mx;
02583   p->hy = p->ly / p->my;
02584   p->hz = p->lz / p->mz;
02585 
02586   p->cpu_only = 0;  // attempt to use CUDA
02587 
02588   // set expected map points per bin length
02589   // note: we want CUDA thread blocks to calculate 4^3 map points
02590   //       we want each thread block contained inside a single bin
02591   //       we want bin volume to be no more than MAX_BIN_VOLUME (27 A^3)
02592   //       and no smaller than MIN_BIN_VOLUME (8 A^3) due to fixed bin depth
02593   //       we expect map spacing to be about 0.25 A but depends on caller
02594 
02595   // start with trying to pack 3^3 thread blocks per atom bin
02596   p->mpblx = 12;
02597   p->mpbly = 12;
02598   p->mpblz = 12;
02599 
02600   // starting bin lengths
02601   p->bx = p->mpblx * p->hx;
02602   p->by = p->mpbly * p->hy;
02603   p->bz = p->mpblz * p->hz;
02604 
02605   // refine bin side lengths if volume of bin is too large
02606   while (p->bx * p->by * p->bz > MAX_BIN_VOLUME) {
02607 
02608     // find longest bin side and reduce its length
02609     if (p->bx > p->by && p->bx > p->bz) {
02610       p->mpblx -= 4;
02611       p->bx = p->mpblx * p->hx;
02612     }
02613     else if (p->by >= p->bx && p->by > p->bz) {
02614       p->mpbly -= 4;
02615       p->by = p->mpbly * p->hy;
02616     }
02617     else {
02618       p->mpblz -= 4;
02619       p->bz = p->mpblz * p->hz;
02620     }
02621 
02622   } // end refinement of bins
02623 
02624   if (p->bx * p->by * p->bz < MIN_BIN_VOLUME) {
02625     // refinement failed due to some hx, hy, hz being too large
02626     // now there is no known correspondence between map points and bins
02627     p->bx = p->by = p->bz = DEFAULT_BIN_LENGTH;
02628     p->mpblx = p->mpbly = p->mpblz = 0;
02629     p->cpu_only = 1;  // CUDA can't be used, map too coarse
02630   }
02631 
02632   p->bx_1 = 1.f / p->bx;
02633   p->by_1 = 1.f / p->by;
02634   p->bz_1 = 1.f / p->bz;
02635 
02636   if (fill_atom_bins(p)) {
02637     return -1;  // failed due to too many extra atoms for bin size
02638   }
02639 
02640   tighten_bin_neighborhood(p);
02641 
02642   return 0;
02643 } // ComputeOccupancyMap_setup()
02644 
02645 
02646 int ComputeOccupancyMap_calculate_slab(ComputeOccupancyMap *p) {
02647 #ifdef TIMING
02648   char report[128] = { 0 };
02649   wkf_timerhandle timer = wkf_timer_create();
02650 #endif
02651 
02652   // each of these routines operates on the slab
02653   // designated by kstart through kstop (z-axis indices)
02654   //
02655   // XXX we are planning CUDA kernels for each of the following routines
02656 
02657 #if 1
02658 #ifdef TIMING
02659   wkf_timer_start(timer);
02660 #endif
02661   find_distance_exclusions(p);
02662   int i, numexcl=0;
02663   for (i=0; i<p->mx * p->my * p->mz; i++) {
02664     if (p->exclusions[i]) numexcl++;
02665   }
02666 #ifdef TIMING
02667   sprintf(report, "ComputeOccupancyMap: find_distance_exclusions()        "
02668       "%f s\n", wkf_timer_timenow(timer));
02669   msgInfo << report << sendmsg;
02670 #endif
02671 #endif
02672 
02673   if (1 == p->num_probes) {
02674 #ifdef TIMING
02675     wkf_timer_start(timer);
02676 #endif
02677     compute_occupancy_monoatom(p);
02678 #ifdef TIMING
02679     sprintf(report, "ComputeOccupancyMap: compute_occupancy_monoatom()      "
02680         "%f s\n", wkf_timer_timenow(timer));
02681     msgInfo << report << sendmsg;
02682 #endif
02683 
02684   }
02685   else {
02686 
02687 #if 1
02688 #ifdef TIMING
02689     wkf_timer_start(timer);
02690 #endif
02691     find_energy_exclusions(p);
02692     int i, numexcl=0;
02693     for (i=0; i<p->mx * p->my * p->mz; i++) {
02694       if (p->exclusions[i]) numexcl++;
02695     }
02696 #ifdef TIMING
02697     sprintf(report, "ComputeOccupancyMap: find_energy_exclusions()          "
02698         "%f s  -> %d exclusions\n", wkf_timer_timenow(timer), numexcl);
02699     msgInfo << report << sendmsg;
02700 #endif
02701 #endif
02702 
02703 #ifdef TIMING
02704     wkf_timer_start(timer);
02705 #endif
02706     compute_occupancy_multiatom(p);
02707 #ifdef TIMING
02708     sprintf(report, "ComputeOccupancyMap: compute_occupancy_multiatom()     "
02709         "%f s\n", wkf_timer_timenow(timer));
02710     msgInfo << report << sendmsg;
02711 #endif
02712 
02713   }
02714 
02715 #ifdef TIMING
02716   wkf_timer_destroy(timer);
02717 #endif
02718 
02719   return 0;
02720 } // ComputeOccupancyMap_calculate_slab()
02721 
02722 
02723 void ComputeOccupancyMap_cleanup(ComputeOccupancyMap *p) {
02724   delete[] p->bin_offsets;
02725   delete[] p->extra;
02726   delete[] p->bincnt;
02727   delete[] p->bin;
02728   delete[] p->exclusions;
02729 } // ComputeOccupancyMap_cleanup()
02730 
02731 
02732 // XXX the CUDA kernels can handle up to 50 extra atoms, set this as
02733 // the bound on "max_extra_atoms" rather than the heuristic below
02734 #define MAX_EXTRA_ATOMS  50
02735 
02736 int fill_atom_bins(ComputeOccupancyMap *p) {
02737   int too_many_extra_atoms = 0;  // be optimistic
02738   int max_extra_atoms = MAX_EXTRA_ATOMS;
02739   //int max_extra_atoms = (int) ceilf(p->num_coords / 10000.f);
02740       // assume no more than 1 over full bin per 10000 atoms
02741   int count_extras = 0;
02742   int n, i, j, k;
02743 
02744   const int *vdw_type = p->vdw_type;
02745   const float *coords = p->coords;
02746   const int num_coords = p->num_coords;
02747   const float lx = p->lx;
02748   const float ly = p->ly;
02749   const float lz = p->lz;
02750   const float bx_1 = p->bx_1;
02751   const float by_1 = p->by_1;
02752   const float bz_1 = p->bz_1;
02753   const float x0 = p->x0;
02754   const float y0 = p->y0;
02755   const float z0 = p->z0;
02756   const float extcutoff = p->extcutoff;
02757 
02758   // padding is based on extended cutoff distance
02759   p->padx = (int) ceilf(extcutoff * bx_1);
02760   p->pady = (int) ceilf(extcutoff * by_1);
02761   p->padz = (int) ceilf(extcutoff * bz_1);
02762 
02763   const int nbx = p->nbx = (int) ceilf(lx * bx_1) + 2*p->padx;
02764   const int nby = p->nby = (int) ceilf(ly * by_1) + 2*p->pady;
02765   p->nbz = (int) ceilf(lz * bz_1) + 2*p->padz;
02766 
02767   int nbins = nbx * nby * p->nbz;
02768 
02769   BinOfAtoms *bin = p->bin = new BinOfAtoms[nbins];
02770   char *bincnt = p->bincnt = new char[nbins];
02771   AtomPosType *extra = p->extra = new AtomPosType[max_extra_atoms];
02772 
02773   memset(bin, 0, nbins * sizeof(BinOfAtoms));
02774   memset(bincnt, 0, nbins * sizeof(char));
02775 
02776   // shift array pointer to the (0,0,0)-bin, which will correspond to
02777   // the map origin
02778   BinOfAtoms *bin_zero
02779     = p->bin_zero = bin + ((p->padz*nby + p->pady)*nbx + p->padx);
02780   char *bincnt_zero
02781     = p->bincnt_zero = bincnt + ((p->padz*nby + p->pady)*nbx + p->padx);
02782 
02783   for (n = 0;  n < num_coords;  n++) {
02784     float x = coords[3*n    ];  // atom coordinates
02785     float y = coords[3*n + 1];
02786     float z = coords[3*n + 2];
02787 
02788     float sx = x - x0;  // translate relative to map origin
02789     float sy = y - y0;
02790     float sz = z - z0;
02791 
02792     if (sx < -extcutoff || sx > lx + extcutoff ||
02793         sy < -extcutoff || sy > ly + extcutoff ||
02794         sz < -extcutoff || sz > lz + extcutoff) {
02795       continue;  // atom is beyond influence of lattice
02796     }
02797 
02798     i = (int) floorf(sx * bx_1);  // bin number
02799     j = (int) floorf(sy * by_1);
02800     k = (int) floorf(sz * bz_1);
02801 
02802     /*
02803     // XXX this test should never be true after passing previous test
02804     if (i < -p->padx || i >= p->nbx + p->padx ||
02805         j < -p->pady || j >= p->nby + p->pady ||
02806         k < -p->padz || k >= p->nbz + p->padz) {
02807       continue;  // atom is outside bin array
02808     }
02809     */
02810 
02811     int index = (k*nby + j)*nbx + i;  // flat index into bin array
02812     int slot = bincnt_zero[index];  // slot within bin to place atom
02813 
02814     if (slot < BIN_DEPTH) {
02815       AtomPosType *atom = bin_zero[index].atom;  // place atom in next slot
02816       atom[slot].x = x;
02817       atom[slot].y = y;
02818       atom[slot].z = z;
02819       atom[slot].vdwtype = vdw_type[n];
02820     }
02821     else if (count_extras < max_extra_atoms) {
02822       extra[count_extras].x = x;
02823       extra[count_extras].y = y;
02824       extra[count_extras].z = z;
02825       extra[count_extras].vdwtype = vdw_type[n];
02826       count_extras++;
02827     }
02828     else {
02829       // XXX debugging
02830       printf("*** too many extras, atom index %d\n", n);
02831       too_many_extra_atoms = 1;
02832     }
02833 
02834     bincnt_zero[index]++;  // increase count of atoms in bin
02835   }
02836   p->num_extras = count_extras;
02837 
02838   // mark unused atom slots
02839   // XXX set vdwtype to -1
02840   for (n = 0;  n < nbins;  n++) {
02841     for (k = bincnt[n];  k < BIN_DEPTH;  k++) {
02842       bin[n].atom[k].vdwtype = -1;
02843     }
02844   }
02845 
02846   return (too_many_extra_atoms ? -1 : 0);
02847 } // fill_atom_bins()
02848 
02849 
02850 // setup tightened bin index offset array of 3-tuples
02851 void tighten_bin_neighborhood(ComputeOccupancyMap *p) {
02852   const int padx = p->padx;
02853   const int pady = p->pady;
02854   const int padz = p->padz;
02855   const float bx2 = p->bx * p->bx;
02856   const float by2 = p->by * p->by;
02857   const float bz2 = p->bz * p->bz;
02858   const float r = p->extcutoff + sqrtf(bx2 + by2 + bz2);  // add bin diagonal
02859   const float r2 = r*r;
02860   int n = 0, i, j, k;
02861   char *bin_offsets
02862     = p->bin_offsets = new char[3 * (2*padx+1)*(2*pady+1)*(2*padz+1)];
02863   for (k = -padz;  k <= padz;  k++) {
02864     for (j = -pady;  j <= pady;  j++) {
02865       for (i = -padx;  i <= padx;  i++) {
02866         if (i*i*bx2 + j*j*by2 + k*k*bz2 >= r2) continue;
02867         bin_offsets[3*n    ] = (char) i;
02868         bin_offsets[3*n + 1] = (char) j;
02869         bin_offsets[3*n + 2] = (char) k;
02870         n++;
02871       }
02872     }
02873   }
02874   p->num_bin_offsets = n;
02875 } // tighten_bin_neighborhood()
02876 
02877 
02878 // For each grid point loop over the close atoms and 
02879 // determine if one of them is closer than excl_dist
02880 // away. If so we assume the clash with the probe will
02881 // result in a very high interaction energy and we can 
02882 // exclude this point from calcultation.
02883 void find_distance_exclusions(ComputeOccupancyMap *p) {
02884   const AtomPosType *extra = p->extra;
02885   const BinOfAtoms *bin_zero = p->bin_zero;
02886   char *excl = p->exclusions;
02887 
02888   int i, j, k, n, index;
02889   int ic, jc, kc;
02890 
02891   const int mx = p->mx;
02892   const int my = p->my;
02893   const int kstart = p->kstart;
02894   const int kstop = p->kstop;
02895   const int nbx = p->nbx;
02896   const int nby = p->nby;
02897   const float excl_dist = p->excl_dist;
02898   const float bx_1 = p->bx_1;
02899   const float by_1 = p->by_1;
02900   const float bz_1 = p->bz_1;
02901   const float hx = p->hx;
02902   const float hy = p->hy;
02903   const float hz = p->hz;
02904   const float x0 = p->x0;
02905   const float y0 = p->y0;
02906   const float z0 = p->z0;
02907   const int num_extras = p->num_extras;
02908   const int bdx = (int) ceilf(excl_dist * bx_1);  // width of nearby bins
02909   const int bdy = (int) ceilf(excl_dist * by_1);  // width of nearby bins
02910   const int bdz = (int) ceilf(excl_dist * bz_1);  // width of nearby bins
02911   const float excldist2 = excl_dist * excl_dist;
02912 
02913   for (k = kstart;  k < kstop;  k++) {  // k index loops over slab
02914     for (j = 0;  j < my;  j++) {
02915       for (i = 0;  i < mx;  i++) {  // loop over map points
02916 
02917         float px = i*hx;
02918         float py = j*hy;
02919         float pz = k*hz;  // translated coordinates of map point
02920 
02921         int ib = (int) floorf(px * bx_1);
02922         int jb = (int) floorf(py * by_1);
02923         int kb = (int) floorf(pz * bz_1);  // zero-based bin index
02924 
02925         px += x0;
02926         py += y0;
02927         pz += z0;  // absolute position
02928 
02929         for (kc = kb - bdz;  kc <= kb + bdz;  kc++) {
02930           for (jc = jb - bdy;  jc <= jb + bdy;  jc++) {
02931             for (ic = ib - bdx;  ic <= ib + bdx;  ic++) {
02932 
02933               const AtomPosType *atom
02934                 = bin_zero[(kc*nby + jc)*nbx + ic].atom;
02935 
02936               for (n = 0;  n < BIN_DEPTH;  n++) {  // atoms in bin
02937                 if (-1 == atom[n].vdwtype) break;  // finished atoms in bin
02938                 float dx = px - atom[n].x;
02939                 float dy = py - atom[n].y;
02940                 float dz = pz - atom[n].z;
02941                 float r2 = dx*dx + dy*dy + dz*dz;
02942                 if (r2 <= excldist2) {
02943                   index = (k*my + j)*mx + i;
02944                   excl[index] = 1;
02945                   goto NEXT_MAP_POINT;  // don't have to look at more atoms
02946                 }
02947               } // end loop over atoms in bin
02948 
02949             }
02950           }
02951         } // end loop over nearby bins
02952 
02953         for (n = 0;  n < num_extras;  n++) {  // extra atoms
02954           float dx = px - extra[n].x;
02955           float dy = py - extra[n].y;
02956           float dz = pz - extra[n].z;
02957           float r2 = dx*dx + dy*dy + dz*dz;
02958           if (r2 <= excldist2) {
02959             index = (k*my + j)*mx + i;
02960             excl[index] = 1;
02961             goto NEXT_MAP_POINT;  // don't have to look at more atoms
02962           }
02963         } // end loop over extra atoms
02964 
02965 NEXT_MAP_POINT:
02966         ; // continue loop over lattice points
02967 
02968       }
02969     }
02970   } // end loop over lattice points
02971 
02972 } // find_distance_exclusions()
02973 
02974 
02975 
02976 // For each grid point sum up the energetic contribution
02977 // of all close atoms. If that interaction energy is above
02978 // the excl_energy cutoff value we don't have to consider
02979 // this grid point in the subsequent calculation.
02980 // Hence, we save the computation of the different probe
02981 // orientation for multiatom probes.
02982 void find_energy_exclusions(ComputeOccupancyMap *p) {
02983   const char *bin_offsets = p->bin_offsets;
02984   const float *vdw_params = p->vdw_params;
02985   const AtomPosType *extra = p->extra;
02986   const BinOfAtoms *bin_zero = p->bin_zero;
02987   char *excl = p->exclusions;
02988 
02989   const float probe_vdweps = p->probe_vdw_params[0];   // use first probe param
02990   const float probe_vdwrmin = p->probe_vdw_params[1];  // for epsilon and rmin
02991 
02992   const int mx = p->mx;
02993   const int my = p->my;
02994   const int kstart = p->kstart;
02995   const int kstop = p->kstop;
02996   const int nbx = p->nbx;
02997   const int nby = p->nby;
02998   const float hx = p->hx;
02999   const float hy = p->hy;
03000   const float hz = p->hz;
03001   const float bx_1 = p->bx_1;
03002   const float by_1 = p->by_1;
03003   const float bz_1 = p->bz_1;
03004   const float x0 = p->x0;
03005   const float y0 = p->y0;
03006   const float z0 = p->z0;
03007   const int num_bin_offsets = p->num_bin_offsets;
03008   const int num_extras = p->num_extras;
03009   const float excl_energy = p->excl_energy;
03010 
03011   int i, j, k, n, index;
03012   const float cutoff2 = p->cutoff * p->cutoff;
03013 
03014   for (k = kstart;  k < kstop;  k++) {  // k index loops over slab
03015     for (j = 0;  j < my;  j++) {
03016       for (i = 0;  i < mx;  i++) {  // loop over map points
03017 
03018         int lindex = (k*my + j)*mx + i;  // map index
03019         if (excl[lindex]) continue;  // already excluded based on distance
03020 
03021         float px = i*hx;
03022         float py = j*hy;
03023         float pz = k*hz;  // translated coordinates of map point
03024 
03025         int ib = (int) floorf(px * bx_1);
03026         int jb = (int) floorf(py * by_1);
03027         int kb = (int) floorf(pz * bz_1);  // zero-based bin index
03028 
03029         px += x0;
03030         py += y0;
03031         pz += z0;  // absolute position
03032 
03033         float u = 0.f;
03034 
03035         for (index = 0;  index < num_bin_offsets;  index++) { // neighborhood
03036           int ic = ib + (int) bin_offsets[3*index    ];
03037           int jc = jb + (int) bin_offsets[3*index + 1];
03038           int kc = kb + (int) bin_offsets[3*index + 2];
03039 
03040           const AtomPosType *atom
03041             = bin_zero[(kc*nby + jc)*nbx + ic].atom;
03042 
03043           for (n = 0;  n < BIN_DEPTH;  n++) {  // atoms in bin
03044             if (-1 == atom[n].vdwtype) break;  // finished atoms in bin
03045             float dx = px - atom[n].x;
03046             float dy = py - atom[n].y;
03047             float dz = pz - atom[n].z;
03048             float r2 = dx*dx + dy*dy + dz*dz;
03049             if (r2 >= cutoff2) continue;
03050             int pindex = 2 * atom[n].vdwtype;
03051             float epsilon = vdw_params[pindex] * probe_vdweps;
03052             float rmin = vdw_params[pindex + 1] + probe_vdwrmin;
03053             float rm6 = rmin*rmin / r2;
03054             rm6 = rm6 * rm6 * rm6;
03055             u += epsilon * rm6 * (rm6 - 2.f);  // sum vdw contribution
03056           } // end loop atoms in bin
03057 
03058         } // end loop bin neighborhood
03059 
03060         for (n = 0;  n < num_extras;  n++) {  // extra atoms
03061           float dx = px - extra[n].x;
03062           float dy = py - extra[n].y;
03063           float dz = pz - extra[n].z;
03064           float r2 = dx*dx + dy*dy + dz*dz;
03065           if (r2 >= cutoff2) continue;
03066           int pindex = 2 * extra[n].vdwtype;
03067           float epsilon = vdw_params[pindex] * probe_vdweps;
03068           float rmin = vdw_params[pindex + 1] + probe_vdwrmin;
03069           float rm6 = rmin*rmin / r2;
03070           rm6 = rm6 * rm6 * rm6;
03071           u += epsilon * rm6 * (rm6 - 2.f);  // sum vdw contribution
03072         } // end loop over extra atoms
03073 
03074         if (u >= excl_energy) excl[lindex] = 1;
03075 
03076       }
03077     }
03078   } // end loop over lattice
03079 
03080 } // find_energy_exclusions()
03081 
03082 
03083 // For a monoatomic probe compute the occupancy rho
03084 // (probability of finding the probe)
03085 //
03086 // For each map point the occupancy is computed as
03087 //
03088 //   rho = exp(-U)
03089 //
03090 // where U is the interaction energy of the probe with the system
03091 // due to the VDW force field.
03092 //
03093 void compute_occupancy_monoatom(ComputeOccupancyMap *p) {
03094   const char *bin_offsets = p->bin_offsets;
03095   const float *vdw_params = p->vdw_params;
03096   const AtomPosType *extra = p->extra;
03097   const BinOfAtoms *bin_zero = p->bin_zero;
03098   const char *excl = p->exclusions;
03099   float *map = p->map;
03100 
03101   const int mx = p->mx;
03102   const int my = p->my;
03103   const int kstart = p->kstart;
03104   const int kstop = p->kstop;
03105   const int nbx = p->nbx;
03106   const int nby = p->nby;
03107   const float hx = p->hx;
03108   const float hy = p->hy;
03109   const float hz = p->hz;
03110   const float bx_1 = p->bx_1;
03111   const float by_1 = p->by_1;
03112   const float bz_1 = p->bz_1;
03113   const float x0 = p->x0;
03114   const float y0 = p->y0;
03115   const float z0 = p->z0;
03116   const int num_bin_offsets = p->num_bin_offsets;
03117   const int num_extras = p->num_extras;
03118 
03119   float probe_vdweps = p->probe_vdw_params[0];   // use first probe param
03120   float probe_vdwrmin = p->probe_vdw_params[1];  //   for epsilon and rmin
03121 
03122   int i, j, k, n, index;
03123   float max_energy = p->excl_energy;
03124   float cutoff2 = p->cutoff * p->cutoff;
03125 
03126   for (k = kstart;  k < kstop;  k++) {  // k index loops over slab
03127     for (j = 0;  j < my;  j++) {
03128       for (i = 0;  i < mx;  i++) {  // loop over lattice points
03129 
03130         int lindex = (k*my + j)*mx + i;  // map index
03131 #if 1
03132         if (excl[lindex]) {  // is map point excluded?
03133           map[lindex] = 0.f;  // clamp occupancy to zero
03134           continue;
03135         }
03136 #endif
03137 
03138         float px = i*hx;
03139         float py = j*hy;
03140         float pz = k*hz;  // translated coordinates of lattice point
03141 
03142         int ib = (int) floorf(px * bx_1);
03143         int jb = (int) floorf(py * by_1);
03144         int kb = (int) floorf(pz * bz_1);  // zero-based bin index
03145 
03146         px += x0;
03147         py += y0;
03148         pz += z0;  // absolute position
03149 
03150         float u = 0.f;
03151 
03152         for (index = 0;  index < num_bin_offsets;  index++) { // neighborhood
03153           int ic = ib + (int) bin_offsets[3*index    ];
03154           int jc = jb + (int) bin_offsets[3*index + 1];
03155           int kc = kb + (int) bin_offsets[3*index + 2];
03156 
03157           const AtomPosType *atom
03158             = bin_zero[(kc*nby + jc)*nbx + ic].atom;
03159 
03160           for (n = 0;  n < BIN_DEPTH;  n++) {  // atoms in bin
03161             if (-1 == atom[n].vdwtype) break;  // finished atoms in bin
03162             float dx = px - atom[n].x;
03163             float dy = py - atom[n].y;
03164             float dz = pz - atom[n].z;
03165             float r2 = dx*dx + dy*dy + dz*dz;
03166             if (r2 >= cutoff2) continue;
03167             int pindex = 2 * atom[n].vdwtype;
03168             float epsilon = vdw_params[pindex] * probe_vdweps;
03169             float rmin = vdw_params[pindex + 1] + probe_vdwrmin;
03170             float rm6 = rmin*rmin / r2;
03171             rm6 = rm6 * rm6 * rm6;
03172             u += epsilon * rm6 * (rm6 - 2.f);  // sum vdw contribution
03173           } // end loop atoms in bin
03174 
03175         } // end loop bin neighborhood
03176 
03177         for (n = 0;  n < num_extras;  n++) {  // extra atoms
03178           float dx = px - extra[n].x;
03179           float dy = py - extra[n].y;
03180           float dz = pz - extra[n].z;
03181           float r2 = dx*dx + dy*dy + dz*dz;
03182           if (r2 >= cutoff2) continue;
03183           int pindex = 2 * extra[n].vdwtype;
03184           float epsilon = vdw_params[pindex] * probe_vdweps;
03185           float rmin = vdw_params[pindex + 1] + probe_vdwrmin;
03186           float rm6 = rmin*rmin / r2;
03187           rm6 = rm6 * rm6 * rm6;
03188           u += epsilon * rm6 * (rm6 - 2.f);  // sum vdw contribution
03189         } // end loop over extra atoms
03190 
03191         float occ = 0.f;
03192         if (u < max_energy) {
03193           occ = expf(-u);
03194         }
03195         map[lindex] = occ;  // the occupancy
03196 
03197       }
03198     }
03199   } // end loop over lattice
03200 
03201 } // compute_occupancy_monoatom()
03202 
03203 
03204 // For a multiatom probe compute the occupancy rho
03205 // (probability of finding the probe)
03206 //
03207 // Calculate occupancy rho at each map point,
03208 // where rho = (1/m) sum_m ( -exp(u[i]) ) over m conformers,
03209 // u[i] is the potential energy of the i-th conformer.
03210 //
03211 void compute_occupancy_multiatom(ComputeOccupancyMap *p) {
03212   const char *bin_offsets = p->bin_offsets;
03213   const float *vdw_params = p->vdw_params;
03214   const float *probe_vdw_params = p->probe_vdw_params;
03215   const float *conformers = p->conformers;
03216   const AtomPosType *extra = p->extra;
03217   const float hx = p->hx;
03218   const float hy = p->hy;
03219   const float hz = p->hz;
03220   const float x0 = p->x0;
03221   const float y0 = p->y0;
03222   const float z0 = p->z0;
03223   const float bx_1 = p->bx_1;
03224   const float by_1 = p->by_1;
03225   const float bz_1 = p->bz_1;
03226   const float inv_num_conformers = 1.f / (float) p->num_conformers;
03227   const int num_bin_offsets = p->num_bin_offsets;
03228   const int num_extras = p->num_extras;
03229   const int num_probes = p->num_probes;
03230   const int num_conformers = p->num_conformers;
03231   const int nbx = p->nbx;
03232   const int nby = p->nby;
03233   const int mx = p->mx;
03234   const int my = p->my;
03235   const int kstart = p->kstart;
03236   const int kstop = p->kstop;
03237 
03238   const BinOfAtoms *bin_zero = p->bin_zero;
03239   const char *excl = p->exclusions;
03240   float *map = p->map;
03241 
03242   int i, j, k, n, nb;
03243 
03244   const float minocc = expf(-p->excl_energy);
03245   const float cutoff2 = p->cutoff * p->cutoff;
03246 
03247   float *u = new float[p->num_conformers];  // cal potential for each conformer
03248 
03249   for (k = kstart;  k < kstop;  k++) {  // k index loops over slab
03250     for (j = 0;  j < my;  j++) {
03251       for (i = 0;  i < mx;  i++) {  // loop over lattice points
03252 
03253         int lindex = (k*my + j)*mx + i;  // map index
03254         if (excl[lindex]) {  // is map point excluded?
03255           map[lindex] = 0.f;  // clamp occupancy to zero
03256           continue;
03257         }
03258 
03259         float px = i*hx;
03260         float py = j*hy;
03261         float pz = k*hz;  // translated coordinates of lattice point
03262 
03263         int ib = (int) floorf(px * bx_1);
03264         int jb = (int) floorf(py * by_1);
03265         int kb = (int) floorf(pz * bz_1);  // zero-based bin index
03266 
03267         int m, ma;
03268 
03269         px += x0;
03270         py += y0;
03271         pz += z0;  // absolute position
03272 
03273         memset(u, 0, num_conformers * sizeof(float));
03274 
03275         for (nb = 0;  nb < num_bin_offsets;  nb++) { // bin neighborhood
03276           int ic = ib + (int) bin_offsets[3*nb    ];
03277           int jc = jb + (int) bin_offsets[3*nb + 1];
03278           int kc = kb + (int) bin_offsets[3*nb + 2];
03279 
03280           const AtomPosType *atom = bin_zero[(kc*nby + jc)*nbx + ic].atom;
03281 
03282           for (n = 0;  n < BIN_DEPTH;  n++) {  // atoms in bin
03283             if (-1 == atom[n].vdwtype) break;  // finished atoms in bin
03284 
03285             for (m = 0;  m < num_conformers;  m++) {  // conformers
03286               float v = 0.f;
03287               for (ma = 0;  ma < num_probes;  ma++) {  // probe
03288                 int index = m*num_probes + ma;
03289                 float dx = conformers[3*index    ] + px - atom[n].x;
03290                 float dy = conformers[3*index + 1] + py - atom[n].y;
03291                 float dz = conformers[3*index + 2] + pz - atom[n].z;
03292                 float r2 = dx*dx + dy*dy + dz*dz;
03293                 if (r2 >= cutoff2) continue;
03294                 int pindex = 2 * atom[n].vdwtype;
03295                 float epsilon = vdw_params[pindex] * probe_vdw_params[2*ma];
03296                 float rmin = vdw_params[pindex + 1] + probe_vdw_params[2*ma + 1];
03297                 float rm6 = rmin*rmin / r2;
03298                 rm6 = rm6 * rm6 * rm6;
03299                 v += epsilon * rm6 * (rm6 - 2.f);  // sum vdw contribution
03300               } // end loop probe
03301 
03302               u[m] += v;  // contribution of one system atom to conformer
03303 
03304             } // end loop conformers
03305 
03306           } // end loop atoms in bin
03307 
03308         } // end loop bin neighborhood
03309 
03310         for (n = 0;  n < num_extras;  n++) {  // extra atoms
03311           for (m = 0;  m < num_conformers;  m++) {  // conformers
03312             float v = 0.f;
03313             for (ma = 0;  ma < num_probes;  ma++) {  // probe
03314               int index = m*num_probes + ma;
03315               float dx = conformers[3*index    ] + px - extra[n].x;
03316               float dy = conformers[3*index + 1] + py - extra[n].y;
03317               float dz = conformers[3*index + 2] + pz - extra[n].z;
03318               float r2 = dx*dx + dy*dy + dz*dz;
03319               if (r2 >= cutoff2) continue;
03320               int pindex = 2 *extra[n].vdwtype;
03321               float epsilon = vdw_params[pindex] * probe_vdw_params[2*ma];
03322               float rmin = vdw_params[pindex + 1] + probe_vdw_params[2*ma + 1];
03323               float rm6 = rmin*rmin / r2;
03324               rm6 = rm6 * rm6 * rm6;
03325               v += epsilon * rm6 * (rm6 - 2.f);  // sum vdw contribution
03326             } // end loop probe
03327 
03328             u[m] += v;  // contribution of one system atom to conformer
03329 
03330           } // end loop conformers
03331         } // end loop over extra atoms
03332 
03333         // now we have energies of all conformers u[i], i=0..m-1
03334 
03335         float z = 0.f;
03336         for (m = 0;  m < num_conformers;  m++) {  // average over conformers
03337           z += expf(-u[m]);
03338         }
03339 
03340         float occ = z * inv_num_conformers;  // the occupency
03341         map[lindex] = (occ > minocc ? occ : 0.f);
03342       }
03343     }
03344   } // end loop over lattice
03345 
03346   delete[] u;  // free extra memory
03347 } // compute_occupancy_multiatom()
03348 
03349 
03350 // Write bin histogram into a dx map
03351 static void write_bin_histogram_map(
03352     const ComputeOccupancyMap *p,
03353     const char *filename
03354     ) {
03355   float xaxis[3], yaxis[3], zaxis[3];
03356   memset(xaxis, 0, 3*sizeof(float));
03357   memset(yaxis, 0, 3*sizeof(float));
03358   memset(zaxis, 0, 3*sizeof(float));
03359   xaxis[0] = p->nbx * p->bx;
03360   yaxis[1] = p->nby * p->by;
03361   zaxis[2] = p->nbz * p->bz;
03362 
03363   int gridsize = p->nbx * p->nby * p->nbz;
03364   float *data = new float[gridsize];
03365 
03366   int i;
03367   for (i=0; i<gridsize; i++) {
03368     data[i] = (float) p->bincnt[i];
03369   }
03370 
03371   float ori[3];
03372   ori[0] = p->x0 - p->padx * p->bx;
03373   ori[1] = p->y0 - p->pady * p->by;
03374   ori[2] = p->z0 - p->padz * p->bz;
03375  
03376   VolumetricData *volhist;
03377   volhist = new VolumetricData("atom binning histogram",
03378                                ori, xaxis, yaxis, zaxis,
03379                                p->nbx, p->nby, p->nbz, data);
03380 
03381   // Call the file writer in the VolMapCreate.C:
03382   volmap_write_dx_file(volhist, filename);
03383 
03384   delete volhist;  // XXX does data get deleted as part of volhist?
03385 }
03386 
03387 
03388 // XXX print out histogram of atom bins
03389 void atom_bin_stats(const ComputeOccupancyMap *p) {
03390   int histogram[10] = { 0 };
03391   int i, j, k;
03392 
03393   printf("*** origin = %g %g %g\n", p->x0, p->y0, p->z0);
03394   printf("*** lenx = %g  leny = %g  lenz = %g\n", p->lx, p->ly, p->lz);
03395   printf("*** bin lengths = %g %g %g\n", p->bx, p->by, p->bz);
03396   printf("*** inverse bin lengths = %g %g %g\n", p->bx_1, p->by_1, p->bz_1);
03397   printf("*** bin array:  %d X %d X %d\n", p->nbx, p->nby, p->nbz);
03398   printf("*** bin padding:  %d %d %d\n", p->padx, p->pady, p->padz);
03399   printf("*** cutoff = %g\n", p->cutoff);
03400   printf("*** extcutoff = %g\n", p->extcutoff);
03401   printf("*** size of tight neighborhood = %d\n", p->num_bin_offsets);
03402 
03403   for (k = 0;  k < p->nbz;  k++) {
03404     for (j = 0;  j < p->nby;  j++) {
03405       for (i = 0;  i < p->nbx;  i++) {
03406         int index = (k*p->nby + j)*p->nbx + i;
03407         int count = p->bincnt[index];
03408         histogram[(count <= 9 ? count : 9)]++;
03409       }
03410     }
03411   }
03412 
03413   printf("*** histogram of bin fill:\n");
03414   for (i = 0;  i < (int) (sizeof(histogram) / sizeof(int));  i++) {
03415     if (i < 9) {
03416       printf("***     atom count %d     number of bins %d\n",
03417           i, histogram[i]);
03418     }
03419     else {
03420       printf("***     atom count > 8   number of bins %d\n", histogram[i]);
03421     }
03422   }
03423   printf("*** number of extra atoms = %d\n", p->num_extras);
03424 
03425   write_bin_histogram_map(p, "binning_histogram.dx");
03426 }
03427 
03428 
03429 // XXX slow quadratic complexity algorithm for checking correctness
03430 //     for every map point, for every probe atom in all conformers,
03431 //     iterate over all atoms
03432 void compute_allatoms(
03433     float *map,                    // return calculated occupancy map
03434     int mx, int my, int mz,        // dimensions of map
03435     float lx, float ly, float lz,  // lengths of map
03436     const float origin[3],         // origin of map
03437     const float axes[9],           // basis vectors of aligned map
03438     const float alignmat[16],      // 4x4 alignment matrix
03439     int num_coords,                // number of atoms
03440     const float *coords,           // atom coordinates, length 3*num_coords
03441     const int *vdw_type,           // vdw type numbers, length num_coords
03442     const float *vdw_params,       // scaled vdw parameters for atoms
03443     float cutoff,                  // cutoff distance
03444     const float *probe_vdw_params, // scaled vdw parameters for probe atoms
03445     int num_probe_atoms,           // number of atoms in probe
03446     int num_conformers,            // number of conformers
03447     const float *conformers,       // length 3*num_probe_atoms*num_conformers
03448     float excl_energy              // exclusion energy threshold
03449     ) {
03450   const float theTrivialConformer[] = { 0.f, 0.f, 0.f };
03451 
03452   if (0 == num_conformers) {  // fix to handle trivial case consistently
03453     num_conformers = 1;
03454     conformers = theTrivialConformer;
03455   }
03456 
03457   float hx = lx / mx;  // lattice spacing
03458   float hy = ly / my;
03459   float hz = lz / mz;
03460 
03461   int i, j, k, n, m, ma;
03462 
03463   float cutoff2 = cutoff * cutoff;
03464   float minocc = expf(-excl_energy);  // minimum nonzero occupancy
03465 
03466   float *u = new float[num_conformers];  // calc potential for each conformer
03467 
03468 #if 1
03469   printf("*** All atoms calculation (quadratic complexity)\n");
03470   printf("***   number of map points = %d\n", mx*my*mz);
03471   printf("***   number of atoms = %d\n", num_coords);
03472   printf("***   number of conformers = %d\n", num_conformers);
03473   printf("***   number of probe atoms = %d\n", num_probe_atoms);
03474 #endif
03475 
03476   for (k = 0;  k < mz;  k++) {
03477     for (j = 0;  j < my;  j++) {
03478       for (i = 0;  i < mx;  i++) {  // loop over lattice points
03479 
03480         int mindex = (k*my + j)*mx + i;  // map flat index
03481 
03482         float px = i*hx + origin[0];
03483         float py = j*hy + origin[1];
03484         float pz = k*hz + origin[2];  // coordinates of lattice point
03485 
03486         memset(u, 0, num_conformers * sizeof(float));
03487 
03488         for (n = 0;  n < num_coords;  n++) {  // all atoms
03489 
03490           for (m = 0;  m < num_conformers;  m++) {  // conformers
03491             float v = 0.f;
03492             for (ma = 0;  ma < num_probe_atoms;  ma++) {  // probe atoms
03493               int index = m*num_probe_atoms + ma;
03494               float dx = conformers[3*index    ] + px - coords[3*n    ];
03495               float dy = conformers[3*index + 1] + py - coords[3*n + 1];
03496               float dz = conformers[3*index + 2] + pz - coords[3*n + 2];
03497               float r2 = dx*dx + dy*dy + dz*dz;
03498               if (r2 >= cutoff2) continue;
03499               int pindex = 2 * vdw_type[n];
03500               float epsilon = vdw_params[pindex] * probe_vdw_params[2*ma];
03501               float rmin = vdw_params[pindex + 1] + probe_vdw_params[2*ma + 1];
03502               float rm6 = rmin*rmin / r2;
03503               rm6 = rm6 * rm6 * rm6;
03504               v += epsilon * rm6 * (rm6 - 2.f);  // sum vdw contribution
03505             } // end loop probe atoms
03506 
03507             u[m] += v;  // contribution of one system atom to conformer
03508 
03509           } // end loop conformers
03510 
03511         } // end loop over all atoms
03512 
03513         float z = 0.f;
03514         for (m = 0;  m < num_conformers;  m++) {  // average conformer energies
03515           z += expf(-u[m]);
03516         }
03517 
03518         map[mindex] = z / num_conformers;  // store average occupancy
03519         if (map[mindex] < minocc) map[mindex] = 0.f;
03520       }
03521     }
03522   } // end loop over map
03523 
03524   delete[] u;
03525 }
03526 

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