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
00023 #include "Debug.h"
00024
00025 #include "GridForceGrid.inl"
00026
00027
00028
00029
00030
00031
00032
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
00053
00054
00055
00056 msg->put(grid->get_grid_type());
00057 grid->pack(msg);
00058 }
00059
00060 GridforceGrid * GridforceGrid::unpack_grid(int gridnum, MIStream *msg)
00061 {
00062
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
00087
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
00104
00105
00106
00107
00108
00109
00110
00111
00112 Position pos;
00113 if (idx >= 8 || idx < 0) {
00114
00115 pos = Vector();
00116 } else if (corners[idx] != Vector()) {
00117
00118 pos = corners[idx];
00119 } else {
00120
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;
00127 DebugM(4, "corner " << idx << " = " << pos << "\n" << endi);
00128 }
00129 return pos;
00130 }
00131
00132
00133
00134
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*)¢er);
00172 msg->put(sizeof(Tensor), (char*)&e);
00173 msg->put(sizeof(Tensor), (char*)&inv);
00174
00175
00176
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
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*)¢er);
00227 msg->get(sizeof(Tensor), (char*)&e);
00228 msg->get(sizeof(Tensor), (char*)&inv);
00229
00230
00231
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);
00263 DebugM(4, "Read line: " << line << endi);
00264 } while (line[0] == '#');
00265 fseek(poten_fp, poten_offset, SEEK_SET);
00266
00267
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
00273 fscanf(poten_fp, "origin %lf %lf %lf\n",
00274 &origin.x, &origin.y, &origin.z);
00275
00276
00277
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
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
00342
00343
00344
00345 subgrids[i]->readSubgridHierarchy(poten, totalGrids);
00346 }
00347 }
00348
00349
00350
00351
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
00374
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
00390
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
00435
00436 poten_fp = Fopen(potfilename, "r");
00437 if (!poten_fp) {
00438 NAMD_die("Problem reading grid force potential file");
00439 }
00440
00441
00442 strcpy(filename, potfilename);
00443
00444
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);
00454
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
00470 readHeader(simParams, mgridParams);
00471
00472 factor = 1.0;
00473 if (mgridParams->gridforceVolts)
00474 {
00475 factor /= 0.0434;
00476 }
00477 scale = mgridParams->gridforceScale;
00478 checksize = mgridParams->gridforceCheckSize;
00479
00480
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;
00490 }
00491 fscanf(poten_fp, "\n");
00492
00493
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
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
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
00536
00537 }
00538 }
00539
00540
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
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
00578
00579
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
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
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
00648 return;
00649 }
00650
00651
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
00674 set_grid(i0, i1, i2, pad_n[dir]);
00675 } else if (var[dir] >= k[dir]-border) {
00676
00677 set_grid(i0, i1, i2, pad_p[dir]);
00678 }
00679 }
00680
00681
00682
00683
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
00697 DebugM(3, "clean up\n" << endi);
00698 delete[] grid_nopad;
00699
00700
00701 for (int i = 0; i < numSubgrids; i++) {
00702 subgrids[i]->poten_fp = poten_fp;
00703 subgrids[i]->initialize(simParams, mgridParams);
00704 }
00705
00706
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
00721
00722
00723
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;
00790 for (int i1 = 0; i1 < 3; i1++) {
00791 inds2[i1] = (inds[i1] + ((i0 & bit) ? 1 : 0)) % k[i1];
00792
00793
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;
00800 }
00801
00802 DebugM(1, "inds2 = " << inds2[0] << " " << inds2[1] << " " << inds2[2] << "\n" << endi);
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
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
00843
00844
00845
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
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]);
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]);
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]);
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]));
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]));
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]));
00874
00875 b[56+i0] = dscales[0] * dscales[1] * dscales[2]
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
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
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
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
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
00954 readHeader(simParams, mgridParams);
00955
00956 factor = 1.0;
00957 if (mgridParams->gridforceVolts)
00958 {
00959 factor /= 0.0434;
00960 }
00961 scale = mgridParams->gridforceScale;
00962 checksize = mgridParams->gridforceCheckSize;
00963
00964 for (int i = 0; i < 3; i++) {
00965 k[i] = k_nopad[i];
00966 }
00967
00968
00969
00970
00971
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
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
01007
01008 double TOL2 = 1e-4;
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
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
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
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
01068
01069
01070
01071
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
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
01103 delete[] grid_tmp;
01104
01105
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;
01168 for (int i1 = 0; i1 < 3; i1++) {
01169 inds2[i1] = (inds[i1] + ((i0 & bit) ? 1 : 0)) % k[i1];
01170
01171
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;
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
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
01218 Position pos = e * Vector(inds2[0], inds2[1], inds2[2]) + origin;
01219 Vector g = parent->inv * (pos - parent->origin);
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
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
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;
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]);
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]);
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]);
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]));
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]));
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]));
01277
01278 b[56+i0] = dscales[0] * dscales[1] * dscales[2]
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
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
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
01331 strcpy(filename, potfilename);
01332
01333
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;
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
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
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
01373
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
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*)¢er);
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*)¢er);
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
01453
01454
01455
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
01491
01492 DebugM(4, "GridforceLiteGrid::set_all_gridvals finished\n" << endi);
01493 }