Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

ComputePme Class Reference

#include <ComputePme.h>

Inheritance diagram for ComputePme:

ComputeHomePatches Compute List of all members.

Public Member Functions

 ComputePme (ComputeID c)
virtual ~ComputePme ()
void doWork ()
void sendData (int, int *, int *, int *)
void sendPencils ()
void copyResults (PmeGridMsg *)
void copyPencils (PmeGridMsg *)
void ungridForces ()
void setMgr (ComputePmeMgr *mgr)

Constructor & Destructor Documentation

ComputePme::ComputePme ComputeID  c  ) 
 

Definition at line 1345 of file ComputePme.C.

References SimParameters::alchDecouple, SimParameters::alchElecLambdaStart, SimParameters::alchFepOn, SimParameters::alchThermIntOn, DebugM, PmeGrid::dim2, PmeGrid::dim3, PmeGrid::K1, PmeGrid::K2, PmeGrid::K3, SimParameters::lesFactor, SimParameters::lesOn, Node::Object(), ReductionMgr::Object(), PmeGrid::order, SimParameters::pairInteractionOn, SimParameters::pairInteractionSelf, SimParameters::PMEGridSizeX, SimParameters::PMEGridSizeY, SimParameters::PMEGridSizeZ, SimParameters::PMEInterpOrder, REDUCTIONS_BASIC, Node::simParameters, simParams, and ReductionMgr::willSubmit().

01345                                   :
01346   ComputeHomePatches(c)
01347 {
01348   DebugM(4,"ComputePme created.\n");
01349 
01350   basePriority = PME_PRIORITY;
01351 
01352   CProxy_ComputePmeMgr::ckLocalBranch(
01353         CkpvAccess(BOCclass_group).computePmeMgr)->setCompute(this);
01354 
01355   useAvgPositions = 1;
01356 
01357   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
01358 
01359   SimParameters *simParams = Node::Object()->simParameters;
01360 
01361   alchFepOn = simParams->alchFepOn;
01362   alchThermIntOn = simParams->alchThermIntOn;
01363   alchDecouple = (alchFepOn || alchThermIntOn) && (simParams->alchDecouple);
01364   alchElecLambdaStart = (alchFepOn || alchThermIntOn) ? 
01365     simParams->alchElecLambdaStart : 0;
01366             
01367   if (alchFepOn || alchThermIntOn) {
01368     numGrids = 2;
01369     if (alchDecouple) numGrids += 2;
01370     if (alchElecLambdaStart || alchThermIntOn) numGrids ++;
01371   }
01372   else numGrids = 1;
01373   lesOn = simParams->lesOn;
01374   if ( lesOn ) {
01375     lesFactor = simParams->lesFactor;
01376     numGrids = lesFactor;
01377   }
01378   selfOn = 0;
01379   pairOn = simParams->pairInteractionOn;
01380   if ( pairOn ) {
01381     selfOn = simParams->pairInteractionSelf;
01382     if ( selfOn ) pairOn = 0;  // make pairOn and selfOn exclusive
01383     numGrids = selfOn ? 1 : 3;
01384   }
01385 
01386   myGrid.K1 = simParams->PMEGridSizeX;
01387   myGrid.K2 = simParams->PMEGridSizeY;
01388   myGrid.K3 = simParams->PMEGridSizeZ;
01389   myGrid.order = simParams->PMEInterpOrder;
01390   myGrid.dim2 = myGrid.K2;
01391   myGrid.dim3 = 2 * (myGrid.K3/2 + 1);
01392   qsize = myGrid.K1 * myGrid.dim2 * myGrid.dim3;
01393   fsize = myGrid.K1 * myGrid.dim2;
01394   q_arr = new double*[fsize*numGrids];
01395   memset( (void*) q_arr, 0, fsize*numGrids * sizeof(double*) );
01396   f_arr = new char[fsize*numGrids];
01397   fz_arr = new char[myGrid.K3];
01398 }

ComputePme::~ComputePme  )  [virtual]
 

Definition at line 1400 of file ComputePme.C.

01401 {
01402   for (int i=0; i<fsize*numGrids; ++i) {
01403     if ( q_arr[i] ) {
01404       delete [] q_arr[i];
01405     }
01406   }
01407   delete [] q_arr;
01408   delete [] f_arr;
01409   delete [] fz_arr;
01410 }


Member Function Documentation

void ComputePme::copyPencils PmeGridMsg  ) 
 

Definition at line 1780 of file ComputePme.C.

References PmeGrid::block1, PmeGrid::block2, PmeGrid::dim2, PmeGrid::dim3, PmeGridMsg::evir, j, PmeGrid::K1, PmeGrid::K2, PmeGridMsg::qgrid, PmeGridMsg::sourceNode, ComputePmeMgr::xBlocks, ComputePmeMgr::yBlocks, ComputePmeMgr::zBlocks, PmeGridMsg::zlist, and PmeGridMsg::zlistlen.

Referenced by ComputePmeMgr::recvUngrid().

