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     poten_fp = parent->poten_fp;
00904     generation = parent->generation + 1;
00905     GridforceFullBaseGrid *tmp = parent;
00906     while (tmp->generation > 0) {
00907         tmp = ((GridforceFullSubGrid *)tmp)->parent;
00908     }
00909     maingrid = (GridforceFullMainGrid *)tmp;
00910     DebugM(4, "generation = " << generation << "\n" << endi);
00911 }
00912 
00913 
00914 void GridforceFullSubGrid::initialize(SimParameters *simParams, MGridforceParams *mgridParams)
00915 {
00916     int tmp;
00917     char line[256];
00918     long int poten_offset;
00919     
00920     // Skip 'attribute's
00921     DebugM(3, "Skipping 'attribute' keywords...\n" << endi);
00922     char str[256];
00923     do {
00924         poten_offset = ftell(poten_fp);
00925         fscanf(poten_fp, "%s", str);
00926         fgets(line, 256, poten_fp);
00927         DebugM(4, "Read line " << str << " " << line << endi);
00928     } while (strcmp(str, "attribute") == 0);
00929     fseek(poten_fp, poten_offset, SEEK_SET);
00930     
00931     // Skip 'field' object
00932     DebugM(3, "Skipping 'field' object\n" << endi);
00933     fscanf(poten_fp, "object");
00934     int n;
00935     n = fscanf(poten_fp, "\"%[^\"]\" class field\n", str);
00936     if (n == 0) {
00937         n = fscanf(poten_fp, "%d class field\n", &tmp);
00938     }
00939     
00940     if (n == 0) {
00941         NAMD_die("Error reading gridforce grid! Could not find field object!\n");
00942     }
00943     
00944     // Skip 'component's
00945     DebugM(3, "Skipping 'component' keywords\n" << endi);
00946     do {
00947         poten_offset = ftell(poten_fp);
00948         fscanf(poten_fp, "%s", str);
00949         fgets(line, 256, poten_fp);
00950     } while (strcmp(str, "component") == 0);
00951     fseek(poten_fp, poten_offset, SEEK_SET);
00952     
00953     // Read header
00954     readHeader(simParams, mgridParams);
00955     
00956     factor = 1.0;
00957     if (mgridParams->gridforceVolts)
00958     {
00959         factor /= 0.0434;  // convert V -> kcal/mol*e
00960     }
00961     scale = mgridParams->gridforceScale;
00962     checksize = mgridParams->gridforceCheckSize;
00963     
00964     for (int i = 0; i < 3; i++) {
00965         k[i] = k_nopad[i];      // subgrids aren't padded
00966     }
00967     
00968     // Make sure that each subgrid dimension is an integral
00969     // number of spanned supergrid cells. This is to ensure that no
00970     // supergrid nodes land in the middle of a subgrid, because in
00971     // this case forces could not be matched properly.
00972     for (int i = 0; i < 3; i++) {
00973         if ((k[i] - 1) % (pmax[i] - pmin[i] + 1) != 0) {
00974             iout << (k[i] - 1) << " % " << (pmax[i] - pmin[i] + 1) << " != 0\n" << endi;
00975             NAMD_die("Error reading gridforce grid! Subgrid dimensions must be an integral number spanned parent cells!");
00976         }
00977     }
00978     
00979     for (int i = 0; i < 3; i++) {
00980         if (parent->cont[i]) {
00981             cont[i] = (pmin[i] == 0 && pmax[i] == parent->k[i]-2) ? TRUE : FALSE;
00982             DebugM(3, "pmin[" << i << "] = " << pmin[i] << " pmax[" << i << "] = " << pmax[i] << " parent->k[" << i << "] = " << parent->k[i] << " cont[" << i << "] = " << cont[i] << "\n" << endi);
00983         } else {
00984             cont[i] = false;
00985             if (parent->generation == 0) {
00986                 // Add to pmin, pmax since parent got extra gridpoint layer(s) (maybe)
00987                 int brd = parent->get_border();
00988                 pmin[i] += brd;
00989                 pmax[i] += brd;
00990             }
00991         }               
00992     }
00993     
00994     DebugM(4, "pmin = " << pmin[0] << " " << pmin[1] << " " << pmin[2] << "\n");
00995     DebugM(4, "pmax = " << pmax[0] << " " << pmax[1] << " " << pmax[2] << "\n" << endi);
00996     
00997     Vector origin2 = parent->origin + parent->e * Position(pmin[0], pmin[1], pmin[2]);
00998     Vector escale, invscale;
00999     for (int i = 0; i < 3; i++) {
01000         escale[i] = double(pmax[i] - pmin[i] + 1)/(k[i]-1);
01001         invscale[i] = 1.0/escale[i];
01002         if (cont[i]) { pmax[i]++; }
01003     }
01004     Tensor e2 = tensorMult(parent->e, Tensor::diagonal(escale));
01005     
01006     // Check that lattice parameters agree with min and max numbers
01007     // from subgrid hierarchy.
01008     double TOL2 = 1e-4; // Totally arbitrary
01009     if (pow(origin2.x-origin.x, 2) > TOL2 ||
01010         pow(origin2.y-origin.y, 2) > TOL2 ||
01011         pow(origin2.z-origin.z, 2) > TOL2 ||
01012         pow(e2.xx-e.xx, 2) > TOL2 ||
01013         pow(e2.xy-e.xy, 2) > TOL2 ||
01014         pow(e2.xz-e.xz, 2) > TOL2 ||
01015         pow(e2.yx-e.yx, 2) > TOL2 ||
01016         pow(e2.yy-e.yy, 2) > TOL2 ||
01017         pow(e2.yz-e.yz, 2) > TOL2 ||
01018         pow(e2.zx-e.zx, 2) > TOL2 ||
01019         pow(e2.zy-e.zy, 2) > TOL2 ||
01020         pow(e2.zz-e.zz, 2) > TOL2)
01021     {
01022         NAMD_die("Error reading gridforce grid! Subgrid lattice does not match!");
01023     }
01024     
01025     // Overwrite what was read from the header
01026     origin = origin2;
01027     e = e2;
01028     
01029     inv = tensorMult(Tensor::diagonal(invscale), parent->inv);
01030     for (int i = 0; i < 3; i++) {
01031         gap[i] = escale[i] * parent->gap[i];
01032         gapinv[i] = invscale[i] * parent->gapinv[i];
01033         offset[i] = parent->offset[i];
01034     }
01035     center = origin + e * 0.5 * Position(k[0], k[1], k[2]);
01036     
01037     DebugM(4, "origin = " << origin << "\n");
01038     DebugM(4, "e = " << e << "\n");
01039     DebugM(4, "inv = " << inv << "\n");
01040     DebugM(4, "gap = " << gap[0] << " " << gap[2] << " " << gap[2] << " " << "\n");
01041     DebugM(4, "gapinv = " << gapinv[0] << " " << gapinv[2] << " " << gapinv[2] << " " << "\n");
01042     DebugM(4, "numSubgrids = " << numSubgrids << "\n");
01043     DebugM(4, "k = " << k[0] << " " << k[1] << " " << k[2] << "\n");
01044     DebugM(4, "escale = " << escale << "\n");
01045     DebugM(4, "invscale = " << invscale << "\n" << endi);
01046     
01047     /*** Set members ***/
01048     size = k[0] * k[1] * k[2];
01049     dk[0] = k[1] * k[2];
01050     dk[1] = k[2];
01051     dk[2] = 1;
01052     
01053     scale_dV = Tensor::diagonal(escale);
01054     scale_d2V = Tensor::diagonal(Vector(escale.x*escale.y, escale.x*escale.z, escale.y*escale.z));
01055     scale_d3V = escale.x * escale.y * escale.z;
01056     
01057     DebugM(4, "scale_dV = " << scale_dV << "\n");
01058     DebugM(4, "scale_d2V = " << scale_d2V << "\n");
01059     DebugM(4, "scale_d3V = " << scale_d3V << "\n" << endi);
01060     
01061     // Allocate storage for potential and read it
01062     float *grid_tmp = new float[size];
01063     
01064     float tmp2;
01065     DebugM(3, "size_nopad = " << size_nopad << "\n");
01066     for (long int count = 0; count < size_nopad; count++) {
01067 //      poten_offset = ftell(poten_fp);
01068 //      fscanf(poten_fp, "%s", str);
01069 //      fgets(line, 256, poten_fp);
01070 //      DebugM(4, "Read line " << str << " " << line << endi);
01071 //      fseek(poten_fp, poten_offset, SEEK_SET);
01072         
01073         int err = fscanf(poten_fp, "%f", &tmp2);
01074         if (err == EOF || err == 0) {
01075             NAMD_die("Grid force potential file incorrectly formatted");
01076         }
01077         grid_tmp[count] = tmp2 * factor;
01078     }
01079     fscanf(poten_fp, "\n");
01080     
01081     // Set real grid
01082     DebugM(3, "allocating grid\n" << endi);
01083     delete[] grid;
01084     grid = new float[size];
01085     for (int i0 = 0; i0 < k_nopad[0]; i0++) {
01086         for (int i1 = 0; i1 < k_nopad[1]; i1++) {
01087             for (int i2 = 0; i2 < k_nopad[2]; i2++) {
01088                 long int ind = i0*dk[0] + i1*dk[1] + i2*dk[2];
01089                 set_grid(i0, i1, i2, grid_tmp[ind]);
01090             }
01091         }
01092     }
01093     
01094     for (int i0 = 0; i0 < k[0]; i0++) {
01095         for (int i1 = 0; i1 < k[1]; i1++) {
01096             for (int i2 = 0; i2 < k[2]; i2++) {
01097                 DebugM(1, "grid[" << i0 << ", " << i1 << ", " << i2 << "] = " << get_grid(i0,i1,i2) << "\n" << endi);
01098             }
01099         }
01100     }
01101     
01102     // Clean up
01103     delete[] grid_tmp;
01104     
01105     // Call initialize for each subgrid
01106     for (int i = 0; i < numSubgrids; i++) {
01107         subgrids[i]->initialize(simParams, mgridParams);
01108     }
01109 }
01110 
01111 
01112 void GridforceFullSubGrid::pack(MOStream *msg) const
01113 {
01114     DebugM(4, "Packing subgrid\n" << endi);
01115     
01116     msg->put(sizeof(Tensor), (char*)&scale_dV);
01117     msg->put(sizeof(Tensor), (char*)&scale_d2V);
01118     msg->put(sizeof(float), (char*)&scale_d3V);
01119     
01120     msg->put(3*sizeof(int), (char*)pmin);
01121     msg->put(3*sizeof(int), (char*)pmax);
01122     msg->put(subgridIdx);
01123     
01124     DebugM(3, "calling GridforceFullBaseGrid::pack\n" << endi);
01125     
01126     GridforceFullBaseGrid::pack(msg);
01127 }
01128 
01129 
01130 void GridforceFullSubGrid::unpack(MIStream *msg)
01131 {
01132     DebugM(4, "Unpacking subgrid\n" << endi);
01133     
01134     msg->get(sizeof(Tensor), (char*)&scale_dV);
01135     msg->get(sizeof(Tensor), (char*)&scale_d2V);
01136     msg->get(sizeof(float), (char*)&scale_d3V);
01137     
01138     msg->get(3*sizeof(int), (char*)pmin);
01139     msg->get(3*sizeof(int), (char*)pmax);
01140     msg->get(subgridIdx);
01141     
01142     GridforceFullBaseGrid::unpack(msg);
01143     
01144     DebugM(4, "size  = " << size << "\n");
01145     DebugM(4, "numSubgrids = " << numSubgrids << "\n");
01146     DebugM(4, "gapinv = " << gapinv[0] << " " << gapinv[2] << " " << gapinv[2] << " " << "\n");
01147     DebugM(4, "generation = " << generation << "\n" << endi);
01148 }
01149 
01150 
01151 void GridforceFullSubGrid::addToSubgridsFlat(void)
01152 {
01153     DebugM(4, "addToSubgridsFlat() called, subgridIdx = " << subgridIdx << ", maingrid->numSubgrids = " << maingrid->numSubgrids << "\n" << endi);
01154     maingrid->subgrids_flat[subgridIdx-1] = this;
01155     for (int i = 0; i < numSubgrids; i++) {
01156         subgrids[i]->addToSubgridsFlat();
01157     }
01158 }
01159 
01160 
01161 void GridforceFullSubGrid::compute_b(float *b, int *inds, Vector gapscale) const
01162 {
01163     for (int i0 = 0; i0 < 8; i0++) {
01164         int inds2[3];
01165         
01166         float voff = 0.0;
01167         int bit = 1;    // bit = 2^i1 in the below loop
01168         for (int i1 = 0; i1 < 3; i1++) {
01169             inds2[i1] = (inds[i1] + ((i0 & bit) ? 1 : 0)) % k[i1];
01170             
01171             // Deal with voltage offsets
01172             if (cont[i1] && inds[i1] == (k[i1]-1) && inds2[i1] == 0) {
01173                 voff += offset[i1];
01174                 DebugM(3, "offset[" << i1 << "] = " << offset[i1] << "\n" << endi);
01175             }
01176             
01177             bit <<= 1;  // i.e. multiply by 2
01178         }
01179         
01180         DebugM(3, "inds2 = " << inds2[0] << " " << inds2[1] << " " << inds2[2] << "\n" << endi);
01181         
01182         int d_hi[3] = {1, 1, 1}; 
01183         int d_lo[3] = {1, 1, 1};
01184         float voffs[3];
01185         float dscales[3] = {0.5, 0.5, 0.5};
01186         for (int i1 = 0; i1 < 3; i1++) {
01187             if (inds2[i1] == 0 && cont[i1]) {
01188                 d_lo[i1] = -(k[i1]-1);
01189                 voffs[i1] = offset[i1];
01190                 dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
01191             }
01192             else if (inds2[i1] == k[i1]-1 && cont[i1]) {
01193                 d_hi[i1] = -(k[i1]-1);
01194                 voffs[i1] = offset[i1];
01195                 dscales[i1] = 1.0/(1.0 + gap[i1]) * 1.0/gapscale[i1];
01196             }
01197             else {
01198                 voffs[i1] = 0.0;
01199             }
01200         }
01201         
01202         bool edge = false;
01203         for (int i1 = 0; i1 < 3; i1++) {
01204             if (!cont[i1] && (inds2[i1] == 0 || inds2[i1] == k[i1]-1)) {
01205                 edge = true;
01206             }
01207         }
01208         
01209         if (inds2[2] == 0) {
01210 //          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);
01211             DebugM(3, "cont = " << cont[0] << " " << cont[1] << " " << cont[2] << "\n" << endi);
01212         }
01213         
01214         if (edge) {
01215             DebugM(2, "Edge!\n" << endi);
01216             
01217             // Must get derivatives from parent
01218             Position pos = e * Vector(inds2[0], inds2[1], inds2[2]) + origin;   // Gridpoint position in realspace
01219             Vector g = parent->inv * (pos - parent->origin);    // Gridpoint position in parent's gridspace
01220             Vector dg;
01221             int inds3[3];
01222             
01223             DebugM(2, "g = " << g << "\n" << endi);
01224             
01225             for (int i = 0; i < 3; i++) {
01226                 inds3[i] = (int)floor(g[i]);
01227                 dg[i] = g[i] - inds3[i];
01228             }
01229     
01230             float x[4], y[4], z[4];
01231             x[0] = 1; y[0] = 1; z[0] = 1;
01232             for (int j = 1; j < 4; j++) {
01233                 x[j] = x[j-1] * dg.x;
01234                 y[j] = y[j-1] * dg.y;
01235                 z[j] = z[j-1] * dg.z;
01236                 DebugM(1, "x[" << j << "] = " << x[j] << "\n");
01237                 DebugM(1, "y[" << j << "] = " << y[j] << "\n");
01238                 DebugM(1, "z[" << j << "] = " << z[j] << "\n" << endi);
01239             }
01240             
01241             // Compute parent matrices
01242             float b_parent[64];
01243             parent->compute_b(b_parent, inds3, gapscale);
01244             
01245             float a_parent[64];
01246             parent->compute_a(a_parent, b_parent);
01247             
01248             // Compute parent derivatives
01249             float V = parent->compute_V(a_parent, x, y, z);
01250             Vector dV = scale_dV * parent->compute_dV(a_parent, x, y, z);
01251             Vector d2V = scale_d2V * parent->compute_d2V(a_parent, x, y, z);
01252             float d3V = scale_d3V * parent->compute_d3V(a_parent, x, y, z);
01253             
01254             b[i0] = V;
01255             b[8+i0] = dV[0];
01256             b[16+i0] = dV[1];
01257             b[24+i0] = dV[2];
01258             b[32+i0] = d2V[0];
01259             b[40+i0] = d2V[1];
01260             b[48+i0] = d2V[2];
01261             b[56+i0] = d3V;
01262         } else {
01263             b[i0] = get_grid(inds2[0],inds2[1],inds2[2]) + voff;        // V
01264             
01265             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
01266             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
01267             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
01268             b[32+i0] = dscales[0] * dscales[1]
01269                 * (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]) -
01270                    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
01271             b[40+i0] = dscales[0] * dscales[2]
01272                 * (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]) -
01273                    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
01274             b[48+i0] = dscales[1] * dscales[2]
01275                 * (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]) -
01276                    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
01277         
01278             b[56+i0] = dscales[0] * dscales[1] * dscales[2]                                     // d3V/dxdydz
01279                 * (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]) -
01280                    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]) +
01281                    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]) +
01282                    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]));
01283         }
01284         
01285         if (inds2[0] == 1 && inds2[1] == 1 && inds2[2] == 0) {
01286             DebugM(1, "Sub V = " << b[i0] << "\n");
01287         
01288             DebugM(1, "Sub dV/dx = " << b[8+i0] << "\n");
01289             DebugM(1, "Sub dV/dy = " << b[16+i0] << "\n");
01290             DebugM(1, "Sub dV/dz = " << b[24+i0] << "\n");
01291         
01292             DebugM(1, "Sub d2V/dxdy = " << b[32+i0] << "\n");
01293             DebugM(1, "Sub d2V/dxdz = " << b[40+i0] << "\n");
01294             DebugM(1, "Sub d2V/dydz = " << b[48+i0] << "\n");
01295         
01296             DebugM(1, "Sub d3V/dxdydz = " << b[56+i0] << "\n" << endi);
01297         }
01298     }
01299 }
01300 
01301 
01302 /*********************/
01303 /* GRIDFORCELITEGRID */
01304 /*********************/
01305 
01306 GridforceLiteGrid::GridforceLiteGrid(int gridnum)
01307 {
01308     mygridnum = gridnum;
01309     grid = NULL;
01310     type = GridforceGridTypeLite;
01311 }
01312 
01313 
01314 GridforceLiteGrid::~GridforceLiteGrid()
01315 {
01316     delete[] grid;
01317 }
01318 
01319 
01320 void GridforceLiteGrid::initialize(char *potfilename, SimParameters *simParams, MGridforceParams *mgridParams)
01321 {
01322     // cheat and use GridforceFullMainGrid to read the file
01323     GridforceFullMainGrid *tmp_grid = new GridforceFullMainGrid(mygridnum);
01324     tmp_grid->initialize(potfilename, simParams, mgridParams, 1);
01325     
01326     if (tmp_grid->get_total_grids() != 1) {
01327         NAMD_die("Cannot use gridforcelite option with multi-resolution grid!");
01328     }
01329     
01330     // save file name so that grid can be re-read via Tcl
01331     strcpy(filename, potfilename);
01332     
01333     // copy parameters
01334     k[0] = tmp_grid->get_k0();
01335     k[1] = tmp_grid->get_k1();
01336     k[2] = tmp_grid->get_k2();
01337     k[3] = 4;   // for V, dV/dx, dV/dy, dV/dz grids
01338     origin = tmp_grid->get_origin();
01339     center = tmp_grid->get_center();
01340     e = tmp_grid->get_e();
01341     inv = tmp_grid->get_inv();
01342     scale = tmp_grid->get_scale();
01343     
01344     // calculate rest of parameters
01345     size = k[0] * k[1] * k[2] * k[3];
01346     dk[0] = k[1] * k[2] * k[3];
01347     dk[1] = k[2] * k[3];
01348     dk[2] = k[3];
01349     dk[3] = 1;
01350     
01351     // copy the potential grid
01352     delete[] grid;
01353     grid = new float[size];
01354     for (int i0 = 0; i0 < k[0]; i0++) {
01355         for (int i1 = 0; i1 < k[1]; i1++) {
01356             for (int i2 = 0; i2 < k[2]; i2++) {
01357                 float V = tmp_grid->get_grid(i0, i1, i2);
01358                 set_grid(i0, i1, i2, 0, V);
01359                 DebugM(1, "V[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 0) << "(" << V << ")\n" << endi);
01360             }
01361         }
01362     }
01363     
01364     delete tmp_grid;
01365     
01366     compute_derivative_grids();
01367 }
01368 
01369 
01370 void GridforceLiteGrid::compute_derivative_grids(void)
01371 {
01372     // calculate derivative grids
01373     // separate loop so all V values have been set already
01374     for (int i0 = 0; i0 < k[0]; i0++) {
01375         for (int i1 = 0; i1 < k[1]; i1++) {
01376             for (int i2 = 0; i2 < k[2]; i2++) {
01377                 float dx, dy, dz;
01378                 if (i0 == 0 || i0 == k[0]-1 || i1 == 0 || i1 == k[1]-1 || i2 == 0 || i2 == k[2]-1) {
01379                     // on edge, set ALL derivatives to zero (make up for lack of padding)
01380                     dx = 0;
01381                     dy = 0;
01382                     dz = 0;
01383                 } else {
01384                     dx = 0.5 * (get_grid_d(i0+1,i1,i2,0) - get_grid_d(i0-1,i1,i2,0));
01385                     dy = 0.5 * (get_grid_d(i0,i1+1,i2,0) - get_grid_d(i0,i1-1,i2,0));
01386                     dz = 0.5 * (get_grid_d(i0,i1,i2+1,0) - get_grid_d(i0,i1,i2-1,0));
01387                 }
01388                 set_grid(i0, i1, i2, 1, dx);
01389                 set_grid(i0, i1, i2, 2, dy);
01390                 set_grid(i0, i1, i2, 3, dz);
01391                 DebugM(1, "dx[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 1) << "(" << dx << ")\n" << endi);
01392                 DebugM(1, "dy[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 2) << "(" << dy << ")\n" << endi);
01393                 DebugM(1, "dz[" << i0 << "," << i1 << "," << i2 << "] = " << get_grid(i0, i1, i2, 3) << "(" << dz << ")\n" << endi);
01394             }
01395         }
01396     }
01397 }
01398 
01399 
01400 void GridforceLiteGrid::reinitialize(SimParameters *simParams, MGridforceParams *mgridParams)
01401 {
01402     initialize(filename, simParams, mgridParams);
01403 }
01404 
01405 
01406 void GridforceLiteGrid::pack(MOStream *msg) const
01407 {
01408     msg->put(4*sizeof(int), (char*)k);
01409     msg->put(size);
01410     msg->put(4*sizeof(long int), (char*)dk);
01411     
01412     msg->put(sizeof(Vector), (char*)&origin);
01413     msg->put(sizeof(Vector), (char*)&center);
01414     msg->put(sizeof(Tensor), (char*)&e);
01415     msg->put(sizeof(Tensor), (char*)&inv);
01416     msg->put(sizeof(Vector), (char*)&scale);
01417     msg->put(sizeof(Bool), (char*)&checksize);
01418     
01419     msg->put(129*sizeof(char), (char*)filename);
01420     
01421     msg->put(size*sizeof(float), (char*)grid);
01422 }
01423 
01424 
01425 void GridforceLiteGrid::unpack(MIStream *msg)
01426 {
01427     delete[] grid;
01428     grid = NULL;
01429 
01430     msg->get(4*sizeof(int), (char*)k);
01431     msg->get(size);
01432     msg->get(4*sizeof(long int), (char*)dk);
01433     
01434     msg->get(sizeof(Vector), (char*)&origin);
01435     msg->get(sizeof(Vector), (char*)&center);
01436     msg->get(sizeof(Tensor), (char*)&e);
01437     msg->get(sizeof(Tensor), (char*)&inv);
01438     msg->get(sizeof(Vector), (char*)&scale);
01439     msg->get(sizeof(Bool), (char*)&checksize);
01440     
01441     msg->get(129*sizeof(char), (char*)filename);
01442     
01443     if (size) {
01444         grid = new float[size];
01445         msg->get(size*sizeof(float), (char*)grid);
01446     }
01447 }
01448 
01449 
01450 long int GridforceLiteGrid::get_all_gridvals(float** all_gridvals) const
01451 {
01452     // Creates a flat array of all grid values and puts it in the
01453     // value pointed to by the 'all_gridvals' argument. Returns the
01454     // resulting array size. Caller is responsible for destroying the
01455     // array via 'delete[]'
01456     
01457     DebugM(4, "GridforceLiteGrid::get_all_gridvals called\n" << endi);
01458     
01459     long int sz = size;
01460     DebugM(4, "size = " << sz << "\n" << endi);
01461     
01462     float *grid_vals = new float[sz];
01463     long int idx = 0;
01464     for (long int i = 0; i < size; i++) {
01465         grid_vals[idx++] = grid[i];
01466     }
01467     CmiAssert(idx == sz);
01468     
01469     *all_gridvals = grid_vals;
01470     
01471     DebugM(4, "GridforceLiteGrid::get_all_gridvals finished\n" << endi);
01472     
01473     return sz;
01474 }
01475 
01476 
01477 void GridforceLiteGrid::set_all_gridvals(float* all_gridvals, long int sz)
01478 {
01479     DebugM(4, "GridforceLiteGrid::set_all_gridvals called\n" << endi);
01480     
01481     long int sz_calc = size;
01482     CmiAssert(sz == sz_calc);
01483     
01484     long int idx = 0;
01485     for (long int i = 0; i < size; i++) {
01486         grid[i] = all_gridvals[idx++];
01487     }
01488     CmiAssert(idx == sz);
01489     
01490     //compute_derivative_grids();       // not needed if we're sending all 4 grids
01491     
01492     DebugM(4, "GridforceLiteGrid::set_all_gridvals finished\n" << endi);
01493 }

Generated on Wed Nov 22 01:17:15 2017 for NAMD by  doxygen 1.4.7