GridForceGrid.h

Go to the documentation of this file.
00001 
00007 #ifndef GRIDFORCEGRID_H
00008 #define GRIDFORCEGRID_H
00009 
00010 #include <set>
00011 
00012 #include "Vector.h"
00013 #include "Tensor.h"
00014 #include "SimParameters.h"
00015 #include "NamdTypes.h"
00016 #include "MStream.h"
00017 #include "charm++.h"
00018 //#include "ComputeGridForce.h"
00019 
00020 #include "MGridforceParams.h"
00021 
00022 
00023 class GridforceFullMainGrid;
00024 class GridforceFullSubGrid;
00025 
00026 
00027 // GridforceGrid is now an abstract class to act as an interface to both GridforceFullMainGrid's and GridforceLiteGrid's
00028 class GridforceGrid {
00029 public:
00030     static GridforceGrid * new_grid(int gridnum, char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams);
00031     virtual ~GridforceGrid();
00032     
00033     virtual void initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams) = 0;
00034     virtual void reinitialize(SimParameters *simParams, MGridforceParams *mgridParams) = 0;
00035     
00036     virtual Position get_center(void) const = 0;
00037     virtual Position get_origin(void) const = 0;
00038     virtual Tensor get_e (void) const = 0;
00039     virtual Tensor get_inv(void) const = 0;
00040     virtual Vector get_scale(void) const = 0;
00041     virtual Bool get_checksize(void) const = 0;
00042     virtual int get_k0(void) const = 0;
00043     virtual int get_k1(void) const = 0; 
00044     virtual int get_k2(void) const = 0;
00045     virtual int get_total_grids(void) const = 0;
00046     
00047     virtual long int get_all_gridvals(float** all_gridvals) const = 0;
00048     virtual void set_all_gridvals(float* all_gridvals, long int sz) = 0;
00049     virtual void set_scale(Vector s) = 0;
00050     
00051     Position wrap_position(const Position &pos, const Lattice &lattice);
00052     bool fits_lattice(const Lattice &lattice);
00053     
00054     inline int compute_VdV(Position pos, float &V, Vector &dV) const { return -1; }
00055 
00056     static void pack_grid(GridforceGrid *grid, MOStream *msg);
00057     static GridforceGrid * unpack_grid(int gridnum, MIStream *msg);
00058     
00059     typedef enum {
00060         GridforceGridTypeUndefined = 0,
00061         GridforceGridTypeFull,
00062         GridforceGridTypeLite
00063     } GridforceGridType;
00064     
00065     inline GridforceGridType get_grid_type(void) { return type; }
00066 
00067 protected:    
00068     virtual void pack(MOStream *msg) const = 0;
00069     virtual void unpack(MIStream *msg) = 0;
00070     
00071     Position get_corner(int idx);
00072     
00073     GridforceGrid() { type = GridforceGridTypeUndefined; }
00074     GridforceGridType type;
00075     int mygridnum;
00076     
00077 private:
00078     Vector corners[8];
00079 };
00080 
00081 
00082 class GridforceFullBaseGrid {
00083     friend class GridforceFullMainGrid;
00084     friend class GridforceFullSubGrid;
00085     
00086 public:
00087     GridforceFullBaseGrid(void);
00088     virtual ~GridforceFullBaseGrid();
00089     
00090 //    int request_box(Vector pos);
00091 //    int get_box(Box *box, Vector pos) const;
00092   
00093     inline Position get_center(void) const { return center; }
00094     inline Position get_origin(void) const { return origin; }
00095     inline Tensor get_e (void) const { return e; }
00096     inline Tensor get_inv(void) const { return inv; }
00097     inline Vector get_scale(void) const { return scale; }
00098     inline Bool get_checksize(void) const { return checksize; }
00099     virtual int get_border(void) const = 0;
00100     
00101     inline float get_grid(int i0, int i1, int i2) const {
00102         return grid[grid_index(i0, i1, i2)];
00103     }
00104     inline double get_grid_d(int i0, int i1, int i2) const {
00105         return double(get_grid(i0, i1, i2));
00106     }
00107     inline void set_grid(int i0, int i1, int i2, float V) {
00108         grid[grid_index(i0, i1, i2)] = V;
00109     }
00110     
00111     inline void set_scale(Vector s) { scale = s; }
00112     
00113     int compute_VdV(Position pos, float &V, Vector &dV) const;
00114     
00115     inline int get_k0(void) const { return k[0]; }
00116     inline int get_k1(void) const { return k[1]; }
00117     inline int get_k2(void) const { return k[2]; }
00118     
00119 protected:
00120     virtual void pack(MOStream *msg) const;
00121     virtual void unpack(MIStream *msg);
00122     
00123     struct GridIndices {
00124       int inds2;
00125       int dk_hi;
00126       int dk_lo;
00127       Bool zero_derivs;
00128     };
00129    
00130     // Utility functions
00131     void readHeader(SimParameters *simParams, MGridforceParams *mgridParams);
00132     
00133     inline long int grid_index(int i0, int i1, int i2) const {
00134         register int inds[3] = {i0, i1, i2};
00135 #ifdef DEBUGM
00136         if (i0 < 0 || i0 >= k[0] || i1 < 0 || i1 >= k[1] || i2 < 0 || i2 >= k[2]) {
00137             char buffer[256];
00138             sprintf(buffer, "Bad grid index! (%d %d %d)", i0, i1, i2);
00139             NAMD_bug(buffer);
00140         }
00141 #endif
00142         return inds[0]*dk[0] + inds[1]*dk[1] + inds[2]*dk[2];
00143     }
00144     
00145     //virtual int get_inds(Position pos, int *inds, Vector &dg, Vector &gapscale) const = 0;
00146     int get_inds(Position pos, int *inds, Vector &dg, Vector &gapscale) const;
00147     void compute_a(float *a, float *b) const;
00148     virtual void compute_b(float *b, int *inds, Vector gapscale)  const = 0;
00149     float compute_V(float *a, float *x, float *y, float *z) const;
00150     Vector compute_dV(float *a, float *x, float *y, float *z) const;
00151     Vector compute_d2V(float *a, float *x, float *y, float *z) const;
00152     float compute_d3V(float *a, float *x, float *y, float *z) const;
00153     
00154     void readSubgridHierarchy(FILE *poten, int &totalGrids);
00155     
00156     FILE *poten_fp;
00157     float *grid;        // Actual grid
00158     
00159     GridforceFullSubGrid **subgrids;
00160     int numSubgrids;
00161     int generation;     // Subgrid level (0 = main grid)
00162     
00163     // should move 'nopad' versions to maingrid only ... or not have them as ivars at all, why are they here?
00164     int k[3];           // Grid dimensions
00165     int k_nopad[3];     // Grid dimensions
00166     long int size;
00167     long int size_nopad;
00168     long int dk[3];
00169     long int dk_nopad[3];
00170     float factor;
00171     
00172     Position origin;    // Grid origin
00173     Position center;    // Center of grid (for wrapping)
00174     Tensor e;           // Grid unit vectors
00175     Tensor inv;         // Inverse of unit vectors
00176     
00177     double p_sum[3];     // Accumulators for sums
00178     double n_sum[3];
00179     double pad_p[3];    // Pad values (p = positive side, n = negative side) for each dimension
00180     double pad_n[3];
00181     Bool cont[3];       // Whether grid is continuous in each dimension
00182     float offset[3];    // Potential offset in each dimension
00183     float gap[3];       // Gap between images of grid in grid units for each dimension
00184     float gapinv[3];    // 1.0/gap
00185 
00186     Vector scale;
00187     Bool checksize;
00188 };
00189 
00190 
00191 class GridforceFullMainGrid : public GridforceGrid, public GridforceFullBaseGrid {
00192     friend class GridforceFullBaseGrid;
00193     friend class GridforceFullSubGrid;
00194 
00195 public:
00196     explicit GridforceFullMainGrid(int gridnum);
00197     virtual ~GridforceFullMainGrid();
00198     
00199     void initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams, int border);
00200     inline void initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams) {
00201         initialize(potfilename, simParams, mgridParams, default_border);
00202     }
00203     void reinitialize(SimParameters *simParams, MGridforceParams *mgridParams);
00204     
00205     inline Position get_center(void) const { return GridforceFullBaseGrid::get_center(); };
00206     inline Position get_origin(void) const { return GridforceFullBaseGrid::get_origin(); };
00207     inline Tensor get_e (void) const { return GridforceFullBaseGrid::get_e(); };
00208     inline Tensor get_inv(void) const { return GridforceFullBaseGrid::get_inv(); };
00209     inline Vector get_scale(void) const { return GridforceFullBaseGrid::get_scale(); };
00210     inline Bool get_checksize(void) const { return GridforceFullBaseGrid::get_checksize(); };
00211     inline int get_k0(void) const { return GridforceFullBaseGrid::get_k0(); };
00212     inline int get_k1(void) const { return GridforceFullBaseGrid::get_k1(); };
00213     inline int get_k2(void) const { return GridforceFullBaseGrid::get_k2(); };
00214     inline int get_border(void) const { return border; }
00215     
00216     inline int compute_VdV(Position pos, float &V, Vector &dV) const { return GridforceFullBaseGrid::compute_VdV(pos, V, dV); };
00217     
00218     inline int get_total_grids(void) const { return totalGrids; }    
00219     inline void set_scale(Vector s) { scale = s; }
00220 
00221 protected:
00222     void pack(MOStream *msg) const;
00223     void unpack(MIStream *msg);
00224     
00225     long int get_all_gridvals(float **all_gridvals) const;
00226     void set_all_gridvals(float *all_gridvals, long int sz);
00227     
00228     //int get_inds(Position pos, int *inds, Vector &dg, Vector &gapscale) const;
00229     void compute_b(float *b, int *inds, Vector gapscale)  const;
00230     void buildSubgridsFlat(void);
00231     
00232     char filename[129];
00233     int totalGrids;
00234     GridforceFullSubGrid **subgrids_flat;
00235     //int mygridnum;
00236     
00237     static const int default_border = 1;
00238     int border;
00239 };
00240 
00241 
00242 class GridforceFullSubGrid : public GridforceFullBaseGrid {
00243     friend class GridforceFullBaseGrid;
00244     friend class GridforceFullMainGrid;
00245 
00246 public:
00247     GridforceFullSubGrid(GridforceFullBaseGrid *parent_in);
00248     
00249     void initialize(SimParameters *simParams, MGridforceParams *mgridParams);
00250     
00251     inline int get_border(void) const { return 0; }
00252     
00253     inline Tensor tensorMult (const Tensor &t1, const Tensor &t2) {
00254         Tensor tmp;
00255         tmp.xx = t1.xx * t2.xx + t1.xy * t2.yx + t1.xz * t2.zx;
00256         tmp.xy = t1.xx * t2.xy + t1.xy * t2.yy + t1.xz * t2.zy;
00257         tmp.xz = t1.xx * t2.xz + t1.xy * t2.yz + t1.xz * t2.zz;
00258         tmp.yx = t1.yx * t2.xx + t1.yy * t2.yx + t1.yz * t2.zx;
00259         tmp.yy = t1.yx * t2.xy + t1.yy * t2.yy + t1.yz * t2.zy;
00260         tmp.yz = t1.yx * t2.xz + t1.yy * t2.yz + t1.yz * t2.zz;
00261         tmp.zx = t1.zx * t2.xx + t1.zy * t2.yx + t1.zz * t2.zx;
00262         tmp.zy = t1.zx * t2.xy + t1.zy * t2.yy + t1.zz * t2.zy;
00263         tmp.zz = t1.zx * t2.xz + t1.zy * t2.yz + t1.zz * t2.zz;
00264         return tmp;
00265     }
00266     
00267 protected:
00268     void pack(MOStream *msg) const;
00269     void unpack(MIStream *msg);
00270     
00271     //int get_inds(Position pos, int *inds, Vector &dg, Vector &gapscale) const;
00272     void compute_b(float *b, int *inds, Vector gapscale) const;
00273     void addToSubgridsFlat(void);
00274     
00275     // Utility numbers
00276     Tensor scale_dV;
00277     Tensor scale_d2V;
00278     float scale_d3V;
00279     
00280     GridforceFullBaseGrid *parent;
00281     int pmin[3], pmax[3];
00282     GridforceFullMainGrid *maingrid;
00283     int subgridIdx;
00284 };
00285 
00286 
00287 class GridforceLiteGrid : public GridforceGrid {
00288 public:
00289     explicit GridforceLiteGrid(int gridnum);
00290     virtual ~GridforceLiteGrid();
00291     
00292     void initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams);
00293     void reinitialize(SimParameters *simParams, MGridforceParams *mgridParams);
00294     
00295     inline Position get_center(void) const { return center; }
00296     inline Position get_origin(void) const { return origin; }
00297     inline Tensor get_e (void) const { return e; }
00298     inline Tensor get_inv(void) const { return inv; }
00299     inline Vector get_scale(void) const { return scale; }
00300     inline Bool get_checksize(void) const { return checksize; }
00301     inline int get_k0(void) const { return k[0]; }
00302     inline int get_k1(void) const { return k[1]; }
00303     inline int get_k2(void) const { return k[2]; }
00304     inline int get_total_grids(void) const { return 1; }
00305     inline void set_scale(Vector s) { scale = s; }
00306     
00307     inline float get_grid(int i0, int i1, int i2, int i3) const {
00308         return grid[grid_index(i0, i1, i2, i3)];
00309     }
00310     inline double get_grid_d(int i0, int i1, int i2, int i3) const {
00311         return double(grid[grid_index(i0, i1, i2, i3)]);
00312     }
00313     inline void set_grid(int i0, int i1, int i2, int i3, float V) {
00314         grid[grid_index(i0, i1, i2, i3)] = V;
00315     }
00316     
00317     long int get_all_gridvals(float** all_gridvals) const;
00318     void set_all_gridvals(float* all_gridvals, long int sz);
00319     
00320     int compute_VdV(Position pos, float &V, Vector &dV) const;
00321     
00322 protected:
00323     void compute_derivative_grids(void);
00324     void compute_wts(float *wts, const Vector &dg) const;
00325     int get_inds(Position pos, int *inds, Vector &dg) const;
00326     float linear_interpolate(int i0, int i1, int i2, int i3, const float *wts) const;
00327     
00328     void pack(MOStream *msg) const;
00329     void unpack(MIStream *msg);
00330     
00331     inline long int grid_index(int i0, int i1, int i2, int i3) const {
00332         // 'i3' is an index for the grid itself (0=V, 1=dV/dx, 2=dV/dy, 3=dV/dz)
00333         register int inds[4] = {i0, i1, i2, i3};
00334         return inds[0]*dk[0] + inds[1]*dk[1] + inds[2]*dk[2] + inds[3]*dk[3];
00335     }
00336     
00337     float *grid;
00338     
00339     int k[4];           // Grid dimensions ... 4th is always 4, for the different grid types
00340     long int size;
00341     long int dk[4];
00342     
00343     Position origin;    // Grid origin
00344     Position center;    // Center of grid (for wrapping)
00345     Tensor e;           // Grid unit vectors
00346     Tensor inv;         // Inverse of unit vectors
00347     
00348     Vector scale;
00349     Bool checksize;
00350     
00351     char filename[129];
00352 };
00353 
00354 #endif

Generated on Thu Sep 21 01:17:12 2017 for NAMD by  doxygen 1.4.7