01780                                             {
01781 
01782   int xBlocks = myMgr->xBlocks;
01783   int yBlocks = myMgr->yBlocks;
01784   int zBlocks = myMgr->zBlocks;
01785   myGrid.block1 = ( myGrid.K1 + xBlocks - 1 ) / xBlocks;
01786   myGrid.block2 = ( myGrid.K2 + yBlocks - 1 ) / yBlocks;
01787   int K1 = myGrid.K1;
01788   int K2 = myGrid.K2;
01789   int dim2 = myGrid.dim2;
01790   int dim3 = myGrid.dim3;
01791   int block1 = myGrid.block1;
01792   int block2 = myGrid.block2;
01793 
01794   // msg->sourceNode = thisIndex.x * initdata.yBlocks + thisIndex.y;
01795   int ib = msg->sourceNode / yBlocks;
01796   int jb = msg->sourceNode % yBlocks;
01797 
01798   int ibegin = ib*block1;
01799   int iend = ibegin + block1;  if ( iend > K1 ) iend = K1;
01800   int jbegin = jb*block2;
01801   int jend = jbegin + block2;  if ( jend > K2 ) jend = K2;
01802 
01803   int zlistlen = msg->zlistlen;
01804   int *zlist = msg->zlist;
01805   float *qmsg = msg->qgrid;
01806   int g;
01807   for ( g=0; g<numGrids; ++g ) {
01808     evir[g] += msg->evir[g];
01809     char *f = f_arr + g*fsize;
01810     double **q = q_arr + g*fsize;
01811     for ( int i=ibegin; i<iend; ++i ) {
01812      for ( int j=jbegin; j<jend; ++j ) {
01813       if( f[i*dim2+j] ) {
01814         for ( int k=0; k<zlistlen; ++k ) {
01815           q[i*dim2+j][zlist[k]] = *(qmsg++);
01816         }
01817       }
01818      }
01819     }
01820   }
01821 }

void ComputePme::copyResults PmeGridMsg  ) 
 

Definition at line 1927 of file ComputePme.C.

References PmeGrid::dim3, PmeGridMsg::evir, PmeGridMsg::fgrid, PmeGridMsg::len, PmeGridMsg::qgrid, PmeGridMsg::start, PmeGridMsg::zlist, and PmeGridMsg::zlistlen.

Referenced by ComputePmeMgr::recvUngrid().

01927                                             {
01928 
01929   int zdim = myGrid.dim3;
01930   int flen = msg->len;
01931   int fstart = msg->start;
01932   int zlistlen = msg->zlistlen;
01933   int *zlist = msg->zlist;
01934   float *qmsg = msg->qgrid;
01935   int g;
01936   for ( g=0; g<numGrids; ++g ) {
01937     evir[g] += msg->evir[g];
01938     char *f = msg->fgrid + g*flen;
01939     double **q = q_arr + fstart + g*fsize;
01940     for ( int i=0; i<flen; ++i ) {
01941       if ( f[i] ) {
01942         for ( int k=0; k<zlistlen; ++k ) {
01943           q[i][zlist[k]] = *(qmsg++);
01944         }
01945       }
01946     }
01947   }
01948 }

void ComputePme::doWork  )  [virtual]
 

Reimplemented from Compute.

Definition at line 1412 of file ComputePme.C.

References ResizeArrayIter< T >::begin(), BigReal, PmeParticle::cg, CompAtom::charge, COLOUMB, DebugM, PmeGrid::dim3, ResizeArrayIter< T >::end(), PmeRealSpace::fill_charges(), ComputePmeMgr::gridPeMap, ComputePmeMgr::gridPeOrder, PmeGrid::K3, NAMD_bug(), ComputePmeMgr::numGridPes, CompAtom::partition, ComputePmeMgr::pmeProxyDir, CompAtom::position, ComputePmeMgr::recipPeDest, scale_coordinates(), sendData(), sendPencils(), SubmitReduction::submit(), TRACE_COMPOBJ_IDOFFSET, ComputePmeMgr::ungrid_count, ComputePmeMgr::ungridCalc(), ComputePmeMgr::usePencils, Vector::x, PmeParticle::x, Vector::y, PmeParticle::y, Vector::z, and PmeParticle::z.

