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 1449 of file ComputePme.C.

References DebugM, SimParameters::decouple, PmeGrid::dim2, PmeGrid::dim3, SimParameters::fepElecLambdaStart, SimParameters::fepOn, 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, SimParameters::thermInt, and ReductionMgr::willSubmit().

01449                                   :
01450   ComputeHomePatches(c)
01451 {
01452   DebugM(4,"ComputePme created.\n");
01453 
01454   CProxy_ComputePmeMgr::ckLocalBranch(
01455         CkpvAccess(BOCclass_group).computePmeMgr)->setCompute(this);
01456 
01457   useAvgPositions = 1;
01458 
01459   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
01460 
01461   SimParameters *simParams = Node::Object()->simParameters;
01462 
01463   fepOn = simParams->fepOn;
01464   thermInt = simParams->thermInt;
01465   decouple = (fepOn || thermInt) && (simParams->decouple);
01466   fepElecLambdaStart = (fepOn || thermInt) ? simParams->fepElecLambdaStart : 0;
01467             
01468   if (fepOn || thermInt) {
01469     numGrids = 2;
01470     if (decouple) numGrids += 2;
01471     if (fepElecLambdaStart || thermInt) numGrids ++;
01472   }
01473   else numGrids = 1;
01474   lesOn = simParams->lesOn;
01475   if ( lesOn ) {
01476     lesFactor = simParams->lesFactor;
01477     numGrids = lesFactor;
01478   }
01479   selfOn = 0;
01480   pairOn = simParams->pairInteractionOn;
01481   if ( pairOn ) {
01482     selfOn = simParams->pairInteractionSelf;
01483     if ( selfOn ) pairOn = 0;  // make pairOn and selfOn exclusive
01484     numGrids = selfOn ? 1 : 3;
01485   }
01486 
01487   myGrid.K1 = simParams->PMEGridSizeX;
01488   myGrid.K2 = simParams->PMEGridSizeY;
01489   myGrid.K3 = simParams->PMEGridSizeZ;
01490   myGrid.order = simParams->PMEInterpOrder;
01491   myGrid.dim2 = myGrid.K2;
01492   myGrid.dim3 = 2 * (myGrid.K3/2 + 1);
01493   qsize = myGrid.K1 * myGrid.dim2 * myGrid.dim3;
01494   fsize = myGrid.K1 * myGrid.dim2;
01495   q_arr = new double*[fsize*numGrids];
01496   memset( (void*) q_arr, 0, fsize*numGrids * sizeof(double*) );
01497   f_arr = new char[fsize*numGrids];
01498   fz_arr = new char[myGrid.K3];
01499 }

ComputePme::~ComputePme  )  [virtual]
 

Definition at line 1501 of file ComputePme.C.

01502 {
01503   for (int i=0; i<fsize*numGrids; ++i) {
01504     if ( q_arr[i] ) {
01505       delete [] q_arr[i];
01506     }
01507   }
01508   delete [] q_arr;
01509   delete [] f_arr;
01510   delete [] fz_arr;
01511 }


Member Function Documentation

void ComputePme::copyPencils PmeGridMsg  ) 
 

Definition at line 1884 of file ComputePme.C.

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

Referenced by ComputePmeMgr::recvUngrid().

01884                                             {
01885 
01886   int xBlocks = myMgr->xBlocks;
01887   int yBlocks = myMgr->yBlocks;
01888   int zBlocks = myMgr->zBlocks;
01889   myGrid.block1 = ( myGrid.K1 + xBlocks - 1 ) / xBlocks;
01890   myGrid.block2 = ( myGrid.K2 + yBlocks - 1 ) / yBlocks;
01891   int K1 = myGrid.K1;
01892   int K2 = myGrid.K2;
01893   int dim2 = myGrid.dim2;
01894   int dim3 = myGrid.dim3;
01895   int block1 = myGrid.block1;
01896   int block2 = myGrid.block2;
01897 
01898   // msg->sourceNode = thisIndex.x * initdata.yBlocks + thisIndex.y;
01899   int ib = msg->sourceNode / yBlocks;
01900   int jb = msg->sourceNode % yBlocks;
01901 
01902   int ibegin = ib*block1;
01903   int iend = ibegin + block1;  if ( iend > K1 ) iend = K1;
01904   int jbegin = jb*block2;
01905   int jend = jbegin + block2;  if ( jend > K2 ) jend = K2;
01906 
01907   int zlistlen = msg->zlistlen;
01908   int *zlist = msg->zlist;
01909   float *qmsg = msg->qgrid;
01910   int g;
01911   for ( g=0; g<numGrids; ++g ) {
01912     evir[g] += msg->evir[g];
01913     char *f = f_arr + g*fsize;
01914     double **q = q_arr + g*fsize;
01915     for ( int i=ibegin; i<iend; ++i ) {
01916      for ( int j=jbegin; j<jend; ++j ) {
01917       if( f[i*dim2+j] ) {
01918         for ( int k=0; k<zlistlen; ++k ) {
01919           q[i*dim2+j][zlist[k]] = *(qmsg++);
01920         }
01921       }
01922      }
01923     }
01924   }
01925 }

