00001 /*************************************************************************** 00002 *cr 00003 *cr (C) Copyright 1995-2019 The Board of Trustees of the 00004 *cr University of Illinois 00005 *cr All Rights Reserved 00006 *cr 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * RCS INFORMATION: 00011 * 00012 * $RCSfile: VolMapCreate.h,v $ 00013 * $Author: johns $ $Locker: $ $State: Exp $ 00014 * $Revision: 1.93 $ $Date: 2019年01月23日 21:33:54 $ 00015 * 00016 **************************************************************************/ 00017 00018 #include "Matrix4.h" 00019 00020 // enable multilevel summation by default 00021 #define VMDUSEMSMPOT 1 00022 00023 class VMDApp; 00024 class VolumetricData; 00025 class AtomSel; 00026 class Molecule; 00027 00033 00034 class VolMapCreate { 00035 public: 00036 typedef enum {COMBINE_AVG, COMBINE_MIN, COMBINE_MAX, COMBINE_STDEV, COMBINE_PMF} CombineType; 00037 00038 protected: 00039 VMDApp *app; 00040 AtomSel *sel; 00041 float delta; // resolution (same along x, y and z) 00042 int computed_frames; // frame counter 00043 int checkpoint_freq; // write checkpoint file every xxx steps 00044 char *checkpoint_name; // checkpoint file name 00045 bool user_minmax; // true = user specified a minmax box, false = compute default minmax 00046 float min_coord[3], max_coord[3]; // used to pass user defaults, avoid using for computations! 00047 00048 protected: 00049 virtual int compute_frame(int frame, float *voldata) = 0; 00050 int compute_init(float padding); 00051 00053 virtual int compute_init() {return compute_init(0.);} 00054 00056 int calculate_minmax (float *min_coord, float *max_coord); 00057 00059 int calculate_max_radius (float &radius); 00060 00062 void combo_begin(CombineType method, void **customptr, void *params); 00063 void combo_addframe(CombineType method, float *voldata, void *customptr, float *framedata); 00064 void combo_export(CombineType method, float *voldata, void *customptr); 00065 void combo_end(CombineType method, void *customptr); 00066 00067 00068 public: 00069 VolumetricData *volmap; 00070 00071 VolMapCreate(VMDApp *app, AtomSel *sel, float resolution); 00072 virtual ~VolMapCreate(); 00073 00074 void set_minmax (float minx, float miny, float minz, float maxx, float maxy, float maxz); 00075 00076 void set_checkpoint (int checkpointfreq, char *checkpointname); 00077 00078 int compute_all(bool allframes, CombineType method, void *params); 00079 00084 virtual void write_map(const char *filename); 00085 00086 // We temporarily need our own file writer until we use molfile plugin 00087 int write_dx_file (const char *filename); 00088 00089 }; 00090 00091 00092 class VolMapCreateMask: public VolMapCreate { 00093 protected: 00094 int compute_init(); 00095 int compute_frame(int frame, float *voldata); 00096 private: 00097 float atomradius; 00098 00099 public: 00100 VolMapCreateMask(VMDApp *app, AtomSel *sel, float res, float the_atomradius) : VolMapCreate(app, sel, res) { 00101 atomradius = the_atomradius; 00102 } 00103 }; 00104 00105 00106 class VolMapCreateDensity : public VolMapCreate { 00107 protected: 00108 float *weight; 00109 char const *weight_string; 00110 int weight_mutable; 00111 int compute_init(); 00112 int compute_frame(int frame, float *voldata); 00113 float radius_scale; // mult. factor for atomic radii 00114 00115 public: 00116 VolMapCreateDensity(VMDApp *app, AtomSel *sel, float res, float *the_weight, char const *the_weight_string, int the_weight_mutable, float the_radscale) : VolMapCreate(app, sel, res) { 00117 weight = the_weight; 00118 weight_string = the_weight_string; 00119 weight_mutable = the_weight_mutable; 00120 // number of random points to use for each atom's gaussian distr. 00121 radius_scale = the_radscale; 00122 } 00123 }; 00124 00125 00126 class VolMapCreateInterp : public VolMapCreate { 00127 protected: 00128 float *weight; 00129 char const *weight_string; 00130 int weight_mutable; 00131 int compute_init(); 00132 int compute_frame(int frame, float *voldata); 00133 00134 public: 00135 VolMapCreateInterp(VMDApp *app, AtomSel *sel, float res, float *the_weight, char const *the_weight_string, int the_weight_mutable) : VolMapCreate(app, sel, res) { 00136 weight = the_weight; 00137 weight_string = the_weight_string; 00138 weight_mutable = the_weight_mutable; 00139 } 00140 }; 00141 00142 00143 class VolMapCreateOccupancy : public VolMapCreate { 00144 private: 00145 bool use_points; 00146 protected: 00147 int compute_init(); 00148 int compute_frame(int frame, float *voldata); 00149 public: 00150 VolMapCreateOccupancy(VMDApp *app, AtomSel *sel, float res, bool use_point_particles) : VolMapCreate(app, sel, res) { 00151 use_points = use_point_particles; 00152 } 00153 }; 00154 00155 00156 class VolMapCreateDistance : public VolMapCreate { 00157 protected: 00158 float max_dist; 00159 int compute_init(); 00160 int compute_frame(int frame, float *voldata); 00161 public: 00162 VolMapCreateDistance(VMDApp *app, AtomSel *sel, float res, float the_max_dist) : VolMapCreate(app, sel, res) { 00163 max_dist = the_max_dist; 00164 } 00165 }; 00166 00167 00168 class VolMapCreateCoulombPotential : public VolMapCreate { 00169 protected: 00170 int compute_init(); 00171 int compute_frame(int frame, float *voldata); 00172 00173 public: 00174 VolMapCreateCoulombPotential(VMDApp *app, AtomSel *sel, float res) : VolMapCreate(app, sel, res) { 00175 } 00176 }; 00177 00178 00179 #if defined(VMDUSEMSMPOT) 00180 class VolMapCreateCoulombPotentialMSM : public VolMapCreate { 00181 protected: 00182 int compute_init(); 00183 int compute_frame(int frame, float *voldata); 00184 00185 public: 00186 VolMapCreateCoulombPotentialMSM(VMDApp *app, AtomSel *sel, float res) : VolMapCreate(app, sel, res) { 00187 } 00188 }; 00189 #endif 00190 00191 00195 class VolMapCreateILS { 00196 private: 00197 VMDApp *app; 00198 int molid; // the molecule we are operating on 00199 00200 int num_atoms; // # atoms in the system 00201 00202 VolumetricData *volmap; // our result: the free energy map 00203 VolumetricData *volmask; // mask defining valid gridpoints in volmap 00204 00205 float delta; // distance of samples for ILS computation 00206 int nsubsamp; // # samples in each dim. downsampled into 00207 // each gridpoint of the final map. 00208 00209 // Number of samples used during computation. 00210 int nsampx, nsampy, nsampz; 00211 00212 float minmax[6]; // minmax coords of bounding box 00213 float gridorigin[3]; // center of the first grid cell 00214 00215 float cutoff; // max interaction dist between any 2 atoms 00216 float extcutoff; // cutoff corrected for the probe size 00217 float excl_dist; // cutoff for the atom clash pre-scanning 00218 00219 bool compute_elec; // compute electrostatics? (currently unused) 00220 00221 // Control of the angular spacing of probe orientation vectors: 00222 // 1 means using 1 orientation only 00223 // 2 corresponds to 6 orientations (vertices of octahedron) 00224 // 3 corresponds to 8 orientations (vertices of hexahedron) 00225 // 4 corresponds to 12 orientations (faces of dodecahedron) 00226 // 5 corresponds to 20 orientations (vertices of dodecahedron) 00227 // 6 corresponds to 32 orientations (faces+vert. of dodecah.) 00228 // 7 and above: geodesic subdivisions of icosahedral faces 00229 // with frequency 1, 2, ... 00230 // Probes with tetrahedral symmetry: 00231 // # number of rotamers for each of the 8 orientations 00232 // (vertices of tetrahedron and its dual tetrahedron). 00233 // 00234 // Note that the angular spacing of the rotations around 00235 // the orientation vectors is chosen to be about the same 00236 // as the angular spacing of the orientation vector 00237 // itself. 00238 int conformer_freq; 00239 00240 int num_conformers; // # probe symmetry unique orientations and 00241 // rotations sampled per grid point 00242 float *conformers; // Stores the precomputed atom positions 00243 // (relative to the center of mass) 00244 // of the different probe orientations and 00245 // rotations. 00246 int num_orientations; // # probe symmetry unique orientations 00247 int num_rotations; // # symmetry unique rotations sampled 00248 // per orientation 00249 00250 // We store the VDW parameters once for each type: 00251 float *vdwparams; // VDW well depths and radii for all types 00252 int *atomtypes; // index list for vdw parameter types 00253 00254 int num_unique_types; // # unique atom types 00255 00256 float temperature; // Temp. in Kelvin at which the MD sim. was performed 00257 00258 int num_probe_atoms; // # atoms in the probe (the ligand) 00259 float probe_effsize; // effective probe radius 00260 float *probe_coords; // probe coordinates 00261 00262 // The two highest symmetry axes of the probe 00263 float probe_symmaxis1[3]; 00264 float probe_symmaxis2[3]; 00265 int probe_axisorder1, probe_axisorder2; 00266 int probe_tetrahedralsymm; // probe has tetrahedral symmetry flag 00267 00268 // VDW parameters for the probe: 00269 // A tuple of eps and rmin is stored for each atom. 00270 // Actually we store beta*sqrt(eps) and rmin/2 00271 // (see function set_probe()). 00272 float *probe_vdw; 00273 float *probe_charge; // charge for each probe atom 00274 00275 int first, last; // trajectory frame range 00276 int computed_frames; // # frames processed 00277 00278 float max_energy; // max energy considered in map, all higher energies 00279 // will be clamped to this value. 00280 float min_occup; // occupancies below this value will be treated 00281 // as zero. 00282 00283 bool pbc; // If flag is set then periodic boundaries are taken 00284 // into account. 00285 bool pbcbox; // If flag is set then the grid dimensions will be chosen 00286 // as the orthogonalized bounding box for the PBC cell. 00287 float pbccenter[3]; // User provided PBC cell center. 00288 00289 AtomSel *alignsel; // Selection to be used for alignment 00290 const float *alignrefpos; // Stores the alignment reference position 00291 00292 Matrix4 transform; // Transformation matrix that was used for the 00293 // alignment of the first frame. 00294 00295 int maskonly; // If set, compute only a mask map telling for which 00296 // gridpoints we expect valid energies, i.e. the points 00297 // for which the maps overlap for all frames. 00298 00299 00300 // Check if the box given by the minmax coordinates is located 00301 // entirely inside the PBC unit cell of the given frame and in 00302 // this case return 1, otherwise return 0. 00303 int box_inside_pbccell(int frame, float *minmax); 00304 00305 // Check if the entire volmap grid is located entirely inside 00306 // the PBC unit cell of the given frame (taking the alignment 00307 // into account) and in this case return 1, otherwise return 0. 00308 int grid_inside_pbccell(int frame, float *voldata, 00309 const Matrix4 &alignment); 00310 00311 // Set grid dimensions to the minmax coordinates and 00312 // align grid with integer coordinates. 00313 int set_grid(); 00314 00315 // Initialize the ILS calculation 00316 int initialize(); 00317 00318 // ILS calculation for the given frame 00319 int compute_frame(int frame, float *voldata); 00320 00321 // Align current frame to the reference 00322 void align_frame(Molecule *mol, int frame, float *coords, 00323 Matrix4 &alignment); 00324 00325 // Get array of coordinates of selected atoms and their 00326 // neighbors (within a cutoff) in the PBC images. 00327 int get_atom_coordinates(int frame, Matrix4 &alignment, 00328 int *(&vdwtypes), 00329 float *(&coords)); 00330 00331 // Check if probe is a linear molecule and returns 00332 // the Cinf axis. 00333 int is_probe_linear(float *axis); 00334 00335 // Simple probe symmetry check 00336 void check_probe_symmetry(); 00337 00338 // Determine probe symmetry and generate probe orientations 00339 // and rotations. 00340 void initialize_probe(); 00341 void get_eff_proberadius(); 00342 00343 00344 // Generate conformers for tetrahedral symmetry 00345 int gen_conf_tetrahedral(float *(&conform), int freq, 00346 int &numorient, int &numrot); 00347 00348 // Generate conformers for all other symmetries 00349 int gen_conf(float *(&conform), int freq, 00350 int &numorient, int &numrot); 00351 00352 float dimple_depth(float phi); 00353 00354 // Create list of unique VDW parameters which can be accessed 00355 // through the atomtypes index list. 00356 int create_unique_paramlist(); 00357 00358 00359 public: 00360 VolMapCreateILS(VMDApp *_app, int molid, int firstframe, 00361 int lastframe, float T, float res, 00362 int subr, float cut, int maskonly); 00363 ~VolMapCreateILS(); 00364 00365 VolumetricData* get_volmap() { return volmap; }; 00366 00367 // Perform ILS calculation for all specified frames. 00368 int compute(); 00369 00371 int add_map_to_molecule(); 00372 00375 int write_map(const char *filename); 00376 00377 // Set probe coordinates,charges and VDW parameters 00378 void set_probe(int num_probe_atoms, int num_conf, 00379 const float *probe_coords, 00380 const float *vdwrmin, const float *vdweps, 00381 const float *charge); 00382 00383 // Set the two highest symmetry axes for the probe and a flag 00384 // telling if we have a tetrahedral symmetry. 00385 // If the axes are not orthogonal the lower axis will be ignored. 00386 void set_probe_symmetry(int order1, const float *axis1, 00387 int order2, const float *axis2, 00388 int tetrahedral); 00389 00390 // Set minmax coordinates of rectangular molecule bounding box 00391 void set_minmax (float minx, float miny, float minz, float maxx, float maxy, float maxz); 00392 00393 // Request PBC aware computation. 00394 void set_pbc(float center[3], int bbox); 00395 00396 // Set maximum energy considered in the calculation. 00397 void set_maxenergy(float maxenergy); 00398 00399 // Set selection to be used for alignment. 00400 void set_alignsel(AtomSel *asel); 00401 00402 // Set transformation matrix that was used for the 00403 // alignment of the first frame. 00404 void set_transform(const Matrix4 *mat); 00405 00406 int get_conformers(float *&conform) const { 00407 conform = conformers; 00408 return num_conformers; 00409 } 00410 00411 void get_statistics(int &numconf, int &numorient, 00412 int &numrot) { 00413 numconf = num_conformers; 00414 numorient = num_orientations; 00415 numrot = num_rotations; 00416 } 00417 }; 00418 00419 00420 // Write given map as a DX file. 00421 int volmap_write_dx_file (VolumetricData *volmap, const char *filename);