00001
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #include <math.h>
00018 #include "charm++.h"
00019
00020 #include "SimParameters.h"
00021 #include "HomePatch.h"
00022 #include "AtomMap.h"
00023 #include "Node.h"
00024 #include "PatchMap.inl"
00025 #include "main.h"
00026 #include "ProxyMgr.decl.h"
00027 #include "ProxyMgr.h"
00028 #include "Migration.h"
00029 #include "Molecule.h"
00030 #include "PatchMgr.h"
00031 #include "Sequencer.h"
00032 #include "LdbCoordinator.h"
00033 #include "Settle.h"
00034 #include "ReductionMgr.h"
00035 #include "Sync.h"
00036 #include "Random.h"
00037 #include "Priorities.h"
00038
00039 #define TINY 1.0e-20;
00040 #define MAXHGS 10
00041 #define MIN_DEBUG_LEVEL 2
00042
00043 #include "Debug.h"
00044
00045 typedef int HGArrayInt[MAXHGS];
00046 typedef BigReal HGArrayBigReal[MAXHGS];
00047 typedef zVector HGArrayVector[MAXHGS];
00048 typedef BigReal HGMatrixBigReal[MAXHGS][MAXHGS];
00049 typedef zVector HGMatrixVector[MAXHGS][MAXHGS];
00050
00051 int average(CompAtom *qtilde,const HGArrayVector &q,BigReal *lambda,const int n,const int m, const HGArrayBigReal &imass, const HGArrayBigReal &length2, const HGArrayInt &ial, const HGArrayInt &ibl, const HGArrayVector &refab, const BigReal tolf, const int ntrial);
00052
00053 void mollify(CompAtom *qtilde,const HGArrayVector &q0,const BigReal *lambda, HGArrayVector &force,const int n, const int m, const HGArrayBigReal &imass,const HGArrayInt &ial,const HGArrayInt &ibl,const HGArrayVector &refab);
00054
00055
00056
00057 #if NAMD_SeparateWaters != 0
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068 #define IS_HYDROGEN_GROUP_WATER(hgs, mass) \
00069 ((hgs >= 3) && ((mass >= 14.0) && (mass <= 18.0)))
00070
00071 #endif
00072
00073
00074 HomePatch::HomePatch(PatchID pd, int atomCnt) : Patch(pd)
00075
00076 #if NAMD_SeparateWaters != 0
00077 ,tempAtom()
00078 #endif
00079 {
00080 min.x = PatchMap::Object()->min_a(patchID);
00081 min.y = PatchMap::Object()->min_b(patchID);
00082 min.z = PatchMap::Object()->min_c(patchID);
00083 max.x = PatchMap::Object()->max_a(patchID);
00084 max.y = PatchMap::Object()->max_b(patchID);
00085 max.z = PatchMap::Object()->max_c(patchID);
00086 center = 0.5*(min+max);
00087
00088 aAway = PatchMap::Object()->numaway_a();
00089 bAway = PatchMap::Object()->numaway_b();
00090 cAway = PatchMap::Object()->numaway_c();
00091
00092 migrationSuspended = false;
00093 allMigrationIn = false;
00094 marginViolations = 0;
00095 patchMapRead = 0;
00096
00097 inMigration = false;
00098 numMlBuf = 0;
00099 flags.sequence = -1;
00100
00101 numAtoms = atomCnt;
00102 replacementForces = 0;
00103
00104 SimParameters *simParams = Node::Object()->simParameters;
00105 doPairlistCheck_newTolerance =
00106 0.5 * ( simParams->pairlistDist - simParams->cutoff );
00107
00108
00109 numFixedAtoms = 0;
00110
00111
00112
00113
00114
00115
00116 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00117 ptnTree.resize(0);
00118
00119
00120 #else
00121 child = new int[proxySpanDim];
00122 nChild = 0;
00123 #endif
00124
00125 #if CMK_PERSISTENT_COMM
00126 phsReady = 0;
00127 nphs = 0;
00128 localphs = new PersistentHandle[CkNumPes()];
00129 for (int i=0; i<CkNumPes(); i++) localphs[i] = 0;
00130 #endif
00131
00132
00133
00134 #if NAMD_SeparateWaters != 0
00135
00136
00137 tempAtom.resize(numAtoms);
00138 numWaterAtoms = -1;
00139
00140 #endif
00141
00142 if (simParams->watmodel == WAT_TIP4) init_tip4();
00143 else if (simParams->watmodel == WAT_SWM4) init_swm4();
00144
00145 isNewProxyAdded = 0;
00146 }
00147
00148 HomePatch::HomePatch(PatchID pd, FullAtomList al) : Patch(pd), atom(al)
00149
00150 #if NAMD_SeparateWaters != 0
00151 ,tempAtom()
00152 #endif
00153 {
00154 min.x = PatchMap::Object()->min_a(patchID);
00155 min.y = PatchMap::Object()->min_b(patchID);
00156 min.z = PatchMap::Object()->min_c(patchID);
00157 max.x = PatchMap::Object()->max_a(patchID);
00158 max.y = PatchMap::Object()->max_b(patchID);
00159 max.z = PatchMap::Object()->max_c(patchID);
00160 center = 0.5*(min+max);
00161
00162 aAway = PatchMap::Object()->numaway_a();
00163 bAway = PatchMap::Object()->numaway_b();
00164 cAway = PatchMap::Object()->numaway_c();
00165
00166 migrationSuspended = false;
00167 allMigrationIn = false;
00168 marginViolations = 0;
00169 patchMapRead = 0;
00170
00171 inMigration = false;
00172 numMlBuf = 0;
00173 flags.sequence = -1;
00174
00175 numAtoms = atom.size();
00176 replacementForces = 0;
00177
00178 SimParameters *simParams = Node::Object()->simParameters;
00179 doPairlistCheck_newTolerance =
00180 0.5 * ( simParams->pairlistDist - simParams->cutoff );
00181
00182
00183 numFixedAtoms = 0;
00184 if ( simParams->fixedAtomsOn ) {
00185 for ( int i = 0; i < numAtoms; ++i ) {
00186 numFixedAtoms += ( atom[i].atomFixed ? 1 : 0 );
00187 }
00188 }
00189
00190 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00191 ptnTree.resize(0);
00192
00193
00194 #else
00195 child = new int[proxySpanDim];
00196 nChild = 0;
00197 #endif
00198
00199 #if CMK_PERSISTENT_COMM
00200 phsReady = 0;
00201 nphs = 0;
00202 localphs = new PersistentHandle[CkNumPes()];
00203 for (int i=0; i<CkNumPes(); i++) localphs[i] = 0;
00204 #endif
00205
00206
00207
00208 #if NAMD_SeparateWaters != 0
00209
00210
00211 tempAtom.resize(numAtoms);
00212 numWaterAtoms = -1;
00213
00214
00215 separateAtoms();
00216
00217 #endif
00218
00219
00220 if (simParams->watmodel == WAT_TIP4) init_tip4();
00221 else if (simParams->watmodel == WAT_SWM4) init_swm4();
00222
00223 isNewProxyAdded = 0;
00224
00225 }
00226
00227 void HomePatch::write_tip4_props() {
00228 printf("Writing r_om and r_ohc: %f | %f\n", r_om, r_ohc);
00229 }
00230
00231 void HomePatch::init_tip4() {
00232
00233
00234 Molecule *mol = Node::Object()->molecule;
00235 r_om = mol->r_om;
00236 r_ohc = mol->r_ohc;
00237 }
00238
00239
00240 void ::HomePatch::init_swm4() {
00241
00242 SimParameters *simParams = Node::Object()->simParameters;
00243 Molecule *mol = Node::Object()->molecule;
00244 int ig;
00245 if (RIGID_NONE == simParams->rigidBonds) return;
00246 if (WAT_SWM4 != simParams->watmodel) return;
00247 for (ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
00248
00249 if (mol->rigid_bond_length(atom[ig].id) > 0) {
00250
00251 BigReal r_hh = mol->rigid_bond_length(atom[ig].id);
00252 BigReal r_oh = mol->rigid_bond_length(atom[ig+3].id);
00253 r_om = mol->rigid_bond_length(atom[ig+2].id);
00254 r_ohc = sqrt(r_oh * r_oh - 0.25 * r_hh * r_hh);
00255
00256 break;
00257 }
00258 }
00259 }
00260
00261 void HomePatch::reinitAtoms(FullAtomList al) {
00262 atom = al;
00263 numAtoms = atom.size();
00264
00265
00266 #if NAMD_SeparateWaters != 0
00267
00268
00269 numWaterAtoms = -1;
00270
00271
00272 separateAtoms();
00273
00274 #endif
00275 }
00276
00277
00278 void HomePatch::useSequencer(Sequencer *sequencerPtr)
00279 { sequencer=sequencerPtr; }
00280
00281
00282 void HomePatch::runSequencer(void)
00283 { sequencer->run(); }
00284
00285 void HomePatch::readPatchMap() {
00286
00287 PatchMap *p = PatchMap::Object();
00288 PatchID nnPatchID[PatchMap::MaxOneAway];
00289
00290 patchMigrationCounter = numNeighbors
00291 = PatchMap::Object()->oneAwayNeighbors(patchID, nnPatchID);
00292 DebugM( 1, "NumNeighbors for pid " <<patchID<<" is "<< numNeighbors << "\n");
00293 int n;
00294 for (n=0; n<numNeighbors; n++) {
00295 realInfo[n].destNodeID = p->node(realInfo[n].destPatchID = nnPatchID[n]);
00296 DebugM( 1, " nnPatchID=" <<nnPatchID[n]<<" nnNodeID="<< realInfo[n].destNodeID<< "\n");
00297 realInfo[n].mList.resize(0);
00298 }
00299
00300
00301 for (int i=0; i<3; i++)
00302 for (int j=0; j<3; j++)
00303 for (int k=0; k<3; k++)
00304 {
00305 int pid = p->pid(p->index_a(patchID)+i-1,
00306 p->index_b(patchID)+j-1, p->index_c(patchID)+k-1);
00307 if (pid < 0) {
00308 DebugM(5, "ERROR, for patchID " << patchID <<" I got neigh pid = " << pid << "\n");
00309 }
00310 if (pid == patchID && ! (
00311 ( (i-1) && p->periodic_a() ) ||
00312 ( (j-1) && p->periodic_b() ) ||
00313 ( (k-1) && p->periodic_c() ) )) {
00314 mInfo[i][j][k] = NULL;
00315 }
00316 else {
00317
00318
00319 for (n = 0; n<numNeighbors; n++) {
00320 if (pid == realInfo[n].destPatchID) {
00321 mInfo[i][j][k] = &realInfo[n];
00322 break;
00323 }
00324 }
00325 if (n == numNeighbors) {
00326 DebugM(4,"BAD News, I could not find PID " << pid << "\n");
00327 }
00328 }
00329 }
00330
00331 DebugM(1,"Patch("<<patchID<<") # of neighbors = " << numNeighbors << "\n");
00332 }
00333
00334 HomePatch::~HomePatch()
00335 {
00336 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00337 ptnTree.resize(0);
00338 delete [] children;
00339 #ifdef USE_NODEPATCHMGR
00340 delete [] nodeChildren;
00341 #endif
00342 #else
00343 delete [] child;
00344 #endif
00345 }
00346
00347
00348 void HomePatch::boxClosed(int)
00349 {
00350 if ( ! --boxesOpen )
00351 {
00352 if ( replacementForces ) {
00353 for ( int i = 0; i < numAtoms; ++i ) {
00354 if ( replacementForces[i].replace ) {
00355 for ( int j = 0; j < Results::maxNumForces; ++j ) { f[j][i] = 0; }
00356 f[Results::normal][i] = replacementForces[i].force;
00357 }
00358 }
00359 replacementForces = 0;
00360 }
00361 DebugM(1,patchID << ": " << CthSelf() << " awakening sequencer "
00362 << sequencer->thread << "(" << patchID << ") @" << CmiTimer() << "\n");
00363
00364 sequencer->awaken();
00365 return;
00366 }
00367 else
00368 {
00369 DebugM(1,patchID << ": " << boxesOpen << " boxes left to close.\n");
00370 }
00371 }
00372
00373 void HomePatch::registerProxy(RegisterProxyMsg *msg) {
00374 DebugM(4, "registerProxy("<<patchID<<") - adding node " <<msg->node<<"\n");
00375 proxy.add(ProxyListElem(msg->node,forceBox.checkOut()));
00376
00377 isNewProxyAdded = 1;
00378
00379 Random((patchID + 37) * 137).reorder(proxy.begin(),proxy.size());
00380 delete msg;
00381 }
00382
00383 void HomePatch::unregisterProxy(UnregisterProxyMsg *msg) {
00384 int n = msg->node;
00385 ProxyListElem *pe = proxy.begin();
00386 for ( ; pe->node != n; ++pe );
00387 forceBox.checkIn(pe->forceBox);
00388 proxy.del(pe - proxy.begin());
00389 delete msg;
00390 }
00391
00392 #if USE_TOPOMAP && USE_SPANNING_TREE
00393
00394 int HomePatch::findSubroots(int dim, int* subroots, int psize, int* pidscopy){
00395 int nChild = 0;
00396 int cones[6][proxySpanDim*proxySpanDim+proxySpanDim];
00397 int conesizes[6] = {0,0,0,0,0,0};
00398 int conecounters[6] = {0,0,0,0,0,0};
00399 int childcounter = 0;
00400 nChild = (psize>proxySpanDim)?proxySpanDim:psize;
00401 TopoManager tmgr;
00402 for(int i=0;i<psize;i++){
00403 int cone = tmgr.getConeNumberForRank(pidscopy[i]);
00404 cones[cone][conesizes[cone]++] = pidscopy[i];
00405 }
00406
00407 while(childcounter<nChild){
00408 for(int i=0;i<6;i++){
00409 if(conecounters[i]<conesizes[i]){
00410 subroots[childcounter++] = cones[i][conecounters[i]++];
00411 }
00412 }
00413 }
00414 for(int i=nChild;i<proxySpanDim;i++)
00415 subroots[i] = -1;
00416 return nChild;
00417 }
00418 #endif // USE_TOPOMAP
00419
00420 static int compDistance(const void *a, const void *b)
00421 {
00422 int d1 = abs(*(int *)a - CkMyPe());
00423 int d2 = abs(*(int *)b - CkMyPe());
00424 if (d1 < d2)
00425 return -1;
00426 else if (d1 == d2)
00427 return 0;
00428 else
00429 return 1;
00430 }
00431
00432 void HomePatch::sendProxies()
00433 {
00434 ProxyListIter pli(proxy);
00435 NodeIDList list;
00436 for ( pli = pli.begin(); pli != pli.end(); ++pli )
00437 {
00438 list.add(pli->node);
00439 }
00440 ProxyMgr::Object()->sendProxies(patchID, &list[0], list.size());
00441 }
00442
00443 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00444 void HomePatch::buildNodeAwareSpanningTree(void){
00445
00446 int *proxyNodeMap = new int[CkNumNodes()];
00447 NodeIDList proxyPeList;
00448 int proxyCnt = proxy.size();
00449 if(proxyCnt==0) {
00450
00451
00452
00453 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
00454 DebugFileTrace *dft = DebugFileTrace::Object();
00455 dft->openTrace();
00456 dft->writeTrace("HomePatch[%d] has 0 proxy on proc[%d] node[%d]\n", patchID, CkMyPe(), CkMyNode());
00457 dft->closeTrace();
00458 #endif
00459 return;
00460 }
00461 proxyPeList.resize(proxyCnt);
00462 ProxyListIter pli(proxy);
00463 int i=0;
00464 for ( pli = pli.begin(); pli != pli.end(); ++pli, i++ )
00465 proxyPeList[i] = pli->node;
00466 ProxyMgr::buildSinglePatchNodeAwareSpanningTree(patchID, proxyPeList, ptnTree, proxyNodeMap);
00467 delete [] proxyNodeMap;
00468 proxyPeList.resize(0);
00469
00470
00471
00472
00473 setupChildrenFromProxySpanningTree();
00474
00475 sendNodeAwareSpanningTree();
00476 }
00477
00478 void HomePatch::setupChildrenFromProxySpanningTree(){
00479 if(ptnTree.size()==0) {
00480 numChild = 0;
00481 delete [] children;
00482 children = NULL;
00483 #ifdef USE_NODEPATCHMGR
00484 numNodeChild = 0;
00485 delete [] nodeChildren;
00486 nodeChildren = NULL;
00487 #endif
00488 return;
00489 }
00490 proxyTreeNode *rootnode = &ptnTree.item(0);
00491 CmiAssert(rootnode->peIDs[0] == CkMyPe());
00492
00493
00494
00495 int internalChild;
00496 int externalChild;
00497 internalChild = rootnode->numPes-1;
00498 numChild = internalChild;
00499 if(numChild > inNodeProxySpanDim) {
00500
00501 CmiAbort("Enabling in-node spanning tree construction has not been implemented yet!\n");
00502 }else{
00503
00504 int treesize = ptnTree.size();
00505 externalChild = (proxySpanDim>(treesize-1))?(treesize-1):proxySpanDim;
00506 numChild += externalChild;
00507
00508 delete [] children;
00509 #ifdef USE_NODEPATCHMGR
00510 delete [] nodeChildren;
00511 #endif
00512 if(numChild==0){
00513 children = NULL;
00514 #ifdef USE_NODEPATCHMGR
00515 nodeChildren = NULL;
00516 numNodeChild = 0;
00517 #endif
00518 return;
00519 }
00520 children = new int[numChild];
00521 for(int i=0; i<externalChild; i++) {
00522 children[i] = ptnTree.item(i+1).peIDs[0];
00523 }
00524 for(int i=externalChild, j=1; i<numChild; i++, j++) {
00525 children[i] = rootnode->peIDs[j];
00526 }
00527 }
00528
00529 #ifdef USE_NODEPATCHMGR
00530
00531
00532 CProxy_NodeProxyMgr pm(CkpvAccess(BOCclass_group).nodeProxyMgr);
00533 NodeProxyMgr *npm = pm[CkMyNode()].ckLocalBranch();
00534 if(rootnode->numPes==1){
00535 npm->registerPatch(patchID, 0, NULL);
00536 }
00537 else{
00538 npm->registerPatch(patchID, rootnode->numPes-1, &rootnode->peIDs[1]);
00539 }
00540
00541
00542 numNodeChild = externalChild;
00543 if(internalChild) numNodeChild++;
00544 nodeChildren = new int[numNodeChild];
00545 for(int i=0; i<externalChild; i++) {
00546 nodeChildren[i] = ptnTree.item(i+1).nodeID;
00547 }
00548
00549
00550 if(internalChild)
00551 nodeChildren[numNodeChild-1] = rootnode->nodeID;
00552 #endif
00553
00554 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
00555 DebugFileTrace *dft = DebugFileTrace::Object();
00556 dft->openTrace();
00557 dft->writeTrace("HomePatch[%d] has %d children: ", patchID, numChild);
00558 for(int i=0; i<numChild; i++)
00559 dft->writeTrace("%d ", children[i]);
00560 dft->writeTrace("\n");
00561 dft->closeTrace();
00562 #endif
00563
00564 }
00565 #endif
00566
00567 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00568
00569 void HomePatch::recvNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *msg){
00570
00571 int treesize = msg->numNodesWithProxies;
00572 ptnTree.resize(treesize);
00573 int *pAllPes = msg->allPes;
00574 for(int i=0; i<treesize; i++) {
00575 proxyTreeNode *oneNode = &ptnTree.item(i);
00576 delete oneNode->peIDs;
00577 oneNode->numPes = msg->numPesOfNode[i];
00578 oneNode->nodeID = CkNodeOf(*pAllPes);
00579 oneNode->peIDs = new int[oneNode->numPes];
00580 for(int j=0; j<oneNode->numPes; j++) {
00581 oneNode->peIDs[j] = *pAllPes;
00582 pAllPes++;
00583 }
00584 }
00585
00586 setupChildrenFromProxySpanningTree();
00587
00588 sendNodeAwareSpanningTree();
00589 }
00590
00591 void HomePatch::sendNodeAwareSpanningTree(){
00592 if(ptnTree.size()==0) return;
00593 ProxyNodeAwareSpanningTreeMsg *msg =
00594 ProxyNodeAwareSpanningTreeMsg::getANewMsg(patchID, CkMyPe(), ptnTree.begin(), ptnTree.size());
00595
00596 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
00597 msg->printOut("HP::sendST");
00598 #endif
00599
00600 CmiAssert(CkMyPe() == msg->allPes[0]);
00601 ProxyMgr::Object()->sendNodeAwareSpanningTree(msg);
00602
00603 }
00604 #else
00605 void HomePatch::recvNodeAwareSpanningTree(ProxyNodeAwareSpanningTreeMsg *msg){}
00606 void HomePatch::sendNodeAwareSpanningTree(){}
00607 #endif
00608
00609 #ifndef NODEAWARE_PROXY_SPANNINGTREE
00610
00611 void HomePatch::recvSpanningTree(int *t, int n)
00612 {
00613 int i;
00614 nChild=0;
00615 tree.resize(n);
00616 for (i=0; i<n; i++) {
00617 tree[i] = t[i];
00618 }
00619
00620 for (i=1; i<=proxySpanDim; i++) {
00621 if (tree.size() <= i) break;
00622 child[i-1] = tree[i];
00623 nChild++;
00624 }
00625
00626
00627 sendSpanningTree();
00628 }
00629
00630 void HomePatch::sendSpanningTree()
00631 {
00632 if (tree.size() == 0) return;
00633 ProxySpanningTreeMsg *msg = new ProxySpanningTreeMsg;
00634 msg->patch = patchID;
00635 msg->node = CkMyPe();
00636 msg->tree = tree;
00637 ProxyMgr::Object()->sendSpanningTree(msg);
00638 }
00639 #else
00640 void HomePatch::recvSpanningTree(int *t, int n){}
00641 void HomePatch::sendSpanningTree(){}
00642 #endif
00643
00644 #ifndef NODEAWARE_PROXY_SPANNINGTREE
00645 void HomePatch::buildSpanningTree(void)
00646 {
00647 nChild = 0;
00648 int psize = proxy.size();
00649 if (psize == 0) return;
00650 NodeIDList oldtree = tree;
00651 int oldsize = oldtree.size();
00652 tree.resize(psize + 1);
00653 tree.setall(-1);
00654 tree[0] = CkMyPe();
00655 int s=1, e=psize+1;
00656 ProxyListIter pli(proxy);
00657 int patchNodesLast =
00658 ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
00659 int nNonPatch = 0;
00660 #if 0
00661 int* pelists = new int[psize];
00662 for (int i=0; i<psize; i++) pelists[i] = -1;
00663 for ( pli = pli.begin(); pli != pli.end(); ++pli ) {
00664 int idx = rand()%psize;
00665 while (pelists[idx] != -1) { idx++; if (idx == psize) idx = 0; }
00666 pelists[idx] = pli->node;
00667 }
00668 for ( int i=0; i<psize; i++)
00669 {
00670 if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pelists[i]) ) {
00671 tree[--e] = pelists[i];
00672 } else {
00673 tree[s++] = pelists[i];
00674 nNonPatch++;
00675 }
00676 }
00677 delete [] pelists;
00678 #else
00679
00680 for ( pli = pli.begin(); pli != pli.end(); ++pli )
00681 {
00682 int oldindex = oldtree.find(pli->node);
00683 if (oldindex != -1 && oldindex < psize) {
00684 tree[oldindex] = pli->node;
00685 }
00686 }
00687 s=1; e=psize;
00688 for ( pli = pli.begin(); pli != pli.end(); ++pli )
00689 {
00690 if (tree.find(pli->node) != -1) continue;
00691 if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pli->node) ) {
00692 while (tree[e] != -1) { e--; if (e==-1) e = psize; }
00693 tree[e] = pli->node;
00694 } else {
00695 while (tree[s] != -1) { s++; if (s==psize+1) s = 1; }
00696 tree[s] = pli->node;
00697 nNonPatch++;
00698 }
00699 }
00700 #if 1
00701 if (oldsize==0 && nNonPatch) {
00702
00703 qsort(&tree[1], nNonPatch, sizeof(int), compDistance);
00704 }
00705 #endif
00706
00707
00708
00709 #if USE_TOPOMAP && USE_SPANNING_TREE
00710
00711
00712 int *treecopy = new int [psize];
00713 int subroots[proxySpanDim];
00714 int subsizes[proxySpanDim];
00715 int subtrees[proxySpanDim][proxySpanDim];
00716 int idxes[proxySpanDim];
00717 int i = 0;
00718
00719 for(i=0;i<proxySpanDim;i++){
00720 subsizes[i] = 0;
00721 idxes[i] = i;
00722 }
00723
00724 for(i=0;i<psize;i++){
00725 treecopy[i] = tree[i+1];
00726 }
00727
00728 TopoManager tmgr;
00729 tmgr.sortRanksByHops(treecopy,nNonPatch);
00730 tmgr.sortRanksByHops(treecopy+nNonPatch,
00731 psize-nNonPatch);
00732
00733
00734 nChild = findSubroots(proxySpanDim,subroots,psize,treecopy);
00735 delete [] treecopy;
00736
00737 for(int i=1;i<psize+1;i++){
00738 int isSubroot=0;
00739 for(int j=0;j<nChild;j++)
00740 if(tree[i]==subroots[j]){
00741 isSubroot=1;
00742 break;
00743 }
00744 if(isSubroot) continue;
00745
00746 int bAdded = 0;
00747 tmgr.sortIndexByHops(tree[i], subroots,
00748 idxes, proxySpanDim);
00749 for(int j=0;j<proxySpanDim;j++){
00750 if(subsizes[idxes[j]]<proxySpanDim){
00751 subtrees[idxes[j]][(subsizes[idxes[j]])++] = tree[i];
00752 bAdded = 1;
00753 break;
00754 }
00755 }
00756 if( psize > proxySpanDim && ! bAdded ) {
00757 NAMD_bug("HomePatch BGL Spanning Tree error: Couldn't find subtree for leaf\n");
00758 }
00759 }
00760
00761 #else
00762
00763 for (int i=1; i<=proxySpanDim; i++) {
00764 if (tree.size() <= i) break;
00765 child[i-1] = tree[i];
00766 nChild++;
00767 }
00768 #endif
00769 #endif
00770
00771 #if 0
00772
00773 CkPrintf("[%d] Spanning tree for %d with %d children %d nNonPatch %d\n", CkMyPe(), patchID, psize, nNonPatch);
00774 for (int i=0; i<psize+1; i++) {
00775 CkPrintf("%d ", tree[i]);
00776 }
00777 CkPrintf("\n");
00778 #endif
00779
00780 sendSpanningTree();
00781 }
00782 #endif
00783
00784
00785 void HomePatch::receiveResults(ProxyResultVarsizeMsg *msg){
00786 DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00787 int n = msg->node;
00788 ProxyListElem *pe = proxy.begin();
00789 for ( ; pe->node != n; ++pe );
00790 Results *r = pe->forceBox->open();
00791
00792 char *iszeroPtr = msg->isZero;
00793 Force *msgFPtr = msg->forceArr;
00794
00795 for ( int k = 0; k < Results::maxNumForces; ++k )
00796 {
00797 Force *rfPtr = r->f[k];
00798 for(int i=0; i<msg->flLen[k]; i++, rfPtr++, iszeroPtr++) {
00799 if((*iszeroPtr)!=1) {
00800 *rfPtr += *msgFPtr;
00801 msgFPtr++;
00802 }
00803 }
00804 }
00805 pe->forceBox->close(&r);
00806 delete msg;
00807 }
00808
00809 void HomePatch::receiveResults(ProxyResultMsg *msg)
00810 {
00811 DebugM(4, "patchID("<<patchID<<") receiveRes() nodeID("<<msg->node<<")\n");
00812 int n = msg->node;
00813 ProxyListElem *pe = proxy.begin();
00814 for ( ; pe->node != n; ++pe );
00815 Results *r = pe->forceBox->open();
00816 for ( int k = 0; k < Results::maxNumForces; ++k )
00817 {
00818 Force *f = r->f[k];
00819 register ForceList::iterator f_i, f_e;
00820 f_i = msg->forceList[k].begin();
00821 f_e = msg->forceList[k].end();
00822 for ( ; f_i != f_e; ++f_i, ++f ) *f += *f_i;
00823 }
00824 pe->forceBox->close(&r);
00825 delete msg;
00826 }
00827
00828 void HomePatch::receiveResults(ProxyCombinedResultMsg *msg)
00829 {
00830 DebugM(4, "patchID("<<patchID<<") receiveRes() #nodes("<<msg->nodes.size()<<")\n");
00831
00832 for (int i=0; i<msg->nodes.size(); i++) {
00833 int node = msg->nodes[i];
00834 ProxyListElem *pe = proxy.begin();
00835 for ( ; pe->node != node; ++pe );
00836 Results *r = pe->forceBox->open();
00837 if (i==0) {
00838 for ( int k = 0; k < Results::maxNumForces; ++k )
00839 {
00840 Force *f = r->f[k];
00841 register ForceList::iterator f_i, f_e;
00842 f_i = msg->forceList[k].begin();
00843 f_e = msg->forceList[k].end();
00844
00845
00846 int nf = f_e - f_i;
00847 #ifdef ARCH_POWERPC
00848 #pragma disjoint (*f_i, *f)
00849 #pragma unroll(4)
00850 #endif
00851 for (int count = 0; count < nf; count++) {
00852 f[count].x += f_i[count].x;
00853 f[count].y += f_i[count].y;
00854 f[count].z += f_i[count].z;
00855 }
00856 }
00857 }
00858 pe->forceBox->close(&r);
00859 }
00860
00861 delete msg;
00862 }
00863
00864 void HomePatch::positionsReady(int doMigration)
00865 {
00866 flags.sequence++;
00867
00868 if (!patchMapRead) { readPatchMap(); }
00869
00870 if (numNeighbors) {
00871 if (doMigration) {
00872 doAtomMigration();
00873 } else {
00874 doMarginCheck();
00875 }
00876 }
00877 doMigration = (doMigration && numNeighbors) || ! patchMapRead;
00878
00879
00880 doGroupSizeCheck();
00881
00882
00883 p.resize(numAtoms);
00884 CompAtom *p_i = p.begin();
00885 pExt.resize(numAtoms);
00886 CompAtomExt *pExt_i = pExt.begin();
00887 FullAtom *a_i = atom.begin();
00888 int i; int n = numAtoms;
00889 for ( i=0; i<n; ++i ) {
00890 p_i[i] = a_i[i];
00891 pExt_i[i] = a_i[i];
00892 }
00893
00894
00895 doPairlistCheck();
00896
00897 if (flags.doMolly) mollyAverage();
00898
00899
00900 ProxyListIter pli(proxy);
00901 int *pids = NULL;
00902 int npid;
00903 if (proxySendSpanning == 0) {
00904 npid = proxy.size();
00905 pids = new int[npid];
00906 int *pidi = pids;
00907 int *pide = pids + proxy.size();
00908 int patchNodesLast =
00909 ( PatchMap::Object()->numNodesWithPatches() < ( 0.7 * CkNumPes() ) );
00910 for ( pli = pli.begin(); pli != pli.end(); ++pli )
00911 {
00912 if ( patchNodesLast && PatchMap::Object()->numPatchesOnNode(pli->node) ) {
00913 *(--pide) = pli->node;
00914 } else {
00915 *(pidi++) = pli->node;
00916 }
00917 }
00918 }
00919 else {
00920 #ifdef NODEAWARE_PROXY_SPANNINGTREE
00921 #ifdef USE_NODEPATCHMGR
00922 npid = numNodeChild;
00923 if(numNodeChild>0) {
00924 pids = new int[npid];
00925 memcpy(pids, nodeChildren, sizeof(int)*npid);
00926 }
00927 #else
00928 npid = numChild;
00929 if(numChild>0) {
00930 pids = new int[numChild];
00931 memcpy(pids, children, sizeof(int)*numChild);
00932 }
00933 #endif
00934 #else
00935 npid = nChild;
00936 pids = new int[proxySpanDim];
00937 for (int i=0; i<nChild; i++) pids[i] = child[i];
00938 #endif
00939 }
00940 if (npid) {
00941 #if CMK_PERSISTENT_COMM
00942 if (phsReady == 0)
00943 {
00944
00945 for (int i=0; i<npid; i++) {
00946 localphs[i] = CmiCreatePersistent(pids[i], 30000);
00947 }
00948 nphs = npid;
00949 phsReady = 1;
00950 }
00951 #endif
00952 int seq = flags.sequence;
00953 int priority = PROXY_DATA_PRIORITY + PATCH_PRIORITY(patchID);
00954
00955 int pdMsgPLLen = p.size();
00956 int pdMsgAvgPLLen = 0;
00957 if(flags.doMolly) {
00958 pdMsgAvgPLLen = p_avg.size();
00959 }
00960 int pdMsgPLExtLen = 0;
00961 if(doMigration || isNewProxyAdded) {
00962 pdMsgPLExtLen = pExt.size();
00963 }
00964 ProxyDataMsg *nmsg = new (pdMsgPLLen, pdMsgAvgPLLen, pdMsgPLExtLen, PRIORITY_SIZE) ProxyDataMsg;
00965 SET_PRIORITY(nmsg,seq,priority);
00966 nmsg->patch = patchID;
00967 nmsg->flags = flags;
00968 nmsg->plLen = pdMsgPLLen;
00969
00970 memcpy(nmsg->positionList, p.begin(), sizeof(CompAtom)*pdMsgPLLen);
00971 nmsg->avgPlLen = pdMsgAvgPLLen;
00972 if(flags.doMolly) {
00973 memcpy(nmsg->avgPositionList, p_avg.begin(), sizeof(CompAtom)*pdMsgAvgPLLen);
00974 }
00975 nmsg->plExtLen = pdMsgPLExtLen;
00976 if(doMigration || isNewProxyAdded){
00977 memcpy(nmsg->positionExtList, pExt.begin(), sizeof(CompAtomExt)*pdMsgPLExtLen);
00978 }
00979
00980 #if NAMD_SeparateWaters != 0
00981
00982 nmsg->numWaterAtoms = numWaterAtoms;
00983 #endif
00984
00985 #if defined(NODEAWARE_PROXY_SPANNINGTREE) && defined(USE_NODEPATCHMGR) && (CMK_SMP) && defined(NAMDSRC_IMMQD_HACK)
00986 nmsg->isFromImmMsgCall = 0;
00987 #endif
00988
00989 #if defined(PROCTRACE_DEBUG) && defined(NAST_DEBUG)
00990 DebugFileTrace *dft = DebugFileTrace::Object();
00991 dft->openTrace();
00992 dft->writeTrace("HP::posReady: for HomePatch[%d], sending proxy msg to: ", patchID);
00993 for(int i=0; i<npid; i++) {
00994 dft->writeTrace("%d ", pids[i]);
00995 }
00996 dft->writeTrace("\n");
00997 dft->closeTrace();
00998 #endif
00999
01000 if(doMigration) {
01001 ProxyMgr::Object()->sendProxyAll(nmsg,npid,pids);
01002 }else{
01003 ProxyMgr::Object()->sendProxyData(nmsg,npid,pids);
01004 }
01005 #if CMK_PERSISTENT_COMM
01006 CmiUsePersistentHandle(NULL, 0);
01007 #endif
01008 isNewProxyAdded = 0;
01009 }
01010 delete [] pids;
01011 DebugM(4, "patchID("<<patchID<<") doing positions Ready\n");
01012
01013 #ifdef REMOVE_PROXYDATAMSG_EXTRACOPY
01014 positionPtrBegin = p.begin();
01015 positionPtrEnd = p.end();
01016 #endif
01017
01018 if(flags.doMolly) {
01019 avgPositionPtrBegin = p_avg.begin();
01020 avgPositionPtrEnd = p_avg.end();
01021 }
01022 Patch::positionsReady(doMigration);
01023
01024 patchMapRead = 1;
01025
01026
01027 if (useSync) Sync::Object()->PatchReady();
01028 }
01029
01030 void HomePatch::replaceForces(ExtForce *f)
01031 {
01032 replacementForces = f;
01033 }
01034
01035 void HomePatch::saveForce(const int ftag)
01036 {
01037 f_saved[ftag].resize(numAtoms);
01038 for ( int i = 0; i < numAtoms; ++i )
01039 {
01040 f_saved[ftag][i] = f[ftag][i];
01041 }
01042 }
01043
01044 #undef DEBUG_REDISTRIB_FORCE
01045 #undef DEBUG_REDISTRIB_FORCE_VERBOSE
01046
01047
01048
01049
01050
01051
01052
01053
01054
01055 void HomePatch::redistrib_lp_force(
01056 Vector& f_ox, Vector& f_h1, Vector& f_h2, Vector& f_lp,
01057 const Vector& p_ox, const Vector& p_h1, const Vector& p_h2,
01058 const Vector& p_lp, Tensor *virial) {
01059
01060 Tensor wc;
01061
01062 #ifdef DEBUG_REDISTRIB_FORCE
01063
01064
01065
01066 Vector totforce(0.0, 0.0, 0.0);
01067 Vector tottorque(0.0, 0.0, 0.0);
01068
01069 totforce = f_ox + f_h1 + f_h2 + f_lp;
01070
01071 tottorque += cross(f_ox, p_ox);
01072 tottorque += cross(f_h1, p_h1);
01073 tottorque += cross(f_h2, p_h2);
01074
01075
01076 tottorque += cross(f_lp, p_lp);
01077
01078
01079 #endif
01080
01081
01082 Vector r_ox_lp = p_lp - p_ox;
01083 BigReal rad_factor = (f_lp * r_ox_lp) * r_ox_lp.rlength() * r_ox_lp.rlength();
01084 Vector f_rad = r_ox_lp * rad_factor;
01085
01086 Tensor vir = outer(f_rad, p_ox);
01087 wc += vir;
01088
01089 f_ox = f_ox + f_rad;
01090
01091
01092 Vector r_hcom_ox = p_ox - ( (p_h1 + p_h2) * 0.5 );
01093 Vector r_h2_h1_2 = (p_h1 - p_h2) * 0.5;
01094
01095
01096 Vector r_oop = cross(r_ox_lp, r_hcom_ox);
01097
01098
01099 Vector r_perp = cross(r_hcom_ox, r_h2_h1_2);
01100
01101
01102
01103
01104
01105
01106 Vector f_ang = f_lp - f_rad;
01107
01108
01109 BigReal oxcomp = (r_hcom_ox.length() - r_ox_lp.length()) *
01110 r_hcom_ox.rlength();
01111 BigReal hydcomp = 0.5 * r_ox_lp.length() * r_hcom_ox.rlength();
01112
01113 f_ox = f_ox + (f_ang * oxcomp);
01114 f_h1 = f_h1 + (f_ang * hydcomp);
01115 f_h2 = f_h2 + (f_ang * hydcomp);
01116
01117
01118 vir = outer(f_ang * oxcomp, p_ox);
01119 wc += vir;
01120 vir = outer(f_ang * hydcomp, p_h1);
01121 wc += vir;
01122 vir = outer(f_ang * hydcomp, p_h2);
01123 wc += vir;
01124 vir = outer(-1.0 * f_lp, p_lp);
01125 wc += vir;
01126
01127 if ( virial ) *virial += wc;
01128
01129
01130 f_lp = Vector(0.0, 0.0, 0.0);
01131
01132 #ifdef DEBUG_REDISTRIB_FORCE
01133
01134 Vector newforce(0.0, 0.0, 0.0);
01135 Vector newtorque(0.0, 0.0, 0.0);
01136 BigReal error = 0.0;
01137
01138 newforce = f_ox + f_h1 + f_h2;
01139
01140 newtorque += cross(f_ox, p_ox);
01141 newtorque += cross(f_h1, p_h1);
01142 newtorque += cross(f_h2, p_h2);
01143
01144 error = fabs(newforce.length() - totforce.length());
01145 if (error > 0.0001) {
01146 printf("Error: Force redistribution for water "
01147 "exceeded force tolerance (%f vs. %f)\n",
01148 newforce.length(), totforce.length());
01149 }
01150 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
01151 printf("Error in force length: %f\n", error);
01152 #endif
01153
01154 error = fabs(newtorque.length() - tottorque.length());
01155 if (error > 0.0001) {
01156 printf("Error: Force redistribution for water "
01157 "exceeded torque tolerance (%f vs. %f)\n",
01158 newtorque.length(), tottorque.length());
01159 }
01160 #ifdef DEBUG_REDISTRIB_FORCE_VERBOSE
01161 printf("Error in torque: %f\n", error);
01162 #endif
01163 #endif
01164 }
01165
01166 void HomePatch::swm4_omrepos(Vector *ref, Vector *pos, Vector *vel,
01167 BigReal invdt) {
01168
01169
01170
01171 pos[2] = pos[0] + (0.5 * (pos[3] + pos[4]) - pos[0]) * (r_om / r_ohc);
01172
01173 if (invdt != 0) {
01174 vel[2] = (pos[2] - ref[2]) * invdt;
01175 }
01176
01177 }
01178
01179 void HomePatch::tip4_omrepos(Vector* ref, Vector* pos, Vector* vel, BigReal invdt) {
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192 pos[3] = pos[0] + (0.5 * (pos[1] + pos[2]) - pos[0]) * (r_om / r_ohc);
01193
01194
01195
01196 if (invdt != 0) {
01197 vel[3] = (pos[3] - ref[3]) * invdt;
01198 }
01199
01200
01201 return;
01202 }
01203
01204 void HomePatch::redistrib_swm4_forces(const int ftag, Tensor *virial) {
01205
01206
01207 ForceList *f_mod = f;
01208 for (int i = 0; i < numAtoms; i++) {
01209 if (atom[i].mass < 0.01) {
01210
01211 redistrib_lp_force(f_mod[ftag][i-2], f_mod[ftag][i+1],
01212 f_mod[ftag][i+2], f_mod[ftag][i],
01213 atom[i-2].position, atom[i+1].position,
01214 atom[i+2].position, atom[i].position, virial);
01215 }
01216 }
01217 }
01218
01219 void HomePatch::redistrib_tip4p_forces(const int ftag, Tensor* virial) {
01220
01221
01222
01223
01224 ForceList *f_mod =f;
01225 for (int i=0; i<numAtoms; i++) {
01226 if (atom[i].mass < 0.01) {
01227
01228 redistrib_lp_force(f_mod[ftag][i-3], f_mod[ftag][i-2],
01229 f_mod[ftag][i-1], f_mod[ftag][i],
01230 atom[i-3].position, atom[i-2].position,
01231 atom[i-1].position, atom[i].position, virial);
01232 }
01233 }
01234 }
01235
01236
01237 void HomePatch::addForceToMomentum(const BigReal timestep, const int ftag,
01238 const int useSaved)
01239 {
01240 SimParameters *simParams = Node::Object()->simParameters;
01241 const BigReal dt = timestep / TIMEFACTOR;
01242 ForceList *f_use = (useSaved ? f_saved : f);
01243
01244 if ( simParams->fixedAtomsOn ) {
01245 for ( int i = 0; i < numAtoms; ++i ) {
01246 if ( atom[i].atomFixed ) {
01247 atom[i].velocity = 0;
01248 } else {
01249 BigReal recip_val = ( atom[i].mass > 0. ? dt * namd_reciprocal( atom[i].mass ) : 0.);
01250 atom[i].velocity += f_use[ftag][i] * recip_val;
01251 }
01252 }
01253 } else {
01254 FullAtom *atom_arr = atom.begin();
01255 const Force *force_arr = f_use[ftag].const_begin();
01256 #ifdef ARCH_POWERPC
01257 #pragma disjoint (*force_arr, *atom_arr)
01258 #endif
01259 for ( int i = 0; i < numAtoms; ++i ) {
01260 if (atom[i].mass == 0.) continue;
01261 BigReal recip_val = ( atom[i].mass > 0. ? dt * namd_reciprocal( atom[i].mass ) : 0.);
01262
01263 atom_arr[i].velocity.x += force_arr[i].x * recip_val;
01264 atom_arr[i].velocity.y += force_arr[i].y * recip_val;
01265 atom_arr[i].velocity.z += force_arr[i].z * recip_val;
01266 }
01267 }
01268 }
01269
01270 void HomePatch::addVelocityToPosition(const BigReal timestep)
01271 {
01272 SimParameters *simParams = Node::Object()->simParameters;
01273 const BigReal dt = timestep / TIMEFACTOR;
01274 if ( simParams->fixedAtomsOn ) {
01275 for ( int i = 0; i < numAtoms; ++i ) {
01276 if ( ! atom[i].atomFixed ) atom[i].position += atom[i].velocity * dt;
01277 }
01278 } else {
01279 FullAtom *atom_arr = atom.begin();
01280 for ( int i = 0; i < numAtoms; ++i ) {
01281 atom_arr[i].position.x += atom_arr[i].velocity.x * dt;
01282 atom_arr[i].position.y += atom_arr[i].velocity.y * dt;
01283 atom_arr[i].position.z += atom_arr[i].velocity.z * dt;
01284 }
01285 }
01286 }
01287
01288
01289 int HomePatch::rattle1(const BigReal timestep, Tensor *virial,
01290 SubmitReduction *ppreduction)
01291 {
01292 Molecule *mol = Node::Object()->molecule;
01293 SimParameters *simParams = Node::Object()->simParameters;
01294 const int fixedAtomsOn = simParams->fixedAtomsOn;
01295 const int useSettle = simParams->useSettle;
01296 const BigReal dt = timestep / TIMEFACTOR;
01297 const BigReal invdt = (dt == 0.) ? 0. : 1.0 / dt;
01298 BigReal tol2 = 2.0 * simParams->rigidTol;
01299 int maxiter = simParams->rigidIter;
01300 int dieOnError = simParams->rigidDie;
01301 int i, iter;
01302 BigReal dsq[10], tmp;
01303 int ial[10], ibl[10];
01304 Vector ref[10];
01305 Vector refab[10];
01306 Vector pos[10];
01307 Vector vel[10];
01308 Vector netdp[10];
01309 BigReal rmass[10];
01310 int fixed[10];
01311 Tensor wc;
01312 BigReal idz, zmin;
01313 int nslabs;
01314
01315
01316
01317
01318
01319
01320 if (!settle1isinitted()) {
01321 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01322
01323 if (mol->rigid_bond_length(atom[ig].id) > 0) {
01324 int oatm;
01325 if (simParams->watmodel == WAT_SWM4) {
01326 oatm = ig+3;
01327
01328
01329 }
01330 else {
01331 oatm = ig+1;
01332
01333 if (atom[ig].mass < 0.5 || atom[ig+1].mass < 0.5) {
01334 oatm += 1;
01335 }
01336 }
01337
01338
01339 settle1init(atom[ig].mass, atom[oatm].mass,
01340 mol->rigid_bond_length(atom[ig].id),
01341 mol->rigid_bond_length(atom[oatm].id));
01342 break;
01343 }
01344 }
01345 }
01346
01347 if (ppreduction) {
01348 nslabs = simParams->pressureProfileSlabs;
01349 idz = nslabs/lattice.c().z;
01350 zmin = lattice.origin().z - 0.5*lattice.c().z;
01351 }
01352
01353
01354 int wathgsize = 3;
01355 int watmodel = simParams->watmodel;
01356 if (watmodel == WAT_TIP4) wathgsize = 4;
01357 else if (watmodel == WAT_SWM4) wathgsize = 5;
01358
01359 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01360 int hgs = atom[ig].hydrogenGroupSize;
01361 if ( hgs == 1 ) continue;
01362
01363 for ( i = 0; i < hgs; ++i ) {
01364 ref[i] = atom[ig+i].position;
01365 pos[i] = atom[ig+i].position;
01366 vel[i] = atom[ig+i].velocity;
01367 rmass[i] = (atom[ig+1].mass > 0. ? 1. / atom[ig+i].mass : 0.);
01368
01369 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
01370
01371
01372 if ( fixed[i] ) rmass[i] = 0.; else pos[i] += vel[i] * dt;
01373 }
01374 int icnt = 0;
01375 if ( ( tmp = mol->rigid_bond_length(atom[ig].id) ) > 0 ) {
01376 if (hgs != wathgsize) {
01377 NAMD_bug("Hydrogen group error caught in rattle1().");
01378 }
01379
01380
01381 if (useSettle && !atom[ig].groupFixed) {
01382 if (simParams->watmodel == WAT_SWM4) {
01383
01384
01385
01386 Vector lp_ref = ref[2];
01387 Vector lp_pos = pos[2];
01388 Vector lp_vel = vel[2];
01389 ref[2] = ref[0];
01390 pos[2] = pos[0];
01391 vel[2] = vel[0];
01392 settle1(ref+2, pos+2, vel+2, invdt);
01393 ref[0] = ref[2];
01394 pos[0] = pos[2];
01395 vel[0] = vel[2];
01396 ref[2] = lp_ref;
01397 pos[2] = lp_pos;
01398 vel[2] = lp_vel;
01399
01400 swm4_omrepos(ref, pos, vel, invdt);
01401 }
01402 else {
01403 settle1(ref, pos, vel, invdt);
01404 if (simParams->watmodel == WAT_TIP4) {
01405 tip4_omrepos(ref, pos, vel, invdt);
01406 }
01407 }
01408
01409
01410
01411 int ppoffset, partition;
01412 if ( invdt == 0 ) for ( i = 0; i < wathgsize; ++i ) {
01413 atom[ig+i].position = pos[i];
01414 } else if ( virial == 0 ) for ( i = 0; i < wathgsize; ++i ) {
01415 atom[ig+i].velocity = vel[i];
01416 } else for ( i = 0; i < wathgsize; ++i ) {
01417 Force df = (vel[i] - atom[ig+i].velocity) * ( atom[ig+i].mass * invdt );
01418 Tensor vir = outer(df, ref[i]);
01419 wc += vir;
01420 f[Results::normal][ig+i] += df;
01421 atom[ig+i].velocity = vel[i];
01422 if (ppreduction) {
01423
01424
01425 if (!i) {
01426 BigReal z = pos[i].z;
01427 partition = atom[ig].partition;
01428 int slab = (int)floor((z-zmin)*idz);
01429 if (slab < 0) slab += nslabs;
01430 else if (slab >= nslabs) slab -= nslabs;
01431 ppoffset = 3*(slab + nslabs*partition);
01432 }
01433 ppreduction->item(ppoffset ) += vir.xx;
01434 ppreduction->item(ppoffset+1) += vir.yy;
01435 ppreduction->item(ppoffset+2) += vir.zz;
01436 }
01437 }
01438 continue;
01439 }
01440 if ( !(fixed[1] && fixed[2]) ) {
01441 dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
01442 }
01443 }
01444 for ( i = 1; i < hgs; ++i ) {
01445 if ( ( tmp = mol->rigid_bond_length(atom[ig+i].id) ) > 0 ) {
01446 if ( !(fixed[0] && fixed[i]) ) {
01447 dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
01448 }
01449 }
01450 }
01451 if ( icnt == 0 ) continue;
01452 for ( i = 0; i < icnt; ++i ) {
01453 refab[i] = ref[ial[i]] - ref[ibl[i]];
01454 }
01455 for ( i = 0; i < hgs; ++i ) {
01456 netdp[i] = 0.;
01457 }
01458 int done;
01459 int consFailure;
01460 for ( iter = 0; iter < maxiter; ++iter ) {
01461
01462 done = 1;
01463 consFailure = 0;
01464 for ( i = 0; i < icnt; ++i ) {
01465 int a = ial[i]; int b = ibl[i];
01466 Vector pab = pos[a] - pos[b];
01467 BigReal pabsq = pab.x*pab.x + pab.y*pab.y + pab.z*pab.z;
01468 BigReal rabsq = dsq[i];
01469 BigReal diffsq = rabsq - pabsq;
01470 if ( fabs(diffsq) > (rabsq * tol2) ) {
01471 Vector &rab = refab[i];
01472 BigReal rpab = rab.x*pab.x + rab.y*pab.y + rab.z*pab.z;
01473 if ( rpab < ( rabsq * 1.0e-6 ) ) {
01474 done = 0;
01475 consFailure = 1;
01476 continue;
01477 }
01478 BigReal rma = rmass[a];
01479 BigReal rmb = rmass[b];
01480 BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
01481 Vector dp = rab * gab;
01482 pos[a] += rma * dp;
01483 pos[b] -= rmb * dp;
01484 if ( invdt != 0. ) {
01485 dp *= invdt;
01486 netdp[a] += dp;
01487 netdp[b] -= dp;
01488 }
01489 done = 0;
01490 }
01491 }
01492 if ( done ) break;
01493 }
01494
01495 if ( consFailure ) {
01496 if ( dieOnError ) {
01497 iout << iERROR << "Constraint failure in RATTLE algorithm for atom "
01498 << (atom[ig].id + 1) << "!\n" << endi;
01499 return -1;
01500 } else {
01501 iout << iWARN << "Constraint failure in RATTLE algorithm for atom "
01502 << (atom[ig].id + 1) << "!\n" << endi;
01503 }
01504 } else if ( ! done ) {
01505 if ( dieOnError ) {
01506 iout << iERROR << "Exceeded RATTLE iteration limit for atom "
01507 << (atom[ig].id + 1) << "!\n" << endi;
01508 return -1;
01509 } else {
01510 iout << iWARN << "Exceeded RATTLE iteration limit for atom "
01511 << (atom[ig].id + 1) << "!\n" << endi;
01512 }
01513 }
01514
01515
01516 int ppoffset, partition;
01517 if ( invdt == 0 ) for ( i = 0; i < hgs; ++i ) {
01518 atom[ig+i].position = pos[i];
01519 } else if ( virial == 0 ) for ( i = 0; i < hgs; ++i ) {
01520 atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
01521 } else for ( i = 0; i < hgs; ++i ) {
01522 Force df = netdp[i] * invdt;
01523 Tensor vir = outer(df, ref[i]);
01524 wc += vir;
01525 f[Results::normal][ig+i] += df;
01526 atom[ig+i].velocity = vel[i] + rmass[i] * netdp[i];
01527 if (ppreduction) {
01528 if (!i) {
01529 BigReal z = pos[i].z;
01530 int partition = atom[ig].partition;
01531 int slab = (int)floor((z-zmin)*idz);
01532 if (slab < 0) slab += nslabs;
01533 else if (slab >= nslabs) slab -= nslabs;
01534 ppoffset = 3*(slab + nslabs*partition);
01535 }
01536 ppreduction->item(ppoffset ) += vir.xx;
01537 ppreduction->item(ppoffset+1) += vir.yy;
01538 ppreduction->item(ppoffset+2) += vir.zz;
01539 }
01540 }
01541 }
01542 if ( dt && virial ) *virial += wc;
01543
01544 return 0;
01545 }
01546
01547
01548 void HomePatch::rattle2(const BigReal timestep, Tensor *virial)
01549 {
01550 Molecule *mol = Node::Object()->molecule;
01551 SimParameters *simParams = Node::Object()->simParameters;
01552 const int fixedAtomsOn = simParams->fixedAtomsOn;
01553 const int useSettle = simParams->useSettle;
01554 const BigReal dt = timestep / TIMEFACTOR;
01555 Tensor wc;
01556 BigReal tol = simParams->rigidTol;
01557 int maxiter = simParams->rigidIter;
01558 int dieOnError = simParams->rigidDie;
01559 int i, iter;
01560 BigReal dsqi[10], tmp;
01561 int ial[10], ibl[10];
01562 Vector ref[10];
01563 Vector refab[10];
01564 Vector vel[10];
01565 BigReal rmass[10];
01566 BigReal redmass[10];
01567 int fixed[10];
01568
01569
01570 int wathgsize = 3;
01571 if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
01572 else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
01573
01574
01575 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01576
01577 int hgs = atom[ig].hydrogenGroupSize;
01578 if ( hgs == 1 ) continue;
01579
01580 for ( i = 0; i < hgs; ++i ) {
01581 ref[i] = atom[ig+i].position;
01582 vel[i] = atom[ig+i].velocity;
01583 rmass[i] = atom[ig+i].mass > 0. ? 1. / atom[ig+i].mass : 0.;
01584 fixed[i] = ( fixedAtomsOn && atom[ig+i].atomFixed );
01585 if ( fixed[i] ) rmass[i] = 0.;
01586 }
01587 int icnt = 0;
01588 if ( ( tmp = mol->rigid_bond_length(atom[ig].id) ) > 0 ) {
01589 if ( wathgsize != 4 ) {
01590 NAMD_bug("Hydrogen group error caught in rattle2().");
01591 }
01592
01593 if (useSettle && !fixed[0] && !fixed[1] && !fixed[2]) {
01594 settle2(atom[ig].mass, atom[ig+1].mass, ref, vel, dt, virial);
01595 for (i=0; i<3; i++) {
01596 atom[ig+i].velocity = vel[i];
01597 }
01598 continue;
01599 }
01600 if ( !(fixed[1] && fixed[2]) ) {
01601 redmass[icnt] = 1. / (rmass[1] + rmass[2]);
01602 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
01603 }
01604 }
01605
01606 for ( i = 1; i < hgs; ++i ) {
01607 if ( ( tmp = mol->rigid_bond_length(atom[ig+i].id) ) > 0 ) {
01608 if ( !(fixed[0] && fixed[i]) ) {
01609 redmass[icnt] = 1. / (rmass[0] + rmass[i]);
01610 dsqi[icnt] = 1. / (tmp * tmp); ial[icnt] = 0;
01611 ibl[icnt] = i; ++icnt;
01612 }
01613 }
01614 }
01615 if ( icnt == 0 ) continue;
01616
01617 for ( i = 0; i < icnt; ++i ) {
01618 refab[i] = ref[ial[i]] - ref[ibl[i]];
01619 }
01620
01621 int done;
01622 for ( iter = 0; iter < maxiter; ++iter ) {
01623 done = 1;
01624 for ( i = 0; i < icnt; ++i ) {
01625 int a = ial[i]; int b = ibl[i];
01626 Vector vab = vel[a] - vel[b];
01627 Vector &rab = refab[i];
01628 BigReal rabsqi = dsqi[i];
01629 BigReal rvab = rab.x*vab.x + rab.y*vab.y + rab.z*vab.z;
01630 if ( (fabs(rvab) * dt * rabsqi) > tol ) {
01631 Vector dp = rab * (-rvab * redmass[i] * rabsqi);
01632 wc += outer(dp,rab);
01633 vel[a] += rmass[a] * dp;
01634 vel[b] -= rmass[b] * dp;
01635 done = 0;
01636 }
01637 }
01638 if ( done ) break;
01639
01640 }
01641 if ( ! done ) {
01642 if ( dieOnError ) {
01643 NAMD_die("Exceeded maximum number of iterations in rattle2().");
01644 } else {
01645 iout << iWARN <<
01646 "Exceeded maximum number of iterations in rattle2().\n" << endi;
01647 }
01648 }
01649
01650 for ( i = 0; i < hgs; ++i ) {
01651 atom[ig+i].velocity = vel[i];
01652 }
01653 }
01654
01655
01656 *virial += wc / ( 0.5 * dt );
01657
01658 }
01659
01660
01661
01662 void HomePatch::mollyAverage()
01663 {
01664 Molecule *mol = Node::Object()->molecule;
01665 SimParameters *simParams = Node::Object()->simParameters;
01666 BigReal tol = simParams->mollyTol;
01667 int maxiter = simParams->mollyIter;
01668 int i, iter;
01669 HGArrayBigReal dsq;
01670 BigReal tmp;
01671 HGArrayInt ial, ibl;
01672 HGArrayVector ref;
01673 HGArrayVector refab;
01674 HGArrayBigReal rmass;
01675 BigReal *lambda;
01676 CompAtom *avg;
01677 int numLambdas = 0;
01678 HGArrayInt fixed;
01679
01680
01681 p_avg.resize(numAtoms);
01682 for ( i=0; i<numAtoms; ++i ) p_avg[i] = p[i];
01683
01684 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01685 int hgs = atom[ig].hydrogenGroupSize;
01686 if ( hgs == 1 ) continue;
01687 for ( i = 0; i < hgs; ++i ) {
01688 ref[i] = atom[ig+i].position;
01689 rmass[i] = 1. / atom[ig+i].mass;
01690 fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
01691 if ( fixed[i] ) rmass[i] = 0.;
01692 }
01693 avg = &(p_avg[ig]);
01694 int icnt = 0;
01695
01696 if ( ( tmp = mol->rigid_bond_length(atom[ig].id) ) ) {
01697 if ( hgs != 3 ) {
01698 NAMD_die("Hydrogen group error caught in mollyAverage(). It's a bug!\n");
01699 }
01700 if ( !(fixed[1] && fixed[2]) ) {
01701 dsq[icnt] = tmp * tmp; ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
01702 }
01703 }
01704 for ( i = 1; i < hgs; ++i ) {
01705 if ( ( tmp = mol->rigid_bond_length(atom[ig+i].id) ) ) {
01706 if ( !(fixed[0] && fixed[i]) ) {
01707 dsq[icnt] = tmp * tmp; ial[icnt] = 0; ibl[icnt] = i; ++icnt;
01708 }
01709 }
01710 }
01711 if ( icnt == 0 ) continue;
01712 numLambdas += icnt;
01713 molly_lambda.resize(numLambdas);
01714 lambda = &(molly_lambda[numLambdas - icnt]);
01715 for ( i = 0; i < icnt; ++i ) {
01716 refab[i] = ref[ial[i]] - ref[ibl[i]];
01717 }
01718
01719 iter=average(avg,ref,lambda,hgs,icnt,rmass,dsq,ial,ibl,refab,tol,maxiter);
01720 if ( iter == maxiter ) {
01721 iout << iWARN << "Exceeded maximum number of iterations in mollyAverage().\n"<<endi;
01722 }
01723 }
01724
01725
01726
01727
01728
01729
01730
01731
01732 }
01733
01734
01735
01736 void HomePatch::mollyMollify(Tensor *virial)
01737 {
01738 Molecule *mol = Node::Object()->molecule;
01739 SimParameters *simParams = Node::Object()->simParameters;
01740 Tensor wc;
01741 int i;
01742 HGArrayInt ial, ibl;
01743 HGArrayVector ref;
01744 CompAtom *avg;
01745 HGArrayVector refab;
01746 HGArrayVector force;
01747 HGArrayBigReal rmass;
01748 BigReal *lambda;
01749 int numLambdas = 0;
01750 HGArrayInt fixed;
01751
01752 for ( int ig = 0; ig < numAtoms; ig += atom[ig].hydrogenGroupSize ) {
01753 int hgs = atom[ig].hydrogenGroupSize;
01754 if (hgs == 1 ) continue;
01755 for ( i = 0; i < hgs; ++i ) {
01756 ref[i] = atom[ig+i].position;
01757 force[i] = f[Results::slow][ig+i];
01758 rmass[i] = 1. / atom[ig+i].mass;
01759 fixed[i] = ( simParams->fixedAtomsOn && atom[ig+i].atomFixed );
01760 if ( fixed[i] ) rmass[i] = 0.;
01761 }
01762 int icnt = 0;
01763
01764 if ( ( mol->rigid_bond_length(atom[ig].id) ) ) {
01765 if ( hgs != 3 ) {
01766 NAMD_die("Hydrogen group error caught in mollyMollify(). It's a bug!\n");
01767 }
01768 if ( !(fixed[1] && fixed[2]) ) {
01769 ial[icnt] = 1; ibl[icnt] = 2; ++icnt;
01770 }
01771 }
01772 for ( i = 1; i < hgs; ++i ) {
01773 if ( ( mol->rigid_bond_length(atom[ig+i].id) ) ) {
01774 if ( !(fixed[0] && fixed[i]) ) {
01775 ial[icnt] = 0; ibl[icnt] = i; ++icnt;
01776 }
01777 }
01778 }
01779
01780 if ( icnt == 0 ) continue;
01781 lambda = &(molly_lambda[numLambdas]);
01782 numLambdas += icnt;
01783 for ( i = 0; i < icnt; ++i ) {
01784 refab[i] = ref[ial[i]] - ref[ibl[i]];
01785 }
01786 avg = &(p_avg[ig]);
01787 mollify(avg,ref,lambda,force,hgs,icnt,rmass,ial,ibl,refab);
01788
01789 for ( i = 0; i < hgs; ++i ) {
01790 wc += outer(force[i]-f[Results::slow][ig+i],ref[i]);
01791 f[Results::slow][ig+i] = force[i];
01792 }
01793 }
01794
01795 *virial += wc;
01796 p_avg.resize(0);
01797 }
01798
01799 void HomePatch::checkpoint(void) {
01800 FullAtomList tmp_a(&atom); checkpoint_atom = tmp_a;
01801 checkpoint_lattice = lattice;
01802
01803
01804 #if NAMD_SeparateWaters != 0
01805 checkpoint_numWaterAtoms = numWaterAtoms;
01806 #endif
01807 }
01808
01809 void HomePatch::revert(void) {
01810 FullAtomList tmp_a(&checkpoint_atom); atom = tmp_a;
01811 numAtoms = atom.size();
01812 lattice = checkpoint_lattice;
01813
01814
01815 #if NAMD_SeparateWaters != 0
01816 numWaterAtoms = checkpoint_numWaterAtoms;
01817 #endif
01818 }
01819
01820 void HomePatch::submitLoadStats(int timestep)
01821 {
01822 LdbCoordinator::Object()->patchLoad(patchID,numAtoms,timestep);
01823 }
01824
01825
01826 void HomePatch::doPairlistCheck()
01827 {
01828 SimParameters *simParams = Node::Object()->simParameters;
01829
01830 if ( ! flags.usePairlists ) {
01831 flags.pairlistTolerance = 0.;
01832 flags.maxAtomMovement = 99999.;
01833 return;
01834 }
01835
01836 int i; int n = numAtoms;
01837 CompAtom *p_i = p.begin();
01838
01839 if ( flags.savePairlists ) {
01840 flags.pairlistTolerance = doPairlistCheck_newTolerance;
01841 flags.maxAtomMovement = 0.;
01842 doPairlistCheck_newTolerance *= (1. - simParams->pairlistShrink);
01843 doPairlistCheck_lattice = lattice;
01844 doPairlistCheck_positions.resize(numAtoms);
01845 CompAtom *psave_i = doPairlistCheck_positions.begin();
01846 for ( i=0; i<n; ++i ) { psave_i[i] = p_i[i]; }
01847 return;
01848 }
01849
01850 Lattice &lattice_old = doPairlistCheck_lattice;
01851 Position center_cur = lattice.unscale(center);
01852 Position center_old = lattice_old.unscale(center);
01853 Vector center_delta = center_cur - center_old;
01854
01855
01856 BigReal max_cd = 0.;
01857 for ( i=0; i<2; ++i ) {
01858 for ( int j=0; j<2; ++j ) {
01859 for ( int k=0; k<2; ++k ) {
01860 ScaledPosition corner( i ? min.x : max.x ,
01861 j ? min.y : max.y ,
01862 k ? min.z : max.z );
01863 Vector corner_delta =
01864 lattice.unscale(corner) - lattice_old.unscale(corner);
01865 corner_delta -= center_delta;
01866 BigReal cd = corner_delta.length2();
01867 if ( cd > max_cd ) max_cd = cd;
01868 }
01869 }
01870 }
01871 max_cd = sqrt(max_cd);
01872
01873
01874 BigReal max_pd = 0.;
01875 CompAtom *p_old_i = doPairlistCheck_positions.begin();
01876 for ( i=0; i<n; ++i ) {
01877 Vector p_delta = p_i[i].position - p_old_i[i].position;
01878 p_delta -= center_delta;
01879 BigReal pd = p_delta.length2();
01880 if ( pd > max_pd ) max_pd = pd;
01881 }
01882 max_pd = sqrt(max_pd);
01883
01884 BigReal max_tol = max_pd + max_cd;
01885
01886 flags.maxAtomMovement = max_tol;
01887
01888
01889
01890 if ( max_tol > ( (1. - simParams->pairlistTrigger) *
01891 doPairlistCheck_newTolerance ) ) {
01892 doPairlistCheck_newTolerance *= (1. + simParams->pairlistGrow);
01893 }
01894
01895 if ( max_tol > doPairlistCheck_newTolerance ) {
01896 doPairlistCheck_newTolerance = max_tol / (1. - simParams->pairlistTrigger);
01897 }
01898
01899 }
01900
01901 void HomePatch::doGroupSizeCheck()
01902 {
01903 if ( ! flags.doNonbonded ) return;
01904
01905 SimParameters *simParams = Node::Object()->simParameters;
01906 BigReal hgcut = 0.5 * simParams->hgroupCutoff; hgcut *= hgcut;
01907 BigReal maxrad2 = 0.;
01908
01909 FullAtomList::iterator p_i = atom.begin();
01910 FullAtomList::iterator p_e = atom.end();
01911
01912 while ( p_i != p_e ) {
01913 int hgs = p_i->hydrogenGroupSize;
01914 p_i->nonbondedGroupSize = hgs;
01915 BigReal x = p_i->position.x;
01916 BigReal y = p_i->position.y;
01917 BigReal z = p_i->position.z;
01918 ++p_i;
01919 int oversize = 0;
01920
01921 for ( int i = 1; i < hgs; ++i ) {
01922 p_i->nonbondedGroupSize = 0;
01923 BigReal dx = p_i->position.x - x;
01924 BigReal dy = p_i->position.y - y;
01925 BigReal dz = p_i->position.z - z;
01926 BigReal r2 = dx * dx + dy * dy + dz * dz;
01927 ++p_i;
01928 if ( r2 > hgcut ) oversize = 1;
01929 else if ( r2 > maxrad2 ) maxrad2 = r2;
01930 }
01931
01932 if ( oversize || hgs > 4 ) {
01933 p_i -= hgs;
01934 for ( int i = 0; i < hgs; ++i ) {
01935 p_i->nonbondedGroupSize = 1;
01936 ++p_i;
01937 }
01938 }
01939 }
01940
01941 flags.maxGroupRadius = sqrt(maxrad2);
01942
01943 }
01944
01945 void HomePatch::doMarginCheck()
01946 {
01947 SimParameters *simParams = Node::Object()->simParameters;
01948
01949 BigReal sysdima = lattice.a_r().unit() * lattice.a();
01950 BigReal sysdimb = lattice.b_r().unit() * lattice.b();
01951 BigReal sysdimc = lattice.c_r().unit() * lattice.c();
01952
01953 BigReal minSize = simParams->patchDimension - simParams->margin;
01954
01955 if ( ( (max.x - min.x)*aAway*sysdima < minSize*0.9999 ) ||
01956 ( (max.y - min.y)*bAway*sysdimb < minSize*0.9999 ) ||
01957 ( (max.z - min.z)*cAway*sysdimc < minSize*0.9999 ) ) {
01958
01959 NAMD_die("Periodic cell has become too small for original patch grid!\n"
01960 "Possible solutions are to restart from a recent checkpoint,\n"
01961 "increase margin, or disable useFlexibleCell for liquid simulation.");
01962 }
01963
01964 BigReal cutoff = simParams->cutoff;
01965
01966 BigReal margina = 0.5 * ( (max.x - min.x) * aAway - cutoff / sysdima );
01967 BigReal marginb = 0.5 * ( (max.y - min.y) * bAway - cutoff / sysdimb );
01968 BigReal marginc = 0.5 * ( (max.z - min.z) * cAway - cutoff / sysdimc );
01969
01970 if ( (margina < -0.0001) || (marginb < -0.0001) || (marginc < -0.0001) ) {
01971 NAMD_die("Periodic cell has become too small for original patch grid!\n"
01972 "There are probably many margin violations already reported.\n"
01973 "Possible solutions are to restart from a recent checkpoint,\n"
01974 "increase margin, or disable useFlexibleCell for liquid simulation.");
01975 }
01976
01977 BigReal minx = min.x - margina;
01978 BigReal miny = min.y - marginb;
01979 BigReal minz = min.z - marginc;
01980 BigReal maxx = max.x + margina;
01981 BigReal maxy = max.y + marginb;
01982 BigReal maxz = max.z + marginc;
01983
01984 int xdev, ydev, zdev;
01985 int problemCount = 0;
01986
01987 FullAtomList::iterator p_i = atom.begin();
01988 FullAtomList::iterator p_e = atom.end();
01989 for ( ; p_i != p_e; ++p_i ) {
01990
01991 ScaledPosition s = lattice.scale(p_i->position);
01992
01993
01994 if (s.x < minx) xdev = 0;
01995 else if (maxx <= s.x) xdev = 2;
01996 else xdev = 1;
01997
01998 if (s.y < miny) ydev = 0;
01999 else if (maxy <= s.y) ydev = 2;
02000 else ydev = 1;
02001
02002 if (s.z < minz) zdev = 0;
02003 else if (maxz <= s.z) zdev = 2;
02004 else zdev = 1;
02005
02006 if (mInfo[xdev][ydev][zdev]) {
02007 ++problemCount;
02008 }
02009
02010 }
02011
02012 marginViolations = problemCount;
02013
02014
02015
02016
02017
02018 }
02019
02020
02021 void
02022 HomePatch::doAtomMigration()
02023 {
02024 int i;
02025
02026 for (i=0; i<numNeighbors; i++) {
02027 realInfo[i].mList.resize(0);
02028 }
02029
02030
02031 AtomMap::Object()->unregisterIDs(patchID,pExt.begin(),pExt.end());
02032
02033
02034
02035 BigReal minx = min.x;
02036 BigReal miny = min.y;
02037 BigReal minz = min.z;
02038 BigReal maxx = max.x;
02039 BigReal maxy = max.y;
02040 BigReal maxz = max.z;
02041
02042 int xdev, ydev, zdev;
02043 int delnum = 0;
02044
02045 FullAtomList::iterator atom_i = atom.begin();
02046 FullAtomList::iterator atom_e = atom.end();
02047
02048
02049 #if NAMD_SeparateWaters != 0
02050 FullAtomList::iterator atom_first = atom_i;
02051 int numLostWaterAtoms = 0;
02052 #endif
02053
02054 while ( atom_i != atom_e ) {
02055 if ( atom_i->hydrogenGroupSize ) {
02056
02057 ScaledPosition s = lattice.scale(atom_i->position);
02058
02059
02060 if (s.x < minx) xdev = 0;
02061 else if (maxx <= s.x) xdev = 2;
02062 else xdev = 1;
02063
02064 if (s.y < miny) ydev = 0;
02065 else if (maxy <= s.y) ydev = 2;
02066 else ydev = 1;
02067
02068 if (s.z < minz) zdev = 0;
02069 else if (maxz <= s.z) zdev = 2;
02070 else zdev = 1;
02071
02072 }
02073
02074 if (mInfo[xdev][ydev][zdev]) {
02075
02076
02077
02078 MigrationList &mCur = mInfo[xdev][ydev][zdev]->mList;
02079 DebugM(3,"Migrating atom " << atom_i->id << " from patch "
02080 << patchID << " with position " << atom_i->position << "\n");
02081 mCur.add(*atom_i);
02082
02083 ++delnum;
02084
02085
02086
02087 #if NAMD_SeparateWaters != 0
02088
02089
02090
02091
02092
02093
02094
02095 int atomIndex = atom_i - atom_first;
02096 if (atomIndex < numWaterAtoms)
02097 numLostWaterAtoms++;
02098 #endif
02099
02100
02101 } else {
02102
02103 if ( delnum ) { *(atom_i-delnum) = *atom_i; }
02104
02105 }
02106
02107 ++atom_i;
02108 }
02109
02110
02111 #if NAMD_SeparateWaters != 0
02112 numWaterAtoms -= numLostWaterAtoms;
02113 #endif
02114
02115
02116 int delpos = numAtoms - delnum;
02117 DebugM(4,"numAtoms " << numAtoms << " deleted " << delnum << "\n");
02118 atom.del(delpos,delnum);
02119
02120 numAtoms = atom.size();
02121
02122 PatchMgr::Object()->sendMigrationMsgs(patchID, realInfo, numNeighbors);
02123
02124
02125 inMigration = true;
02126
02127
02128 for (i=0; i<numMlBuf; i++) {
02129 DebugM(1, "Draining migration buffer ("<<i<<","<<numMlBuf<<")\n");
02130 depositMigration(msgbuf[i]);
02131 }
02132 numMlBuf = 0;
02133
02134 if (!allMigrationIn) {
02135 DebugM(3,"All Migrations NOT in, we are suspending patch "<<patchID<<"\n");
02136 migrationSuspended = true;
02137 sequencer->suspend();
02138 migrationSuspended = false;
02139 }
02140 allMigrationIn = false;
02141
02142 inMigration = false;
02143 marginViolations = 0;
02144 }
02145
02146 void
02147 HomePatch::depositMigration(MigrateAtomsMsg *msg)
02148 {
02149
02150 if (!inMigration) {
02151
02152 msgbuf[numMlBuf++] = msg;
02153 return;
02154 }
02155
02156
02157
02158 #if NAMD_SeparateWaters != 0
02159
02160
02161
02162
02163
02164
02165 mergeAtomList(msg->migrationList);
02166
02167
02168 #else
02169
02170
02171
02172
02173 {
02174 MigrationList &migrationList = msg->migrationList;
02175 MigrationList::iterator mi;
02176 Transform mother_transform;
02177 for (mi = migrationList.begin(); mi != migrationList.end(); mi++) {
02178 DebugM(1,"Migrating atom " << mi->id << " to patch "
02179 << patchID << " with position " << mi->position << "\n");
02180 if ( mi->hydrogenGroupSize ) {
02181 mi->position = lattice.nearest(mi->position,center,&(mi->transform));
02182 mother_transform = mi->transform;
02183 } else {
02184 mi->position = lattice.reverse_transform(mi->position,mi->transform);
02185 mi->position = lattice.apply_transform(mi->position,mother_transform);
02186 mi->transform = mother_transform;
02187 }
02188 atom.add(*mi);
02189 }
02190 }
02191
02192
02193 #endif // if (NAMD_SeparateWaters != 0)
02194
02195
02196 numAtoms = atom.size();
02197 delete msg;
02198
02199 DebugM(3,"Counter on " << patchID << " = " << patchMigrationCounter << "\n");
02200 if (!--patchMigrationCounter) {
02201 DebugM(3,"All Migrations are in for patch "<<patchID<<"\n");
02202 allMigrationIn = true;
02203 patchMigrationCounter = numNeighbors;
02204 if (migrationSuspended) {
02205 DebugM(3,"patch "<<patchID<<" is being awakened\n");
02206 migrationSuspended = false;
02207 sequencer->awaken();
02208 return;
02209 }
02210 else {
02211 DebugM(3,"patch "<<patchID<<" was not suspended\n");
02212 }
02213 }
02214 }
02215
02216
02217
02218
02219 #if NAMD_SeparateWaters != 0
02220
02221
02222
02223
02224 void HomePatch::separateAtoms() {
02225
02226
02227
02228
02229
02230
02231
02232
02233
02234
02235 if (atom.size() < 0) return;
02236
02237
02238 tempAtom.resize(numAtoms);
02239
02240
02241 int wathgsize = 3;
02242 if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
02243 else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
02244
02245
02246 int i = 0;
02247 int waterIndex = 0;
02248 int nonWaterIndex = 0;
02249 while (i < numAtoms) {
02250
02251 FullAtom &atom_i = atom[i];
02252 Mass mass = atom_i.mass;
02253 int hgs = atom_i.hydrogenGroupSize;
02254
02255 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
02256
02257
02258 if (waterIndex != i) {
02259 atom[waterIndex ] = atom[i ];
02260 atom[waterIndex + 1] = atom[i + 1];
02261 atom[waterIndex + 2] = atom[i + 2];
02262 if (wathgsize > 3) atom[waterIndex + 3] = atom[i + 3];
02263 if (wathgsize > 4) atom[waterIndex + 4] = atom[i + 4];
02264
02265 }
02266
02267 waterIndex += wathgsize;
02268 i += wathgsize;
02269
02270 } else {
02271
02272
02273 for (int j = 0; j < hgs; j++)
02274 tempAtom[nonWaterIndex + j] = atom[i + j];
02275
02276 nonWaterIndex += hgs;
02277 i += hgs;
02278 }
02279
02280 }
02281
02282
02283
02284
02285
02286
02287
02288 for (i = 0; i < nonWaterIndex; i++)
02289 atom[waterIndex + i] = tempAtom[i];
02290
02291
02292 numWaterAtoms = waterIndex;
02293 }
02294
02295
02296
02297
02298
02299
02300
02301 void HomePatch::mergeAtomList(FullAtomList &al) {
02302
02303
02304 if (al.size() <= 0) return;
02305
02306 const int orig_atomSize = atom.size();
02307 const int orig_alSize = al.size();
02308
02309
02310 atom.resize(orig_atomSize + orig_alSize);
02311
02312
02313 #if 0 // version where non-waters are moved to scratch first
02314
02315
02316
02317
02318
02319
02320
02321
02322
02323
02324
02325
02326
02327
02328
02329
02330
02331
02332
02333
02334
02335 const int orig_atom_numNonWaters = orig_atomSize - numWaterAtoms;
02336 tempAtom.resize(orig_atom_numNonWaters + al.size());
02337 for (int i = 0; i < orig_atom_numNonWaters; i++)
02338 tempAtom[i] = atom[numWaterAtoms + i];
02339
02340
02341
02342 int atom_waterIndex = numWaterAtoms;
02343 int atom_nonWaterIndex = orig_atom_numNonWaters;
02344 int i = 0;
02345 while (i < orig_alSize) {
02346
02347 FullAtom &atom_i = al[i];
02348 int hgs = atom_i.hydrogenGroupSize;
02349 Mass mass = atom_i.mass;
02350
02351 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
02352
02353
02354
02355
02356 al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
02357 Transform mother_transform = al[i].transform;
02358
02359
02360 al[i+1].position = lattice.reverse_transform(al[i+1].position, al[i+1].transform);
02361 al[i+1].position = lattice.apply_transform(al[i+1].position, mother_transform);
02362 al[i+1].transform = mother_transform;
02363
02364
02365 al[i+2].position = lattice.reverse_transform(al[i+2].position, al[i+2].transform);
02366 al[i+2].position = lattice.apply_transform(al[i+2].position, mother_transform);
02367 al[i+2].transform = mother_transform;
02368
02369
02370 atom[atom_waterIndex ] = al[i ];
02371 atom[atom_waterIndex + 1] = al[i + 1];
02372 atom[atom_waterIndex + 2] = al[i + 2];
02373
02374 atom_waterIndex += 3;
02375 i += 3;
02376
02377 } else {
02378
02379
02380
02381
02382 al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
02383 Transform mother_transform = al[i].transform;
02384
02385
02386 for (int j = 1; j < hgs; j++) {
02387 al[i+j].position = lattice.reverse_transform(al[i+j].position, al[i+j].transform);
02388 al[i+j].position = lattice.apply_transform(al[i+j].position, mother_transform);
02389 al[i+j].transform = mother_transform;
02390 }
02391
02392
02393 for (int j = 0; j < hgs; j++)
02394 tempAtom[atom_nonWaterIndex + j] = al[i + j];
02395
02396 atom_nonWaterIndex += hgs;
02397 i += hgs;
02398 }
02399
02400 }
02401
02402
02403 for (int i = 0; i < atom_nonWaterIndex; i++)
02404 atom[atom_waterIndex + i] = tempAtom[i];
02405
02406
02407 numWaterAtoms = atom_waterIndex;
02408 numAtoms = atom.size();
02409
02410
02411 #else
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421
02422
02423
02424
02425
02426 int wathgsize = 3;
02427 if (simParams->watmodel == WAT_TIP4) wathgsize = 4;
02428 else if (simParams->watmodel == WAT_SWM4) wathgsize = 5;
02429
02430
02431 int al_numWaterAtoms = 0;
02432 int i = 0;
02433 while (i < orig_alSize) {
02434
02435 FullAtom &atom_i = al[i];
02436 int hgs = atom_i.hydrogenGroupSize;
02437 Mass mass = atom_i.mass;
02438
02439 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
02440 al_numWaterAtoms += wathgsize;
02441 }
02442
02443 i += hgs;
02444 }
02445
02446
02447
02448 if (al_numWaterAtoms > 0) {
02449 for (i = orig_atomSize - 1; i >= numWaterAtoms; i--) {
02450 atom[i + al_numWaterAtoms] = atom[i];
02451 }
02452 }
02453
02454
02455
02456
02457 int atom_waterIndex = numWaterAtoms;
02458 int atom_nonWaterIndex = orig_atomSize + al_numWaterAtoms;
02459 i = 0;
02460 while (i < orig_alSize) {
02461
02462 FullAtom &atom_i = al[i];
02463 int hgs = atom_i.hydrogenGroupSize;
02464 Mass mass = atom_i.mass;
02465
02466 if (IS_HYDROGEN_GROUP_WATER(hgs, mass)) {
02467
02468
02469
02470
02471 al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
02472 Transform mother_transform = al[i].transform;
02473
02474
02475 al[i+1].position = lattice.reverse_transform(al[i+1].position, al[i+1].transform);
02476 al[i+1].position = lattice.apply_transform(al[i+1].position, mother_transform);
02477 al[i+1].transform = mother_transform;
02478
02479
02480 al[i+2].position = lattice.reverse_transform(al[i+2].position, al[i+2].transform);
02481 al[i+2].position = lattice.apply_transform(al[i+2].position, mother_transform);
02482 al[i+2].transform = mother_transform;
02483
02484
02485 atom[atom_waterIndex ] = al[i ];
02486 atom[atom_waterIndex + 1] = al[i + 1];
02487 atom[atom_waterIndex + 2] = al[i + 2];
02488
02489 if (wathgsize > 3) atom[atom_waterIndex + 3] = al[i + 3];
02490
02491 atom_waterIndex += wathgsize;
02492 i += wathgsize;
02493
02494 } else {
02495
02496
02497
02498
02499 al[i].position = lattice.nearest(al[i].position, center, &(al[i].transform));
02500 Transform mother_transform = al[i].transform;
02501
02502
02503 for (int j = 1; j < hgs; j++) {
02504 al[i+j].position = lattice.reverse_transform(al[i+j].position, al[i+j].transform);
02505 al[i+j].position = lattice.apply_transform(al[i+j].position, mother_transform);
02506 al[i+j].transform = mother_transform;
02507 }
02508
02509
02510 for (int j = 0; j < hgs; j++)
02511 atom[atom_nonWaterIndex + j] = al[i + j];
02512
02513 atom_nonWaterIndex += hgs;
02514 i += hgs;
02515 }
02516
02517 }
02518
02519
02520 numWaterAtoms = atom_waterIndex;
02521 numAtoms = atom_nonWaterIndex;
02522
02523 #endif
02524 }
02525
02526 #endif
02527
02528
02529
02530 inline void lubksb(HGMatrixBigReal &a, int n, HGArrayInt &indx,
02531 HGArrayBigReal &b)
02532 {
02533 int i,ii=-1,ip,j;
02534 double sum;
02535
02536 for (i=0;i<n;i++) {
02537 ip=indx[i];
02538 sum=b[ip];
02539 b[ip]=b[i];
02540 if (ii >= 0)
02541 for (j=ii;j<i;j++) sum -= a[i][j]*b[j];
02542 else if (sum) ii=i;
02543 b[i]=sum;
02544 }
02545 for (i=n-1;i>=0;i--) {
02546 sum=b[i];
02547 for (j=i+1;j<n;j++) sum -= a[i][j]*b[j];
02548 b[i]=sum/a[i][i];
02549 }
02550 }
02551
02552
02553 inline void ludcmp(HGMatrixBigReal &a, int n, HGArrayInt &indx, BigReal *d)
02554 {
02555
02556 int i,imax,j,k;
02557 double big,dum,sum,temp;
02558 HGArrayBigReal vv;
02559 *d=1.0;
02560 for (i=0;i<n;i++) {
02561 big=0.0;
02562 for (j=0;j<n;j++)
02563 if ((temp=fabs(a[i][j])) > big) big=temp;
02564 if (big == 0.0) NAMD_die("Singular matrix in routine ludcmp\n");
02565 vv[i]=1.0/big;
02566 }
02567 for (j=0;j<n;j++) {
02568 for (i=0;i<j;i++) {
02569 sum=a[i][j];
02570 for (k=0;k<i;k++) sum -= a[i][k]*a[k][j];
02571 a[i][j]=sum;
02572 }
02573 big=0.0;
02574 for (i=j;i<n;i++) {
02575 sum=a[i][j];
02576 for (k=0;k<j;k++)
02577 sum -= a[i][k]*a[k][j];
02578 a[i][j]=sum;
02579 if ( (dum=vv[i]*fabs(sum)) >= big) {
02580 big=dum;
02581 imax=i;
02582 }
02583 }
02584 if (j != imax) {
02585 for (k=0;k<n;k++) {
02586 dum=a[imax][k];
02587 a[imax][k]=a[j][k];
02588 a[j][k]=dum;
02589 }
02590 *d = -(*d);
02591 vv[imax]=vv[j];
02592 }
02593 indx[j]=imax;
02594 if (a[j][j] == 0.0) a[j][j]=TINY;
02595 if (j != n-1) {
02596 dum=1.0/(a[j][j]);
02597 for (i=j+1;i<n;i++) a[i][j] *= dum;
02598 }
02599 }
02600 }
02601
02602
02603 inline void G_q(const HGArrayVector &refab, HGMatrixVector &gqij,
02604 const int n, const int m, const HGArrayInt &ial, const HGArrayInt &ibl) {
02605 int i;
02606
02607 for(i=0;i<m;i++) {
02608 gqij[i][ial[i]]=2.0*refab[i];
02609 gqij[i][ibl[i]]=-gqij[i][ial[i]];
02610 }
02611 };
02612
02613
02614
02615 int average(CompAtom *qtilde,const HGArrayVector &q,BigReal *lambda,const int n,const int m, const HGArrayBigReal &imass, const HGArrayBigReal &length2, const HGArrayInt &ial, const HGArrayInt &ibl, const HGArrayVector &refab, const BigReal tolf, const int ntrial) {
02616
02617
02618
02619
02620
02621
02622
02623
02624
02625
02626
02627
02628
02629 int k,k1,i,j;
02630 BigReal errx,errf,d,tolx;
02631
02632 HGArrayInt indx;
02633 HGArrayBigReal p;
02634 HGArrayBigReal fvec;
02635 HGMatrixBigReal fjac;
02636 HGArrayVector avgab;
02637 HGMatrixVector grhs;
02638 HGMatrixVector auxrhs;
02639 HGMatrixVector glhs;
02640
02641
02642 tolx=tolf;
02643
02644
02645
02646 for (i=0;i<m;i++) {
02647 lambda[i]=0.0;
02648 }
02649
02650
02651
02652
02653 G_q(refab,grhs,n,m,ial,ibl);
02654 for (k=1;k<=ntrial;k++) {
02655
02656 HGArrayBigReal gij;
02657
02658
02659 {
02660 BigReal multiplier;
02661 HGArrayVector tmp;
02662 for (i=0;i<m;i++) {
02663 multiplier = lambda[i];
02664
02665 for (j=0;j<n;j++) {
02666 auxrhs[i][j]=multiplier*imass[j]*grhs[i][j];
02667 }
02668 }
02669 for (j=0;j<n;j++) {
02670
02671 for (i=0;i<m;i++) {
02672 tmp[j]+=auxrhs[i][j];
02673 }
02674 }
02675
02676 for (j=0;j<n;j++) {
02677 qtilde[j].position=q[j]+tmp[j];
02678 }
02679
02680 }
02681
02682 for ( i = 0; i < m; i++ ) {
02683 avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
02684 }
02685
02686
02687
02688 {
02689
02690
02691 HGMatrixVector grhs2;
02692
02693 G_q(avgab,glhs,n,m,ial,ibl);
02694 #ifdef DEBUG0
02695 iout<<iINFO << "G_q:" << std::endl<<endi;
02696 for (i=0;i<m;i++) {
02697 iout<<iINFO << glhs[i*n+0] << " " << glhs[i*n+1] << " " << glhs[i*n+2] << std::endl<<endi;
02698 }
02699 #endif
02700
02701
02702 for (j=0; j<n; j++) {
02703 for (i=0; i<m;i++) {
02704 grhs2[i][j] = grhs[i][j]*imass[j];
02705 }
02706 }
02707
02708
02709
02710 for (i=0;i<m;i++) {
02711 for (j=0;j<m;j++) {
02712 fjac[i][j] = 0;
02713 for (k1=0;k1<n;k1++) {
02714 fjac[i][j] += glhs[i][k1]*grhs2[j][k1];
02715 }
02716 }
02717 }
02718 #ifdef DEBUG0
02719 iout<<iINFO << "glhs" <<endi;
02720 for(i=0;i<9;i++) {
02721 iout<<iINFO << glhs[i] << ","<<endi;
02722 }
02723 iout<<iINFO << std::endl<<endi;
02724 for(i=0;i<9;i++) {
02725 iout<<iINFO << grhs2[i] << ","<<endi;
02726 }
02727 iout<<iINFO << std::endl<<endi;
02728 #endif
02729
02730 }
02731
02732 #ifdef DEBUG0
02733 iout<<iINFO << "Jac" << std::endl<<endi;
02734 for (i=0;i<m;i++)
02735 for (j=0;j<m;j++)
02736 iout<<iINFO << fjac[i][j] << " "<<endi;
02737 iout<< std::endl<<endi;
02738 #endif
02739
02740
02741 for (i=0;i<m;i++) {
02742 gij[i]=avgab[i]*avgab[i]-length2[i];
02743 }
02744 #ifdef DEBUG0
02745 iout<<iINFO << "G" << std::endl<<endi;
02746 iout<<iINFO << "( "<<endi;
02747 for(i=0;i<m-1;i++) {
02748 iout<<iINFO << gij[i] << ", "<<endi;
02749 }
02750 iout<<iINFO << gij[m-1] << ")" << std::endl<<endi;
02751 #endif
02752
02753 for(i=0;i<m;i++) {
02754 fvec[i] = gij[i];
02755 }
02756
02757
02758
02759 errf=0.0;
02760 for (i=0;i<m;i++) errf += fabs(fvec[i]);
02761 #ifdef DEBUG0
02762 iout<<iINFO << "errf: " << errf << std::endl<<endi;
02763 #endif
02764 if (errf <= tolf) {
02765 break;
02766 }
02767 for (i=0;i<m;i++) p[i] = -fvec[i];
02768
02769 ludcmp(fjac,m,indx,&d);
02770 lubksb(fjac,m,indx,p);
02771
02772 errx=0.0;
02773 for (i=0;i<m;i++) {
02774 errx += fabs(p[i]);
02775 }
02776 for (i=0;i<m;i++)
02777 lambda[i] += p[i];
02778
02779 #ifdef DEBUG0
02780 iout<<iINFO << "lambda:" << lambda[0]
02781 << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
02782 iout<<iINFO << "errx: " << errx << std::endl<<endi;
02783 #endif
02784 if (errx <= tolx) break;
02785 #ifdef DEBUG0
02786 iout<<iINFO << "Qtilde:" << std::endl<<endi;
02787 iout<<iINFO << qtilde[0].position << " " << qtilde[1].position << " " << qtilde[2].position << std::endl<<endi;
02788 #endif
02789 }
02790 #ifdef DEBUG
02791 iout<<iINFO << "LAMBDA:" << lambda[0] << " " << lambda[1] << " " << lambda[2] << std::endl<<endi;
02792 #endif
02793
02794 return k;
02795 }
02796
02797 void mollify(CompAtom *qtilde,const HGArrayVector &q0,const BigReal *lambda, HGArrayVector &force,const int n, const int m, const HGArrayBigReal &imass,const HGArrayInt &ial,const HGArrayInt &ibl,const HGArrayVector &refab) {
02798 int i,j,k;
02799 BigReal d;
02800 HGMatrixBigReal fjac;
02801 Vector zero(0.0,0.0,0.0);
02802
02803 HGArrayVector tmpforce;
02804 HGArrayVector tmpforce2;
02805 HGArrayVector y;
02806 HGMatrixVector grhs;
02807 HGMatrixVector glhs;
02808 HGArrayBigReal aux;
02809 HGArrayInt indx;
02810
02811 for(i=0;i<n;i++) {
02812 tmpforce[i]=imass[i]*force[i];
02813 }
02814
02815 HGMatrixVector grhs2;
02816 HGArrayVector avgab;
02817
02818 for ( i = 0; i < m; i++ ) {
02819 avgab[i] = qtilde[ial[i]].position - qtilde[ibl[i]].position;
02820 }
02821
02822 G_q(avgab,glhs,n,m,ial,ibl);
02823 G_q(refab,grhs,n,m,ial,ibl);
02824
02825 for (j=0; j<n; j++) {
02826 for (i=0; i<m;i++) {
02827 grhs2[i][j] = grhs[i][j]*imass[j];
02828 }
02829 }
02830
02831
02832
02833 for (i=0;i<m;i++) {
02834 for (j=0;j<m;j++) {
02835 fjac[j][i] = 0;
02836 for (k=0;k<n;k++) {
02837 fjac[j][i] += glhs[i][k]*grhs2[j][k];
02838 }
02839 }
02840 }
02841
02842
02843
02844
02845 for(i=0;i<m;i++) {
02846 aux[i]=0.0;
02847 for(j=0;j<n;j++) {
02848 aux[i]+=grhs[i][j]*tmpforce[j];
02849 }
02850 }
02851
02852 ludcmp(fjac,m,indx,&d);
02853 lubksb(fjac,m,indx,aux);
02854
02855 for(j=0;j<n;j++) {
02856 y[j] = zero;
02857 for(i=0;i<m;i++) {
02858 y[j] += aux[i]*glhs[i][j];
02859 }
02860 }
02861 for(i=0;i<n;i++) {
02862 y[i]=force[i]-y[i];
02863 }
02864
02865
02866 for(i=0;i<n;i++) {
02867 tmpforce2[i]=imass[i]*y[i];
02868 }
02869
02870
02871 for (i=0;i<n;i++) {
02872 tmpforce[i]=zero;
02873 }
02874
02875 for (j=0;j<m;j++) {
02876 Vector tmpf = 2.0*lambda[j]*(tmpforce2[ial[j]]-tmpforce2[ibl[j]]);
02877 tmpforce[ial[j]] += tmpf;
02878 tmpforce[ibl[j]] -= tmpf;
02879 }
02880
02881
02882 for(i=0;i<n;i++) {
02883 force[i]=tmpforce[i]+y[i];
02884 }
02885
02886 }
02887
02888 #if CMK_PERSISTENT_COMM
02889 void HomePatch::destoryPersistComm()
02890 {
02891 for (int i=0; i<nphs; i++) {
02892 CmiDestoryPersistent(localphs[i]);
02893 }
02894 phsReady = 0;
02895 }
02896 #endif