void ComputePme::copyResults PmeGridMsg  ) 
 

Definition at line 2035 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().

02035                                             {
02036 
02037   int zdim = myGrid.dim3;
02038   int flen = msg->len;
02039   int fstart = msg->start;
02040   int zlistlen = msg->zlistlen;
02041   int *zlist = msg->zlist;
02042   float *qmsg = msg->qgrid;
02043   int g;
02044   for ( g=0; g<numGrids; ++g ) {
02045     evir[g] += msg->evir[g];
02046     char *f = msg->fgrid + g*flen;
02047     double **q = q_arr + fstart + g*fsize;
02048     for ( int i=0; i<flen; ++i ) {
02049       if ( f[i] ) {
02050         for ( int k=0; k<zlistlen; ++k ) {
02051           q[i][zlist[k]] = *(qmsg++);
02052         }
02053       }
02054     }
02055   }
02056 }

void ComputePme::doWork  )  [virtual]
 

Reimplemented from Compute.

Definition at line 1513 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.

01514 {
01515   DebugM(4,"Entering ComputePme::doWork().\n");
01516 
01517 #ifdef TRACE_COMPUTE_OBJECTS
01518     double traceObjStartTime = CmiWallTimer();
01519 #endif
01520 
01521   ResizeArrayIter<PatchElem> ap(patchList);
01522 
01523   // Skip computations if nothing to do.
01524   if ( ! patchList[0].p->flags.doFullElectrostatics )
01525   {
01526     for (ap = ap.begin(); ap != ap.end(); ap++) {
01527       CompAtom *x = (*ap).positionBox->open();
01528       Results *r = (*ap).forceBox->open();
01529       (*ap).positionBox->close(&x);
01530       (*ap).forceBox->close(&r);
01531     }
01532     reduction->submit();
01533     return;
01534   }
01535 
01536   // allocate storage
01537   numLocalAtoms = 0;
01538   for (ap = ap.begin(); ap != ap.end(); ap++) {
01539     numLocalAtoms += (*ap).p->getNumAtoms();
01540   }
01541 
01542   Lattice lattice = patchList[0].p->flags.lattice;
01543 
01544   localData = new PmeParticle[numLocalAtoms*(numGrids+
01545                                         ((numGrids>1 || selfOn)?1:0))];
01546   localPartition = new char[numLocalAtoms];
01547 
01548   int g;
01549   for ( g=0; g<numGrids; ++g ) {
01550     localGridData[g] = localData + numLocalAtoms*(g+1);
01551   }
01552 
01553   // get positions and charges
01554   PmeParticle * data_ptr = localData;
01555   char * part_ptr = localPartition;
01556   const BigReal coloumb_sqrt = sqrt( COLOUMB * ComputeNonbondedUtil::scaling
01557                                 * ComputeNonbondedUtil::dielectric_1 );
01558 
01559   for (ap = ap.begin(); ap != ap.end(); ap++) {
01560 #ifdef NETWORK_PROGRESS
01561     CmiNetworkProgress();
01562 #endif
01563 
01564     CompAtom *x = (*ap).positionBox->open();
01565     if ( patchList[0].p->flags.doMolly ) {
01566       (*ap).positionBox->close(&x);
01567       x = (*ap).avgPositionBox->open();
01568     }
01569     int numAtoms = (*ap).p->getNumAtoms();
01570 
01571     for(int i=0; i<numAtoms; ++i)
01572     {
01573       data_ptr->x = x[i].position.x;
01574       data_ptr->y = x[i].position.y;
01575       data_ptr->z = x[i].position.z;
01576       data_ptr->cg = coloumb_sqrt * x[i].charge;
01577       ++data_ptr;
01578       *part_ptr = x[i].partition;
01579       ++part_ptr;
01580     }
01581 
01582     if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); }
01583     else { (*ap).positionBox->close(&x); }
01584   }
01585 
01586   // copy to other grids if needed
01587   if ( ((fepOn || thermInt) && (!decouple)) || lesOn ) {
01588     for ( g=0; g<numGrids; ++g ) {
01589 #ifdef NETWORK_PROGRESS
01590       CmiNetworkProgress();
01591 #endif
01592 
01593       PmeParticle *lgd = localGridData[g];
01594       int nga = 0;
01595       for(int i=0; i<numLocalAtoms; ++i) {
01596         if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
01597           // for FEP/TI: grid 0 gets non-alch + partition 1;
01598           // grid 1 gets non-alch + partition 2;
01599           // grid 2 (only if called for with numGrids=3) gets only non-alch
01600           lgd[nga++] = localData[i];
01601         }
01602       }
01603       numGridAtoms[g] = nga;
01604     }
01605   } else if ( (fepOn || thermInt) && decouple) {
01606     // alchemical decoupling: four grids
01607     // g=0: partition 0 and partition 1
01608     // g=1: partition 0 and partition 2
01609     // g=2: only partition 1 atoms
01610     // g=3: only partition 2 atoms
01611     // plus one grid g=4, only partition 0, if numGrids=5
01612     for ( g=0; g<2; ++g ) {  // same as before for first 2
01613 #ifdef NETWORK_PROGRESS
01614       CmiNetworkProgress();
01615 #endif
01616 
01617       PmeParticle *lgd = localGridData[g];
01618       int nga = 0;
01619       for(int i=0; i<numLocalAtoms; ++i) {
01620         if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
01621           lgd[nga++] = localData[i];
01622         }
01623       }
01624       numGridAtoms[g] = nga;
01625     }
01626     for (g=2 ; g<4 ; ++g ) {  // only alchemical atoms for these 2
01627 #ifdef NETWORK_PROGRESS
01628       CmiNetworkProgress();
01629 #endif
01630 
01631       PmeParticle *lgd = localGridData[g];
01632       int nga = 0;
01633       for(int i=0; i<numLocalAtoms; ++i) {
01634         if ( localPartition[i] == (g-1) ) {
01635           lgd[nga++] = localData[i];
01636         }
01637       }
01638       numGridAtoms[g] = nga;
01639     }
01640     for (g=4 ; g<numGrids ; ++g ) {  // only non-alchemical atoms 
01641       // numGrids=5 only if fepElecLambdaStart > 0
01642 #ifdef NETWORK_PROGRESS
01643       CmiNetworkProgress();
01644 #endif
01645 
01646       PmeParticle *lgd = localGridData[g];
01647       int nga = 0;
01648       for(int i=0; i<numLocalAtoms; ++i) {
01649         if ( localPartition[i] == 0 ) {
01650           lgd[nga++] = localData[i];
01651         }
01652       }
01653       numGridAtoms[g] = nga;
01654     }
01655   } else if ( selfOn ) {
01656     if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 1 failed");
01657     g = 0;
01658     PmeParticle *lgd = localGridData[g];
01659     int nga = 0;
01660     for(int i=0; i<numLocalAtoms; ++i) {
01661       if ( localPartition[i] == 1 ) {
01662         lgd[nga++] = localData[i];
01663       }
01664     }
01665     numGridAtoms[g] = nga;
01666   } else if ( pairOn ) {
01667     if ( numGrids != 3 ) NAMD_bug("ComputePme::doWork assertion 2 failed");
01668     g = 0;
01669     PmeParticle *lgd = localGridData[g];
01670     int nga = 0;
01671     for(int i=0; i<numLocalAtoms; ++i) {
01672       if ( localPartition[i] == 1 || localPartition[i] == 2 ) {
01673         lgd[nga++] = localData[i];
01674       }
01675     }
01676     numGridAtoms[g] = nga;
01677     for ( g=1; g<3; ++g ) {
01678       PmeParticle *lgd = localGridData[g];
01679       int nga = 0;
01680       for(int i=0; i<numLocalAtoms; ++i) {
01681         if ( localPartition[i] == g ) {
01682           lgd[nga++] = localData[i];
01683         }
01684       }
01685       numGridAtoms[g] = nga;
01686     }
01687   } else {
01688     if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 3 failed");
01689     localGridData[0] = localData;
01690     numGridAtoms[0] = numLocalAtoms;
01691   }
01692 
01693   memset( (void*) fz_arr, 0, myGrid.K3 * sizeof(char) );
01694 
01695   // calculate self energy
01696   BigReal ewaldcof = ComputeNonbondedUtil::ewaldcof;
01697   for ( g=0; g<numGrids; ++g ) {
01698 #ifdef NETWORK_PROGRESS
01699     CmiNetworkProgress();
01700 #endif
01701 
01702     evir[g] = 0;
01703     BigReal selfEnergy = 0;
01704     data_ptr = localGridData[g];
01705     int i;
01706     for(i=0; i<numGridAtoms[g]; ++i)
01707     {
01708       selfEnergy += data_ptr->cg * data_ptr->cg;
01709       ++data_ptr;
01710     }
01711     selfEnergy *= -1. * ewaldcof / SQRT_PI;
01712     evir[g][0] += selfEnergy;
01713 
01714     double **q = q_arr + g*fsize;
01715     for (i=0; i<fsize; ++i) {
01716       if ( q[i] ) {
01717         memset( (void*) (q[i]), 0, myGrid.dim3 * sizeof(double) );
01718       }
01719     }
01720 
01721     char *f = f_arr + g*fsize;
01722     memset( (void*) f, 0, fsize * sizeof(char) );
01723     myRealSpace[g] = new PmeRealSpace(myGrid,numGridAtoms[g]);
01724     scale_coordinates(localGridData[g], numGridAtoms[g], lattice, myGrid);
01725     myRealSpace[g]->fill_charges(q, f, fz_arr, localGridData[g]);
01726   }
01727 
01728 #ifdef TRACE_COMPUTE_OBJECTS
01729     traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+this->cid, traceObjStartTime, CmiWallTimer());
01730 #endif
01731 
01732   if ( myMgr->usePencils ) {
01733     sendPencils();
01734   } else {
01735 #if 0
01736   CProxy_ComputePmeMgr pmeProxy(CkpvAccess(BOCclass_group).computePmeMgr);
01737 #if CHARM_VERSION > 050402
01738   pmeProxy[CkMyPe()].sendGrid();
01739 #else
01740   pmeProxy.sendGrid(CkMyPe());
01741 #endif
01742 #else
01743   sendData(myMgr->numGridPes,myMgr->gridPeOrder,
01744                 myMgr->recipPeDest,myMgr->gridPeMap);
01745 #endif
01746   }
01747 
01748   if ( myMgr->ungrid_count == 0 ) {
01749     myMgr->pmeProxyDir[CkMyPe()].ungridCalc();
01750   }
01751 
01752 }

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

