GridForceGrid.C

Go to the documentation of this file.
00001 
00007 #include <iostream>
00008 #include <typeinfo>
00009 
00010 #include "GridForceGrid.h"
00011 #include "Vector.h"
00012 #include "Node.h"
00013 #include "SimParameters.h"
00014 #include "Molecule.h"
00015 #include "InfoStream.h"
00016 #include "common.h"
00017 #include "ComputeGridForce.h"
00018 
00019 #include "MGridforceParams.h"
00020 
00021 #define MIN_DEBUG_LEVEL 4
00022 //#define DEBUGM
00023 #include "Debug.h"
00024 
00025 #include "GridForceGrid.inl"
00026 
00027 
00028 /*****************/
00029 /* GRIDFORCEGRID */
00030 /*****************/
00031 
00032 // Class methods
00033 
00034 GridforceGrid * GridforceGrid::new_grid(int gridnum, char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
00035 {
00036     GridforceGrid *grid = NULL;
00037     if (mgridParams->gridforceLite) {
00038         grid = new GridforceLiteGrid(gridnum);
00039     } else {
00040         grid = new GridforceFullMainGrid(gridnum);
00041     }
00042     
00043     grid->initialize(potfilename, simParams, mgridParams);
00044     
00045     return grid;
00046 }
00047 
00048 GridforceGrid::~GridforceGrid() { ; }
00049 
00050 void GridforceGrid::pack_grid(GridforceGrid *grid, MOStream *msg)
00051 {
00052     // Abstract interface for packing a grid into a message.  This
00053     // could easily be a non-static function as it was before, but is
00054     // static to make it similar to unpack_grid below, which does need
00055     // to be static since it creates a GridforceGrid object.
00056     msg->put(grid->get_grid_type());
00057     grid->pack(msg);
00058 }
00059 
00060 GridforceGrid * GridforceGrid::unpack_grid(int gridnum, MIStream *msg)
00061 {
00062     // Abstract interface for unpacking a grid from a message.
00063     GridforceGrid *grid = NULL;
00064     int type;
00065     
00066     msg->get(type);
00067     
00068     switch (type) {
00069     case GridforceGridTypeFull:
00070         grid = new GridforceFullMainGrid(gridnum);
00071         break;
00072     case GridforceGridTypeLite:
00073         grid = new GridforceLiteGrid(gridnum);
00074         break;
00075     default:
00076         NAMD_bug("GridforceGrid::unpack_grid called with unknown grid type!");
00077     }
00078     
00079     grid->unpack(msg);
00080     
00081     return grid;
00082 }
00083 
00084 bool GridforceGrid::fits_lattice(const Lattice &lattice)
00085 {
00086     // Idea: Go through each grid corner and wrap it to the grid center;
00087     // if the position moves, then the grid is too big and we return false
00088     DebugM(4, "Checking whether grid fits periodic cell\n");
00089     Position center = get_center();
00090     for (int i = 0; i < 8; i++) {
00091         Position pos = get_corner(i);
00092         Position pos_wrapped = wrap_position(pos, lattice);
00093         if ((pos - pos_wrapped).length() > 1.) {
00094             DebugM(5, "(" << pos << ") != (" << pos_wrapped << ")\n" << endi);
00095             return false;
00096         }
00097     }
00098     return true;
00099 }
00100 
00101 Position GridforceGrid::get_corner(int idx)
00102 {
00103     // idx -> (x,y,z) (cell basis coordinates)
00104     // 0 -> (0,0,0)
00105     // 1 -> (1,0,0)
00106     // 2 -> (0,1,0)
00107     // 3 -> (1,1,0)
00108     // 4 -> (0,0,1)
00109     // 5 -> (1,0,1)
00110     // 6 -> (0,1,1)
00111     // 7 -> (1,1,1)
00112     Position pos;
00113     if (idx >= 8 || idx < 0) {
00114         // invalid index
00115         pos = Vector(); // default value of Vector() is strange enough to be a decent "undefined" value (-99999, -99999, -999999)
00116     } else if (corners[idx] != Vector()) {
00117         // use cached value if possible
00118         pos = corners[idx];
00119     } else {
00120         // must calculate
00121         Tensor e = get_e();
00122         pos = get_origin();
00123         if (idx & (1 << 0)) pos += e * Vector(get_k0()-1, 0, 0);
00124         if (idx & (1 << 1)) pos += e * Vector(0, get_k1()-1, 0);
00125         if (idx & (1 << 2)) pos += e * Vector(0, 0, get_k2()-1);
00126         corners[idx] = pos;     // cache for future use
00127         DebugM(4, "corner " << idx << " = " << pos << "\n" << endi);
00128     }
00129     return pos;
00130 }
00131 
00132 
00133 /*************************/
00134 /* GRIDFORCEFULLBASEGRID */
00135 /*************************/
00136 
00137 GridforceFullBaseGrid::GridforceFullBaseGrid(void)
00138 {
00139     cont[0] = cont[1] = cont[2] = FALSE;
00140     grid = NULL;
00141     numSubgrids = 0;
00142     subgrids = NULL;
00143 }
00144 
00145 GridforceFullBaseGrid::~GridforceFullBaseGrid()
00146 {
00147     delete[] grid;
00148     for (int i = 0; i < numSubgrids; i++) {
00149         delete subgrids[i];
00150     }
00151     delete[] subgrids;
00152 }
00153 
00154 
00155 void GridforceFullBaseGrid::pack(MOStream *msg) const
00156 {
00157     DebugM(2, "Packing message\n" << endi);
00158     
00159     msg->put(numSubgrids);
00160     msg->put(generation);
00161     
00162     msg->put(3*sizeof(int), (char*)k);
00163     msg->put(3*sizeof(int), (char*)k_nopad);
00164     msg->put(size);
00165     msg->put(size_nopad);
00166     msg->put(3*sizeof(long int), (char*)dk);
00167     msg->put(3*sizeof(long int), (char*)dk_nopad);
00168     msg->put(factor);
00169     
00170     msg->put(sizeof(Vector), (char*)&origin);
00171     msg->put(sizeof(Vector), (char*)&center);
00172     msg->put(sizeof(Tensor), (char*)&e);
00173     msg->put(sizeof(Tensor), (char*)&inv);
00174              
00175 //    msg->put(3*sizeof(float), (char*)pad_p);
00176 //    msg->put(3*sizeof(float), (char*)pad_n);
00177     msg->put(3*sizeof(Bool), (char*)cont);
00178     msg->put(3*sizeof(float), (char*)offset);
00179     msg->put(3*sizeof(float), (char*)gap);
00180     msg->put(3*sizeof(float), (char*)gapinv);
00181     msg->put(sizeof(Vector), (char*)&scale);
00182     msg->put(sizeof(Bool), (char*)&checksize);
00183     
00184     DebugM(2, "Packing grid, size = " << size << "\n" << endi);
00185     
00186     msg->put(size*sizeof(float), (char*)grid);
00187     
00188     DebugM(2, "Packing subgrids\n" << endi);
00189     
00190     for (int i = 0; i < numSubgrids; i++) {
00191         subgrids[i]->pack(msg);
00192     }
00193 }
00194 
00195 void GridforceFullBaseGrid::unpack(MIStream *msg)
00196 {
00197     DebugM(3, "Unpacking message\n" << endi);
00198 //    iout << iINFO << CkMyPe() << " Unpacking message\n" << endi;
00199 
00200     delete[] grid;
00201     grid = NULL;
00202     for (int i = 0; i < numSubgrids; i++) {
00203         delete subgrids[i];
00204     }
00205     numSubgrids = 0;
00206     delete[] subgrids;
00207     subgrids = NULL;
00208     
00209     msg->get(numSubgrids);
00210     msg->get(generation);
00211     
00212     DebugM(3, "numSubgrids = " << numSubgrids << "\n");
00213     DebugM(3, "generation = " << generation << "\n" << endi);
00214     
00215     msg->get(3*sizeof(int), (char*)k);
00216     msg->get(3*sizeof(int), (char*)k_nopad);
00217     msg->get(size);
00218     msg->get(size_nopad);
00219     msg->get(3*sizeof(long int), (char*)dk);
00220     msg->get(3*sizeof(long int), (char*)dk_nopad);
00221     msg->get(factor);
00222     
00223     DebugM(3, "size = " << size << "\n" << endi);
00224     
00225     msg->get(sizeof(Vector), (char*)&origin);
00226     msg->get(sizeof(Vector), (char*)&center);
00227     msg->get(sizeof(Tensor), (char*)&e);
00228     msg->get(sizeof(Tensor), (char*)&inv);
00229              
00230 //    msg->get(3*sizeof(float), (char*)pad_p);
00231 //    msg->get(3*sizeof(float), (char*)pad_n);
00232     msg->get(3*sizeof(Bool), (char*)cont);
00233     msg->get(3*sizeof(float), (char*)offset);
00234     msg->get(3*sizeof(float), (char*)gap);
00235     msg->get(3*sizeof(float), (char*)gapinv);
00236     msg->get(sizeof(Vector), (char*)&scale);
00237     msg->get(sizeof(Bool), (char*)&checksize);
00238     
00239     if (size) {
00240         DebugM(3, "allocating grid, size = " << size << "\n" << endi);
00241         grid = new float[size];
00242         msg->get(size*sizeof(float), (char*)grid);
00243     }
00244     
00245     if (numSubgrids) {
00246         DebugM(3, "Creating subgrids array, size " << numSubgrids << "\n" << endi);
00247         subgrids = new GridforceFullSubGrid *[numSubgrids];
00248         for (int i = 0; i < numSubgrids; i++) {
00249             subgrids[i] = new GridforceFullSubGrid(this);
00250             subgrids[i]->unpack(msg);
00251         }
00252     }
00253 }
00254 
00255 
00256 void GridforceFullBaseGrid::readHeader(SimParameters *simParams, MGridforceParams *mgridParams)
00257 {
00258   char line[256];
00259   long int poten_offset;
00260   do {
00261       poten_offset = ftell(poten_fp);
00262       fgets(line, 256, poten_fp);       // Read comment lines
00263       DebugM(4, "Read line: " << line << endi);
00264   } while (line[0] == '#');
00265   fseek(poten_fp, poten_offset, SEEK_SET);
00266   
00267   // read grid dimensions
00268   fscanf(poten_fp, "object %*d class gridpositions counts %d %d %d\n",
00269          &k_nopad[0], &k_nopad[1], &k_nopad[2]);
00270   size_nopad = k_nopad[0] * k_nopad[1] * k_nopad[2];
00271   
00272   // Read origin
00273   fscanf(poten_fp, "origin %lf %lf %lf\n",
00274          &origin.x, &origin.y, &origin.z);
00275         
00276   // Read delta (unit vectors)
00277   // These are column vectors, so must fill gridfrcE tensor to reflect this
00278   fscanf(poten_fp, "delta %lf %lf %lf\n", &e.xx, &e.yx, &e.zx);
00279   fscanf(poten_fp, "delta %lf %lf %lf\n", &e.xy, &e.yy, &e.zy);
00280   fscanf(poten_fp, "delta %lf %lf %lf\n", &e.xz, &e.yz, &e.zz);
00281   
00282   center = origin + e * 0.5 
00283            * Position(k_nopad[0]-1, k_nopad[1]-1, k_nopad[2]-1);
00284   
00285   fscanf(poten_fp, "object %*d class gridconnections counts %*lf %*lf %*lf\n");
00286   fscanf(poten_fp, "object %*d class array type double rank 0 items %*d data follows\n");
00287     
00288   // Calculate inverse tensor
00289   BigReal det;
00290   det = e.xx*(e.yy*e.zz - e.yz*e.zy) - e.xy*(e.yx*e.zz - e.yz*e.zx) 
00291         + e.xz*(e.yx*e.zy - e.yy*e.zx);
00292   inv.xx =  (e.yy*e.zz - e.yz*e.zy)/det;
00293   inv.xy = -(e.xy*e.zz - e.xz*e.zy)/det;
00294   inv.xz =  (e.xy*e.yz - e.xz*e.yy)/det;
00295   inv.yx = -(e.yx*e.zz - e.yz*e.zx)/det;
00296   inv.yy =  (e.xx*e.zz - e.xz*e.zx)/det;
00297   inv.yz = -(e.xx*e.yz - e.xz*e.yx)/det;
00298   inv.zx =  (e.yx*e.zy - e.yy*e.zx)/det;
00299   inv.zy = -(e.xx*e.zy - e.xy*e.zx)/det;
00300   inv.zz =  (e.xx*e.yy - e.xy*e.yx)/det;
00301   
00302   DebugM(4, "origin = " << origin << "\n");
00303   DebugM(4, "e = " << e << "\n");
00304   DebugM(4, "inv = " << inv << "\n" << endi);
00305 }
00306 
00307 void GridforceFullBaseGrid::readSubgridHierarchy(FILE *poten, int &totalGrids)
00308 {
00309     DebugM(4, "Beginning of readSubgridHierarchy, generation = " << generation << ", totalGrids = " << totalGrids << "\n" << endi);
00310     
00311     int elems, generation_in;
00312     
00313     subgrids = new GridforceFullSubGrid *[numSubgrids];
00314     
00315     for (int i = 0; i < numSubgrids; i++) {
00316         subgrids[i] = new GridforceFullSubGrid(this);
00317         elems = fscanf(poten_fp, "# namdnugrid subgrid %d generation %d min %d %d %d max %d %d %d subgrids count %d\n",
00318                &subgrids[i]->subgridIdx, &generation_in,
00319                &subgrids[i]->pmin[0], &subgrids[i]->pmin[1], &subgrids[i]->pmin[2],
00320                &subgrids[i]->pmax[0], &subgrids[i]->pmax[1], &subgrids[i]->pmax[2],
00321                &subgrids[i]->numSubgrids);
00322         if (elems < 9) {
00323             char msg[256];
00324             sprintf(msg, "Problem reading Gridforce potential file! (%d < 9)", elems);
00325             NAMD_die(msg);
00326         }
00327         
00328         totalGrids++;
00329         
00330         if (subgrids[i]->subgridIdx != (totalGrids - 1)) {
00331             char msg[256];
00332             sprintf(msg, "Problem reading Gridforce potential file! (%d != %d)", subgrids[i]->subgridIdx, totalGrids - 1);
00333             NAMD_die(msg);
00334         }
00335         if (subgrids[i]->generation != generation_in) {
00336             char msg[256];
00337             sprintf(msg, "Problem reading Gridforce potential file! (%d != %d)", subgrids[i]->generation, generation_in);
00338             NAMD_die(msg);
00339         }
00340         
00341 //      DebugM(3, "setting maingrid\n");
00342 //      subgrids[i]->maingrid->subgrids_flat[subgrids[i]->subgridIdx] = subgrids[i];
00343 //      DebugM(3, "reading subgrid hierarchy\n");
00344         
00345         subgrids[i]->readSubgridHierarchy(poten, totalGrids);
00346     }
00347 }
00348 
00349 
00350 /*************************/
00351 /* GRIDFORCEFULLMAINGRID */
00352 /*************************/
00353 
00354 GridforceFullMainGrid::GridforceFullMainGrid(int gridnum)
00355 {
00356     mygridnum = gridnum;
00357     generation = 0;
00358     subgrids_flat = NULL;
00359     type = GridforceGridTypeFull;
00360 }
00361 
00362 
00363 GridforceFullMainGrid::~GridforceFullMainGrid()
00364 {
00365     delete[] subgrids_flat;
00366 }
00367 
00368 
00369 void GridforceFullMainGrid::pack(MOStream *msg) const
00370 {
00371     DebugM(4, "Packing maingrid\n" << endi);
00372     
00373 //     msg->put(3*sizeof(float), (char*)pad_p);
00374 //     msg->put(3*sizeof(float), (char*)pad_n);
00375     msg->put(totalGrids);
00376     msg->put(mygridnum);
00377     msg->put(129*sizeof(char), (char*)filename);
00378     
00379     DebugM(3, "calling GridforceFullBaseGrid::pack\n" << endi);
00380     
00381     GridforceFullBaseGrid::pack(msg);
00382 }
00383 
00384 
00385 void GridforceFullMainGrid::unpack(MIStream *msg)
00386 {
00387     DebugM(4, "Unpacking maingrid\n" << endi);
00388     
00389 //     msg->get(3*sizeof(float), (char*)pad_p);
00390 //     msg->get(3*sizeof(float), (char*)pad_n);
00391     msg->get(totalGrids);
00392     msg->get(mygridnum);
00393     msg->get(129*sizeof(char), (char*)filename);
00394     
00395     GridforceFullBaseGrid::unpack(msg);
00396     
00397     DebugM(4, "size  = " << size << "\n");
00398     DebugM(4, "numSubgrids = " << numSubgrids << "\n");
00399     DebugM(4, "gapinv = " << gapinv[0] << " " << gapinv[2] << " " << gapinv[2] << " " << "\n");
00400     DebugM(4, "generation = " << generation << "\n" << endi);
00401     DebugM(4, "filename = " << filename << "\n" << endi);
00402     
00403     buildSubgridsFlat();
00404 }
00405 
00406 
00407 void GridforceFullMainGrid::buildSubgridsFlat(void)
00408 {
00409     DebugM(4, "buildSubgridsFlat() called, totalGrids-1 = " << totalGrids-1 << "\n" << endi);
00410     delete[] subgrids_flat;
00411     subgrids_flat = new GridforceFullSubGrid *[totalGrids-1];
00412     for (int i = 0; i < numSubgrids; i++) {
00413         DebugM(3, "adding to subgridsFlat\n" << endi);
00414         subgrids[i]->addToSubgridsFlat();
00415         DebugM(3, "success!\n" << endi);
00416     }
00417     for (int i = 0; i < totalGrids-1; i++) {
00418         DebugM(4, "subgrids_flat[" << i << "]->numSubgrids = " << subgrids_flat[i]->numSubgrids << "\n" << endi);
00419     }
00420     for (int i = 0; i < numSubgrids; i++) {
00421         DebugM(4, "subgrids[" << i << "]->numSubgrids = " << subgrids[i]->numSubgrids << "\n" << endi);
00422     }
00423 }
00424 
00425 
00426 void GridforceFullMainGrid::initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams, int brd)
00427 {
00428     if (brd >= 0) {
00429         border = brd;
00430     } else {
00431         border = default_border;
00432     }
00433     
00434     // FROM init1
00435     //FILE *poten = Fopen(potfilename, "r");
00436     poten_fp = Fopen(potfilename, "r");
00437     if (!poten_fp) {
00438         NAMD_die("Problem reading grid force potential file");
00439     }
00440     
00441     // save file name so that grid can be re-read via Tcl
00442     strcpy(filename, potfilename);
00443     
00444     // Read special comment fields and create subgrid objects
00445     totalGrids = 1;
00446     char line[256];
00447     Bool flag = FALSE;
00448     numSubgrids = 0;
00449     float version;
00450     long int poten_offset;
00451     do {
00452         poten_offset = ftell(poten_fp);
00453         fgets(line, 256, poten_fp);     // Read comment lines
00454         //flag = sscanf(line, "# maingrid subgrids count %d\n", &numSubgrids);
00455         flag = sscanf(line, "# namdnugrid version %f\n", &version);
00456     } while (line[0] == '#' && !flag);
00457     
00458     if (flag) {
00459         if (version != 1.0) {
00460             NAMD_die("Unsupported version of non-uniform grid file format!");
00461         }
00462         fscanf(poten_fp, "# namdnugrid maingrid subgrids count %d\n", &numSubgrids);
00463         readSubgridHierarchy(poten_fp, totalGrids);
00464         buildSubgridsFlat();
00465     } else {
00466         fseek(poten_fp, poten_offset, SEEK_SET);
00467     }
00468     
00469     // Read header
00470     readHeader(simParams, mgridParams);
00471     
00472     factor = 1.0;
00473     if (mgridParams->gridforceVolts)
00474     {
00475         factor /= 0.0434;  // convert V -> kcal/mol*e
00476     }
00477     scale = mgridParams->gridforceScale;
00478     checksize = mgridParams->gridforceCheckSize;
00479 
00480     // Allocate storage for potential and read it
00481     float *grid_nopad = new float[size_nopad];
00482     
00483     float tmp2;
00484     for (long int count = 0; count < size_nopad; count++) {
00485         int err = fscanf(poten_fp, "%f", &tmp2);
00486         if (err == EOF || err == 0) {
00487             NAMD_die("Grid force potential file incorrectly formatted");
00488         }
00489         grid_nopad[count] = tmp2 * factor;      // temporary, so just store flat
00490     }
00491     fscanf(poten_fp, "\n");
00492     
00493     // Shortcuts for accessing 1-D array with four indices
00494     dk_nopad[0] = k_nopad[1] * k_nopad[2];
00495     dk_nopad[1] = k_nopad[2];
00496     dk_nopad[2] = 1;
00497     
00498     Vector Kvec[3];
00499     Kvec[0] = e * Position(k_nopad[0]-1, 0, 0);
00500     Kvec[1] = e * Position(0, k_nopad[1]-1, 0);
00501     Kvec[2] = e * Position(0, 0, k_nopad[2]-1);
00502     Vector Avec[3];
00503     Avec[0] = simParams->lattice.a();
00504     Avec[1] = simParams->lattice.b();
00505     Avec[2] = simParams->lattice.c();
00506     
00507     // Decide whether we're wrapping
00508     for (int i0 = 0; i0 < 3; i0++) {
00509         if (mgridParams->gridforceCont[i0])
00510         {
00511             Bool found = FALSE;
00512             for (int i1 = 0; i1 < 3; i1++) {
00513                 if (cross(Avec[i0].unit(), Kvec[i1].unit()).length() < 1e-4) {
00514                     found = TRUE;
00515                     cont[i1] = TRUE;
00516                     offset[i1] = mgridParams->gridforceVOffset[i0] * factor;
00517                     // want in grid-point units (normal = 1)
00518                     gap[i1] = (inv * (Avec[i0] - Kvec[i1])).length();   
00519                     gapinv[i1] = 1.0/gap[i1];
00520                     
00521                     if (gap[i1] < 0) {
00522                         NAMD_die("Gridforce Grid overlap!");
00523                     }
00524                     
00525                     DebugM(4, "cont[" << i1 << "] = " << cont[i1] << "\n");
00526                     DebugM(4, "gap[" << i1 << "] = " << gap[i1] << "\n");
00527                     DebugM(4, "gapinv[" << i1 << "] = " << gapinv[i1] << "\n" << endi);
00528                 }
00529             }
00530             
00531             if (!found) {
00532                 NAMD_die("No Gridforce unit vector found parallel to requested continuous grid direction!");
00533             }
00534         } else {
00535             // check for grid overlap in non-wrapping dimensions
00536             // occurs below
00537         }
00538     }
00539     
00540     // Figure out size of true grid (padded on non-periodic sides)
00541     Vector delta = 0;
00542     for (int i = 0; i < 3; i++) {
00543         if (cont[i]) {
00544             k[i] = k_nopad[i];
00545         } else {
00546             k[i] = k_nopad[i] + 2*border;
00547             delta[i] -= border;
00548         }
00549     }
00550     DebugM(4, "delta = " << e * delta << " (" << delta << ")\n" << endi);
00551     origin += e * delta;
00552     
00553     // Check for grid overlap
00554     if (!fits_lattice(simParams->lattice)) {
00555       char errmsg[512];
00556       if (checksize) {
00557         sprintf(errmsg, "Warning: Periodic cell basis too small for Gridforce grid %d.  Set gridforcechecksize off in configuration file to ignore.\n", mygridnum);
00558         NAMD_die(errmsg);      
00559       }
00560     }
00561     
00562     size = k[0] * k[1] * k[2];
00563     dk[0] = k[1] * k[2];
00564     dk[1] = k[2];
00565     dk[2] = 1;
00566     
00567     DebugM(3, "size = " << size << ", size_nopad = " << size_nopad << "\n" << endi);
00568     
00569     delete[] grid;
00570     grid = new float[size];
00571     
00572     n_sum[0] = n_sum[1] = n_sum[2] = 0;
00573     p_sum[0] = p_sum[1] = p_sum[2] = 0;
00574     for (int i0 = 0; i0 < k_nopad[0]; i0++) {
00575         for (int i1 = 0; i1 < k_nopad[1]; i1++) {
00576             for (int i2 = 0; i2 < k_nopad[2]; i2++) {
00577                 // Edges are special cases -- take force there to be
00578                 // zero for smooth transition across potential
00579                 // boundary
00580                 
00581                 long int ind_nopad = i0*dk_nopad[0] + i1*dk_nopad[1] + i2*dk_nopad[2];
00582                 int j0 = (cont[0]) ? i0 : i0 + border;
00583                 int j1 = (cont[1]) ? i1 : i1 + border;
00584                 int j2 = (cont[2]) ? i2 : i2 + border;
00585                 long int ind = j0*dk[0] + j1*dk[1] + j2*dk[2];
00586                 
00587                 if (i0 == 0)                    n_sum[0] += grid_nopad[ind_nopad];
00588                 else if (i0 == k_nopad[0]-1)    p_sum[0] += grid_nopad[ind_nopad];
00589                 if (i1 == 0)                    n_sum[1] += grid_nopad[ind_nopad];
00590                 else if (i1 == k_nopad[1]-1)    p_sum[1] += grid_nopad[ind_nopad];
00591                 if (i2 == 0)                    n_sum[2] += grid_nopad[ind_nopad];
00592                 else if (i2 == k_nopad[2]-1)    p_sum[2] += grid_nopad[ind_nopad];
00593                 
00594                 //grid[ind] = grid_nopad[ind_nopad];
00595                 set_grid(j0, j1, j2, grid_nopad[ind_nopad]);
00596             }
00597         }
00598     }
00599     
00600     const BigReal modThresh = 1.0;
00601     
00602     BigReal n_avg[3], p_avg[3];
00603     int i0;
00604     for (int i0 = 0; i0 < 3; i0++) {
00605         int i1 = (i0 + 1) % 3;
00606         int i2 = (i0 + 2) % 3;
00607         n_avg[i0] = n_sum[i0] / (k_nopad[i1] * k_nopad[i2]);
00608         p_avg[i0] = p_sum[i0] / (k_nopad[i1] * k_nopad[i2]);
00609         
00610         if (cont[i0] && fabs(offset[i0] - (p_avg[i0]-n_avg[i0])) > modThresh) 
00611         {
00612             iout << iWARN << "GRID FORCE POTENTIAL DIFFERENCE IN K" << i0
00613                  << " DIRECTION IS " 
00614                  << offset[i0] - (p_avg[i0]-n_avg[i0]) 
00615                  << " KCAL/MOL*E\n" << endi;
00616         }
00617     }
00618     
00619     Bool twoPadVals = (cont[0] + cont[1] + cont[2] == 2);
00620     double padVal = 0.0;
00621     long int weight = 0;
00622     if (!twoPadVals) {
00623         // Determine pad value (must average)
00624         if (!cont[0]) {
00625             padVal += p_sum[0] + n_sum[0];
00626             weight += 2 * k_nopad[1] * k_nopad[2];
00627         }
00628         if (!cont[1]) {
00629             padVal += p_sum[1] + n_sum[1];
00630             weight += 2 * k_nopad[0] * k_nopad[2];
00631         }
00632         if (!cont[2]) {
00633             padVal += p_sum[2] + n_sum[2];
00634             weight += 2 * k_nopad[0] * k_nopad[1];
00635         }
00636         padVal /= weight;
00637     }
00638     
00639     for (int i = 0; i < 3; i++) {
00640         pad_n[i] = (cont[i]) ? 0.0 : (twoPadVals) ? n_avg[i] : padVal;
00641         pad_p[i] = (cont[i]) ? 0.0 : (twoPadVals) ? p_avg[i] : padVal;
00642         DebugM(4, "pad_n[" << i << "] = " << pad_n[i] << "\n");
00643         DebugM(4, "pad_p[" << i << "] = " << pad_p[i] << "\n" << endi);
00644     }
00645     
00646     if (cont[0] && cont[1] && cont[2]) {
00647         // Nothing to do
00648         return;
00649     }
00650     
00651     // Now fill in rest of new grid
00652     for (int i0 = 0; i0 < k[0]; i0++) {
00653         for (int i1 = 0; i1 < k[1]; i1++) {
00654             for (int i2 = 0; i2 < k[2]; i2++) {
00655                 if ( (cont[0] || (i0 >= border && i0 < k[0]-border)) 
00656                      && (cont[1] || (i1 >= border && i1 < k[1]-border)) 
00657                      && (cont[2] || i2 == border) )
00658                 {
00659                     i2 += k_nopad[2]-1;
00660                     continue;
00661                 }
00662                 
00663                 long int ind = i0*dk[0] + i1*dk[1] + i2*dk[2];
00664 
00665                 Position pos = e * Position(i0, i1, i2);
00666                 int var[3] = {i0, i1, i2};
00667                 
00668                 for (int dir = 0; dir < 3; dir++) {
00669                     if (cont[dir]) 
00670                         continue;
00671                     
00672                     if (var[dir] < border) {
00673                         //grid[ind] = pad_n[dir];
00674                         set_grid(i0, i1, i2, pad_n[dir]);
00675                     } else if (var[dir] >= k[dir]-border) {
00676                         //grid[ind] = pad_p[dir];
00677                         set_grid(i0, i1, i2, pad_p[dir]);
00678                     }
00679                 }
00680                 
00681 //              DebugM(2, "grid[" << ind << "; " << i0 << ", " << i1
00682 //                     << ", " << i2 << "] = " << get_grid(ind)
00683 //                     << "\n" << endi);
00684             }
00685         }
00686     }
00687     
00688     for (int i0 = 0; i0 < k[0]; i0++) {
00689         for (int i1 = 0; i1 < k[1]; i1++) {
00690             for (int i2 = 0; i2 < k[2]; i2++) {
00691                 DebugM(1, "grid[" << i0 << ", " << i1 << ", " << i2 << "] = " << get_grid(i0,i1,i2) << "\n" << endi);
00692             }
00693         }
00694     }
00695     
00696     // Clean up
00697     DebugM(3, "clean up\n" << endi);
00698     delete[] grid_nopad;
00699     
00700     // Call initialize for each subgrid
00701     for (int i = 0; i < numSubgrids; i++) {
00702         subgrids[i]->poten_fp = poten_fp;
00703         subgrids[i]->initialize(simParams, mgridParams);
00704     }
00705     
00706     // close file pointer
00707     fclose(poten_fp);
00708 }
00709 
00710 
00711 void GridforceFullMainGrid::reinitialize(SimParameters *simParams, MGridforceParams *mgridParams)
00712 {
00713     DebugM(4, "reinitializing grid\n" << endi);
00714     initialize(filename, simParams, mgridParams);
00715 }
00716 
00717 
00718 long int GridforceFullMainGrid::get_all_gridvals(float **all_gridvals) const
00719 {
00720     // Creates a flat array of all grid values, including subgrids,
00721     // and puts it in the value pointed to by the 'grids'
00722     // argument. Returns the resulting array size. Caller is
00723     // responsible for destroying the array via 'delete[]'
00724     
00725     DebugM(4, "get_all_gridvals called\n" << endi);
00726     
00727     long int sz = 0;
00728     sz += size;
00729     for (int i = 0; i < totalGrids-1; i++) {
00730         sz += subgrids_flat[i]->size;
00731     }
00732     DebugM(4, "size = " << sz << "\n" << endi);
00733     
00734     float *grid_vals = new float[sz];
00735     long int idx = 0;
00736     for (long int i = 0; i < size; i++) {
00737         grid_vals[idx++] = grid[i];
00738     }
00739     for (int j = 0; j < totalGrids-1; j++) {
00740         for (long int i = 0; i < subgrids_flat[j]->size; i++) {
00741             grid_vals[idx++] = subgrids_flat[j]->grid[i];
00742         }
00743     }
00744     CmiAssert(idx == sz);
00745     
00746     *all_gridvals = grid_vals;
00747     
00748     DebugM(4, "get_all_gridvals finished\n" << endi);
00749     
00750     return sz;
00751 }
00752 
00753 
00754 void GridforceFullMainGrid::set_all_gridvals(float *all_gridvals, long int sz)
00755 {
00756     DebugM(4, "set_all_gridvals called\n" << endi);
00757     
00758     long int sz_calc = 0;
00759     sz_calc += size;
00760     for (int i = 0; i < totalGrids-1; i++) {
00761         sz_calc += subgrids_flat[i]->size;
00762     }
00763     CmiAssert(sz == sz_calc);
00764     
00765     long int idx = 0;
00766     for (long int i = 0; i < size; i++) {
00767         DebugM(1, "all_gridvals[" << idx << "] = " << all_gridvals[idx] << "\n" << endi);
00768         grid[i] = all_gridvals[idx++];
00769     }
00770     for (int j = 0; j < totalGrids-1; j++) {
00771         for (long int i = 0; i < subgrids_flat[j]->size; i++) {
00772             DebugM(1, "all_gridvals[" << idx << "] = " << all_gridvals[idx] << "\n" << endi);
00773             subgrids_flat[j]->grid[i] = all_gridvals[idx++];
00774         }
00775     }
00776     CmiAssert(idx == sz);
00777 
00778     DebugM(4, "set_all_gridvals finished\n" << endi);
00779 }
00780 
00781 
00782 void GridforceFullMainGrid::compute_b(float *b, int *inds, Vector gapscale) const
00783 {
00784     for (int i0 = 0; i0 < 8; i0++) {
00785         int inds2[3];
00786         int zero_derivs = FALSE;
00787         
00788         float voff = 0.0;
00789         int bit = 1;    // bit = 2^i1 in the below loop
00790         for (int i1 = 0; i1 < 3; i1++) {
00791             inds2[i1] = (inds[i1] + ((i0 & bit) ? 1 : 0)) % k[i1];
00792             
00793             // Deal with voltage offsets
00794             if (cont[i1] && inds[i1] == (k[i1]-1) && inds2[i1] == 0) {
00795                 voff += offset[i1];
00796                 DebugM(3, "offset[" << i1 << "] = " << offset[i1] << "\n" << endi);
00797             }
00798             
00799             bit <<= 1;  // i.e. multiply by 2
00800         }
00801         
00802         DebugM(1, "inds2 = " << inds2[0] << " " << inds2[1] << " " << inds2[2] << "\n" << endi);
00803         
00804         // NOTE: leaving everything in terms of unit cell coordinates for now,
00805         // eventually will multiply by inv tensor when applying the force
00806         
00807         // First set variables 'dk_{hi,lo}' (glob notation). The 'hi'
00808         // ('lo') variable in a given dimension is the number added (subtracted)
00809         // to go up (down) one grid point in that dimension; both are normally
00810         // just the corresponding 'dk[i]'. However, if we are sitting on a
00811         // boundary and we are using a continuous grid, then we want to map the
00812         // next point off the grid back around to the other side. e.g. point
00813         // (k[0], i1, k) maps to point (0, i1, k), which would be
00814         // accomplished by changing 'dk1_hi' to -(k[0]-1)*dk1.
00815         
00816         int d_hi[3] = {1, 1, 1};
00817         int d_lo[3] = {1, 1, 1};
00818         float voffs[3];
00819         float dscales[3] = {0.5, 0.5, 0.5};
00820         for (int i1 = 0; i1 < 3; i1++) {
00821             if (inds2[i1] == 0) {
00822                 if (cont[i1]) {
00823                     d_lo[i1] = -(k[i1]-1);
00824                     voffs[i1] = offset[i1];
00825                     dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
00826                 }
00827                 else zero_derivs = TRUE;
00828             }
00829             else if (inds2[i1] == k[i1]-1) {
00830                 if (cont[i1]) {
00831                     d_hi[i1] = -(k[i1]-1);
00832                     voffs[i1] = offset[i1];
00833                     dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
00834                 }
00835                 else zero_derivs = TRUE;
00836             }
00837             else {
00838                 voffs[i1] = 0.0;
00839             }
00840         }
00841         
00842 //      DebugM(2, "cont = " << cont[0] << " " << cont[1] << " " << cont[2] << "\n" << endi);
00843 //      DebugM(2, "zero_derivs = " << zero_derivs << "\n" << endi);
00844 //      DebugM(2, "d_hi = " << d_hi[0] << " " << d_hi[1] << " " << d_hi[2] << "\n" << endi);
00845 //      DebugM(2, "d_lo = " << d_lo[0] << " " << d_lo[1] << " " << d_lo[2] << "\n" << endi);
00846         DebugM(1, "dscales = " << dscales[0] << " " << dscales[1] << " " << dscales[2] << "\n" << endi);
00847         DebugM(1, "voffs = " << voffs[0] << " " << voffs[1] << " " << voffs[2] << "\n" << endi);
00848         
00849         // V
00850         b[i0] = get_grid(inds2[0],inds2[1],inds2[2]) + voff;
00851         
00852         if (zero_derivs) {
00853             DebugM(2, "zero_derivs\n" << endi);
00854             b[8+i0] = 0.0;
00855             b[16+i0] = 0.0;
00856             b[24+i0] = 0.0;
00857             b[32+i0] = 0.0;
00858             b[40+i0] = 0.0;
00859             b[48+i0] = 0.0;
00860             b[56+i0] = 0.0;
00861         } else {
00862             b[8+i0]  = dscales[0] * (get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]) + voffs[0]);       //  dV/dx
00863             b[16+i0] = dscales[1] * (get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]) - get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]) + voffs[1]);       //  dV/dy
00864             b[24+i0] = dscales[2] * (get_grid_d(inds2[0],inds2[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0],inds2[1],inds2[2]-d_lo[2]) + voffs[2]);       //  dV/dz
00865             b[32+i0] = dscales[0] * dscales[1]
00866                 * (get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]) -
00867                    get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]));    //  d2V/dxdy
00868             b[40+i0] = dscales[0] * dscales[2]
00869                 * (get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]+d_hi[2]) -
00870                    get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]-d_lo[2]));    //  d2V/dxdz
00871             b[48+i0] = dscales[1] * dscales[2]
00872                 * (get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) -
00873                    get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]));    //  d2V/dydz
00874         
00875             b[56+i0] = dscales[0] * dscales[1] * dscales[2]                                     // d3V/dxdydz
00876                 * (get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) -
00877                    get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) +
00878                    get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) +
00879                    get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]));
00880         }
00881         
00882         DebugM(1, "V = " << b[i0] << "\n");
00883         
00884         DebugM(1, "dV/dx = " << b[8+i0] << "\n");
00885         DebugM(1, "dV/dy = " << b[16+i0] << "\n");
00886         DebugM(1, "dV/dz = " << b[24+i0] << "\n");
00887         
00888         DebugM(1, "d2V/dxdy = " << b[32+i0] << "\n");
00889         DebugM(1, "d2V/dxdz = " << b[40+i0] << "\n");
00890         DebugM(1, "d2V/dydz = " << b[48+i0] << "\n");
00891         
00892         DebugM(1, "d3V/dxdydz = " << b[56+i0] << "\n" << endi);
00893     }
00894 }
00895 
00896 
00897 /************************/
00898 /* GRIDFORCEFULLSUBGRID */
00899 /************************/
00900 
00901 GridforceFullSubGrid::GridforceFullSubGrid(GridforceFullBaseGrid *parent_in) {
00902     parent = parent_in;
00903     generation = parent->generation + 1;
00904     GridforceFullBaseGrid *tmp = parent;
00905     while (tmp->generation > 0) {
00906         tmp = ((GridforceFullSubGrid *)tmp)->parent;
00907     }
00908     maingrid = (GridforceFullMainGrid *)tmp;
00909     DebugM(4, "generation = " << generation << "\n" << endi);
00910 }
00911 
00912 
00913 void GridforceFullSubGrid::initialize(SimParameters *simParams, MGridforceParams *mgridParams)
00914 {
00915     int tmp;
00916     char line[256];
00917     long int poten_offset;
00918     
00919     // Skip 'attribute's
00920     DebugM(3, "Skipping 'attribute' keywords...\n" << endi);
00921     char str[256];
00922     do {
00923         poten_offset = ftell(poten_fp);
00924         fscanf(poten_fp, "%s", str);
00925         fgets(line, 256, poten_fp);
00926         DebugM(4, "Read line " << str << " " << line << endi);
00927     } while (strcmp(str, "attribute") == 0);
00928     fseek(poten_fp, poten_offset, SEEK_SET);
00929     
00930     // Skip 'field' object
00931     DebugM(3, "Skipping 'field' object\n" << endi);
00932     fscanf(poten_fp, "object");
00933     int n;
00934     n = fscanf(poten_fp, "\"%[^\"]\" class field\n", str);
00935     if (n == 0) {
00936         n = fscanf(poten_fp, "%d class field\n", &tmp);
00937     }
00938     
00939     if (n == 0) {
00940         NAMD_die("Error reading gridforce grid! Could not find field object!\n");
00941     }
00942     
00943     // Skip 'component's
00944     DebugM(3, "Skipping 'component' keywords\n" << endi);
00945     do {
00946         poten_offset = ftell(poten_fp);
00947         fscanf(poten_fp, "%s", str);
00948         fgets(line, 256, poten_fp);
00949     } while (strcmp(str, "component") == 0);
00950     fseek(poten_fp, poten_offset, SEEK_SET);
00951     
00952     // Read header
00953     readHeader(simParams, mgridParams);
00954     
00955     factor = 1.0;
00956     if (mgridParams->gridforceVolts)
00957     {
00958         factor /= 0.0434;  // convert V -> kcal/mol*e
00959     }
00960     scale = mgridParams->gridforceScale;
00961     checksize = mgridParams->gridforceCheckSize;
00962     
00963     for (int i = 0; i < 3; i++) {
00964         k[i] = k_nopad[i];      // subgrids aren't padded
00965     }
00966     
00967     // Make sure that each subgrid dimension is an integral
00968     // number of spanned supergrid cells. This is to ensure that no
00969     // supergrid nodes land in the middle of a subgrid, because in
00970     // this case forces could not be matched properly.
00971     for (int i = 0; i < 3; i++) {
00972         if ((k[i] - 1) % (pmax[i] - pmin[i] + 1) != 0) {
00973             iout << (k[i] - 1) << " % " << (pmax[i] - pmin[i] + 1) << " != 0\n" << endi;
00974             NAMD_die("Error reading gridforce grid! Subgrid dimensions must be an integral number spanned parent cells!");
00975         }
00976     }
00977     
00978     for (int i = 0; i < 3; i++) {
00979         if (parent->cont[i]) {
00980             cont[i] = (pmin[i] == 0 && pmax[i] == parent->k[i]-2) ? TRUE : FALSE;
00981             DebugM(3, "pmin[" << i << "] = " << pmin[i] << " pmax[" << i << "] = " << pmax[i] << " parent->k[" << i << "] = " << parent->k[i] << " cont[" << i << "] = " << cont[i] << "\n" << endi);
00982         } else {
00983             cont[i] = false;
00984             if (parent->generation == 0) {
00985                 // Add to pmin, pmax since parent got extra gridpoint layer(s) (maybe)
00986                 int brd = parent->get_border();
00987                 pmin[i] += brd;
00988                 pmax[i] += brd;
00989             }
00990         }               
00991     }
00992     
00993     DebugM(4, "pmin = " << pmin[0] << " " << pmin[1] << " " << pmin[2] << "\n");
00994     DebugM(4, "pmax = " << pmax[0] << " " << pmax[1] << " " << pmax[2] << "\n" << endi);
00995     
00996     Vector origin2 = parent->origin + parent->e * Position(pmin[0], pmin[1], pmin[2]);
00997     Vector escale, invscale;
00998     for (int i = 0; i < 3; i++) {
00999         escale[i] = double(pmax[i] - pmin[i] + 1)/(k[i]-1);
01000         invscale[i] = 1.0/escale[i];
01001         if (cont[i]) { pmax[i]++; }
01002     }
01003     Tensor e2 = tensorMult(parent->e, Tensor::diagonal(escale));
01004     
01005     // Check that lattice parameters agree with min and max numbers
01006     // from subgrid hierarchy.
01007     double TOL2 = 1e-4; // Totally arbitrary
01008     if (pow(origin2.x-origin.x, 2) > TOL2 ||
01009         pow(origin2.y-origin.y, 2) > TOL2 ||
01010         pow(origin2.z-origin.z, 2) > TOL2 ||
01011         pow(e2.xx-e.xx, 2) > TOL2 ||
01012         pow(e2.xy-e.xy, 2) > TOL2 ||
01013         pow(e2.xz-e.xz, 2) > TOL2 ||
01014         pow(e2.yx-e.yx, 2) > TOL2 ||
01015         pow(e2.yy-e.yy, 2) > TOL2 ||
01016         pow(e2.yz-e.yz, 2) > TOL2 ||
01017         pow(e2.zx-e.zx, 2) > TOL2 ||
01018         pow(e2.zy-e.zy, 2) > TOL2 ||
01019         pow(e2.zz-e.zz, 2) > TOL2)
01020     {
01021         NAMD_die("Error reading gridforce grid! Subgrid lattice does not match!");
01022     }
01023     
01024     // Overwrite what was read from the header
01025     origin = origin2;
01026     e = e2;
01027     
01028     inv = tensorMult(Tensor::diagonal(invscale), parent->inv);
01029     for (int i = 0; i < 3; i++) {
01030         gap[i] = escale[i] * parent->gap[i];
01031         gapinv[i] = invscale[i] * parent->gapinv[i];
01032         offset[i] = parent->offset[i];
01033     }
01034     center = origin + e * 0.5 * Position(k[0], k[1], k[2]);
01035     
01036     DebugM(4, "origin = " << origin << "\n");
01037     DebugM(4, "e = " << e << "\n");
01038     DebugM(4, "inv = " << inv << "\n");
01039     DebugM(4, "gap = " << gap[0] << " " << gap[2] << " " << gap[2] << " " << "\n");
01040     DebugM(4, "gapinv = " << gapinv[0] << " " << gapinv[2] << " " << gapinv[2] << " " << "\n");
01041     DebugM(4, "numSubgrids = " << numSubgrids << "\n");
01042     DebugM(4, "k = " << k[0] << " " << k[1] << " " << k[2] << "\n");
01043     DebugM(4, "escale = " << escale << "\n");
01044     DebugM(4, "invscale = " << invscale << "\n" << endi);
01045     
01046     /*** Set members ***/
01047     size = k[0] * k[1] * k[2];
01048     dk[0] = k[1] * k[2];
01049     dk[1] = k[2];
01050     dk[2] = 1;
01051     
01052     scale_dV = Tensor::diagonal(escale);
01053     scale_d2V = Tensor::diagonal(Vector(escale.x*escale.y, escale.x*escale.z, escale.y*escale.z));
01054     scale_d3V = escale.x * escale.y * escale.z;
01055     
01056     DebugM(4, "scale_dV = " << scale_dV << "\n");
01057     DebugM(4, "scale_d2V = " << scale_d2V << "\n");
01058     DebugM(4, "scale_d3V = " << scale_d3V << "\n" << endi);
01059     
01060     // Allocate storage for potential and read it
01061     float *grid_tmp = new float[size];
01062     
01063     float tmp2;
01064     DebugM(3, "size_nopad = " << size_nopad << "\n");
01065     for (long int count = 0; count < size_nopad; count++) {
01066 //      poten_offset = ftell(poten_fp);
01067 //      fscanf(poten_fp, "%s", str);
01068 //      fgets(line, 256, poten_fp);
01069 //      DebugM(4, "Read line " << str << " " << line << endi);
01070 //      fseek(poten_fp, poten_offset, SEEK_SET);
01071         
01072         int err = fscanf(poten_fp, "%f", &tmp2);
01073         if (err == EOF || err == 0) {
01074             NAMD_die("Grid force potential file incorrectly formatted");
01075         }
01076         grid_tmp[count] = tmp2 * factor;
01077     }
01078     fscanf(poten_fp, "\n");
01079     
01080     // Set real grid
01081     DebugM(3, "allocating grid\n" << endi);
01082     delete[] grid;
01083     grid = new float[size];
01084     for (int i0 = 0; i0 < k_nopad[0]; i0++) {
01085         for (int i1 = 0; i1 < k_nopad[1]; i1++) {
01086             for (int i2 = 0; i2 < k_nopad[2]; i2++) {
01087                 long int ind = i0*dk[0] + i1*dk[1] + i2*dk[2];
01088                 set_grid(i0, i1, i2, grid_tmp[ind]);
01089             }
01090         }
01091     }
01092     
01093     for (int i0 = 0; i0 < k[0]; i0++) {
01094         for (int i1 = 0; i1 < k[1]; i1++) {
01095             for (int i2 = 0; i2 < k[2]; i2++) {
01096                 DebugM(1, "grid[" << i0 << ", " << i1 << ", " << i2 << "] = " << get_grid(i0,i1,i2) << "\n" << endi);
01097             }
01098         }
01099     }
01100     
01101     // Clean up
01102     delete[] grid_tmp;
01103     
01104     // Call initialize for each subgrid
01105     for (int i = 0; i < numSubgrids; i++) {
01106         subgrids[i]->initialize(simParams, mgridParams);
01107     }
01108 }
01109 
01110 
01111 void GridforceFullSubGrid::pack(MOStream *msg) const
01112 {
01113     DebugM(4, "Packing subgrid\n" << endi);
01114     
01115     msg->put(sizeof(Tensor), (char*)&scale_dV);
01116     msg->put(sizeof(Tensor), (char*)&scale_d2V);
01117     msg->put(sizeof(float), (char*)&scale_d3V);
01118     
01119     msg->put(3*sizeof(int), (char*)pmin);
01120     msg->put(3*sizeof(int), (char*)pmax);
01121     msg->put(subgridIdx);
01122     
01123     DebugM(3, "calling GridforceFullBaseGrid::pack\n" << endi);
01124     
01125     GridforceFullBaseGrid::pack(msg);
01126 }
01127 
01128 
01129 void GridforceFullSubGrid::unpack(MIStream *msg)
01130 {
01131     DebugM(4, "Unpacking subgrid\n" << endi);
01132     
01133     msg->get(sizeof(Tensor), (char*)&scale_dV);
01134     msg->get(sizeof(Tensor), (char*)&scale_d2V);
01135     msg->get(sizeof(float), (char*)&scale_d3V);
01136     
01137     msg->get(3*sizeof(int), (char*)pmin);
01138     msg->get(3*sizeof(int), (char*)pmax);
01139     msg->get(subgridIdx);
01140     
01141     GridforceFullBaseGrid::unpack(msg);
01142     
01143     DebugM(4, "size  = " << size << "\n");
01144     DebugM(4, "numSubgrids = " << numSubgrids << "\n");
01145     DebugM(4, "gapinv = " << gapinv[0] << " " << gapinv[2] << " " << gapinv[2] << " " << "\n");
01146     DebugM(4, "generation = " << generation << "\n" << endi);
01147 }
01148 
01149 
01150 void GridforceFullSubGrid::addToSubgridsFlat(void)
01151 {
01152     DebugM(4, "addToSubgridsFlat() called, subgridIdx = " << subgridIdx << ", maingrid->numSubgrids = " << maingrid->numSubgrids << "\n" << endi);
01153     maingrid->subgrids_flat[subgridIdx-1] = this;
01154     for (int i = 0; i < numSubgrids; i++) {
01155         subgrids[i]->addToSubgridsFlat();
01156     }
01157 }
01158 
01159 
01160 void GridforceFullSubGrid::compute_b(float *b, int *inds, Vector gapscale) const
01161 {
01162     for (int i0 = 0; i0 < 8; i0++) {
01163         int inds2[3];
01164         
01165         float voff = 0.0;
01166         int bit = 1;    // bit = 2^i1 in the below loop
01167         for (int i1 = 0; i1 < 3; i1++) {
01168             inds2[i1] = (inds[i1] + ((i0 & bit) ? 1 : 0)) % k[i1];
01169             
01170             // Deal with voltage offsets
01171             if (cont[i1] && inds[i1] == (k[i1]-1) && inds2[i1] == 0) {
01172                 voff += offset[i1];
01173                 DebugM(3, "offset[" << i1 << "] = " << offset[i1] << "\n" << endi);
01174             }
01175             
01176             bit <<= 1;  // i.e. multiply by 2
01177         }
01178         
01179         DebugM(3, "inds2 = " << inds2[0] << " " << inds2[1] << " " << inds2[2] << "\n" << endi);
01180         
01181         int d_hi[3] = {1, 1, 1}; 
01182         int d_lo[3] = {1, 1, 1};
01183         float voffs[3];
01184         float dscales[3] = {0.5, 0.5, 0.5};
01185         for (int i1 = 0; i1 < 3; i1++) {
01186             if (inds2[i1] == 0 && cont[i1]) {
01187                 d_lo[i1] = -(k[i1]-1);
01188                 voffs[i1] = offset[i1];
01189                 dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
01190             }
01191             else if (inds2[i1] == k[i1]-1 && cont[i1]) {
01192                 d_hi[i1] = -(k[i1]-1);
01193                 voffs[i1] = offset[i1];
01194                 dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
01195             }
01196             else {
01197                 voffs[i1] = 0.0;
01198             }
01199         }
01200         
01201         bool edge = false;
01202         for (int i1 = 0; i1 < 3; i1++) {
01203             if (!cont[i1] && (inds2[i1] == 0 || inds2[i1] == k[i1]-1)) {
01204                 edge = true;
01205             }
01206         }
01207         
01208         if (inds2[2] == 0) {
01209 //          DebugM(3, "cont = " << cont[0] << " " << cont[1] << " " << cont[2] << " d_hi = " << d_hi[0] << " " << d_hi[1] << " " << d_hi[2] << " d_lo = " << d_lo[0] << " " << d_lo[1] << " " << d_lo[2] << " dscales = " << dscales[0] << " " << dscales[1] << " " << dscales[2] << "\n" << endi);
01210             DebugM(3, "cont = " << cont[0] << " " << cont[1] << " " << cont[2] << "\n" << endi);
01211         }
01212         
01213         if (edge) {
01214             DebugM(2, "Edge!\n" << endi);
01215             
01216             // Must get derivatives from parent
01217             Position pos = e * Vector(inds2[0], inds2[1], inds2[2]) + origin;   // Gridpoint position in realspace
01218             Vector g = parent->inv * (pos - parent->origin);    // Gridpoint position in parent's gridspace
01219             Vector dg;
01220             int inds3[3];
01221             
01222             DebugM(2, "g = " << g << "\n" << endi);
01223             
01224             for (int i = 0; i < 3; i++) {
01225                 inds3[i] = (int)floor(g[i]);
01226                 dg[i] = g[i] - inds3[i];
01227             }
01228     
01229             float x[4], y[4], z[4];
01230             x[0] = 1; y[0] = 1; z[0] = 1;
01231             for (int j = 1; j < 4; j++) {
01232                 x[j] = x[j-1] * dg.x;
01233                 y[j] = y[j-1] * dg.y;
01234                 z[j] = z[j-1] * dg.z;
01235                 DebugM(1, "x[" << j << "] = " << x[j] << "\n");
01236                 DebugM(1, "y[" << j << "] = " << y[j] << "\n");
01237                 DebugM(1, "z[" << j << "] = " << z[j] << "\n" << endi);
01238             }
01239             
01240             // Compute parent matrices
01241             float b_parent[64];
01242             parent->compute_b(b_parent, inds3, gapscale);
01243             
01244             float a_parent[64];
01245             parent->compute_a(a_parent, b_parent);
01246             
01247             // Compute parent derivatives
01248             float V = parent->compute_V(a_parent, x, y, z);
01249             Vector dV = scale_dV * parent->compute_dV(a_parent, x, y, z);
01250             Vector d2V = scale_d2V * parent->compute_d2V(a_parent, x, y, z);
01251             float d3V = scale_d3V * parent->compute_d3V(a_parent, x, y, z);
01252             
01253             b[i0] = V;
01254             b[8+i0] = dV[0];
01255             b[16+i0] = dV[1];
01256             b[24+i0] = dV[2];
01257             b[32+i0] = d2V[0];
01258             b[40+i0] = d2V[1];
01259             b[48+i0] = d2V[2];
01260             b[56+i0] = d3V;
01261         } else {
01262             b[i0] = get_grid(inds2[0],inds2[1],inds2[2]) + voff;        // V
01263             
01264             b[8+i0]  = dscales[0] * (get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]) + voffs[0]);       //  dV/dx
01265             b[16+i0] = dscales[1] * (get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]) - get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]) + voffs[1]);       //  dV/dy
01266             b[24+i0] = dscales[2] * (get_grid_d(inds2[0],inds2[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0],inds2[1],inds2[2]-d_lo[2]) + voffs[2]);       //  dV/dz
01267             b[32+i0] = dscales[0] * dscales[1]
01268                 * (get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]) -
01269                    get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]));    //  d2V/dxdy
01270             b[40+i0] = dscales[0] * dscales[2]
01271                 * (get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]+d_hi[2]) -
01272                    get_grid_d(inds2[0]+d_hi[0],inds2[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1],inds2[2]-d_lo[2]));    //  d2V/dxdz
01273             b[48+i0] = dscales[1] * dscales[2]
01274                 * (get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) -
01275                    get_grid_d(inds2[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]));    //  d2V/dydz
01276         
01277             b[56+i0] = dscales[0] * dscales[1] * dscales[2]                                     // d3V/dxdydz
01278                 * (get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]+d_hi[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) -
01279                    get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]+d_hi[2]) +
01280                    get_grid_d(inds2[0]+d_hi[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]) + get_grid_d(inds2[0]-d_lo[0],inds2[1]+d_hi[1],inds2[2]-d_lo[2]) +
01281                    get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]+d_hi[2]) - get_grid_d(inds2[0]-d_lo[0],inds2[1]-d_lo[1],inds2[2]-d_lo[2]));
01282         }
01283         
01284         if (inds2[0] == 1 && inds2[1] == 1 && inds2[2] == 0) {
01285             DebugM(1, "Sub V = " << b[i0] << "\n");
01286         
01287             DebugM(1, "Sub dV/dx = " << b[8+i0] << "\n");
01288             DebugM(1, "Sub dV/dy = " << b[16+i0] << "\n");
01289             DebugM(1, "Sub dV/dz = " << b[24+i0] << "\n");
01290         
01291             DebugM(1, "Sub d2V/dxdy = " << b[32+i0] << "\n");
01292             DebugM(1, "Sub d2V/dxdz = " << b[40+i0] << "\n");
01293             DebugM(1, "Sub d2V/dydz = " << b[48+i0] << "\n");
01294         
01295             DebugM(1, "Sub d3V/dxdydz = " << b[56+i0] << "\n" << endi);
01296         }
01297     }
01298 }
01299 
01300 
01301 /*********************/
01302 /* GRIDFORCELITEGRID */
01303 /*********************/
01304 
01305 GridforceLiteGrid::GridforceLiteGrid(int gridnum)
01306 {
01307     mygridnum = gridnum;
01308     grid = NULL;
01309     type = GridforceGridTypeLite;
01310 }
01311 
01312 
01313 GridforceLiteGrid::~GridforceLiteGrid()
01314 {
01315     delete[] grid;
01316 }
01317 
01318 
01319 void GridforceLiteGrid::initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
01320 {
01321     // cheat and use GridforceFullMainGrid to read the file
01322     GridforceFullMainGrid *tmp_grid = new GridforceFullMainGrid(mygridnum);
01323     tmp_grid->initialize(potfilename, simParams, mgridParams, 1);
01324     
01325     if (tmp_grid->get_total_grids() != 1) {
01326         NAMD_die("Cannot use gridforcelite option with multi-resolution grid!");
01327     }
01328     
01329     // save file name so that grid can be re-read via Tcl
01330     strcpy(filename, potfilename);
01331     
01332     // copy parameters
01333     k[0] = tmp_grid->get_k0();
01334     k[1] = tmp_grid->get_k1();
01335     k[2] = tmp_grid->get_k2();
01336     k[3] = 4;   // for V, dV/dx, dV/dy, dV/dz grids
01337     origin = tmp_grid->get_origin();
01338     center = tmp_grid->get_center();
01339     e = tmp_grid->get_e();
01340     inv = tmp_grid->get_inv();
01341     scale = tmp_grid->get_scale();
01342     
01343     // calculate rest of parameters
01344     size = k[0] * k[1] * k[2] * k[3];
01345     dk[0] = k[1] * k[2] * k[3];
01346     dk[1] = k[2] * k[3];
01347     dk[2] = k[3];
01348     dk[3] = 1;
01349     
01350     // copy the potential grid
01351     delete[] grid;
01352     grid = new float[size];
01353     for (int i0 = 0; i0 < k[0]; i0++) {
01354         for (int i1 = 0; i1 < k[1]; i1++) {
01355             for (int i2 = 0; i2 < k[2]; i2++) {
01356                 float V = tmp_grid->get_grid(i0, i1, i2);
01357                 set_grid(i0, i1, i2, 0, V);
01358                 DebugM(1, "V[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 0) << "(" << V << ")\n" << endi);
01359             }
01360         }
01361     }
01362     
01363     delete tmp_grid;
01364     
01365     compute_derivative_grids();
01366 }
01367 
01368 
01369 void GridforceLiteGrid::compute_derivative_grids(void)
01370 {
01371     // calculate derivative grids
01372     // separate loop so all V values have been set already
01373     for (int i0 = 0; i0 < k[0]; i0++) {
01374         for (int i1 = 0; i1 < k[1]; i1++) {
01375             for (int i2 = 0; i2 < k[2]; i2++) {
01376                 float dx, dy, dz;
01377                 if (i0 == 0 || i0 == k[0]-1 || i1 == 0 || i1 == k[1]-1 || i2 == 0 || i2 == k[2]-1) {
01378                     // on edge, set ALL derivatives to zero (make up for lack of padding)
01379                     dx = 0;
01380                     dy = 0;
01381                     dz = 0;
01382                 } else {
01383                     dx = 0.5 * (get_grid_d(i0+1,i1,i2,0) - get_grid_d(i0-1,i1,i2,0));
01384                     dy = 0.5 * (get_grid_d(i0,i1+1,i2,0) - get_grid_d(i0,i1-1,i2,0));
01385                     dz = 0.5 * (get_grid_d(i0,i1,i2+1,0) - get_grid_d(i0,i1,i2-1,0));
01386                 }
01387                 set_grid(i0, i1, i2, 1, dx);
01388                 set_grid(i0, i1, i2, 2, dy);
01389                 set_grid(i0, i1, i2, 3, dz);
01390                 DebugM(1, "dx[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 1) << "(" << dx << ")\n" << endi);
01391                 DebugM(1, "dy[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 2) << "(" << dy << ")\n" << endi);
01392                 DebugM(1, "dz[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 3) << "(" << dz << ")\n" << endi);
01393             }
01394         }
01395     }
01396 }
01397 
01398 
01399 void GridforceLiteGrid::reinitialize(SimParameters *simParams, MGridforceParams *mgridParams)
01400 {
01401     initialize(filename, simParams, mgridParams);
01402 }
01403 
01404 
01405 void GridforceLiteGrid::pack(MOStream *msg) const
01406 {
01407     msg->put(4*sizeof(int), (char*)k);
01408     msg->put(size);
01409     msg->put(4*sizeof(long int), (char*)dk);
01410     
01411     msg->put(sizeof(Vector), (char*)&origin);
01412     msg->put(sizeof(Vector), (char*)&center);
01413     msg->put(sizeof(Tensor), (char*)&e);
01414     msg->put(sizeof(Tensor), (char*)&inv);
01415     msg->put(sizeof(Vector), (char*)&scale);
01416     msg->put(sizeof(Bool), (char*)&checksize);
01417     
01418     msg->put(129*sizeof(char), (char*)filename);
01419     
01420     msg->put(size*sizeof(float), (char*)grid);
01421 }
01422 
01423 
01424 void GridforceLiteGrid::unpack(MIStream *msg)
01425 {
01426     delete[] grid;
01427     grid = NULL;
01428 
01429     msg->get(4*sizeof(int), (char*)k);
01430     msg->get(size);
01431     msg->get(4*sizeof(long int), (char*)dk);
01432     
01433     msg->get(sizeof(Vector), (char*)&origin);
01434     msg->get(sizeof(Vector), (char*)&center);
01435     msg->get(sizeof(Tensor), (char*)&e);
01436     msg->get(sizeof(Tensor), (char*)&inv);
01437     msg->get(sizeof(Vector), (char*)&scale);
01438     msg->get(sizeof(Bool), (char*)&checksize);
01439     
01440     msg->get(129*sizeof(char), (char*)filename);
01441     
01442     if (size) {
01443         grid = new float[size];
01444         msg->get(size*sizeof(float), (char*)grid);
01445     }
01446 }
01447 
01448 
01449 long int GridforceLiteGrid::get_all_gridvals(float** all_gridvals) const
01450 {
01451     // Creates a flat array of all grid values and puts it in the
01452     // value pointed to by the 'all_gridvals' argument. Returns the
01453     // resulting array size. Caller is responsible for destroying the
01454     // array via 'delete[]'
01455     
01456     DebugM(4, "GridforceLiteGrid::get_all_gridvals called\n" << endi);
01457     
01458     long int sz = size;
01459     DebugM(4, "size = " << sz << "\n" << endi);
01460     
01461     float *grid_vals = new float[sz];
01462     long int idx = 0;
01463     for (long int i = 0; i < size; i++) {
01464         grid_vals[idx++] = grid[i];
01465     }
01466     CmiAssert(idx == sz);
01467     
01468     *all_gridvals = grid_vals;
01469     
01470     DebugM(4, "GridforceLiteGrid::get_all_gridvals finished\n" << endi);
01471     
01472     return sz;
01473 }
01474 
01475 
01476 void GridforceLiteGrid::set_all_gridvals(float* all_gridvals, long int sz)
01477 {
01478     DebugM(4, "GridforceLiteGrid::set_all_gridvals called\n" << endi);
01479     
01480     long int sz_calc = size;
01481     CmiAssert(sz == sz_calc);
01482     
01483     long int idx = 0;
01484     for (long int i = 0; i < size; i++) {
01485         grid[i] = all_gridvals[idx++];
01486     }
01487     CmiAssert(idx == sz);
01488     
01489     //compute_derivative_grids();       // not needed if we're sending all 4 grids
01490     
01491     DebugM(4, "GridforceLiteGrid::set_all_gridvals finished\n" << endi);
01492 }

Generated on Tue Sep 19 01:17:12 2017 for NAMD by  doxygen 1.4.7