01413 {
01414   DebugM(4,"Entering ComputePme::doWork().\n");
01415 
01416 #ifdef TRACE_COMPUTE_OBJECTS
01417     double traceObjStartTime = CmiWallTimer();
01418 #endif
01419 
01420   ResizeArrayIter<PatchElem> ap(patchList);
01421 
01422   // Skip computations if nothing to do.
01423   if ( ! patchList[0].p->flags.doFullElectrostatics )
01424   {
01425     for (ap = ap.begin(); ap != ap.end(); ap++) {
01426       CompAtom *x = (*ap).positionBox->open();
01427       Results *r = (*ap).forceBox->open();
01428       (*ap).positionBox->close(&x);
01429       (*ap).forceBox->close(&r);
01430     }
01431     reduction->submit();
01432     return;
01433   }
01434 
01435   // allocate storage
01436   numLocalAtoms = 0;
01437   for (ap = ap.begin(); ap != ap.end(); ap++) {
01438     numLocalAtoms += (*ap).p->getNumAtoms();
01439   }
01440 
01441   Lattice lattice = patchList[0].p->flags.lattice;
01442 
01443   localData = new PmeParticle[numLocalAtoms*(numGrids+
01444                                         ((numGrids>1 || selfOn)?1:0))];
01445   localPartition = new char[numLocalAtoms];
01446 
01447   int g;
01448   for ( g=0; g<numGrids; ++g ) {
01449     localGridData[g] = localData + numLocalAtoms*(g+1);
01450   }
01451 
01452   // get positions and charges
01453   PmeParticle * data_ptr = localData;
01454   char * part_ptr = localPartition;
01455   const BigReal coloumb_sqrt = sqrt( COLOUMB * ComputeNonbondedUtil::scaling
01456                                 * ComputeNonbondedUtil::dielectric_1 );
01457 
01458   for (ap = ap.begin(); ap != ap.end(); ap++) {
01459 #ifdef NETWORK_PROGRESS
01460     CmiNetworkProgress();
01461 #endif
01462 
01463     CompAtom *x = (*ap).positionBox->open();
01464     // CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo();
01465     if ( patchList[0].p->flags.doMolly ) {
01466       (*ap).positionBox->close(&x);
01467       x = (*ap).avgPositionBox->open();
01468     }
01469     int numAtoms = (*ap).p->getNumAtoms();
01470 
01471     for(int i=0; i<numAtoms; ++i)
01472     {
01473       data_ptr->x = x[i].position.x;
01474       data_ptr->y = x[i].position.y;
01475       data_ptr->z = x[i].position.z;
01476       data_ptr->cg = coloumb_sqrt * x[i].charge;
01477       ++data_ptr;
01478       *part_ptr = x[i].partition;
01479       ++part_ptr;
01480     }
01481 
01482     if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
01483     else { (*ap).positionBox->close(&x); }
01484   }
01485 
01486   // copy to other grids if needed
01487   if ( ((alchFepOn || alchThermIntOn) && (!alchDecouple)) || lesOn ) {
01488     for ( g=0; g<numGrids; ++g ) {
01489 #ifdef NETWORK_PROGRESS
01490       CmiNetworkProgress();
01491 #endif
01492 
01493       PmeParticle *lgd = localGridData[g];
01494       int nga = 0;
01495       for(int i=0; i<numLocalAtoms; ++i) {
01496         if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
01497           // for FEP/TI: grid 0 gets non-alch + partition 1;
01498           // grid 1 gets non-alch + partition 2;
01499           // grid 2 (only if called for with numGrids=3) gets only non-alch
01500           lgd[nga++] = localData[i];
01501         }
01502       }
01503       numGridAtoms[g] = nga;
01504     }
01505   } else if ( (alchFepOn || alchThermIntOn) && alchDecouple) {
01506     // alchemical decoupling: four grids
01507     // g=0: partition 0 and partition 1
01508     // g=1: partition 0 and partition 2
01509     // g=2: only partition 1 atoms
01510     // g=3: only partition 2 atoms
01511     // plus one grid g=4, only partition 0, if numGrids=5
01512     for ( g=0; g<2; ++g ) {  // same as before for first 2
01513 #ifdef NETWORK_PROGRESS
01514       CmiNetworkProgress();
01515 #endif
01516 
01517       PmeParticle *lgd = localGridData[g];
01518       int nga = 0;
01519       for(int i=0; i<numLocalAtoms; ++i) {
01520         if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
01521           lgd[nga++] = localData[i];
01522         }
01523       }
01524       numGridAtoms[g] = nga;
01525     }
01526     for (g=2 ; g<4 ; ++g ) {  // only alchemical atoms for these 2
01527 #ifdef NETWORK_PROGRESS
01528       CmiNetworkProgress();
01529 #endif
01530 
01531       PmeParticle *lgd = localGridData[g];
01532       int nga = 0;
01533       for(int i=0; i<numLocalAtoms; ++i) {
01534         if ( localPartition[i] == (g-1) ) {
01535           lgd[nga++] = localData[i];
01536         }
01537       }
01538       numGridAtoms[g] = nga;
01539     }
01540     for (g=4 ; g<numGrids ; ++g ) {  // only non-alchemical atoms 
01541       // numGrids=5 only if alchElecLambdaStart > 0
01542 #ifdef NETWORK_PROGRESS
01543       CmiNetworkProgress();
01544 #endif
01545 
01546       PmeParticle *lgd = localGridData[g];
01547       int nga = 0;
01548       for(int i=0; i<numLocalAtoms; ++i) {
01549         if ( localPartition[i] == 0 ) {
01550           lgd[nga++] = localData[i];
01551         }
01552       }
01553       numGridAtoms[g] = nga;
01554     }
01555   } else if ( selfOn ) {
01556     if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 1 failed");
01557     g = 0;
01558     PmeParticle *lgd = localGridData[g];
01559     int nga = 0;
01560     for(int i=0; i<numLocalAtoms; ++i) {
01561       if ( localPartition[i] == 1 ) {
01562         lgd[nga++] = localData[i];
01563       }
01564     }
01565     numGridAtoms[g] = nga;
01566   } else if ( pairOn ) {
01567     if ( numGrids != 3 ) NAMD_bug("ComputePme::doWork assertion 2 failed");
01568     g = 0;
01569     PmeParticle *lgd = localGridData[g];
01570     int nga = 0;
01571     for(int i=0; i<numLocalAtoms; ++i) {
01572       if ( localPartition[i] == 1 || localPartition[i] == 2 ) {
01573         lgd[nga++] = localData[i];
01574       }
01575     }
01576     numGridAtoms[g] = nga;
01577     for ( g=1; g<3; ++g ) {
01578       PmeParticle *lgd = localGridData[g];
01579       int nga = 0;
01580       for(int i=0; i<numLocalAtoms; ++i) {
01581         if ( localPartition[i] == g ) {
01582           lgd[nga++] = localData[i];
01583         }
01584       }
01585       numGridAtoms[g] = nga;
01586     }
01587   } else {
01588     if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 3 failed");
01589     localGridData[0] = localData;
01590     numGridAtoms[0] = numLocalAtoms;
01591   }
01592 
01593   memset( (void*) fz_arr, 0, myGrid.K3 * sizeof(char) );
01594 
01595   // calculate self energy
01596   BigReal ewaldcof = ComputeNonbondedUtil::ewaldcof;
01597   for ( g=0; g<numGrids; ++g ) {
01598 #ifdef NETWORK_PROGRESS
01599     CmiNetworkProgress();
01600 #endif
01601 
01602     evir[g] = 0;
01603     BigReal selfEnergy = 0;
01604     data_ptr = localGridData[g];
01605     int i;
01606     for(i=0; i<numGridAtoms[g]; ++i)
01607     {
01608       selfEnergy += data_ptr->cg * data_ptr->cg;
01609       ++data_ptr;
01610     }
01611     selfEnergy *= -1. * ewaldcof / SQRT_PI;
01612     evir[g][0] += selfEnergy;
01613 
01614     double **q = q_arr + g*fsize;
01615     for (i=0; i<fsize; ++i) {
01616       if ( q[i] ) {
01617         memset( (void*) (q[i]), 0, myGrid.dim3 * sizeof(double) );
01618       }
01619     }
01620 
01621     char *f = f_arr + g*fsize;
01622     memset( (void*) f, 0, fsize * sizeof(char) );
01623     myRealSpace[g] = new PmeRealSpace(myGrid,numGridAtoms[g]);
01624     scale_coordinates(localGridData[g], numGridAtoms[g], lattice, myGrid);
01625     myRealSpace[g]->fill_charges(q, f, fz_arr, localGridData[g]);
01626   }
01627 
01628 #ifdef TRACE_COMPUTE_OBJECTS
01629     traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+this->cid, traceObjStartTime, CmiWallTimer());
01630 #endif
01631 
01632   if ( myMgr->usePencils ) {
01633     sendPencils();
01634   } else {
01635 #if 0
01636   CProxy_ComputePmeMgr pmeProxy(CkpvAccess(BOCclass_group).computePmeMgr);
01637   pmeProxy[CkMyPe()].sendGrid();
01638 #else
01639   sendData(myMgr->numGridPes,myMgr->gridPeOrder,
01640                 myMgr->recipPeDest,myMgr->gridPeMap);
01641 #endif
01642   }
01643 
01644   if ( myMgr->ungrid_count == 0 ) {
01645     myMgr->pmeProxyDir[CkMyPe()].ungridCalc();
01646   }
01647 
01648 }