Definition at line 1928 of file ComputePme.C.

References PmeGrid::block1, PmeGrid::block2, CmiMemcpy, PmeGrid::dim2, PmeGrid::dim3, PmeGridMsg::fgrid, iERROR(), iout, 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().

01929                                                                   {
01930 
01931   // iout << "Sending charge grid for " << numLocalAtoms << " atoms to FFT on " << iPE << ".\n" << endi;
01932 
01933   myGrid.block1 = ( myGrid.K1 + numRecipPes - 1 ) / numRecipPes;
01934   myGrid.block2 = ( myGrid.K2 + numRecipPes - 1 ) / numRecipPes;
01935   bsize = myGrid.block1 * myGrid.dim2 * myGrid.dim3;
01936 
01937   Lattice lattice = patchList[0].p->flags.lattice;
01938 
01939   resultsRemaining = numRecipPes;
01940 
01941   strayChargeErrors = 0;
01942 
01943   CProxy_ComputePmeMgr pmeProxy(CkpvAccess(BOCclass_group).computePmeMgr);
01944   for (int j=0; j<numRecipPes; j++) {
01945     int pe = recipPeOrder[j];  // different order
01946     int start = pe * bsize;
01947     int len = bsize;
01948     if ( start >= qsize ) { start = 0; len = 0; }
01949     if ( start + len > qsize ) { len = qsize - start; }
01950     int zdim = myGrid.dim3;
01951     int fstart = start / zdim;
01952     int flen = len / zdim;
01953     int fcount = 0;
01954     int i;
01955 
01956     int g;
01957     for ( g=0; g<numGrids; ++g ) {
01958       char *f = f_arr + fstart + g*fsize;
01959       int fcount_g = 0;
01960       for ( i=0; i<flen; ++i ) {
01961         fcount_g += ( f[i] ? 1 : 0 );
01962       }
01963       fcount += fcount_g;
01964       if ( ! recipPeDest[pe] ) {
01965         if ( fcount_g ) {
01966           ++strayChargeErrors;
01967           iout << iERROR << "Stray PME grid charges detected: "
01968                 << CkMyPe() << " sending to " << gridPeMap[pe] << " for planes";
01969           int iz = -1;
01970           for ( i=0; i<flen; ++i ) {
01971             if ( f[i] ) {
01972               int jz = (i+fstart)/myGrid.K2;
01973               if ( iz != jz ) { iout << " " << jz;  iz = jz; }
01974             }
01975           }
01976           iout << "\n" << endi;
01977         }
01978       }
01979     }
01980 
01981 #ifdef NETWORK_PROGRESS
01982     CmiNetworkProgress();
01983 #endif
01984 
01985     if ( ! recipPeDest[pe] ) continue;
01986 
01987     int zlistlen = 0;
01988     for ( i=0; i<myGrid.K3; ++i ) {
01989       if ( fz_arr[i] ) ++zlistlen;
01990     }
01991 
01992     PmeGridMsg *msg = new (fcount*zlistlen, zlistlen, flen*numGrids,
01993                                 numGrids, PRIORITY_SIZE) PmeGridMsg;
01994     msg->sourceNode = CkMyPe();
01995     msg->lattice = lattice;
01996     msg->start = fstart;
01997     msg->len = flen;
01998     msg->zlistlen = zlistlen;
01999     int *zlist = msg->zlist;
02000     zlistlen = 0;
02001     for ( i=0; i<myGrid.K3; ++i ) {
02002       if ( fz_arr[i] ) zlist[zlistlen++] = i;
02003     }
02004     float *qmsg = msg->qgrid;
02005     for ( g=0; g<numGrids; ++g ) {
02006       char *f = f_arr + fstart + g*fsize;
02007       CmiMemcpy((void*)(msg->fgrid+g*flen),(void*)f,flen*sizeof(char));
02008       double **q = q_arr + fstart + g*fsize;
02009       for ( i=0; i<flen; ++i ) {
02010         if ( f[i] ) {
02011           for ( int k=0; k<zlistlen; ++k ) {
02012             *(qmsg++) = q[i][zlist[k]];
02013           }
02014         }
02015       }
02016     }
02017 
02018     msg->sequence = sequence();
02019     SET_PRIORITY(msg,sequence(),PME_GRID_PRIORITY)
02020 #if CHARM_VERSION > 050402
02021     pmeProxy[gridPeMap[pe]].recvGrid(msg);
02022 #else
02023     pmeProxy.recvGrid(msg,gridPeMap[pe]);
02024 #endif
02025   }
02026 
02027   for (int i=0; i<fsize; ++i) {
02028     if ( q_arr[i] ) {
02029       memset( (void*) (q_arr[i]), -1, myGrid.dim3 * sizeof(double) );
02030     }
02031   }
02032 
02033 }

