00001
00007 #include "ComputeGridForce.h"
00008 #include "GridForceGrid.h"
00009 #include "Node.h"
00010 #include "SimParameters.h"
00011 #include "HomePatch.h"
00012 #include "Molecule.h"
00013
00014 #define MIN_DEBUG_LEVEL 4
00015
00016 #include "Debug.h"
00017
00018 #include "MGridforceParams.h"
00019
00020 #define GF_FORCE_OUTPUT
00021 #define GF_FORCE_OUTPUT_FREQ 100
00022
00023 ComputeGridForceNodeMgr::ComputeGridForceNodeMgr()
00024 {
00025 CkpvAccess(BOCclass_group).computeGridForceNodeMgr = thisgroup;
00026
00027 grids_deposited = 0;
00028 grids = NULL;
00029 requestGrids_count = 0;
00030 myNode = CkMyNode();
00031 nodeSize = CkNodeSize(myNode);
00032 gf = new ComputeGridForce*[nodeSize];
00033 proc_count = 0;
00034 cur_grid = 0;
00035
00036
00037
00038
00039
00040
00041 }
00042
00043 void ComputeGridForceNodeMgr::depositInitialGrid(GridDepositMsg *msg)
00044 {
00045 const int gridnum = msg->gridnum;
00046 GridforceGrid *grid = msg->grid;
00047 num_grids = msg->num_grids;
00048 delete msg;
00049
00050
00051
00052
00053 if (grids_deposited == 0) {
00054
00055
00056
00057 grids = new GridforceGrid*[num_grids];
00058 }
00059 grids[gridnum] = grid;
00060 grids_deposited++;
00061 if (grids_deposited == num_grids) {
00062
00063
00064 int i;
00065 for(i=0;i<num_grids; i++) {
00066 grids[i]->init2();
00067 }
00068
00069
00070
00071
00072 CProxy_ComputeGridForceNodeMgr myproxy(thisgroup);
00073 myproxy[0].requestInitialGridData();
00074 }
00075 }
00076
00077 void ComputeGridForceNodeMgr::requestInitialGridData()
00078 {
00079
00080
00081
00082 requestGrids_count++;
00083
00084
00085
00086
00087
00088 if (requestGrids_count == CkNumNodes()) {
00089 requestGrids_count = 0;
00090 bool sent = false;
00091 while (!sent && cur_grid < num_grids) {
00092 GridSegmentMsg *msg = new GridSegmentMsg;
00093
00094 msg->gridnum = cur_grid;
00095 grids[cur_grid]->init3(msg->grid,&(msg->start_index),
00096 &(msg->count));
00097
00098 if (msg->count > 0) {
00099
00100
00101
00102
00103
00104
00105
00106 CProxy_ComputeGridForceNodeMgr myproxy(thisgroup);
00107
00108
00109 myproxy.receiveInitialGridData(msg);
00110
00111
00112
00113 sent = true;
00114 } else {
00115
00116
00117
00118 cur_grid++;
00119 }
00120 }
00121
00122
00123
00124
00125 }
00126
00127 }
00128
00129 void ComputeGridForceNodeMgr::receiveInitialGridData(GridSegmentMsg *msg)
00130 {
00131
00132
00133
00134
00135
00136
00137 grids[msg->gridnum]->init4(msg->grid,msg->start_index,msg->count);
00138 delete msg;
00139
00140 CProxy_ComputeGridForceNodeMgr myproxy(thisgroup);
00141 myproxy[0].requestInitialGridData();
00142 }
00143
00144
00145
00146 void ComputeGridForceNodeMgr::submitRequest(SubReqMsg *msg)
00147 {
00148
00149
00150 const int msgrank = msg->rank;
00151
00152
00153 gf[msgrank] = msg->gf;
00154 delete msg;
00155 proc_count++;
00156
00157 if (proc_count < nodeSize) {
00158
00159
00160
00161
00162 return;
00163 } else proc_count = 0;
00164
00165
00166 int i,rank;
00167
00168
00169
00170 for(i=0; i < num_grids; i++) {
00171 grids[i]->consolidateGridRequests();
00172 }
00173
00174
00175 int num_indices = 0;
00176 for(i=0; i < num_grids; i++) {
00177 num_indices += grids[i]->getNumGridRequests();
00178 }
00179
00180 if (num_indices == 0) {
00181
00182
00183 resume();
00184 } else {
00185
00186
00187 msgSentCount=0;
00188 int *gridStart = new int[num_grids+1];
00189 int *gridIndices = new int[num_indices];
00190 int *gridNodeList = new int[num_indices];
00191 int *gridNodeCount = new int[CkNumNodes()];
00192 GridforceGrid::
00193 getGridIndices(gridStart,gridIndices,gridNodeList,gridNodeCount);
00194
00195 GridRequestMsg **msgs = new GridRequestMsg*[CkNumNodes()];
00196 int *cur_ind = new int[CkNumNodes()];
00197 int i;
00198 for(i=0; i < CkNumNodes(); i++) {
00199 msgs[i] = new(num_grids+1,gridNodeCount[i]) GridRequestMsg;
00200 msgs[i]->from_node = myNode;
00201 msgs[i]->num_grid_indices = gridNodeCount[i];
00202 cur_ind[i] = 0;
00203 }
00204
00205 int grid;
00206 CProxy_ComputeGridForceNodeMgr me(thisgroup);
00207
00208 for(grid=0; grid < num_grids; grid++) {
00209 for(i=0; i < CkNumNodes(); i++)
00210 msgs[i]->gridStartIndex[grid] = cur_ind[i];
00211
00212 for(i=gridStart[grid]; i < gridStart[grid+1]; i++) {
00213 int node = gridNodeList[i];
00214 msgs[node]->gridIndexList[cur_ind[node]] = gridIndices[i];
00215 cur_ind[node]++;
00216 }
00217 }
00218 for(i=0; i < CkNumNodes(); i++)
00219 msgs[i]->gridStartIndex[num_grids] = cur_ind[i];
00220
00221 for(i=0;i<CkNumNodes();i++) {
00222 if (gridNodeCount[i] > 0) {
00223
00224
00225 me[i].fetchGridValues(msgs[i]);
00226 msgSentCount++;
00227 } else {
00228 delete msgs[i];
00229 }
00230 }
00231 delete [] cur_ind;
00232 delete [] msgs;
00233 delete [] gridNodeCount;
00234 delete [] gridNodeList;
00235 delete [] gridIndices;
00236 delete [] gridStart;
00237 }
00238
00239
00240
00241 }
00242
00243 void ComputeGridForceNodeMgr::fetchGridValues(GridRequestMsg *msg)
00244 {
00245
00246
00247
00248
00249
00250
00251
00252 int *startIndex = new int[num_grids+1];
00253 startIndex[0] = 0;
00254 int gridnum;
00255 int i;
00256 for(gridnum=0;gridnum<num_grids;gridnum++) {
00257 int nvals = 0;
00258 for(i=msg->gridStartIndex[gridnum];
00259 i < msg->gridStartIndex[gridnum+1]; i++) {
00260 const int block = msg->gridIndexList[i];
00261 nvals += grids[gridnum]->getBlockSize(block);
00262
00263
00264
00265 }
00266 startIndex[gridnum+1] = startIndex[gridnum] + nvals;
00267 }
00268 GridValuesMsg *outmsg
00269 = new(num_grids+1,startIndex[num_grids]) GridValuesMsg;
00270 outmsg->num_grid_indices = startIndex[num_grids];
00271
00272 for(gridnum=0;gridnum<=num_grids;gridnum++) {
00273
00274
00275
00276 outmsg->gridStartIndex[gridnum] = startIndex[gridnum];
00277 }
00278 delete [] startIndex;
00279
00280
00281 for (gridnum = 0; gridnum < num_grids; gridnum++) {
00282 GridforceGrid* grid = grids[gridnum];
00283
00284
00285
00286
00287
00288 int outindex = outmsg->gridStartIndex[gridnum];
00289 for(i=msg->gridStartIndex[gridnum];
00290 i < msg->gridStartIndex[gridnum+1]; i++) {
00291 const int block = msg->gridIndexList[i];
00292 int j;
00293 const int jmax = grid->getBlockSize(block);
00294
00295
00296
00297
00298
00299 for(j=0; j < jmax; j++, outindex++) {
00300 const int idx = grid->getIndexForBlock(block,j);
00301
00302
00303 outmsg->gridVals[outindex].index = idx;
00304 outmsg->gridVals[outindex].val = grid->get_grid(idx);
00305
00306
00307
00308
00309 }
00310 }
00311 }
00312
00313 CProxy_ComputeGridForceNodeMgr me(thisgroup);
00314 me[msg->from_node].recvGridValues(outmsg);
00315 delete msg;
00316 }
00317
00318 void ComputeGridForceNodeMgr::recvGridValues(GridValuesMsg *msg)
00319 {
00320
00321
00322
00323 int gridnum;
00324 const int myrank = CkMyRank();
00325 for (gridnum = 0; gridnum < num_grids; gridnum++) {
00326 GridforceGrid *grid = grids[gridnum];
00327 int i;
00328
00329
00330 for(i=msg->gridStartIndex[gridnum];
00331 i < msg->gridStartIndex[gridnum+1]; i++) {
00332 const int idx = msg->gridVals[i].index;
00333 int bi,boffset;
00334 grid->getGridBlockIdxOffset(idx,bi,boffset);
00335
00336
00337 grid->addGridValue(idx,msg->gridVals[i].val);
00338
00339
00340
00341
00342
00343
00344 }
00345 }
00346 msgSentCount--;
00347
00348
00349
00350
00351 if (msgSentCount == 0) {
00352 resume();
00353 }
00354 delete msg;
00355 }
00356
00357 ComputeGridForceMgr::ComputeGridForceMgr()
00358 {
00359 CkpvAccess(BOCclass_group).computeGridForceMgr = thisgroup;
00360 }
00361
00362 void ComputeGridForceMgr::request(ComputeGridForce *mygf) {
00363
00364
00365 gf = mygf;
00366 SubReqMsg *msg = new SubReqMsg;
00367 msg->rank = CkMyRank();
00368 msg->gf = mygf;
00369 CProxy_ComputeGridForceNodeMgr nodemgr(CkpvAccess(BOCclass_group).
00370 computeGridForceNodeMgr);
00371 nodemgr[CkMyNode()].submitRequest(msg);
00372 }
00373
00374 void ComputeGridForceNodeMgr::resume()
00375 {
00376
00377
00378
00379
00380
00381 int rank;
00382 CProxy_ComputeGridForceMgr pemgr(CkpvAccess(BOCclass_group).
00383 computeGridForceMgr);
00384
00385 for(rank=0; rank < nodeSize; rank++) {
00386 pemgr[CkNodeFirst(myNode) + rank].finishWork();
00387 }
00388 }
00389
00390 void ComputeGridForceMgr::finishWork() {
00391 gf->finishWork();
00392 }
00393
00394 ComputeGridForce::ComputeGridForce(ComputeID c)
00395 : ComputeHomePatches(c)
00396 {
00397
00398 reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00399
00400 }
00401
00402
00403
00404 ComputeGridForce::~ComputeGridForce()
00405 {
00406 delete reduction;
00407 }
00408
00409
00410 void ComputeGridForce::doWork() {
00411 Molecule *mol = Node::Object()->molecule;
00412 ResizeArrayIter<PatchElem> ap(patchList);
00413 DebugM(3,patchID << ": doWork() called.\n");
00414
00415
00416 int i = 0;
00417 num_boxes = patchList.size();
00418 box_state = new BoxState[num_boxes];
00419
00420 for (ap = ap.begin(); ap != ap.end(); ap++) {
00421
00422 box_state[i].posBox = (*ap).positionBox;
00423 box_state[i].atoms = box_state[i].posBox->open();
00424 box_state[i].forceBox = (*ap).forceBox;
00425 box_state[i].results = box_state[i].forceBox->open();
00426 box_state[i].homePatch = (*ap).p;
00427 FullAtom* a = box_state[i].homePatch->getAtomList().begin();
00428 doForce(box_state[i].homePatch,a);
00429
00430
00431
00432 i++;
00433 }
00434 ComputeGridForceMgr* mgr = CProxy_ComputeGridForceMgr::
00435 ckLocalBranch(CkpvAccess(BOCclass_group).computeGridForceMgr);
00436 mgr->request(this);
00437
00438
00439
00440
00441
00442
00443
00444
00445 DebugM(2,patchID << ": doWork() completed.\n");
00446 }
00447
00448 void ComputeGridForce::finishWork()
00449 {
00450 Molecule *mol = Node::Object()->molecule;
00451 ResizeArrayIter<PatchElem> ap(patchList);
00452 DebugM(3,patchID << ": doWork() called.\n");
00453
00454
00455
00456 int i;
00457 for (i=0; i < num_boxes; i++) {
00458
00459
00460 CompAtom *x = box_state[i].atoms;
00461 Results *r = box_state[i].results;
00462 HomePatch *p = box_state[i].homePatch;
00463 FullAtom* a = p->getAtomList().begin();
00464 finishForce(p,a,r);
00465 box_state[i].posBox->close(&box_state[i].atoms);
00466 box_state[i].forceBox->close(&box_state[i].results);
00467
00468
00469 }
00470
00471
00472
00473
00474
00475 DebugM(2,patchID << ": finishWork() completed.\n");
00476 reduction->submit();
00477 delete [] box_state;
00478 box_state = 0;
00479 }
00480
00481 void ComputeGridForce::doForce(HomePatch* homePatch,FullAtom* p)
00482 {
00483 SimParameters *simParams = Node::Object()->simParameters;
00484 Molecule *mol = Node::Object()->molecule;
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497 int numAtoms = homePatch->getNumAtoms();
00498
00499 DebugM(4, "numGridforceGrids = " << mol->numGridforceGrids << "\n" << endi);
00500
00501 for (int gridnum = 0; gridnum < mol->numGridforceGrids; gridnum++) {
00502 GridforceGrid *grid = mol->get_gridfrc_grid(gridnum);
00503 DebugM(4, "scale = " << grid->get_scale() << "\n" << endi);
00504
00505 Position center = grid->get_center();
00506 Tensor inv = grid->get_inv();
00507 Vector gfScale = grid->get_scale();
00508
00509
00510 if (CkMyPe() == 0) {
00511 Position origin = grid->get_origin();
00512 Tensor e = grid->get_e();
00513
00514 for (int i0 = 0; i0 < 2; i0++) {
00515 float x = i0 * (grid->get_k0()-1);
00516 for (int i1 = 0; i1 < 2; i1++) {
00517 float y = i1 * (grid->get_k1()-1);
00518 for (int i2 = 0; i2 < 2; i2++) {
00519 float z = i2 * (grid->get_k2()-1);
00520
00521 Position pos_orig = origin + e * Position(x, y, z);
00522 Position pos = pos_orig + homePatch->lattice.wrap_delta(pos_orig);
00523 pos += homePatch->lattice.delta(pos, center) - (pos - center);
00524
00525 if (pos != pos_orig) {
00526 char err_msg[128];
00527
00528 sprintf(err_msg, "Gridforce grid too large for periodic cell: (%f %f %f) != (%f %f %f)",
00529 pos.x, pos.y, pos.z, pos_orig.x, pos_orig.y, pos_orig.z);
00530 NAMD_die(err_msg);
00531 }
00532
00533 }
00534 }
00535 }
00536 }
00537
00538
00539 for (int i = 0; i < numAtoms; i++) {
00540
00541
00542 if (mol->is_atom_gridforced(p[i].id, gridnum))
00543 {
00544
00545
00546
00547
00548 Position pos = p[i].position;
00549 pos += homePatch->lattice.wrap_delta(p[i].position);
00550 pos += homePatch->lattice.delta(pos, center) - (pos - center);
00551
00552 grid->request_box(pos);
00553 }
00554 }
00555 }
00556 }
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574 void ComputeGridForce::finishForce(HomePatch* homePatch,FullAtom* p,
00575 Results* r)
00576 {
00577 SimParameters *simParams = Node::Object()->simParameters;
00578 Molecule *mol = Node::Object()->molecule;
00579
00580
00581 Force *forces = r->f[Results::normal];
00582 BigReal energy = 0;
00583 Force extForce = 0.;
00584 Tensor extVirial;
00585
00586 Real scale;
00587 Charge charge;
00588 GridforceGrid::Box box;
00589 Force f;
00590 float v;
00591 int numAtoms = homePatch->getNumAtoms();
00592
00593 DebugM(4, "numGridforceGrids = " << mol->numGridforceGrids << "\n" << endi);
00594
00595 for (int gridnum = 0; gridnum < mol->numGridforceGrids; gridnum++) {
00596 const GridforceGrid *grid = mol->get_gridfrc_grid(gridnum);
00597 DebugM(4, "scale = " << grid->get_scale() << "\n" << endi);
00598
00599 Position center = grid->get_center();
00600 Tensor inv = grid->get_inv();
00601 Vector gfScale = grid->get_scale();
00602
00603 DebugM(4, "center = " << center << "\n" << endi);
00604
00605
00606 for (int i = 0; i < numAtoms; i++) {
00607
00608
00609 if (mol->is_atom_gridforced(p[i].id, gridnum))
00610 {
00611
00612 mol->get_gridfrc_params(scale, charge, p[i].id, gridnum);
00613
00614
00615 Position pos = p[i].position;
00616 pos += homePatch->lattice.wrap_delta(p[i].position);
00617 pos += homePatch->lattice.delta(pos, center) - (pos - center);
00618
00619 int err = grid->get_box(&box, pos);
00620 if (err) {
00621 DebugM(4, "ts = " << patch->flags.step << " gridnum = " << gridnum << " force4 = 0 0 0\n" << endi);
00622 DebugM(4, "v4 = 0\n" << endi);
00623
00624 continue;
00625 }
00626
00627 float a[64];
00628
00629
00630
00631
00632
00633
00634
00635
00636 a[0] = box.b[0];
00637 a[1] = box.b[8];
00638 a[2] = -3*box.b[0] + 3*box.b[1] - 2*box.b[8] - box.b[9];
00639 a[3] = 2*box.b[0] - 2*box.b[1] + box.b[8] + box.b[9];
00640 a[4] = box.b[16];
00641 a[5] = box.b[32];
00642 a[6] = -3*box.b[16] + 3*box.b[17] - 2*box.b[32] - box.b[33];
00643 a[7] = 2*box.b[16] - 2*box.b[17] + box.b[32] + box.b[33];
00644 a[8] = -3*box.b[0] + 3*box.b[2] - 2*box.b[16] - box.b[18];
00645 a[9] = -3*box.b[8] + 3*box.b[10] - 2*box.b[32] - box.b[34];
00646 a[10] = 9*box.b[0] - 9*box.b[1] - 9*box.b[2] + 9*box.b[3] + 6*box.b[8] + 3*box.b[9] - 6*box.b[10] - 3*box.b[11]
00647 + 6*box.b[16] - 6*box.b[17] + 3*box.b[18] - 3*box.b[19] + 4*box.b[32] + 2*box.b[33] + 2*box.b[34] + box.b[35];
00648 a[11] = -6*box.b[0] + 6*box.b[1] + 6*box.b[2] - 6*box.b[3] - 3*box.b[8] - 3*box.b[9] + 3*box.b[10] + 3*box.b[11]
00649 - 4*box.b[16] + 4*box.b[17] - 2*box.b[18] + 2*box.b[19] - 2*box.b[32] - 2*box.b[33] - box.b[34] - box.b[35];
00650 a[12] = 2*box.b[0] - 2*box.b[2] + box.b[16] + box.b[18];
00651 a[13] = 2*box.b[8] - 2*box.b[10] + box.b[32] + box.b[34];
00652 a[14] = -6*box.b[0] + 6*box.b[1] + 6*box.b[2] - 6*box.b[3] - 4*box.b[8] - 2*box.b[9] + 4*box.b[10] + 2*box.b[11]
00653 - 3*box.b[16] + 3*box.b[17] - 3*box.b[18] + 3*box.b[19] - 2*box.b[32] - box.b[33] - 2*box.b[34] - box.b[35];
00654 a[15] = 4*box.b[0] - 4*box.b[1] - 4*box.b[2] + 4*box.b[3] + 2*box.b[8] + 2*box.b[9] - 2*box.b[10] - 2*box.b[11]
00655 + 2*box.b[16] - 2*box.b[17] + 2*box.b[18] - 2*box.b[19] + box.b[32] + box.b[33] + box.b[34] + box.b[35];
00656 a[16] = box.b[24];
00657 a[17] = box.b[40];
00658 a[18] = -3*box.b[24] + 3*box.b[25] - 2*box.b[40] - box.b[41];
00659 a[19] = 2*box.b[24] - 2*box.b[25] + box.b[40] + box.b[41];
00660 a[20] = box.b[48];
00661 a[21] = box.b[56];
00662 a[22] = -3*box.b[48] + 3*box.b[49] - 2*box.b[56] - box.b[57];
00663 a[23] = 2*box.b[48] - 2*box.b[49] + box.b[56] + box.b[57];
00664 a[24] = -3*box.b[24] + 3*box.b[26] - 2*box.b[48] - box.b[50];
00665 a[25] = -3*box.b[40] + 3*box.b[42] - 2*box.b[56] - box.b[58];
00666 a[26] = 9*box.b[24] - 9*box.b[25] - 9*box.b[26] + 9*box.b[27] + 6*box.b[40] + 3*box.b[41] - 6*box.b[42] - 3*box.b[43]
00667 + 6*box.b[48] - 6*box.b[49] + 3*box.b[50] - 3*box.b[51] + 4*box.b[56] + 2*box.b[57] + 2*box.b[58] + box.b[59];
00668 a[27] = -6*box.b[24] + 6*box.b[25] + 6*box.b[26] - 6*box.b[27] - 3*box.b[40] - 3*box.b[41] + 3*box.b[42] + 3*box.b[43]
00669 - 4*box.b[48] + 4*box.b[49] - 2*box.b[50] + 2*box.b[51] - 2*box.b[56] - 2*box.b[57] - box.b[58] - box.b[59];
00670 a[28] = 2*box.b[24] - 2*box.b[26] + box.b[48] + box.b[50];
00671 a[29] = 2*box.b[40] - 2*box.b[42] + box.b[56] + box.b[58];
00672 a[30] = -6*box.b[24] + 6*box.b[25] + 6*box.b[26] - 6*box.b[27] - 4*box.b[40] - 2*box.b[41] + 4*box.b[42] + 2*box.b[43]
00673 - 3*box.b[48] + 3*box.b[49] - 3*box.b[50] + 3*box.b[51] - 2*box.b[56] - box.b[57] - 2*box.b[58] - box.b[59];
00674 a[31] = 4*box.b[24] - 4*box.b[25] - 4*box.b[26] + 4*box.b[27] + 2*box.b[40] + 2*box.b[41] - 2*box.b[42] - 2*box.b[43]
00675 + 2*box.b[48] - 2*box.b[49] + 2*box.b[50] - 2*box.b[51] + box.b[56] + box.b[57] + box.b[58] + box.b[59];
00676 a[32] = -3*box.b[0] + 3*box.b[4] - 2*box.b[24] - box.b[28];
00677 a[33] = -3*box.b[8] + 3*box.b[12] - 2*box.b[40] - box.b[44];
00678 a[34] = 9*box.b[0] - 9*box.b[1] - 9*box.b[4] + 9*box.b[5] + 6*box.b[8] + 3*box.b[9] - 6*box.b[12] - 3*box.b[13]
00679 + 6*box.b[24] - 6*box.b[25] + 3*box.b[28] - 3*box.b[29] + 4*box.b[40] + 2*box.b[41] + 2*box.b[44] + box.b[45];
00680 a[35] = -6*box.b[0] + 6*box.b[1] + 6*box.b[4] - 6*box.b[5] - 3*box.b[8] - 3*box.b[9] + 3*box.b[12] + 3*box.b[13]
00681 - 4*box.b[24] + 4*box.b[25] - 2*box.b[28] + 2*box.b[29] - 2*box.b[40] - 2*box.b[41] - box.b[44] - box.b[45];
00682 a[36] = -3*box.b[16] + 3*box.b[20] - 2*box.b[48] - box.b[52];
00683 a[37] = -3*box.b[32] + 3*box.b[36] - 2*box.b[56] - box.b[60];
00684 a[38] = 9*box.b[16] - 9*box.b[17] - 9*box.b[20] + 9*box.b[21] + 6*box.b[32] + 3*box.b[33] - 6*box.b[36] - 3*box.b[37]
00685 + 6*box.b[48] - 6*box.b[49] + 3*box.b[52] - 3*box.b[53] + 4*box.b[56] + 2*box.b[57] + 2*box.b[60] + box.b[61];
00686 a[39] = -6*box.b[16] + 6*box.b[17] + 6*box.b[20] - 6*box.b[21] - 3*box.b[32] - 3*box.b[33] + 3*box.b[36] + 3*box.b[37]
00687 - 4*box.b[48] + 4*box.b[49] - 2*box.b[52] + 2*box.b[53] - 2*box.b[56] - 2*box.b[57] - box.b[60] - box.b[61];
00688 a[40] = 9*box.b[0] - 9*box.b[2] - 9*box.b[4] + 9*box.b[6] + 6*box.b[16] + 3*box.b[18] - 6*box.b[20] - 3*box.b[22]
00689 + 6*box.b[24] - 6*box.b[26] + 3*box.b[28] - 3*box.b[30] + 4*box.b[48] + 2*box.b[50] + 2*box.b[52] + box.b[54];
00690 a[41] = 9*box.b[8] - 9*box.b[10] - 9*box.b[12] + 9*box.b[14] + 6*box.b[32] + 3*box.b[34] - 6*box.b[36] - 3*box.b[38]
00691 + 6*box.b[40] - 6*box.b[42] + 3*box.b[44] - 3*box.b[46] + 4*box.b[56] + 2*box.b[58] + 2*box.b[60] + box.b[62];
00692 a[42] = -27*box.b[0] + 27*box.b[1] + 27*box.b[2] - 27*box.b[3] + 27*box.b[4] - 27*box.b[5] - 27*box.b[6] + 27*box.b[7]
00693 - 18*box.b[8] - 9*box.b[9] + 18*box.b[10] + 9*box.b[11] + 18*box.b[12] + 9*box.b[13] - 18*box.b[14] - 9*box.b[15]
00694 - 18*box.b[16] + 18*box.b[17] - 9*box.b[18] + 9*box.b[19] + 18*box.b[20] - 18*box.b[21] + 9*box.b[22] - 9*box.b[23]
00695 - 18*box.b[24] + 18*box.b[25] + 18*box.b[26] - 18*box.b[27] - 9*box.b[28] + 9*box.b[29] + 9*box.b[30] - 9*box.b[31]
00696 - 12*box.b[32] - 6*box.b[33] - 6*box.b[34] - 3*box.b[35] + 12*box.b[36] + 6*box.b[37] + 6*box.b[38] + 3*box.b[39]
00697 - 12*box.b[40] - 6*box.b[41] + 12*box.b[42] + 6*box.b[43] - 6*box.b[44] - 3*box.b[45] + 6*box.b[46] + 3*box.b[47]
00698 - 12*box.b[48] + 12*box.b[49] - 6*box.b[50] + 6*box.b[51] - 6*box.b[52] + 6*box.b[53] - 3*box.b[54] + 3*box.b[55]
00699 - 8*box.b[56] - 4*box.b[57] - 4*box.b[58] - 2*box.b[59] - 4*box.b[60] - 2*box.b[61] - 2*box.b[62] - box.b[63];
00700 a[43] = 18*box.b[0] - 18*box.b[1] - 18*box.b[2] + 18*box.b[3] - 18*box.b[4] + 18*box.b[5] + 18*box.b[6] - 18*box.b[7]
00701 + 9*box.b[8] + 9*box.b[9] - 9*box.b[10] - 9*box.b[11] - 9*box.b[12] - 9*box.b[13] + 9*box.b[14] + 9*box.b[15]
00702 + 12*box.b[16] - 12*box.b[17] + 6*box.b[18] - 6*box.b[19] - 12*box.b[20] + 12*box.b[21] - 6*box.b[22] + 6*box.b[23]
00703 + 12*box.b[24] - 12*box.b[25] - 12*box.b[26] + 12*box.b[27] + 6*box.b[28] - 6*box.b[29] - 6*box.b[30] + 6*box.b[31]
00704 + 6*box.b[32] + 6*box.b[33] + 3*box.b[34] + 3*box.b[35] - 6*box.b[36] - 6*box.b[37] - 3*box.b[38] - 3*box.b[39]
00705 + 6*box.b[40] + 6*box.b[41] - 6*box.b[42] - 6*box.b[43] + 3*box.b[44] + 3*box.b[45] - 3*box.b[46] - 3*box.b[47]
00706 + 8*box.b[48] - 8*box.b[49] + 4*box.b[50] - 4*box.b[51] + 4*box.b[52] - 4*box.b[53] + 2*box.b[54] - 2*box.b[55]
00707 + 4*box.b[56] + 4*box.b[57] + 2*box.b[58] + 2*box.b[59] + 2*box.b[60] + 2*box.b[61] + box.b[62] + box.b[63];
00708 a[44] = -6*box.b[0] + 6*box.b[2] + 6*box.b[4] - 6*box.b[6] - 3*box.b[16] - 3*box.b[18] + 3*box.b[20] + 3*box.b[22]
00709 - 4*box.b[24] + 4*box.b[26] - 2*box.b[28] + 2*box.b[30] - 2*box.b[48] - 2*box.b[50] - box.b[52] - box.b[54];
00710 a[45] = -6*box.b[8] + 6*box.b[10] + 6*box.b[12] - 6*box.b[14] - 3*box.b[32] - 3*box.b[34] + 3*box.b[36] + 3*box.b[38]
00711 - 4*box.b[40] + 4*box.b[42] - 2*box.b[44] + 2*box.b[46] - 2*box.b[56] - 2*box.b[58] - box.b[60] - box.b[62];
00712 a[46] = 18*box.b[0] - 18*box.b[1] - 18*box.b[2] + 18*box.b[3] - 18*box.b[4] + 18*box.b[5] + 18*box.b[6] - 18*box.b[7]
00713 + 12*box.b[8] + 6*box.b[9] - 12*box.b[10] - 6*box.b[11] - 12*box.b[12] - 6*box.b[13] + 12*box.b[14] + 6*box.b[15]
00714 + 9*box.b[16] - 9*box.b[17] + 9*box.b[18] - 9*box.b[19] - 9*box.b[20] + 9*box.b[21] - 9*box.b[22] + 9*box.b[23]
00715 + 12*box.b[24] - 12*box.b[25] - 12*box.b[26] + 12*box.b[27] + 6*box.b[28] - 6*box.b[29] - 6*box.b[30] + 6*box.b[31]
00716 + 6*box.b[32] + 3*box.b[33] + 6*box.b[34] + 3*box.b[35] - 6*box.b[36] - 3*box.b[37] - 6*box.b[38] - 3*box.b[39]
00717 + 8*box.b[40] + 4*box.b[41] - 8*box.b[42] - 4*box.b[43] + 4*box.b[44] + 2*box.b[45] - 4*box.b[46] - 2*box.b[47]
00718 + 6*box.b[48] - 6*box.b[49] + 6*box.b[50] - 6*box.b[51] + 3*box.b[52] - 3*box.b[53] + 3*box.b[54] - 3*box.b[55]
00719 + 4*box.b[56] + 2*box.b[57] + 4*box.b[58] + 2*box.b[59] + 2*box.b[60] + box.b[61] + 2*box.b[62] + box.b[63];
00720 a[47] = -12*box.b[0] + 12*box.b[1] + 12*box.b[2] - 12*box.b[3] + 12*box.b[4] - 12*box.b[5] - 12*box.b[6] + 12*box.b[7]
00721 - 6*box.b[8] - 6*box.b[9] + 6*box.b[10] + 6*box.b[11] + 6*box.b[12] + 6*box.b[13] - 6*box.b[14] - 6*box.b[15]
00722 - 6*box.b[16] + 6*box.b[17] - 6*box.b[18] + 6*box.b[19] + 6*box.b[20] - 6*box.b[21] + 6*box.b[22] - 6*box.b[23]
00723 - 8*box.b[24] + 8*box.b[25] + 8*box.b[26] - 8*box.b[27] - 4*box.b[28] + 4*box.b[29] + 4*box.b[30] - 4*box.b[31]
00724 - 3*box.b[32] - 3*box.b[33] - 3*box.b[34] - 3*box.b[35] + 3*box.b[36] + 3*box.b[37] + 3*box.b[38] + 3*box.b[39]
00725 - 4*box.b[40] - 4*box.b[41] + 4*box.b[42] + 4*box.b[43] - 2*box.b[44] - 2*box.b[45] + 2*box.b[46] + 2*box.b[47]
00726 - 4*box.b[48] + 4*box.b[49] - 4*box.b[50] + 4*box.b[51] - 2*box.b[52] + 2*box.b[53] - 2*box.b[54] + 2*box.b[55]
00727 - 2*box.b[56] - 2*box.b[57] - 2*box.b[58] - 2*box.b[59] - box.b[60] - box.b[61] - box.b[62] - box.b[63];
00728 a[48] = 2*box.b[0] - 2*box.b[4] + box.b[24] + box.b[28];
00729 a[49] = 2*box.b[8] - 2*box.b[12] + box.b[40] + box.b[44];
00730 a[50] = -6*box.b[0] + 6*box.b[1] + 6*box.b[4] - 6*box.b[5] - 4*box.b[8] - 2*box.b[9] + 4*box.b[12] + 2*box.b[13]
00731 - 3*box.b[24] + 3*box.b[25] - 3*box.b[28] + 3*box.b[29] - 2*box.b[40] - box.b[41] - 2*box.b[44] - box.b[45];
00732 a[51] = 4*box.b[0] - 4*box.b[1] - 4*box.b[4] + 4*box.b[5] + 2*box.b[8] + 2*box.b[9] - 2*box.b[12] - 2*box.b[13]
00733 + 2*box.b[24] - 2*box.b[25] + 2*box.b[28] - 2*box.b[29] + box.b[40] + box.b[41] + box.b[44] + box.b[45];
00734 a[52] = 2*box.b[16] - 2*box.b[20] + box.b[48] + box.b[52];
00735 a[53] = 2*box.b[32] - 2*box.b[36] + box.b[56] + box.b[60];
00736 a[54] = -6*box.b[16] + 6*box.b[17] + 6*box.b[20] - 6*box.b[21] - 4*box.b[32] - 2*box.b[33] + 4*box.b[36] + 2*box.b[37]
00737 - 3*box.b[48] + 3*box.b[49] - 3*box.b[52] + 3*box.b[53] - 2*box.b[56] - box.b[57] - 2*box.b[60] - box.b[61];
00738 a[55] = 4*box.b[16] - 4*box.b[17] - 4*box.b[20] + 4*box.b[21] + 2*box.b[32] + 2*box.b[33] - 2*box.b[36] - 2*box.b[37]
00739 + 2*box.b[48] - 2*box.b[49] + 2*box.b[52] - 2*box.b[53] + box.b[56] + box.b[57] + box.b[60] + box.b[61];
00740 a[56] = -6*box.b[0] + 6*box.b[2] + 6*box.b[4] - 6*box.b[6] - 4*box.b[16] - 2*box.b[18] + 4*box.b[20] + 2*box.b[22]
00741 - 3*box.b[24] + 3*box.b[26] - 3*box.b[28] + 3*box.b[30] - 2*box.b[48] - box.b[50] - 2*box.b[52] - box.b[54];
00742 a[57] = -6*box.b[8] + 6*box.b[10] + 6*box.b[12] - 6*box.b[14] - 4*box.b[32] - 2*box.b[34] + 4*box.b[36] + 2*box.b[38]
00743 - 3*box.b[40] + 3*box.b[42] - 3*box.b[44] + 3*box.b[46] - 2*box.b[56] - box.b[58] - 2*box.b[60] - box.b[62];
00744 a[58] = 18*box.b[0] - 18*box.b[1] - 18*box.b[2] + 18*box.b[3] - 18*box.b[4] + 18*box.b[5] + 18*box.b[6] - 18*box.b[7]
00745 + 12*box.b[8] + 6*box.b[9] - 12*box.b[10] - 6*box.b[11] - 12*box.b[12] - 6*box.b[13] + 12*box.b[14] + 6*box.b[15]
00746 + 12*box.b[16] - 12*box.b[17] + 6*box.b[18] - 6*box.b[19] - 12*box.b[20] + 12*box.b[21] - 6*box.b[22] + 6*box.b[23]
00747 + 9*box.b[24] - 9*box.b[25] - 9*box.b[26] + 9*box.b[27] + 9*box.b[28] - 9*box.b[29] - 9*box.b[30] + 9*box.b[31]
00748 + 8*box.b[32] + 4*box.b[33] + 4*box.b[34] + 2*box.b[35] - 8*box.b[36] - 4*box.b[37] - 4*box.b[38] - 2*box.b[39]
00749 + 6*box.b[40] + 3*box.b[41] - 6*box.b[42] - 3*box.b[43] + 6*box.b[44] + 3*box.b[45] - 6*box.b[46] - 3*box.b[47]
00750 + 6*box.b[48] - 6*box.b[49] + 3*box.b[50] - 3*box.b[51] + 6*box.b[52] - 6*box.b[53] + 3*box.b[54] - 3*box.b[55]
00751 + 4*box.b[56] + 2*box.b[57] + 2*box.b[58] + box.b[59] + 4*box.b[60] + 2*box.b[61] + 2*box.b[62] + box.b[63];
00752 a[59] = -12*box.b[0] + 12*box.b[1] + 12*box.b[2] - 12*box.b[3] + 12*box.b[4] - 12*box.b[5] - 12*box.b[6] + 12*box.b[7]
00753 - 6*box.b[8] - 6*box.b[9] + 6*box.b[10] + 6*box.b[11] + 6*box.b[12] + 6*box.b[13] - 6*box.b[14] - 6*box.b[15]
00754 - 8*box.b[16] + 8*box.b[17] - 4*box.b[18] + 4*box.b[19] + 8*box.b[20] - 8*box.b[21] + 4*box.b[22] - 4*box.b[23]
00755 - 6*box.b[24] + 6*box.b[25] + 6*box.b[26] - 6*box.b[27] - 6*box.b[28] + 6*box.b[29] + 6*box.b[30] - 6*box.b[31]
00756 - 4*box.b[32] - 4*box.b[33] - 2*box.b[34] - 2*box.b[35] + 4*box.b[36] + 4*box.b[37] + 2*box.b[38] + 2*box.b[39]
00757 - 3*box.b[40] - 3*box.b[41] + 3*box.b[42] + 3*box.b[43] - 3*box.b[44] - 3*box.b[45] + 3*box.b[46] + 3*box.b[47]
00758 - 4*box.b[48] + 4*box.b[49] - 2*box.b[50] + 2*box.b[51] - 4*box.b[52] + 4*box.b[53] - 2*box.b[54] + 2*box.b[55]
00759 - 2*box.b[56] - 2*box.b[57] - box.b[58] - box.b[59] - 2*box.b[60] - 2*box.b[61] - box.b[62] - box.b[63];
00760 a[60] = 4*box.b[0] - 4*box.b[2] - 4*box.b[4] + 4*box.b[6] + 2*box.b[16] + 2*box.b[18] - 2*box.b[20] - 2*box.b[22]
00761 + 2*box.b[24] - 2*box.b[26] + 2*box.b[28] - 2*box.b[30] + box.b[48] + box.b[50] + box.b[52] + box.b[54];
00762 a[61] = 4*box.b[8] - 4*box.b[10] - 4*box.b[12] + 4*box.b[14] + 2*box.b[32] + 2*box.b[34] - 2*box.b[36] - 2*box.b[38]
00763 + 2*box.b[40] - 2*box.b[42] + 2*box.b[44] - 2*box.b[46] + box.b[56] + box.b[58] + box.b[60] + box.b[62];
00764 a[62] = -12*box.b[0] + 12*box.b[1] + 12*box.b[2] - 12*box.b[3] + 12*box.b[4] - 12*box.b[5] - 12*box.b[6] + 12*box.b[7]
00765 - 8*box.b[8] - 4*box.b[9] + 8*box.b[10] + 4*box.b[11] + 8*box.b[12] + 4*box.b[13] - 8*box.b[14] - 4*box.b[15]
00766 - 6*box.b[16] + 6*box.b[17] - 6*box.b[18] + 6*box.b[19] + 6*box.b[20] - 6*box.b[21] + 6*box.b[22] - 6*box.b[23]
00767 - 6*box.b[24] + 6*box.b[25] + 6*box.b[26] - 6*box.b[27] - 6*box.b[28] + 6*box.b[29] + 6*box.b[30] - 6*box.b[31]
00768 - 4*box.b[32] - 2*box.b[33] - 4*box.b[34] - 2*box.b[35] + 4*box.b[36] + 2*box.b[37] + 4*box.b[38] + 2*box.b[39]
00769 - 4*box.b[40] - 2*box.b[41] + 4*box.b[42] + 2*box.b[43] - 4*box.b[44] - 2*box.b[45] + 4*box.b[46] + 2*box.b[47]
00770 - 3*box.b[48] + 3*box.b[49] - 3*box.b[50] + 3*box.b[51] - 3*box.b[52] + 3*box.b[53] - 3*box.b[54] + 3*box.b[55]
00771 - 2*box.b[56] - box.b[57] - 2*box.b[58] - box.b[59] - 2*box.b[60] - box.b[61] - 2*box.b[62] - box.b[63];
00772 a[63] = 8*box.b[0] - 8*box.b[1] - 8*box.b[2] + 8*box.b[3] - 8*box.b[4] + 8*box.b[5] + 8*box.b[6] - 8*box.b[7]
00773 + 4*box.b[8] + 4*box.b[9] - 4*box.b[10] - 4*box.b[11] - 4*box.b[12] - 4*box.b[13] + 4*box.b[14] + 4*box.b[15]
00774 + 4*box.b[16] - 4*box.b[17] + 4*box.b[18] - 4*box.b[19] - 4*box.b[20] + 4*box.b[21] - 4*box.b[22] + 4*box.b[23]
00775 + 4*box.b[24] - 4*box.b[25] - 4*box.b[26] + 4*box.b[27] + 4*box.b[28] - 4*box.b[29] - 4*box.b[30] + 4*box.b[31]
00776 + 2*box.b[32] + 2*box.b[33] + 2*box.b[34] + 2*box.b[35] - 2*box.b[36] - 2*box.b[37] - 2*box.b[38] - 2*box.b[39]
00777 + 2*box.b[40] + 2*box.b[41] - 2*box.b[42] - 2*box.b[43] + 2*box.b[44] + 2*box.b[45] - 2*box.b[46] - 2*box.b[47]
00778 + 2*box.b[48] - 2*box.b[49] + 2*box.b[50] - 2*box.b[51] + 2*box.b[52] - 2*box.b[53] + 2*box.b[54] - 2*box.b[55]
00779 + box.b[56] + box.b[57] + box.b[58] + box.b[59] + box.b[60] + box.b[61] + box.b[62] + box.b[63];
00780
00781 for (int j = 0; j < 64; j++) DebugM(2, "a[" << j << "] = " << a[j] << "\n" << endi);
00782 for (int j = 0; j < 64; j++) DebugM(2, "b[" << j << "] = " << box.b[j] << "\n" << endi);
00783
00784
00785
00786 float x[4], y[4], z[4];
00787 x[0] = 1; y[0] = 1; z[0] = 1;
00788 for (int j = 1; j < 4; j++) {
00789 x[j] = x[j-1] * box.loc.x;
00790 y[j] = y[j-1] * box.loc.y;
00791 z[j] = z[j-1] * box.loc.z;
00792 }
00793
00794 int ind = 0;
00795 f = 0;
00796 v = 0;
00797 for (int l = 0; l < 4; l++) {
00798 for (int k = 0; k < 4; k++) {
00799 for (int j = 0; j < 4; j++) {
00800 v += a[ind] * x[j] * y[k] * z[l];
00801 if (j > 0) f.x -= a[ind] * j * x[j-1] * y[k] * z[l];
00802 if (k > 0) f.y -= a[ind] * k * x[j] * y[k-1] * z[l];
00803 if (l > 0) f.z -= a[ind] * l * x[j] * y[k] * z[l-1];
00804 ind++;
00805 }
00806 }
00807 }
00808
00809
00810
00811 Force force = scale * Tensor::diagonal(gfScale) * charge * ((box.scale * f) * inv);
00812
00813 DebugM(4, "f4 = " << f << "\n" << endi);
00814 DebugM(4, "ts = " << patch->flags.step << " gridnum = " << gridnum << " force4 = " << force << "\n" << endi);
00815 DebugM(4, "v4 = " << v << "\n" << endi);
00816
00817
00818
00819
00820
00821
00822
00823
00824 forces[i] += force;
00825 extForce += force;
00826 Position vpos = homePatch->lattice.reverse_transform(p[i].position, p[i].transform);
00827
00828 DebugM(4, "transform = " << (int)p[i].transform.i << " "
00829 << (int)p[i].transform.j << " " << (int)p[i].transform.k << "\n" << endi);
00830
00831
00832 if (gfScale.x == gfScale.y && gfScale.x == gfScale.z)
00833 {
00834
00835 energy += v * charge * scale * gfScale.x;
00836 }
00837 extVirial += outer(force,vpos);
00838 }
00839 }
00840
00841 }
00842
00843
00844
00845
00846
00847 reduction->item(REDUCTION_MISC_ENERGY) += energy;
00848 ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,extForce);
00849 ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,extVirial);
00850
00851 }
00852
00853
00854
00855 #include "ComputeGridForceMgr.def.h"