VolMapCreateILS.C

00001 /***************************************************************************
00002 *cr
00003 *cr (C) Copyright 1995-2019 The Board of Trustees of the
00004 *cr University of Illinois
00005 *cr All Rights Reserved
00006 *cr
00007 ***************************************************************************/
00008 /***************************************************************************
00009 * RCS INFORMATION:
00010 *
00011 * $RCSfile: 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