void ComputePme::sendPencils  ) 
 

Definition at line 1755 of file ComputePme.C.

References PmeGrid::block1, PmeGrid::block2, PmeGrid::dim2, PmeGrid::dim3, PmeGridMsg::fgrid, PmeGridMsg::hasData, iERROR(), iout, 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().

01755                              {
01756 
01757   // iout << "Sending charge grid for " << numLocalAtoms << " atoms to FFT on " << iPE << ".\n" << endi;
01758 
01759   int xBlocks = myMgr->xBlocks;
01760   int yBlocks = myMgr->yBlocks;
01761   int zBlocks = myMgr->zBlocks;
01762   myGrid.block1 = ( myGrid.K1 + xBlocks - 1 ) / xBlocks;
01763   myGrid.block2 = ( myGrid.K2 + yBlocks - 1 ) / yBlocks;
01764   int K1 = myGrid.K1;
01765   int K2 = myGrid.K2;
01766   int dim2 = myGrid.dim2;
01767   int dim3 = myGrid.dim3;
01768   int block1 = myGrid.block1;
01769   int block2 = myGrid.block2;
01770 
01771   Lattice lattice = patchList[0].p->flags.lattice;
01772 
01773   resultsRemaining = myMgr->numPencilsActive;
01774   const char *pencilActive = myMgr->pencilActive;
01775 
01776   strayChargeErrors = 0;
01777   // int savedMessages = 0;
01778 
01779   for (int ib=0; ib<xBlocks; ++ib) {
01780    for (int jb=0; jb<yBlocks; ++jb) {
01781     int ibegin = ib*block1;
01782     int iend = ibegin + block1;  if ( iend > K1 ) iend = K1;
01783     int jbegin = jb*block2;
01784     int jend = jbegin + block2;  if ( jend > K2 ) jend = K2;
01785     int flen = numGrids * (iend - ibegin) * (jend - jbegin);
01786     int fcount = 0;
01787 
01788     for ( int g=0; g<numGrids; ++g ) {
01789       char *f = f_arr + g*fsize;
01790       int fcount_g = 0;
01791       for ( int i=ibegin; i<iend; ++i ) {
01792        for ( int j=jbegin; j<jend; ++j ) {
01793         fcount_g += ( f[i*dim2+j] ? 1 : 0 );
01794        }
01795       }
01796       fcount += fcount_g;
01797       if ( ! pencilActive[ib*yBlocks+jb] ) {
01798         if ( fcount_g ) {
01799           ++strayChargeErrors;
01800           iout << iERROR << "Stray PME grid charges detected: "
01801                 << CkMyPe() << " sending to (x,y)";
01802           for ( int i=ibegin; i<iend; ++i ) {
01803            for ( int j=jbegin; j<jend; ++j ) {
01804             if ( f[i*dim2+j] ) { iout << " (" << i << "," << j << ")"; }
01805            }
01806           }
01807           iout << "\n" << endi;
01808         }
01809       }
01810     }
01811 
01812 #ifdef NETWORK_PROGRESS
01813     CmiNetworkProgress();
01814 #endif
01815 
01816     if ( ! pencilActive[ib*yBlocks+jb] ) continue;
01817 
01818     int zlistlen = 0;
01819     for ( int i=0; i<myGrid.K3; ++i ) {
01820       if ( fz_arr[i] ) ++zlistlen;
01821     }
01822 
01823     int hd = ( fcount? 1 : 0 );  // has data?
01824     // if ( ! hd ) ++savedMessages;
01825 
01826     PmeGridMsg *msg = new (hd*fcount*zlistlen, hd*zlistlen, hd*flen,
01827         hd*numGrids, PRIORITY_SIZE) PmeGridMsg;
01828     msg->sourceNode = CkMyPe();
01829     msg->hasData = hd;
01830     msg->lattice = lattice;
01831    if ( hd ) {
01832 #if 0
01833     msg->start = fstart;
01834     msg->len = flen;
01835 #else
01836     msg->start = -1;   // obsolete?
01837     msg->len = -1;   // obsolete?
01838 #endif
01839     msg->zlistlen = zlistlen;
01840     int *zlist = msg->zlist;
01841     zlistlen = 0;
01842     for ( int i=0; i<myGrid.K3; ++i ) {
01843       if ( fz_arr[i] ) zlist[zlistlen++] = i;
01844     }
01845     char *fmsg = msg->fgrid;
01846     float *qmsg = msg->qgrid;
01847     for ( int g=0; g<numGrids; ++g ) {
01848       char *f = f_arr + g*fsize;
01849       double **q = q_arr + g*fsize;
01850       for ( int i=ibegin; i<iend; ++i ) {
01851        for ( int j=jbegin; j<jend; ++j ) {
01852         *(fmsg++) = f[i*dim2+j];
01853         if( f[i*dim2+j] ) {
01854           for ( int k=0; k<zlistlen; ++k ) {
01855             *(qmsg++) = q[i*dim2+j][zlist[k]];
01856           }
01857         }
01858        }
01859       }
01860     }
01861    } else { // no data no reply
01862      myMgr->ungrid_count--;
01863    }
01864 
01865     msg->sequence = sequence();
01866     SET_PRIORITY(msg,sequence(),PME_GRID_PRIORITY)
01867     myMgr->zPencil(ib,jb,0).recvGrid(msg);
01868    }
01869   }
01870 
01871   // if ( savedMessages ) {
01872   //   CkPrintf("Pe %d eliminated %d PME messages\n",CkMyPe(),savedMessages);
01873   // }
01874 
01875   for (int i=0; i<fsize; ++i) {
01876     if ( q_arr[i] ) {
01877       memset( (void*) (q_arr[i]), -1, myGrid.dim3 * sizeof(double) );
01878     }
01879   }
01880 
01881 }

