ComputePme Class Reference

#include <ComputePme.h>

Inheritance diagram for ComputePme:

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(). : ComputeHomePatches(c) { DebugM(4,"ComputePme created.\n"); basePriority = PME_PRIORITY; CProxy_ComputePmeMgr::ckLocalBranch( CkpvAccess(BOCclass_group).computePmeMgr)->setCompute(this); useAvgPositions = 1; reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC); SimParameters *simParams = Node::Object()->simParameters; alchFepOn = simParams->alchFepOn; alchThermIntOn = simParams->alchThermIntOn; alchDecouple = (alchFepOn || alchThermIntOn) && (simParams->alchDecouple); alchElecLambdaStart = (alchFepOn || alchThermIntOn) ? simParams->alchElecLambdaStart : 0; if (alchFepOn || alchThermIntOn) { numGrids = 2; if (alchDecouple) numGrids += 2; if (alchElecLambdaStart || alchThermIntOn) numGrids ++; } else numGrids = 1; lesOn = simParams->lesOn; if ( lesOn ) { lesFactor = simParams->lesFactor; numGrids = lesFactor; } selfOn = 0; pairOn = simParams->pairInteractionOn; if ( pairOn ) { selfOn = simParams->pairInteractionSelf; if ( selfOn ) pairOn = 0; // make pairOn and selfOn exclusive numGrids = selfOn ? 1 : 3; } myGrid.K1 = simParams->PMEGridSizeX; myGrid.K2 = simParams->PMEGridSizeY; myGrid.K3 = simParams->PMEGridSizeZ; myGrid.order = simParams->PMEInterpOrder; myGrid.dim2 = myGrid.K2; myGrid.dim3 = 2 * (myGrid.K3/2 + 1); qsize = myGrid.K1 * myGrid.dim2 * myGrid.dim3; fsize = myGrid.K1 * myGrid.dim2; q_arr = new double*[fsize*numGrids]; memset( (void*) q_arr, 0, fsize*numGrids * sizeof(double*) ); f_arr = new char[fsize*numGrids]; fz_arr = new char[myGrid.K3]; }

ComputePme::~ComputePme (   )  [virtual]

Definition at line 1400 of file ComputePme.C. { for (int i=0; i<fsize*numGrids; ++i) { if ( q_arr[i] ) { delete [] q_arr[i]; } } delete [] q_arr; delete [] f_arr; delete [] fz_arr; }

---

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(). { int xBlocks = myMgr->xBlocks; int yBlocks = myMgr->yBlocks; int zBlocks = myMgr->zBlocks; myGrid.block1 = ( myGrid.K1 + xBlocks - 1 ) / xBlocks; myGrid.block2 = ( myGrid.K2 + yBlocks - 1 ) / yBlocks; int K1 = myGrid.K1; int K2 = myGrid.K2; int dim2 = myGrid.dim2; int dim3 = myGrid.dim3; int block1 = myGrid.block1; int block2 = myGrid.block2; // msg->sourceNode = thisIndex.x * initdata.yBlocks + thisIndex.y; int ib = msg->sourceNode / yBlocks; int jb = msg->sourceNode % yBlocks; int ibegin = ib*block1; int iend = ibegin + block1; if ( iend > K1 ) iend = K1; int jbegin = jb*block2; int jend = jbegin + block2; if ( jend > K2 ) jend = K2; int zlistlen = msg->zlistlen; int *zlist = msg->zlist; float *qmsg = msg->qgrid; int g; for ( g=0; g<numGrids; ++g ) { evir[g] += msg->evir[g]; char *f = f_arr + g*fsize; double **q = q_arr + g*fsize; for ( int i=ibegin; i<iend; ++i ) { for ( int j=jbegin; j<jend; ++j ) { if( f[i*dim2+j] ) { for ( int k=0; k<zlistlen; ++k ) { q[i*dim2+j][zlist[k]] = *(qmsg++); } } } } } }

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(). { int zdim = myGrid.dim3; int flen = msg->len; int fstart = msg->start; int zlistlen = msg->zlistlen; int *zlist = msg->zlist; float *qmsg = msg->qgrid; int g; for ( g=0; g<numGrids; ++g ) { evir[g] += msg->evir[g]; char *f = msg->fgrid + g*flen; double **q = q_arr + fstart + g*fsize; for ( int i=0; i<flen; ++i ) { if ( f[i] ) { for ( int k=0; k<zlistlen; ++k ) { q[i][zlist[k]] = *(qmsg++); } } } } }

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. { DebugM(4,"Entering ComputePme::doWork().\n"); #ifdef TRACE_COMPUTE_OBJECTS double traceObjStartTime = CmiWallTimer(); #endif ResizeArrayIter<PatchElem> ap(patchList); // Skip computations if nothing to do. if ( ! patchList[0].p->flags.doFullElectrostatics ) { for (ap = ap.begin(); ap != ap.end(); ap++) { CompAtom *x = (*ap).positionBox->open(); Results *r = (*ap).forceBox->open(); (*ap).positionBox->close(&x); (*ap).forceBox->close(&r); } reduction->submit(); return; } // allocate storage numLocalAtoms = 0; for (ap = ap.begin(); ap != ap.end(); ap++) { numLocalAtoms += (*ap).p->getNumAtoms(); } Lattice lattice = patchList[0].p->flags.lattice; localData = new PmeParticle[numLocalAtoms*(numGrids+ ((numGrids>1 || selfOn)?1:0))]; localPartition = new char[numLocalAtoms]; int g; for ( g=0; g<numGrids; ++g ) { localGridData[g] = localData + numLocalAtoms*(g+1); } // get positions and charges PmeParticle * data_ptr = localData; char * part_ptr = localPartition; const BigReal coloumb_sqrt = sqrt( COLOUMB * ComputeNonbondedUtil::scaling * ComputeNonbondedUtil::dielectric_1 ); for (ap = ap.begin(); ap != ap.end(); ap++) { #ifdef NETWORK_PROGRESS CmiNetworkProgress(); #endif CompAtom *x = (*ap).positionBox->open(); // CompAtomExt *xExt = (*ap).p->getCompAtomExtInfo(); if ( patchList[0].p->flags.doMolly ) { (*ap).positionBox->close(&x); x = (*ap).avgPositionBox->open(); } int numAtoms = (*ap).p->getNumAtoms(); for(int i=0; i<numAtoms; ++i) { data_ptr->x = x[i].position.x; data_ptr->y = x[i].position.y; data_ptr->z = x[i].position.z; data_ptr->cg = coloumb_sqrt * x[i].charge; ++data_ptr; *part_ptr = x[i].partition; ++part_ptr; } if ( patchList[0].p->flags.doMolly ) { (*ap).avgPositionBox->close(&x); } else { (*ap).positionBox->close(&x); } } // copy to other grids if needed if ( ((alchFepOn || alchThermIntOn) && (!alchDecouple)) || lesOn ) { for ( g=0; g<numGrids; ++g ) { #ifdef NETWORK_PROGRESS CmiNetworkProgress(); #endif PmeParticle *lgd = localGridData[g]; int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) { // for FEP/TI: grid 0 gets non-alch + partition 1; // grid 1 gets non-alch + partition 2; // grid 2 (only if called for with numGrids=3) gets only non-alch lgd[nga++] = localData[i]; } } numGridAtoms[g] = nga; } } else if ( (alchFepOn || alchThermIntOn) && alchDecouple) { // alchemical decoupling: four grids // g=0: partition 0 and partition 1 // g=1: partition 0 and partition 2 // g=2: only partition 1 atoms // g=3: only partition 2 atoms // plus one grid g=4, only partition 0, if numGrids=5 for ( g=0; g<2; ++g ) { // same as before for first 2 #ifdef NETWORK_PROGRESS CmiNetworkProgress(); #endif PmeParticle *lgd = localGridData[g]; int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) { lgd[nga++] = localData[i]; } } numGridAtoms[g] = nga; } for (g=2 ; g<4 ; ++g ) { // only alchemical atoms for these 2 #ifdef NETWORK_PROGRESS CmiNetworkProgress(); #endif PmeParticle *lgd = localGridData[g]; int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == (g-1) ) { lgd[nga++] = localData[i]; } } numGridAtoms[g] = nga; } for (g=4 ; g<numGrids ; ++g ) { // only non-alchemical atoms // numGrids=5 only if alchElecLambdaStart > 0 #ifdef NETWORK_PROGRESS CmiNetworkProgress(); #endif PmeParticle *lgd = localGridData[g]; int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == 0 ) { lgd[nga++] = localData[i]; } } numGridAtoms[g] = nga; } } else if ( selfOn ) { if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 1 failed"); g = 0; PmeParticle *lgd = localGridData[g]; int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == 1 ) { lgd[nga++] = localData[i]; } } numGridAtoms[g] = nga; } else if ( pairOn ) { if ( numGrids != 3 ) NAMD_bug("ComputePme::doWork assertion 2 failed"); g = 0; PmeParticle *lgd = localGridData[g]; int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == 1 || localPartition[i] == 2 ) { lgd[nga++] = localData[i]; } } numGridAtoms[g] = nga; for ( g=1; g<3; ++g ) { PmeParticle *lgd = localGridData[g]; int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == g ) { lgd[nga++] = localData[i]; } } numGridAtoms[g] = nga; } } else { if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 3 failed"); localGridData[0] = localData; numGridAtoms[0] = numLocalAtoms; } memset( (void*) fz_arr, 0, myGrid.K3 * sizeof(char) ); // calculate self energy BigReal ewaldcof = ComputeNonbondedUtil::ewaldcof; for ( g=0; g<numGrids; ++g ) { #ifdef NETWORK_PROGRESS CmiNetworkProgress(); #endif evir[g] = 0; BigReal selfEnergy = 0; data_ptr = localGridData[g]; int i; for(i=0; i<numGridAtoms[g]; ++i) { selfEnergy += data_ptr->cg * data_ptr->cg; ++data_ptr; } selfEnergy *= -1. * ewaldcof / SQRT_PI; evir[g][0] += selfEnergy; double **q = q_arr + g*fsize; for (i=0; i<fsize; ++i) { if ( q[i] ) { memset( (void*) (q[i]), 0, myGrid.dim3 * sizeof(double) ); } } char *f = f_arr + g*fsize; memset( (void*) f, 0, fsize * sizeof(char) ); myRealSpace[g] = new PmeRealSpace(myGrid,numGridAtoms[g]); scale_coordinates(localGridData[g], numGridAtoms[g], lattice, myGrid); myRealSpace[g]->fill_charges(q, f, fz_arr, localGridData[g]); } #ifdef TRACE_COMPUTE_OBJECTS traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+this->cid, traceObjStartTime, CmiWallTimer()); #endif if ( myMgr->usePencils ) { sendPencils(); } else { #if 0 CProxy_ComputePmeMgr pmeProxy(CkpvAccess(BOCclass_group).computePmeMgr); pmeProxy[CkMyPe()].sendGrid(); #else sendData(myMgr->numGridPes,myMgr->gridPeOrder, myMgr->recipPeDest,myMgr->gridPeMap); #endif } if ( myMgr->ungrid_count == 0 ) { myMgr->pmeProxyDir[CkMyPe()].ungridCalc(); } }

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(). { // iout << "Sending charge grid for " << numLocalAtoms << " atoms to FFT on " << iPE << ".\n" << endi; myGrid.block1 = ( myGrid.K1 + numRecipPes - 1 ) / numRecipPes; myGrid.block2 = ( myGrid.K2 + numRecipPes - 1 ) / numRecipPes; bsize = myGrid.block1 * myGrid.dim2 * myGrid.dim3; Lattice lattice = patchList[0].p->flags.lattice; resultsRemaining = numRecipPes; strayChargeErrors = 0; CProxy_ComputePmeMgr pmeProxy(CkpvAccess(BOCclass_group).computePmeMgr); for (int j=0; j<numRecipPes; j++) { int pe = recipPeOrder[j]; // different order int start = pe * bsize; int len = bsize; if ( start >= qsize ) { start = 0; len = 0; } if ( start + len > qsize ) { len = qsize - start; } int zdim = myGrid.dim3; int fstart = start / zdim; int flen = len / zdim; int fcount = 0; int i; int g; for ( g=0; g<numGrids; ++g ) { char *f = f_arr + fstart + g*fsize; int fcount_g = 0; for ( i=0; i<flen; ++i ) { fcount_g += ( f[i] ? 1 : 0 ); } fcount += fcount_g; if ( ! recipPeDest[pe] ) { if ( fcount_g ) { ++strayChargeErrors; iout << iERROR << "Stray PME grid charges detected: " << CkMyPe() << " sending to " << gridPeMap[pe] << " for planes"; int iz = -1; for ( i=0; i<flen; ++i ) { if ( f[i] ) { int jz = (i+fstart)/myGrid.K2; if ( iz != jz ) { iout << " " << jz; iz = jz; } } } iout << "\n" << endi; } } } #ifdef NETWORK_PROGRESS CmiNetworkProgress(); #endif if ( ! recipPeDest[pe] ) continue; int zlistlen = 0; for ( i=0; i<myGrid.K3; ++i ) { if ( fz_arr[i] ) ++zlistlen; } PmeGridMsg *msg = new (fcount*zlistlen, zlistlen, flen*numGrids, numGrids, PRIORITY_SIZE) PmeGridMsg; msg->sourceNode = CkMyPe(); msg->lattice = lattice; msg->start = fstart; msg->len = flen; msg->zlistlen = zlistlen; int *zlist = msg->zlist; zlistlen = 0; for ( i=0; i<myGrid.K3; ++i ) { if ( fz_arr[i] ) zlist[zlistlen++] = i; } float *qmsg = msg->qgrid; for ( g=0; g<numGrids; ++g ) { char *f = f_arr + fstart + g*fsize; CmiMemcpy((void*)(msg->fgrid+g*flen),(void*)f,flen*sizeof(char)); double **q = q_arr + fstart + g*fsize; for ( i=0; i<flen; ++i ) { if ( f[i] ) { for ( int k=0; k<zlistlen; ++k ) { *(qmsg++) = q[i][zlist[k]]; } } } } msg->sequence = sequence(); SET_PRIORITY(msg,sequence(),PME_GRID_PRIORITY) pmeProxy[gridPeMap[pe]].recvGrid(msg); } for (int i=0; i<fsize; ++i) { if ( q_arr[i] ) { memset( (void*) (q_arr[i]), -1, myGrid.dim3 * sizeof(double) ); } } }

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(). { // iout << "Sending charge grid for " << numLocalAtoms << " atoms to FFT on " << iPE << ".\n" << endi; int xBlocks = myMgr->xBlocks; int yBlocks = myMgr->yBlocks; int zBlocks = myMgr->zBlocks; myGrid.block1 = ( myGrid.K1 + xBlocks - 1 ) / xBlocks; myGrid.block2 = ( myGrid.K2 + yBlocks - 1 ) / yBlocks; int K1 = myGrid.K1; int K2 = myGrid.K2; int dim2 = myGrid.dim2; int dim3 = myGrid.dim3; int block1 = myGrid.block1; int block2 = myGrid.block2; Lattice lattice = patchList[0].p->flags.lattice; resultsRemaining = myMgr->numPencilsActive; const char *pencilActive = myMgr->pencilActive; strayChargeErrors = 0; // int savedMessages = 0; for (int ib=0; ib<xBlocks; ++ib) { for (int jb=0; jb<yBlocks; ++jb) { int ibegin = ib*block1; int iend = ibegin + block1; if ( iend > K1 ) iend = K1; int jbegin = jb*block2; int jend = jbegin + block2; if ( jend > K2 ) jend = K2; int flen = numGrids * (iend - ibegin) * (jend - jbegin); int fcount = 0; for ( int g=0; g<numGrids; ++g ) { char *f = f_arr + g*fsize; int fcount_g = 0; for ( int i=ibegin; i<iend; ++i ) { for ( int j=jbegin; j<jend; ++j ) { fcount_g += ( f[i*dim2+j] ? 1 : 0 ); } } fcount += fcount_g; if ( ! pencilActive[ib*yBlocks+jb] ) { if ( fcount_g ) { ++strayChargeErrors; iout << iERROR << "Stray PME grid charges detected: " << CkMyPe() << " sending to (x,y)"; for ( int i=ibegin; i<iend; ++i ) { for ( int j=jbegin; j<jend; ++j ) { if ( f[i*dim2+j] ) { iout << " (" << i << "," << j << ")"; } } } iout << "\n" << endi; } } } #ifdef NETWORK_PROGRESS CmiNetworkProgress(); #endif if ( ! pencilActive[ib*yBlocks+jb] ) continue; int zlistlen = 0; for ( int i=0; i<myGrid.K3; ++i ) { if ( fz_arr[i] ) ++zlistlen; } int hd = ( fcount? 1 : 0 ); // has data? // if ( ! hd ) ++savedMessages; PmeGridMsg *msg = new (hd*fcount*zlistlen, hd*zlistlen, hd*flen, hd*numGrids, PRIORITY_SIZE) PmeGridMsg; msg->sourceNode = CkMyPe(); msg->hasData = hd; msg->lattice = lattice; if ( hd ) { #if 0 msg->start = fstart; msg->len = flen; #else msg->start = -1; // obsolete? msg->len = -1; // obsolete? #endif msg->zlistlen = zlistlen; int *zlist = msg->zlist; zlistlen = 0; for ( int i=0; i<myGrid.K3; ++i ) { if ( fz_arr[i] ) zlist[zlistlen++] = i; } char *fmsg = msg->fgrid; float *qmsg = msg->qgrid; for ( int g=0; g<numGrids; ++g ) { char *f = f_arr + g*fsize; double **q = q_arr + g*fsize; for ( int i=ibegin; i<iend; ++i ) { for ( int j=jbegin; j<jend; ++j ) { *(fmsg++) = f[i*dim2+j]; if( f[i*dim2+j] ) { for ( int k=0; k<zlistlen; ++k ) { *(qmsg++) = q[i*dim2+j][zlist[k]]; } } } } } } else { // no data no reply myMgr->ungrid_count--; } msg->sequence = sequence(); SET_PRIORITY(msg,sequence(),PME_GRID_PRIORITY) myMgr->zPencil(ib,jb,0).recvGrid(msg); } } // if ( savedMessages ) { // CkPrintf("Pe %d eliminated %d PME messages\n",CkMyPe(),savedMessages); // } for (int i=0; i<fsize; ++i) { if ( q_arr[i] ) { memset( (void*) (q_arr[i]), -1, myGrid.dim3 * sizeof(double) ); } } }

void ComputePme::setMgr (  ComputePmeMgr *  mgr  )  [inline]

Definition at line 30 of file ComputePme.h. Referenced by ComputePmeMgr::setCompute(). { 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(). { SimParameters *simParams = Node::Object()->simParameters; Vector *localResults = new Vector[numLocalAtoms* ((numGrids>1 || selfOn)?2:1)]; Vector *gridResults; if ( alchFepOn || alchThermIntOn || lesOn || selfOn || pairOn ) { for(int i=0; i<numLocalAtoms; ++i) { localResults[i] = 0.; } gridResults = localResults + numLocalAtoms; } else { gridResults = localResults; } Vector pairForce = 0.; Lattice lattice = patchList[0].p->flags.lattice; int g = 0; if(!simParams->commOnly) { for ( g=0; g<numGrids; ++g ) { #ifdef NETWORK_PROGRESS CmiNetworkProgress(); #endif myRealSpace[g]->compute_forces(q_arr+g*fsize, localGridData[g], gridResults); delete myRealSpace[g]; scale_forces(gridResults, numGridAtoms[g], lattice); if ( alchFepOn || alchThermIntOn ) { double scale = 1.; elecLambdaUp = (simParams->alchLambda <= alchElecLambdaStart)? 0. : \ (simParams->alchLambda - alchElecLambdaStart) / (1. - alchElecLambdaStart); elecLambdaDown = ((1-simParams->alchLambda) <= alchElecLambdaStart)? 0. : ((1-simParams->alchLambda) - alchElecLambdaStart) / (1. - alchElecLambdaStart); if ( g == 0 ) scale = elecLambdaUp; else if ( g == 1 ) scale = elecLambdaDown; else if ( g == 2 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1); if (alchDecouple) { if ( g == 2 ) scale = 1 - elecLambdaUp; else if ( g == 3 ) scale = 1 - elecLambdaDown; else if ( g == 4 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1); } int nga = 0; if (!alchDecouple) { for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) { // (g=2: only partition 0) localResults[i] += gridResults[nga++] * scale; } } } else { // alchDecouple if ( g < 2 ) { for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) { // g = 0: partition 0 or partition 1 // g = 1: partition 0 or partition 2 localResults[i] += gridResults[nga++] * scale; } } } else { for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == (g-1) || localPartition[i] == (g-4)) { // g = 2: partition 1 only // g = 3: partition 2 only // g = 4: partition 0 only localResults[i] += gridResults[nga++] * scale; } } } } } else if ( lesOn ) { double scale = 1.; if ( alchFepOn ) { if ( g == 0 ) scale = simParams->alchLambda; else if ( g == 1 ) scale = 1. - simParams->alchLambda; } else if ( lesOn ) { scale = 1.0 / (double)lesFactor; } int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) { localResults[i] += gridResults[nga++] * scale; } } } else if ( selfOn ) { PmeParticle *lgd = localGridData[g]; int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == 1 ) { pairForce += gridResults[nga]; // should add up to almost zero localResults[i] += gridResults[nga++]; } } } else if ( pairOn ) { if ( g == 0 ) { int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == 1 ) { pairForce += gridResults[nga]; } if ( localPartition[i] == 1 || localPartition[i] == 2 ) { localResults[i] += gridResults[nga++]; } } } else if ( g == 1 ) { int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == g ) { pairForce -= gridResults[nga]; // should add up to almost zero localResults[i] -= gridResults[nga++]; } } } else { int nga = 0; for(int i=0; i<numLocalAtoms; ++i) { if ( localPartition[i] == g ) { localResults[i] -= gridResults[nga++]; } } } } } } else {// cleanup for commOnly to avoid leaking for ( g=0; g<numGrids; ++g ) {delete myRealSpace[g];} } delete [] localData; delete [] localPartition; Vector *results_ptr = localResults; ResizeArrayIter<PatchElem> ap(patchList); // add in forces for (ap = ap.begin(); ap != ap.end(); ap++) { Results *r = (*ap).forceBox->open(); Force *f = r->f[Results::slow]; int numAtoms = (*ap).p->getNumAtoms(); if ( ! strayChargeErrors && ! simParams->commOnly ) { for(int i=0; i<numAtoms; ++i) { f[i].x += results_ptr->x; f[i].y += results_ptr->y; f[i].z += results_ptr->z; ++results_ptr; } } (*ap).forceBox->close(&r); } delete [] localResults; if ( pairOn || selfOn ) { ADD_VECTOR_OBJECT(reduction,REDUCTION_PAIR_ELECT_FORCE,pairForce); } for ( g=0; g<numGrids; ++g ) { double scale = 1.; if ( alchFepOn || alchThermIntOn ) { if ( g == 0 ) scale = elecLambdaUp; else if ( g == 1 ) scale = elecLambdaDown; else if ( g == 2 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1); if (alchDecouple) { if ( g == 2 ) scale = 1-elecLambdaUp; else if ( g == 3 ) scale = 1-elecLambdaDown; else if ( g == 4 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1); } } else if ( lesOn ) { scale = 1.0 / (double)lesFactor; } else if ( pairOn ) { scale = ( g == 0 ? 1. : -1. ); } reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += evir[g][0] * scale; reduction->item(REDUCTION_VIRIAL_SLOW_XX) += evir[g][1] * scale; reduction->item(REDUCTION_VIRIAL_SLOW_XY) += evir[g][2] * scale; reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += evir[g][3] * scale; reduction->item(REDUCTION_VIRIAL_SLOW_YX) += evir[g][2] * scale; reduction->item(REDUCTION_VIRIAL_SLOW_YY) += evir[g][4] * scale; reduction->item(REDUCTION_VIRIAL_SLOW_YZ) += evir[g][5] * scale; reduction->item(REDUCTION_VIRIAL_SLOW_ZX) += evir[g][3] * scale; reduction->item(REDUCTION_VIRIAL_SLOW_ZY) += evir[g][5] * scale; reduction->item(REDUCTION_VIRIAL_SLOW_ZZ) += evir[g][6] * scale; double scale2 = 0.; //alchElecLambdaStart = (alchFepOn || alchThermIntOn) ? simParams->alchElecLambdaStart : 0; if (alchFepOn) { BigReal elecLambda2Up = (simParams->alchLambda2 <= alchElecLambdaStart)? 0. : \ (simParams->alchLambda2 - alchElecLambdaStart) / (1. - alchElecLambdaStart); BigReal elecLambda2Down = ((1-simParams->alchLambda2) <= alchElecLambdaStart)? 0. : \ ((1-simParams->alchLambda2) - alchElecLambdaStart) / (1. - alchElecLambdaStart); if ( g == 0 ) scale2 = elecLambda2Up; else if ( g == 1 ) scale2 = elecLambda2Down; else if ( g == 2 ) scale2 = (elecLambda2Up + elecLambda2Down - 1)*(-1); if (alchDecouple && g == 2 ) scale2 = 1 - elecLambda2Up; else if (alchDecouple && g == 3 ) scale2 = 1 - elecLambda2Down; else if (alchDecouple && g == 4 ) scale2 = (elecLambda2Up + elecLambda2Down - 1)*(-1); } reduction->item(REDUCTION_ELECT_ENERGY_SLOW_F) += evir[g][0] * scale2; if (alchThermIntOn) { // no decoupling: // part. 1 <-> all of system except partition 2: g[0] - g[2] // (interactions between all atoms [partition 0 OR partition 1], // minus all [within partition 0]) // U = elecLambdaUp * (U[0] - U[2]) // dU/dl = U[0] - U[2]; // part. 2 <-> all of system except partition 1: g[1] - g[2] // (interactions between all atoms [partition 0 OR partition 2], // minus all [within partition 0]) // U = elecLambdaDown * (U[1] - U[2]) // dU/dl = U[1] - U[2]; // alchDecouple: // part. 1 <-> part. 0: g[0] - g[2] - g[4] // (interactions between all atoms [partition 0 OR partition 1] // minus all [within partition 1] minus all [within partition 0] // U = elecLambdaUp * (U[0] - U[4]) + (1-elecLambdaUp)* U[2] // dU/dl = U[0] - U[2] - U[4]; // part. 2 <-> part. 0: g[1] - g[3] - g[4] // (interactions between all atoms [partition 0 OR partition 2] // minus all [within partition 2] minus all [within partition 0] // U = elecLambdaDown * (U[1] - U[4]) + (1-elecLambdaDown)* U[3] // dU/dl = U[1] - U[3] - U[4]; if ( g == 0 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) += evir[g][0]; if ( g == 1 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) += evir[g][0]; if (!alchDecouple) { if ( g == 2 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) -= evir[g][0]; if ( g == 2 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) -= evir[g][0]; } else { // alchDecouple if ( g == 2 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) -= evir[g][0]; if ( g == 3 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) -= evir[g][0]; if ( g == 4 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_1) -= evir[g][0]; if ( g == 4 ) reduction->item(REDUCTION_ELECT_ENERGY_PME_TI_2) -= evir[g][0]; } } } reduction->item(REDUCTION_STRAY_CHARGE_ERRORS) += strayChargeErrors; reduction->submit(); }

---

The documentation for this class was generated from the following files:

• ComputePme.h
• ComputePme.C

---