void ComputePme::sendData int  ,
int *  ,
int *  ,
int * 
 

Definition at line 1824 of file ComputePme.C.

References PmeGrid::block1, PmeGrid::block2, PmeGrid::dim2, PmeGrid::dim3, PmeGridMsg::fgrid, iERROR(), iout, j, PmeGrid::K1, PmeGrid::K2, PmeGrid::K3, PmeGridMsg::lattice, PmeGridMsg::len, PME_GRID_PRIORITY, PmeGridMsg::qgrid, Compute::sequence(), PmeGridMsg::sequence, SET_PRIORITY, PmeGridMsg::sourceNode, PmeGridMsg::start, PmeGridMsg::zlist, and PmeGridMsg::zlistlen.

Referenced by doWork(), and ComputePmeMgr::sendGrid().

01825                                                                   {
01826 
01827   // iout << "Sending charge grid for " << numLocalAtoms << " atoms to FFT on " << iPE << ".\n" << endi;
01828 
01829   myGrid.block1 = ( myGrid.K1 + numRecipPes - 1 ) / numRecipPes;
01830   myGrid.block2 = ( myGrid.K2 + numRecipPes - 1 ) / numRecipPes;
01831   bsize = myGrid.block1 * myGrid.dim2 * myGrid.dim3;
01832 
01833   Lattice lattice = patchList[0].p->flags.lattice;
01834 
01835   resultsRemaining = numRecipPes;
01836 
01837   strayChargeErrors = 0;
01838 
01839   CProxy_ComputePmeMgr pmeProxy(CkpvAccess(BOCclass_group).computePmeMgr);
01840   for (int j=0; j<numRecipPes; j++) {
01841     int pe = recipPeOrder[j];  // different order
01842     int start = pe * bsize;
01843     int len = bsize;
01844     if ( start >= qsize ) { start = 0; len = 0; }
01845     if ( start + len > qsize ) { len = qsize - start; }
01846     int zdim = myGrid.dim3;
01847     int fstart = start / zdim;
01848     int flen = len / zdim;
01849     int fcount = 0;
01850     int i;
01851 
01852     int g;
01853     for ( g=0; g<numGrids; ++g ) {
01854       char *f = f_arr + fstart + g*fsize;
01855       int fcount_g = 0;
01856       for ( i=0; i<flen; ++i ) {
01857         fcount_g += ( f[i] ? 1 : 0 );
01858       }
01859       fcount += fcount_g;
01860       if ( ! recipPeDest[pe] ) {
01861         if ( fcount_g ) {
01862           ++strayChargeErrors;
01863           iout << iERROR << "Stray PME grid charges detected: "
01864                 << CkMyPe() << " sending to " << gridPeMap[pe] << " for planes";
01865           int iz = -1;
01866           for ( i=0; i<flen; ++i ) {
01867             if ( f[i] ) {
01868               int jz = (i+fstart)/myGrid.K2;
01869               if ( iz != jz ) { iout << " " << jz;  iz = jz; }
01870             }
01871           }
01872           iout << "\n" << endi;
01873         }
01874       }
01875     }
01876 
01877 #ifdef NETWORK_PROGRESS
01878     CmiNetworkProgress();
01879 #endif
01880 
01881     if ( ! recipPeDest[pe] ) continue;
01882 
01883     int zlistlen = 0;
01884     for ( i=0; i<myGrid.K3; ++i ) {
01885       if ( fz_arr[i] ) ++zlistlen;
01886     }
01887 
01888     PmeGridMsg *msg = new (fcount*zlistlen, zlistlen, flen*numGrids,
01889                                 numGrids, PRIORITY_SIZE) PmeGridMsg;
01890     msg->sourceNode = CkMyPe();
01891     msg->lattice = lattice;
01892     msg->start = fstart;
01893     msg->len = flen;
01894     msg->zlistlen = zlistlen;
01895     int *zlist = msg->zlist;
01896     zlistlen = 0;
01897     for ( i=0; i<myGrid.K3; ++i ) {
01898       if ( fz_arr[i] ) zlist[zlistlen++] = i;
01899     }
01900     float *qmsg = msg->qgrid;
01901     for ( g=0; g<numGrids; ++g ) {
01902       char *f = f_arr + fstart + g*fsize;
01903       CmiMemcpy((void*)(msg->fgrid+g*flen),(void*)f,flen*sizeof(char));
01904       double **q = q_arr + fstart + g*fsize;
01905       for ( i=0; i<flen; ++i ) {
01906         if ( f[i] ) {
01907           for ( int k=0; k<zlistlen; ++k ) {
01908             *(qmsg++) = q[i][zlist[k]];
01909           }
01910         }
01911       }
01912     }
01913 
01914     msg->sequence = sequence();
01915     SET_PRIORITY(msg,sequence(),PME_GRID_PRIORITY)
01916     pmeProxy[gridPeMap[pe]].recvGrid(msg);
01917   }
01918 
01919   for (int i=0; i<fsize; ++i) {
01920     if ( q_arr[i] ) {
01921       memset( (void*) (q_arr[i]), -1, myGrid.dim3 * sizeof(double) );
01922     }
01923   }
01924 
01925 }