void ComputePme::setMgr ComputePmeMgr mgr  )  [inline]
 

Definition at line 33 of file ComputePme.h.

Referenced by ComputePmeMgr::setCompute().

00033 { myMgr = mgr; }

void ComputePme::ungridForces  ) 
 

Definition at line 2058 of file ComputePme.C.

References ADD_VECTOR_OBJECT, ResizeArrayIter< T >::begin(), BigReal, SimParameters::commOnly, PmeRealSpace::compute_forces(), ResizeArrayIter< T >::end(), Results::f, Force, SubmitReduction::item(), SimParameters::lambda, SimParameters::lambda2, 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().

02058                               {
02059 
02060     SimParameters *simParams = Node::Object()->simParameters;
02061 
02062     Vector *localResults = new Vector[numLocalAtoms*
02063                                         ((numGrids>1 || selfOn)?2:1)];
02064     Vector *gridResults;
02065     if ( fepOn || thermInt || lesOn || selfOn || pairOn ) {
02066       for(int i=0; i<numLocalAtoms; ++i) { localResults[i] = 0.; }
02067       gridResults = localResults + numLocalAtoms;
02068     } else {
02069       gridResults = localResults;
02070     }
02071 
02072     Vector pairForce = 0.;
02073     Lattice lattice = patchList[0].p->flags.lattice;
02074     int g = 0;
02075     if(!simParams->commOnly) {
02076     for ( g=0; g<numGrids; ++g ) {
02077 #ifdef NETWORK_PROGRESS
02078       CmiNetworkProgress();
02079 #endif
02080 
02081       myRealSpace[g]->compute_forces(q_arr+g*fsize, localGridData[g], gridResults);
02082       delete myRealSpace[g];
02083       scale_forces(gridResults, numGridAtoms[g], lattice);
02084 
02085       if ( fepOn || thermInt ) {
02086         double scale = 1.;
02087         elecLambdaUp =  (simParams->lambda <= fepElecLambdaStart)? 0. : \
02088           (simParams->lambda - fepElecLambdaStart) / (1. - fepElecLambdaStart);
02089         elecLambdaDown =  ((1-simParams->lambda) <= fepElecLambdaStart)? 0. : \
02090           ((1-simParams->lambda) - fepElecLambdaStart) / (1. - fepElecLambdaStart);
02091         if ( g == 0 ) scale = elecLambdaUp;
02092         else if ( g == 1 ) scale = elecLambdaDown;
02093         else if ( g == 2 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
02094         if (decouple) {
02095           if ( g == 2 ) scale = 1 - elecLambdaUp;
02096           else if ( g == 3 ) scale = 1 - elecLambdaDown;
02097           else if ( g == 4 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
02098         }
02099         int nga = 0;
02100         if (!decouple) {
02101           for(int i=0; i<numLocalAtoms; ++i) {
02102             if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
02103               // (g=2: only partition 0)
02104               localResults[i] += gridResults[nga++] * scale;
02105             }
02106           }
02107         }
02108         else {  // decouple
02109           if ( g < 2 ) {
02110             for(int i=0; i<numLocalAtoms; ++i) {
02111               if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
02112                 // g = 0: partition 0 or partition 1
02113                 // g = 1: partition 0 or partition 2
02114                 localResults[i] += gridResults[nga++] * scale;
02115               }
02116             }
02117           }
02118           else {
02119             for(int i=0; i<numLocalAtoms; ++i) {
02120               if ( localPartition[i] == (g-1) || localPartition[i] == (g-4)) {
02121                 // g = 2: partition 1 only
02122                 // g = 3: partition 2 only
02123                 // g = 4: partition 0 only
02124                 localResults[i] += gridResults[nga++] * scale;
02125               }
02126             }
02127           }
02128         }
02129       } else if ( lesOn ) {
02130         double scale = 1.;
02131         if ( fepOn ) {
02132           if ( g == 0 ) scale = simParams->lambda;
02133           else if ( g == 1 ) scale = 1. - simParams->lambda;
02134         } else if ( lesOn ) {
02135           scale = 1.0 / (double)lesFactor;
02136         }
02137         int nga = 0;
02138         for(int i=0; i<numLocalAtoms; ++i) {
02139           if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
02140             localResults[i] += gridResults[nga++] * scale;
02141           }
02142         }
02143       } else if ( selfOn ) {
02144         PmeParticle *lgd = localGridData[g];
02145         int nga = 0;
02146         for(int i=0; i<numLocalAtoms; ++i) {
02147           if ( localPartition[i] == 1 ) {
02148             pairForce += gridResults[nga];  // should add up to almost zero
02149             localResults[i] += gridResults[nga++];
02150           }
02151         }
02152       } else if ( pairOn ) {
02153         if ( g == 0 ) {
02154           int nga = 0;
02155           for(int i=0; i<numLocalAtoms; ++i) {
02156             if ( localPartition[i] == 1 ) {
02157               pairForce += gridResults[nga];
02158             }
02159             if ( localPartition[i] == 1 || localPartition[i] == 2 ) {
02160               localResults[i] += gridResults[nga++];
02161             }
02162           }
02163         } else if ( g == 1 ) {
02164           int nga = 0;
02165           for(int i=0; i<numLocalAtoms; ++i) {
02166             if ( localPartition[i] == g ) {
02167               pairForce -= gridResults[nga];  // should add up to almost zero
02168               localResults[i] -= gridResults[nga++];
02169             }
02170           }
02171         } else {
02172           int nga = 0;
02173           for(int i=0; i<numLocalAtoms; ++i) {
02174             if ( localPartition[i] == g ) {
02175               localResults[i] -= gridResults[nga++];
02176             }
02177           }
02178         }
02179       }
02180     }
02181     }
02182 
02183     delete [] localData;
02184     delete [] localPartition;
02185 
02186     Vector *results_ptr = localResults;
02187     ResizeArrayIter<PatchElem> ap(patchList);
02188 
02189     // add in forces
02190     for (ap = ap.begin(); ap != ap.end(); ap++) {
02191       Results *r = (*ap).forceBox->open();
02192       Force *f = r->f[Results::slow];
02193       int numAtoms = (*ap).p->getNumAtoms();
02194 
02195       if ( ! strayChargeErrors && ! simParams->commOnly ) {
02196         for(int i=0; i<numAtoms; ++i) {
02197           f[i].x += results_ptr->x;
02198           f[i].y += results_ptr->y;
02199           f[i].z += results_ptr->z;
02200           ++results_ptr;
02201         }
02202       }
02203   
02204       (*ap).forceBox->close(&r);
02205     }
02206 
02207     delete [] localResults;
02208    
02209     if ( pairOn || selfOn ) {
02210         ADD_VECTOR_OBJECT(reduction,REDUCTION_PAIR_ELECT_FORCE,pairForce);
02211     }
02212 
02213     for ( g=0; g<numGrids; ++g ) {
02214       double scale = 1.;
02215       if ( fepOn || thermInt ) {
02216         if ( g == 0 ) scale = elecLambdaUp;
02217         else if ( g == 1 ) scale = elecLambdaDown;
02218         else if ( g == 2 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
02219         if (decouple) {
02220           if ( g == 2 ) scale = 1-elecLambdaUp;
02221           else if ( g == 3 ) scale = 1-elecLambdaDown;
02222           else if ( g == 4 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
02223         }
02224       } else if ( lesOn ) {
02225         scale = 1.0 / (double)lesFactor;
02226       } else if ( pairOn ) {
02227         scale = ( g == 0 ? 1. : -1. );
02228       }
02229       reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += evir[g][0] * scale;
02230       reduction->item(REDUCTION_VIRIAL_SLOW_XX) += evir[g][1] * scale;
02231       reduction->item(REDUCTION_VIRIAL_SLOW_XY) += evir[g][2] * scale;
02232       reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += evir[g][3] * scale;
02233       reduction->item(REDUCTION_VIRIAL_SLOW_YX) += evir[g][2] * scale;
02234       reduction->item(REDUCTION_VIRIAL_SLOW_YY) += evir[g][4] * scale;
02235       reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += evir[g][5] * scale;
02236       reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += evir[g][3] * scale;
02237       reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += evir[g][5] * scale;
02238       reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += evir[g][6] * scale;
02239 
02240       double scale2 = 0.;
02241 
02242       //fepElecLambdaStart = (fepOn || thermInt) ? simParams->fepElecLambdaStart : 0;
02243 
02244       if (fepOn) {
02245         BigReal elecLambda2Up =  (simParams->lambda2 <= fepElecLambdaStart)? 0. : \
02246               (simParams->lambda2 - fepElecLambdaStart) / (1. - fepElecLambdaStart);
02247         BigReal elecLambda2Down =  ((1-simParams->lambda2) <= fepElecLambdaStart)? 0. : \
02248               ((1-simParams->lambda2) - fepElecLambdaStart) / (1. - fepElecLambdaStart);
02249         
02250         
02251         if ( g == 0 ) scale2 = elecLambda2Up;
02252         else if ( g == 1 ) scale2 = elecLambda2Down;
02253         else if ( g == 2 ) scale2 = (elecLambda2Up + elecLambda2Down - 1)*(-1);
02254         if (decouple && g == 2 ) scale2 = 1 - elecLambda2Up;
02255         else if (decouple && g == 3 ) scale2 = 1 - elecLambda2Down;
02256         else if (decouple && g == 4 ) scale2 = (elecLambda2Up + elecLambda2Down - 1)*(-1);
02257       }
02258       reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F) += evir[g][0] * scale2;
02259       
02260       if (thermInt) {
02261         
02262         // no decoupling:
02263         // part. 1 <-> all of system except partition 2: g[0] - g[2] 
02264         // (interactions between all atoms [partition 0 OR partition 1], 
02265         // minus all [within partition 0])
02266         // U = elecLambdaUp * (U[0] - U[2])
02267         // dU/dl = U[0] - U[2];
02268         
02269         // part. 2 <-> all of system except partition 1: g[1] - g[2] 
02270         // (interactions between all atoms [partition 0 OR partition 2], 
02271         // minus all [within partition 0])
02272         // U = elecLambdaDown * (U[1] - U[2])
02273         // dU/dl = U[1] - U[2];
02274 
02275         // decouple:
02276         // part. 1 <-> part. 0: g[0] - g[2] - g[4] 
02277         // (interactions between all atoms [partition 0 OR partition 1]
02278         // minus all [within partition 1] minus all [within partition 0]
02279         // U = elecLambdaUp * (U[0] - U[4]) + (1-elecLambdaUp)* U[2]
02280         // dU/dl = U[0] - U[2] - U[4];
02281 
02282         // part. 2 <-> part. 0: g[1] - g[3] - g[4] 
02283         // (interactions between all atoms [partition 0 OR partition 2]
02284         // minus all [within partition 2] minus all [within partition 0]
02285         // U = elecLambdaDown * (U[1] - U[4]) + (1-elecLambdaDown)* U[3]
02286         // dU/dl = U[1] - U[3] - U[4];
02287         
02288         
02289         if ( g == 0 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) += evir[g][0];
02290         if ( g == 1 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) += evir[g][0];
02291         if (!decouple) {
02292           if ( g == 2 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) -= evir[g][0];
02293           if ( g == 2 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) -= evir[g][0];
02294         }
02295         else {  // decouple
02296           if ( g == 2 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) -= evir[g][0];
02297           if ( g == 3 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) -= evir[g][0];
02298           if ( g == 4 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) -= evir[g][0];
02299           if ( g == 4 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) -= evir[g][0];
02300         }
02301       }
02302     }
02303     reduction->item(REDUCTION_STRAY_CHARGE_ERRORS) += strayChargeErrors;
02304     reduction->submit();
02305 
02306 }


The documentation for this class was generated from the following files:
Generated on Fri Sep 5 04:07:18 2008 for NAMD by  doxygen 1.3.9.1