void ComputePme::sendPencils  ) 
 

Definition at line 1651 of file ComputePme.C.

References PmeGrid::block1, PmeGrid::block2, PmeGrid::dim2, PmeGrid::dim3, PmeGridMsg::fgrid, PmeGridMsg::hasData, iERROR(), iout, j, PmeGrid::K1, PmeGrid::K2, PmeGrid::K3, PmeGridMsg::lattice, PmeGridMsg::len, ComputePmeMgr::numPencilsActive, ComputePmeMgr::pencilActive, PME_GRID_PRIORITY, PmeGridMsg::qgrid, ComputePmeMgr::recvGrid(), Compute::sequence(), PmeGridMsg::sequence, SET_PRIORITY, PmeGridMsg::sourceNode, PmeGridMsg::start, ComputePmeMgr::ungrid_count, ComputePmeMgr::xBlocks, ComputePmeMgr::yBlocks, ComputePmeMgr::zBlocks, PmeGridMsg::zlist, PmeGridMsg::zlistlen, and ComputePmeMgr::zPencil.

Referenced by doWork().

01651                              {
01652 
01653   // iout << "Sending charge grid for " << numLocalAtoms << " atoms to FFT on " << iPE << ".\n" << endi;
01654 
01655   int xBlocks = myMgr->xBlocks;
01656   int yBlocks = myMgr->yBlocks;
01657   int zBlocks = myMgr->zBlocks;
01658   myGrid.block1 = ( myGrid.K1 + xBlocks - 1 ) / xBlocks;
01659   myGrid.block2 = ( myGrid.K2 + yBlocks - 1 ) / yBlocks;
01660   int K1 = myGrid.K1;
01661   int K2 = myGrid.K2;
01662   int dim2 = myGrid.dim2;
01663   int dim3 = myGrid.dim3;
01664   int block1 = myGrid.block1;
01665   int block2 = myGrid.block2;
01666 
01667   Lattice lattice = patchList[0].p->flags.lattice;
01668 
01669   resultsRemaining = myMgr->numPencilsActive;
01670   const char *pencilActive = myMgr->pencilActive;
01671 
01672   strayChargeErrors = 0;
01673   // int savedMessages = 0;
01674 
01675   for (int ib=0; ib<xBlocks; ++ib) {
01676    for (int jb=0; jb<yBlocks; ++jb) {
01677     int ibegin = ib*block1;
01678     int iend = ibegin + block1;  if ( iend > K1 ) iend = K1;
01679     int jbegin = jb*block2;
01680     int jend = jbegin + block2;  if ( jend > K2 ) jend = K2;
01681     int flen = numGrids * (iend - ibegin) * (jend - jbegin);
01682     int fcount = 0;
01683 
01684     for ( int g=0; g<numGrids; ++g ) {
01685       char *f = f_arr + g*fsize;
01686       int fcount_g = 0;
01687       for ( int i=ibegin; i<iend; ++i ) {
01688        for ( int j=jbegin; j<jend; ++j ) {
01689         fcount_g += ( f[i*dim2+j] ? 1 : 0 );
01690        }
01691       }
01692       fcount += fcount_g;
01693       if ( ! pencilActive[ib*yBlocks+jb] ) {
01694         if ( fcount_g ) {
01695           ++strayChargeErrors;
01696           iout << iERROR << "Stray PME grid charges detected: "
01697                 << CkMyPe() << " sending to (x,y)";
01698           for ( int i=ibegin; i<iend; ++i ) {
01699            for ( int j=jbegin; j<jend; ++j ) {
01700             if ( f[i*dim2+j] ) { iout << " (" << i << "," << j << ")"; }
01701            }
01702           }
01703           iout << "\n" << endi;
01704         }
01705       }
01706     }
01707 
01708 #ifdef NETWORK_PROGRESS
01709     CmiNetworkProgress();
01710 #endif
01711 
01712     if ( ! pencilActive[ib*yBlocks+jb] ) continue;
01713 
01714     int zlistlen = 0;
01715     for ( int i=0; i<myGrid.K3; ++i ) {
01716       if ( fz_arr[i] ) ++zlistlen;
01717     }
01718 
01719     int hd = ( fcount? 1 : 0 );  // has data?
01720     // if ( ! hd ) ++savedMessages;
01721 
01722     PmeGridMsg *msg = new (hd*fcount*zlistlen, hd*zlistlen, hd*flen,
01723         hd*numGrids, PRIORITY_SIZE) PmeGridMsg;
01724     msg->sourceNode = CkMyPe();
01725     msg->hasData = hd;
01726     msg->lattice = lattice;
01727    if ( hd ) {
01728 #if 0
01729     msg->start = fstart;
01730     msg->len = flen;
01731 #else
01732     msg->start = -1;   // obsolete?
01733     msg->len = -1;   // obsolete?
01734 #endif
01735     msg->zlistlen = zlistlen;
01736     int *zlist = msg->zlist;
01737     zlistlen = 0;
01738     for ( int i=0; i<myGrid.K3; ++i ) {
01739       if ( fz_arr[i] ) zlist[zlistlen++] = i;
01740     }
01741     char *fmsg = msg->fgrid;
01742     float *qmsg = msg->qgrid;
01743     for ( int g=0; g<numGrids; ++g ) {
01744       char *f = f_arr + g*fsize;
01745       double **q = q_arr + g*fsize;
01746       for ( int i=ibegin; i<iend; ++i ) {
01747        for ( int j=jbegin; j<jend; ++j ) {
01748         *(fmsg++) = f[i*dim2+j];
01749         if( f[i*dim2+j] ) {
01750           for ( int k=0; k<zlistlen; ++k ) {
01751             *(qmsg++) = q[i*dim2+j][zlist[k]];
01752           }
01753         }
01754        }
01755       }
01756     }
01757    } else { // no data no reply
01758      myMgr->ungrid_count--;
01759    }
01760 
01761     msg->sequence = sequence();
01762     SET_PRIORITY(msg,sequence(),PME_GRID_PRIORITY)
01763     myMgr->zPencil(ib,jb,0).recvGrid(msg);
01764    }
01765   }
01766 
01767   // if ( savedMessages ) {
01768   //   CkPrintf("Pe %d eliminated %d PME messages\n",CkMyPe(),savedMessages);
01769   // }
01770 
01771   for (int i=0; i<fsize; ++i) {
01772     if ( q_arr[i] ) {
01773       memset( (void*) (q_arr[i]), -1, myGrid.dim3 * sizeof(double) );
01774     }
01775   }
01776 
01777 }

void ComputePme::setMgr ComputePmeMgr mgr  )  [inline]
 

Definition at line 30 of file ComputePme.h.

Referenced by ComputePmeMgr::setCompute().

00030 { myMgr = mgr; }

void ComputePme::ungridForces  ) 
 

Definition at line 1950 of file ComputePme.C.

References ADD_VECTOR_OBJECT, SimParameters::alchLambda, SimParameters::alchLambda2, ResizeArrayIter< T >::begin(), BigReal, SimParameters::commOnly, PmeRealSpace::compute_forces(), ResizeArrayIter< T >::end(), Results::f, Force, SubmitReduction::item(), Node::Object(), REDUCTION_ELECT_ENERGY_PME_TI_1, REDUCTION_ELECT_ENERGY_PME_TI_2, REDUCTION_ELECT_ENERGY_SLOW, REDUCTION_ELECT_ENERGY_SLOW_F, REDUCTION_PAIR_ELECT_FORCE, REDUCTION_STRAY_CHARGE_ERRORS, scale_forces(), Node::simParameters, simParams, SubmitReduction::submit(), Vector::x, Vector::y, and Vector::z.

Referenced by ComputePmeMgr::ungridCalc().

01950                               {
01951 
01952     SimParameters *simParams = Node::Object()->simParameters;
01953 
01954     Vector *localResults = new Vector[numLocalAtoms*
01955                                         ((numGrids>1 || selfOn)?2:1)];
01956     Vector *gridResults;
01957     if ( alchFepOn || alchThermIntOn || lesOn || selfOn || pairOn ) {
01958       for(int i=0; i<numLocalAtoms; ++i) { localResults[i] = 0.; }
01959       gridResults = localResults + numLocalAtoms;
01960     } else {
01961       gridResults = localResults;
01962     }
01963 
01964     Vector pairForce = 0.;
01965     Lattice lattice = patchList[0].p->flags.lattice;
01966     int g = 0;
01967     if(!simParams->commOnly) {
01968     for ( g=0; g<numGrids; ++g ) {
01969 #ifdef NETWORK_PROGRESS
01970       CmiNetworkProgress();
01971 #endif
01972 
01973       myRealSpace[g]->compute_forces(q_arr+g*fsize, localGridData[g], gridResults);
01974       delete myRealSpace[g];
01975       scale_forces(gridResults, numGridAtoms[g], lattice);
01976 
01977       if ( alchFepOn || alchThermIntOn ) {
01978         double scale = 1.;
01979         elecLambdaUp =  (simParams->alchLambda <= alchElecLambdaStart)? 0. : \
01980           (simParams->alchLambda - alchElecLambdaStart) / (1. - alchElecLambdaStart);
01981         elecLambdaDown =  ((1-simParams->alchLambda) <= alchElecLambdaStart)? 0. : ((1-simParams->alchLambda) - alchElecLambdaStart) / (1. - alchElecLambdaStart);
01982         if ( g == 0 ) scale = elecLambdaUp;
01983         else if ( g == 1 ) scale = elecLambdaDown;
01984         else if ( g == 2 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
01985         if (alchDecouple) {
01986           if ( g == 2 ) scale = 1 - elecLambdaUp;
01987           else if ( g == 3 ) scale = 1 - elecLambdaDown;
01988           else if ( g == 4 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
01989         }
01990         int nga = 0;
01991         if (!alchDecouple) {
01992           for(int i=0; i<numLocalAtoms; ++i) {
01993             if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
01994               // (g=2: only partition 0)
01995               localResults[i] += gridResults[nga++] * scale;
01996             }
01997           }
01998         }
01999         else {  // alchDecouple
02000           if ( g < 2 ) {
02001             for(int i=0; i<numLocalAtoms; ++i) {
02002               if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
02003                 // g = 0: partition 0 or partition 1
02004                 // g = 1: partition 0 or partition 2
02005                 localResults[i] += gridResults[nga++] * scale;
02006               }
02007             }
02008           }
02009           else {
02010             for(int i=0; i<numLocalAtoms; ++i) {
02011               if ( localPartition[i] == (g-1) || localPartition[i] == (g-4)) {
02012                 // g = 2: partition 1 only
02013                 // g = 3: partition 2 only
02014                 // g = 4: partition 0 only
02015                 localResults[i] += gridResults[nga++] * scale;
02016               }
02017             }
02018           }
02019         }
02020       } else if ( lesOn ) {
02021         double scale = 1.;
02022         if ( alchFepOn ) {
02023           if ( g == 0 ) scale = simParams->alchLambda;
02024           else if ( g == 1 ) scale = 1. - simParams->alchLambda;
02025         } else if ( lesOn ) {
02026           scale = 1.0 / (double)lesFactor;
02027         }
02028         int nga = 0;
02029         for(int i=0; i<numLocalAtoms; ++i) {
02030           if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
02031             localResults[i] += gridResults[nga++] * scale;
02032           }
02033         }
02034       } else if ( selfOn ) {
02035         PmeParticle *lgd = localGridData[g];
02036         int nga = 0;
02037         for(int i=0; i<numLocalAtoms; ++i) {
02038           if ( localPartition[i] == 1 ) {
02039             pairForce += gridResults[nga];  // should add up to almost zero
02040             localResults[i] += gridResults[nga++];
02041           }
02042         }
02043       } else if ( pairOn ) {
02044         if ( g == 0 ) {
02045           int nga = 0;
02046           for(int i=0; i<numLocalAtoms; ++i) {
02047             if ( localPartition[i] == 1 ) {
02048               pairForce += gridResults[nga];
02049             }
02050             if ( localPartition[i] == 1 || localPartition[i] == 2 ) {
02051               localResults[i] += gridResults[nga++];
02052             }
02053           }
02054         } else if ( g == 1 ) {
02055           int nga = 0;
02056           for(int i=0; i<numLocalAtoms; ++i) {
02057             if ( localPartition[i] == g ) {
02058               pairForce -= gridResults[nga];  // should add up to almost zero
02059               localResults[i] -= gridResults[nga++];
02060             }
02061           }
02062         } else {
02063           int nga = 0;
02064           for(int i=0; i<numLocalAtoms; ++i) {
02065             if ( localPartition[i] == g ) {
02066               localResults[i] -= gridResults[nga++];
02067             }
02068           }
02069         }
02070       }
02071     }
02072     }
02073     else
02074     {// cleanup for commOnly to avoid leaking
02075       for ( g=0; g<numGrids; ++g ) {delete myRealSpace[g];}
02076     }
02077     delete [] localData;
02078     delete [] localPartition;
02079 
02080     Vector *results_ptr = localResults;
02081     ResizeArrayIter<PatchElem> ap(patchList);
02082 
02083     // add in forces
02084     for (ap = ap.begin(); ap != ap.end(); ap++) {
02085       Results *r = (*ap).forceBox->open();
02086       Force *f = r->f[Results::slow];
02087       int numAtoms = (*ap).p->getNumAtoms();
02088 
02089       if ( ! strayChargeErrors && ! simParams->commOnly ) {
02090         for(int i=0; i<numAtoms; ++i) {
02091           f[i].x += results_ptr->x;
02092           f[i].y += results_ptr->y;
02093           f[i].z += results_ptr->z;
02094           ++results_ptr;
02095         }
02096       }
02097   
02098       (*ap).forceBox->close(&r);
02099     }
02100 
02101     delete [] localResults;
02102    
02103     if ( pairOn || selfOn ) {
02104         ADD_VECTOR_OBJECT(reduction,REDUCTION_PAIR_ELECT_FORCE,pairForce);
02105     }
02106 
02107     for ( g=0; g<numGrids; ++g ) {
02108       double scale = 1.;
02109       if ( alchFepOn || alchThermIntOn ) {
02110         if ( g == 0 ) scale = elecLambdaUp;
02111         else if ( g == 1 ) scale = elecLambdaDown;
02112         else if ( g == 2 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
02113         if (alchDecouple) {
02114           if ( g == 2 ) scale = 1-elecLambdaUp;
02115           else if ( g == 3 ) scale = 1-elecLambdaDown;
02116           else if ( g == 4 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
02117         }
02118       } else if ( lesOn ) {
02119         scale = 1.0 / (double)lesFactor;
02120       } else if ( pairOn ) {
02121         scale = ( g == 0 ? 1. : -1. );
02122       }
02123       reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += evir[g][0] * scale;
02124       reduction->item(REDUCTION_VIRIAL_SLOW_XX) += evir[g][1] * scale;
02125       reduction->item(REDUCTION_VIRIAL_SLOW_XY) += evir[g][2] * scale;
02126       reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += evir[g][3] * scale;
02127       reduction->item(REDUCTION_VIRIAL_SLOW_YX) += evir[g][2] * scale;
02128       reduction->item(REDUCTION_VIRIAL_SLOW_YY) += evir[g][4] * scale;
02129       reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += evir[g][5] * scale;
02130       reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += evir[g][3] * scale;
02131       reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += evir[g][5] * scale;
02132       reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += evir[g][6] * scale;
02133 
02134       double scale2 = 0.;
02135 
02136       //alchElecLambdaStart = (alchFepOn || alchThermIntOn) ? simParams->alchElecLambdaStart : 0;
02137 
02138       if (alchFepOn) {
02139         BigReal elecLambda2Up =  (simParams->alchLambda2 <= alchElecLambdaStart)? 0. : \
02140               (simParams->alchLambda2 - alchElecLambdaStart) / (1. - alchElecLambdaStart);
02141         BigReal elecLambda2Down =  ((1-simParams->alchLambda2) <= alchElecLambdaStart)? 0. : \
02142               ((1-simParams->alchLambda2) - alchElecLambdaStart) / (1. - alchElecLambdaStart);
02143         
02144         
02145         if ( g == 0 ) scale2 = elecLambda2Up;
02146         else if ( g == 1 ) scale2 = elecLambda2Down;
02147         else if ( g == 2 ) scale2 = (elecLambda2Up + elecLambda2Down - 1)*(-1);
02148         if (alchDecouple && g == 2 ) scale2 = 1 - elecLambda2Up;
02149         else if (alchDecouple && g == 3 ) scale2 = 1 - elecLambda2Down;
02150         else if (alchDecouple && g == 4 ) scale2 = (elecLambda2Up + elecLambda2Down - 1)*(-1);
02151       }
02152       reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F) += evir[g][0] * scale2;
02153       
02154       if (alchThermIntOn) {
02155         
02156         // no decoupling:
02157         // part. 1 <-> all of system except partition 2: g[0] - g[2] 
02158         // (interactions between all atoms [partition 0 OR partition 1], 
02159         // minus all [within partition 0])
02160         // U = elecLambdaUp * (U[0] - U[2])
02161         // dU/dl = U[0] - U[2];
02162         
02163         // part. 2 <-> all of system except partition 1: g[1] - g[2] 
02164         // (interactions between all atoms [partition 0 OR partition 2], 
02165         // minus all [within partition 0])
02166         // U = elecLambdaDown * (U[1] - U[2])
02167         // dU/dl = U[1] - U[2];
02168 
02169         // alchDecouple:
02170         // part. 1 <-> part. 0: g[0] - g[2] - g[4] 
02171         // (interactions between all atoms [partition 0 OR partition 1]
02172         // minus all [within partition 1] minus all [within partition 0]
02173         // U = elecLambdaUp * (U[0] - U[4]) + (1-elecLambdaUp)* U[2]
02174         // dU/dl = U[0] - U[2] - U[4];
02175 
02176         // part. 2 <-> part. 0: g[1] - g[3] - g[4] 
02177         // (interactions between all atoms [partition 0 OR partition 2]
02178         // minus all [within partition 2] minus all [within partition 0]
02179         // U = elecLambdaDown * (U[1] - U[4]) + (1-elecLambdaDown)* U[3]
02180         // dU/dl = U[1] - U[3] - U[4];
02181         
02182         
02183         if ( g == 0 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) += evir[g][0];
02184         if ( g == 1 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) += evir[g][0];
02185         if (!alchDecouple) {
02186           if ( g == 2 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) -= evir[g][0];
02187           if ( g == 2 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) -= evir[g][0];
02188         }
02189         else {  // alchDecouple
02190           if ( g == 2 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) -= evir[g][0];
02191           if ( g == 3 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) -= evir[g][0];
02192           if ( g == 4 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) -= evir[g][0];
02193           if ( g == 4 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) -= evir[g][0];
02194         }
02195       }
02196     }
02197     reduction->item(REDUCTION_STRAY_CHARGE_ERRORS) += strayChargeErrors;
02198     reduction->submit();
02199 
02200 }


The documentation for this class was generated from the following files:
Generated on Mon Nov 23 04:59:36 2009 for NAMD by  doxygen 1.3.9.1