ComputePme Class Reference

#include <ComputePme.h>

Inheritance diagram for ComputePme:

Compute List of all members.

Public Member Functions

 ComputePme (ComputeID c, PatchID pid)
virtual ~ComputePme ()
void initialize ()
void atomUpdate ()
int noWork ()
void doWork ()
void doQMWork ()
void ungridForces ()
void setMgr (ComputePmeMgr *mgr)

Friends

class ComputePmeMgr

Detailed Description

Definition at line 24 of file ComputePme.h.


Constructor & Destructor Documentation

ComputePme::ComputePme ( ComputeID  c,
PatchID  pid 
)

Definition at line 2667 of file ComputePme.C.

References Compute::basePriority, DebugM, PmeGrid::dim2, PmeGrid::dim3, PmeGrid::K1, PmeGrid::K2, PmeGrid::K3, Node::Object(), PmeGrid::order, PME_PRIORITY, Compute::setNumPatches(), Node::simParameters, and simParams.

02667                                                : Compute(c), patchID(pid)
02668 {
02669   DebugM(4,"ComputePme created.\n");
02670   basePriority = PME_PRIORITY;
02671   setNumPatches(1);
02672 
02673   CProxy_ComputePmeMgr::ckLocalBranch(
02674         CkpvAccess(BOCclass_group).computePmeMgr)->addCompute(this);
02675 
02676   SimParameters *simParams = Node::Object()->simParameters;
02677 
02678   qmForcesOn =  simParams->qmForcesOn;
02679   offload = simParams->PMEOffload;
02680 
02681   alchOn = simParams->alchOn;
02682   alchFepOn = simParams->alchFepOn;
02683   alchThermIntOn = simParams->alchThermIntOn;
02684   alchDecouple = alchOn && simParams->alchDecouple;
02685   alchElecLambdaStart = alchOn ? simParams->alchElecLambdaStart : 0;
02686             
02687   if (alchOn) {
02688     numGrids = 2;
02689     if (alchDecouple) numGrids += 2;
02690     if (alchElecLambdaStart || alchThermIntOn) numGrids ++;
02691   }
02692   else numGrids = 1;
02693   lesOn = simParams->lesOn;
02694   if ( lesOn ) {
02695     lesFactor = simParams->lesFactor;
02696     numGrids = lesFactor;
02697   }
02698   selfOn = 0;
02699   pairOn = simParams->pairInteractionOn;
02700   if ( pairOn ) {
02701     selfOn = simParams->pairInteractionSelf;
02702     if ( selfOn ) pairOn = 0;  // make pairOn and selfOn exclusive
02703     numGrids = selfOn ? 1 : 3;
02704   }
02705 
02706   myGrid.K1 = simParams->PMEGridSizeX;
02707   myGrid.K2 = simParams->PMEGridSizeY;
02708   myGrid.K3 = simParams->PMEGridSizeZ;
02709   myGrid.order = simParams->PMEInterpOrder;
02710   myGrid.dim2 = myGrid.K2;
02711   myGrid.dim3 = 2 * (myGrid.K3/2 + 1);
02712 
02713 #ifdef NAMD_CUDA
02714   cuda_atoms_offset = 0;
02715   f_data_host = 0;
02716   f_data_dev = 0;
02717  if ( ! offload )
02718 #endif
02719  {
02720   for ( int g=0; g<numGrids; ++g ) myRealSpace[g] = new PmeRealSpace(myGrid);
02721  }
02722 
02723   atomsChanged = 0;
02724   
02725   qmLoclIndx = 0;
02726   qmLocalCharges = 0;
02727 }

ComputePme::~ComputePme (  )  [virtual]

Definition at line 2963 of file ComputePme.C.

02964 {
02965 #ifdef NAMD_CUDA
02966   if ( ! offload )
02967 #endif
02968   {
02969     for ( int g=0; g<numGrids; ++g ) delete myRealSpace[g];
02970   }
02971 }


Member Function Documentation

void ComputePme::atomUpdate (  )  [virtual]

Reimplemented from Compute.

Definition at line 2665 of file ComputePme.C.

02665 { atomsChanged = 1; }

void ComputePme::doQMWork (  ) 

Definition at line 3056 of file ComputePme.C.

References doWork(), Molecule::get_numQMAtoms(), Molecule::get_qmAtmChrg(), Molecule::get_qmAtmIndx(), Molecule::get_qmAtomGroup(), Patch::getCompAtomExtInfo(), Patch::getNumAtoms(), Node::molecule, and Node::Object().

03056                           {
03057     
03058 //     iout << CkMyPe() << ") ----> PME doQMWork.\n" << endi ;
03059     
03060     
03061     int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
03062     const Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
03063     const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
03064     const Real *qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
03065     
03066     const CompAtomExt *xExt = patch->getCompAtomExtInfo();
03067     
03068     // Determine number of qm atoms in this patch for the current step.
03069     numLocalQMAtoms = 0;
03070     for (int paIter=0; paIter<patch->getNumAtoms(); paIter++) {
03071         if ( qmAtomGroup[xExt[paIter].id] != 0 ) {
03072             numLocalQMAtoms++;
03073         }
03074     }
03075     
03076     // We prepare a charge vector with QM charges for use in the PME calculation.
03077     
03078     // Clears data from last step, if there is any.
03079     if (qmLoclIndx != 0)
03080         delete [] qmLoclIndx;
03081     if (qmLocalCharges != 0)
03082         delete [] qmLocalCharges;
03083     
03084     qmLoclIndx = new int[numLocalQMAtoms] ;
03085     qmLocalCharges = new Real[numLocalQMAtoms] ;
03086     
03087     // I am assuming there will be (in general) more QM atoms among all QM groups
03088     // than MM atoms in a patch.
03089     int procAtms = 0;
03090     
03091     for (int paIter=0; paIter<patch->getNumAtoms(); paIter++) {
03092         
03093         for (int i=0; i<numQMAtms; i++) {
03094             
03095             if (qmAtmIndx[i] == xExt[paIter].id) {
03096                 
03097                 qmLoclIndx[procAtms] = paIter ;
03098                 qmLocalCharges[procAtms] = qmAtmChrg[i];
03099                 
03100                 procAtms++;
03101                 break;
03102             }
03103             
03104         }
03105         
03106         if (procAtms == numLocalQMAtoms)
03107             break;
03108     }
03109     
03110     doWork();
03111     return ;
03112 }

void ComputePme::doWork (  )  [virtual]

Reimplemented from Compute.

Definition at line 3114 of file ComputePme.C.

References Compute::basePriority, ResizeArray< Elem >::begin(), PmeParticle::cg, Box< Owner, Data >::close(), COMPUTE_HOME_PRIORITY, COULOMB, DebugM, deviceCUDA, ComputeNonbondedUtil::dielectric_1, Flags::doMolly, Patch::flags, DeviceCUDA::getDeviceID(), Patch::getNumAtoms(), Flags::lattice, Box< Owner, Data >::open(), PATCH_PRIORITY, PME_OFFLOAD_PRIORITY, PME_PRIORITY, ComputePmeMgr::recipEvirCount, ResizeArray< Elem >::resize(), ComputeNonbondedUtil::scaling, ComputePmeMgr::submitReductions(), ungridForces(), ComputePmeMgr::ungridForcesCount, PmeParticle::x, x, PmeParticle::y, and PmeParticle::z.

Referenced by doQMWork().

03115 {
03116   DebugM(4,"Entering ComputePme::doWork().\n");
03117 
03118   if ( basePriority >= COMPUTE_HOME_PRIORITY ) {
03119 #ifdef NAMD_CUDA
03120     basePriority = ( offload ? PME_OFFLOAD_PRIORITY : PME_PRIORITY );
03121 #else
03122     basePriority = PME_PRIORITY;
03123 #endif
03124     ungridForces();
03125     // CkPrintf("doWork 2 pe %d %d %d\n", CkMyPe(), myMgr->ungridForcesCount, myMgr->recipEvirCount);
03126     if ( ! --(myMgr->ungridForcesCount) && ! myMgr->recipEvirCount ) myMgr->submitReductions();
03127     return;
03128   }
03129   basePriority = COMPUTE_HOME_PRIORITY + PATCH_PRIORITY(patchID);
03130   // CkPrintf("doWork 1 pe %d %d %d\n", CkMyPe(), myMgr->ungridForcesCount, myMgr->recipEvirCount);
03131 
03132 #ifdef TRACE_COMPUTE_OBJECTS
03133     double traceObjStartTime = CmiWallTimer();
03134 #endif
03135 
03136 #ifdef NAMD_CUDA
03137   if ( offload ) cudaSetDevice(deviceCUDA->getDeviceID());
03138 #endif
03139 
03140   // allocate storage
03141   numLocalAtoms = patch->getNumAtoms();
03142 
03143   Lattice &lattice = patch->flags.lattice;
03144 
03145   localData_alloc.resize(numLocalAtoms*(numGrids+ ((numGrids>1 || selfOn)?1:0)));
03146   localData = localData_alloc.begin();
03147   localPartition_alloc.resize(numLocalAtoms);
03148   localPartition = localPartition_alloc.begin();
03149 
03150   int g;
03151   for ( g=0; g<numGrids; ++g ) {
03152     localGridData[g] = localData + numLocalAtoms*(g+1);
03153   }
03154 
03155   // get positions and charges
03156   PmeParticle * data_ptr = localData;
03157   unsigned char * part_ptr = localPartition;
03158   const BigReal coulomb_sqrt = sqrt( COULOMB * ComputeNonbondedUtil::scaling
03159                                 * ComputeNonbondedUtil::dielectric_1 );
03160 
03161   {
03162     CompAtom *x = positionBox->open();
03163     // CompAtomExt *xExt = patch->getCompAtomExtInfo();
03164     if ( patch->flags.doMolly ) {
03165       positionBox->close(&x);
03166       x = avgPositionBox->open();
03167     }
03168     int numAtoms = patch->getNumAtoms();
03169 
03170     for(int i=0; i<numAtoms; ++i)
03171     {
03172       data_ptr->x = x[i].position.x;
03173       data_ptr->y = x[i].position.y;
03174       data_ptr->z = x[i].position.z;
03175       data_ptr->cg = coulomb_sqrt * x[i].charge;
03176       ++data_ptr;
03177       *part_ptr = x[i].partition;
03178       ++part_ptr;
03179     }
03180 
03181     // QM loop to overwrite charges of QM atoms.
03182     // They are zero for NAMD, but are updated in ComputeQM.
03183     if ( qmForcesOn ) {
03184         
03185         for(int i=0; i<numLocalQMAtoms; ++i)
03186         {
03187           localData[qmLoclIndx[i]].cg = coulomb_sqrt * qmLocalCharges[i];
03188         }
03189         
03190     }
03191     
03192     if ( patch->flags.doMolly ) { avgPositionBox->close(&x); }
03193     else { positionBox->close(&x); }
03194   }
03195 
03196   // copy to other grids if needed
03197   if ( (alchOn && (!alchDecouple)) || lesOn ) {
03198     for ( g=0; g<numGrids; ++g ) {
03199       PmeParticle *lgd = localGridData[g];
03200       int nga = 0;
03201       for(int i=0; i<numLocalAtoms; ++i) {
03202         if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
03203           // for FEP/TI: grid 0 gets non-alch + partition 1;
03204           // grid 1 gets non-alch + partition 2;
03205           // grid 2 (only if called for with numGrids=3) gets only non-alch
03206           lgd[nga++] = localData[i];
03207         }
03208       }
03209       numGridAtoms[g] = nga;
03210     }
03211   } else if ( alchOn && alchDecouple) {
03212     // alchemical decoupling: four grids
03213     // g=0: partition 0 and partition 1
03214     // g=1: partition 0 and partition 2
03215     // g=2: only partition 1 atoms
03216     // g=3: only partition 2 atoms
03217     // plus one grid g=4, only partition 0, if numGrids=5
03218     for ( g=0; g<2; ++g ) {  // same as before for first 2
03219       PmeParticle *lgd = localGridData[g];
03220       int nga = 0;
03221       for(int i=0; i<numLocalAtoms; ++i) {
03222         if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
03223           lgd[nga++] = localData[i];
03224         }
03225       }
03226       numGridAtoms[g] = nga;
03227     }
03228     for (g=2 ; g<4 ; ++g ) {  // only alchemical atoms for these 2
03229       PmeParticle *lgd = localGridData[g];
03230       int nga = 0;
03231       for(int i=0; i<numLocalAtoms; ++i) {
03232         if ( localPartition[i] == (g-1) ) {
03233           lgd[nga++] = localData[i];
03234         }
03235       }
03236       numGridAtoms[g] = nga;
03237     }
03238     for (g=4 ; g<numGrids ; ++g ) {  // only non-alchemical atoms 
03239       // numGrids=5 only if alchElecLambdaStart > 0
03240       PmeParticle *lgd = localGridData[g];
03241       int nga = 0;
03242       for(int i=0; i<numLocalAtoms; ++i) {
03243         if ( localPartition[i] == 0 ) {
03244           lgd[nga++] = localData[i];
03245         }
03246       }
03247       numGridAtoms[g] = nga;
03248     }
03249   } else if ( selfOn ) {
03250     if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 1 failed");
03251     g = 0;
03252     PmeParticle *lgd = localGridData[g];
03253     int nga = 0;
03254     for(int i=0; i<numLocalAtoms; ++i) {
03255       if ( localPartition[i] == 1 ) {
03256         lgd[nga++] = localData[i];
03257       }
03258     }
03259     numGridAtoms[g] = nga;
03260   } else if ( pairOn ) {
03261     if ( numGrids != 3 ) NAMD_bug("ComputePme::doWork assertion 2 failed");
03262     g = 0;
03263     PmeParticle *lgd = localGridData[g];
03264     int nga = 0;
03265     for(int i=0; i<numLocalAtoms; ++i) {
03266       if ( localPartition[i] == 1 || localPartition[i] == 2 ) {
03267         lgd[nga++] = localData[i];
03268       }
03269     }
03270     numGridAtoms[g] = nga;
03271     for ( g=1; g<3; ++g ) {
03272       PmeParticle *lgd = localGridData[g];
03273       int nga = 0;
03274       for(int i=0; i<numLocalAtoms; ++i) {
03275         if ( localPartition[i] == g ) {
03276           lgd[nga++] = localData[i];
03277         }
03278       }
03279       numGridAtoms[g] = nga;
03280     }
03281   } else {
03282     if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 3 failed");
03283     localGridData[0] = localData;
03284     numGridAtoms[0] = numLocalAtoms;
03285   }
03286 
03287  if ( ! myMgr->doWorkCount ) {
03288   myMgr->doWorkCount = myMgr->pmeComputes.size();
03289 
03290 #ifdef NAMD_CUDA
03291  if ( !  offload )
03292 #endif // NAMD_CUDA
03293  {
03294   memset( (void*) myMgr->fz_arr, 0, (myGrid.K3+myGrid.order-1) * sizeof(char) );
03295 
03296   for (int i=0; i<myMgr->q_count; ++i) {
03297     memset( (void*) (myMgr->q_list[i]), 0, (myGrid.K3+myGrid.order-1) * sizeof(float) );
03298   }
03299  }
03300 
03301   for ( g=0; g<numGrids; ++g ) {
03302     myMgr->evir[g] = 0;
03303   }
03304 
03305   myMgr->strayChargeErrors = 0;
03306 
03307   myMgr->compute_sequence = sequence();
03308  }
03309 
03310   if ( sequence() != myMgr->compute_sequence ) NAMD_bug("ComputePme sequence mismatch in doWork()");
03311 
03312   int strayChargeErrors = 0;
03313 
03314   // calculate self energy
03315   BigReal ewaldcof = ComputeNonbondedUtil::ewaldcof;
03316   for ( g=0; g<numGrids; ++g ) {
03317     BigReal selfEnergy = 0;
03318     data_ptr = localGridData[g];
03319     int i;
03320     for(i=0; i<numGridAtoms[g]; ++i)
03321     {
03322       selfEnergy += data_ptr->cg * data_ptr->cg;
03323       ++data_ptr;
03324     }
03325     selfEnergy *= -1. * ewaldcof / SQRT_PI;
03326     myMgr->evir[g][0] += selfEnergy;
03327 
03328     float **q = myMgr->q_arr + g*myMgr->fsize;
03329     char *f = myMgr->f_arr + g*myMgr->fsize;
03330 
03331     scale_coordinates(localGridData[g], numGridAtoms[g], lattice, myGrid);
03332 #ifdef NAMD_CUDA
03333    if ( offload ) {
03334     if ( myMgr->cuda_atoms_alloc == 0 ) {  // first call
03335       int na = myMgr->cuda_atoms_alloc = 1.2 * (myMgr->cuda_atoms_count + 1000);
03336       cuda_errcheck("before malloc atom data for pme");
03337       cudaMallocHost((void**) &(myMgr->a_data_host), 7*na*sizeof(float));
03338       cudaMalloc((void**) &(myMgr->a_data_dev), 7*na*sizeof(float));
03339       cuda_errcheck("malloc atom data for pme");
03340       myMgr->cuda_atoms_count = 0;
03341     }
03342     cuda_atoms_offset = myMgr->cuda_atoms_count;
03343     int n = numGridAtoms[g];
03344     myMgr->cuda_atoms_count += n;
03345     if ( myMgr->cuda_atoms_count > myMgr->cuda_atoms_alloc ) {
03346       CkPrintf("Pe %d expanding CUDA PME atoms allocation because %d > %d\n",
03347                         CkMyPe(), myMgr->cuda_atoms_count, myMgr->cuda_atoms_alloc);
03348       cuda_errcheck("before malloc expanded atom data for pme");
03349       int na = myMgr->cuda_atoms_alloc = 1.2 * (myMgr->cuda_atoms_count + 1000);
03350       const float *a_data_host_old = myMgr->a_data_host;
03351       cudaMallocHost((void**) &(myMgr->a_data_host), 7*na*sizeof(float));
03352       cuda_errcheck("malloc expanded host atom data for pme");
03353       memcpy(myMgr->a_data_host, a_data_host_old, 7*cuda_atoms_offset*sizeof(float));
03354       cudaFreeHost((void*) a_data_host_old);
03355       cuda_errcheck("free expanded host atom data for pme");
03356       cudaFree(myMgr->a_data_dev);
03357       cuda_errcheck("free expanded dev atom data for pme");
03358       cudaMalloc((void**) &(myMgr->a_data_dev), 7*na*sizeof(float));
03359       cuda_errcheck("malloc expanded dev atom data for pme");
03360     }
03361     float *a_data_host = myMgr->a_data_host + 7 * cuda_atoms_offset;
03362     data_ptr = localGridData[g];
03363     double order_1 = myGrid.order - 1;
03364     double K1 = myGrid.K1;
03365     double K2 = myGrid.K2;
03366     double K3 = myGrid.K3;
03367     int found_negative = 0;
03368     for ( int i=0; i<n; ++i ) {
03369       if ( data_ptr[i].x < 0 || data_ptr[i].y < 0 || data_ptr[i].z < 0 ) {
03370         found_negative = 1;
03371         // CkPrintf("low coord: %f %f %f\n", data_ptr[i].x, data_ptr[i].y, data_ptr[i].z);
03372       }
03373       double x_int = (int) data_ptr[i].x;
03374       double y_int = (int) data_ptr[i].y;
03375       double z_int = (int) data_ptr[i].z;
03376       a_data_host[7*i  ] = data_ptr[i].x - x_int;  // subtract in double precision
03377       a_data_host[7*i+1] = data_ptr[i].y - y_int;
03378       a_data_host[7*i+2] = data_ptr[i].z - z_int;
03379       a_data_host[7*i+3] = data_ptr[i].cg;
03380       x_int -= order_1;  if ( x_int < 0 ) x_int += K1;
03381       y_int -= order_1;  if ( y_int < 0 ) y_int += K2;
03382       z_int -= order_1;  if ( z_int < 0 ) z_int += K3;
03383       a_data_host[7*i+4] = x_int;
03384       a_data_host[7*i+5] = y_int;
03385       a_data_host[7*i+6] = z_int;
03386     }
03387     if ( found_negative ) NAMD_bug("found negative atom coordinate in ComputePme::doWork");
03388    } else
03389 #endif // NAMD_CUDA
03390    {
03391     myRealSpace[g]->set_num_atoms(numGridAtoms[g]);
03392     myRealSpace[g]->fill_charges(q, myMgr->q_list, myMgr->q_count, strayChargeErrors, f, myMgr->fz_arr, localGridData[g]);
03393    }
03394   }
03395   myMgr->strayChargeErrors += strayChargeErrors;
03396 
03397 #ifdef TRACE_COMPUTE_OBJECTS
03398     traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+this->cid, traceObjStartTime, CmiWallTimer());
03399 #endif
03400 
03401  if ( --(myMgr->doWorkCount) == 0 ) {
03402 // cudaDeviceSynchronize();  // XXXX
03403 #ifdef NAMD_CUDA
03404   if ( offload ) {
03405     ComputePmeMgr::cuda_submit_charges_args args;
03406     args.mgr = myMgr;
03407     args.lattice = &lattice;
03408     args.sequence = sequence();
03409     CmiLock(ComputePmeMgr::cuda_lock);
03410     if ( ComputePmeMgr::cuda_busy ) {
03411       ComputePmeMgr::cuda_submit_charges_deque.push_back(args);
03412     } else if ( CkMyPe() == deviceCUDA->getMasterPe() ) {
03413       // avoid adding work to nonbonded data preparation pe
03414       args.mgr->cuda_submit_charges(*args.lattice, args.sequence);
03415     } else {
03416       ComputePmeMgr::cuda_busy = true;
03417       while ( 1 ) {
03418         CmiUnlock(ComputePmeMgr::cuda_lock);
03419         args.mgr->cuda_submit_charges(*args.lattice, args.sequence);
03420         CmiLock(ComputePmeMgr::cuda_lock);
03421         if ( ComputePmeMgr::cuda_submit_charges_deque.size() ) {
03422           args = ComputePmeMgr::cuda_submit_charges_deque.front();
03423           ComputePmeMgr::cuda_submit_charges_deque.pop_front();
03424         } else {
03425           ComputePmeMgr::cuda_busy = false;
03426           break;
03427         }
03428       }
03429     }
03430     CmiUnlock(ComputePmeMgr::cuda_lock);
03431   } else
03432 #endif // NAMD_CUDA
03433   {
03434     myMgr->chargeGridReady(lattice,sequence());
03435   }
03436  }
03437  atomsChanged = 0;
03438 }

void ComputePme::initialize (  )  [virtual]

Reimplemented from Compute.

Definition at line 2729 of file ComputePme.C.

References ComputePmeMgr::cuda_atoms_count, Patch::getNumAtoms(), NAMD_bug(), PatchMap::Object(), Patch::registerAvgPositionPickup(), Patch::registerForceDeposit(), and Patch::registerPositionPickup().

02729                             {
02730   if (!(patch = PatchMap::Object()->patch(patchID))) {
02731     NAMD_bug("ComputePme used with unknown patch.");
02732   }
02733   positionBox = patch->registerPositionPickup(this);
02734   avgPositionBox = patch->registerAvgPositionPickup(this);
02735   forceBox = patch->registerForceDeposit(this);
02736 #ifdef NAMD_CUDA
02737  if ( offload ) {
02738   myMgr->cuda_atoms_count += patch->getNumAtoms();
02739  }
02740 #endif
02741 }

int ComputePme::noWork (  )  [virtual]

Reimplemented from Compute.

Definition at line 3017 of file ComputePme.C.

References ResizeArray< Elem >::add(), Flags::doFullElectrostatics, Patch::flags, ComputePmeMgr::heldComputes, ComputePmeMgr::noWorkCount, ComputePmeMgr::pmeComputes, ComputePmeMgr::recipEvirCount, ComputePmeMgr::reduction, ResizeArray< Elem >::size(), Box< Owner, Data >::skip(), SubmitReduction::submit(), and ComputePmeMgr::ungridForcesCount.

03017                        {
03018 
03019   if ( patch->flags.doFullElectrostatics ) {
03020     // In QM/MM simulations, atom charges form QM regions need special treatment.
03021     if ( qmForcesOn ) {
03022         return 1;
03023     }
03024     if ( ! myMgr->ungridForcesCount && ! myMgr->recipEvirCount ) return 0;  // work to do, enqueue as usual
03025     myMgr->heldComputes.add(this);
03026     return 1;  // don't enqueue yet
03027   }
03028 
03029   positionBox->skip();
03030   forceBox->skip();
03031 
03032   if ( ++(myMgr->noWorkCount) == myMgr->pmeComputes.size() ) {
03033     myMgr->noWorkCount = 0;
03034     myMgr->reduction->submit();
03035   }
03036 
03037   atomsChanged = 0;
03038 
03039   return 1;  // no work for this step
03040 }

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

Definition at line 34 of file ComputePme.h.

00034 { myMgr = mgr; }

void ComputePme::ungridForces (  ) 

Definition at line 4023 of file ComputePme.C.

References ComputePmeMgr::alchLambda, ResizeArray< Elem >::begin(), PmeRealSpace::compute_forces(), ComputePmeMgr::compute_sequence, endi(), Patch::flags, ComputePmeMgr::fsize, iERROR(), iout, Flags::lattice, NAMD_bug(), Node::Object(), ComputePmeMgr::q_arr, ResizeArray< Elem >::resize(), scale_forces(), Compute::sequence(), Node::simParameters, simParams, Flags::step, Vector::x, Vector::y, and Vector::z.

Referenced by doWork().

04023                               {
04024 
04025     if ( sequence() != myMgr->compute_sequence ) NAMD_bug("ComputePme sequence mismatch in ungridForces()");
04026  
04027     SimParameters *simParams = Node::Object()->simParameters;
04028 
04029     localResults_alloc.resize(numLocalAtoms* ((numGrids>1 || selfOn)?2:1));
04030     Vector *localResults = localResults_alloc.begin();
04031     Vector *gridResults;
04032     if ( alchOn || lesOn || selfOn || pairOn ) {
04033       for(int i=0; i<numLocalAtoms; ++i) { localResults[i] = 0.; }
04034       gridResults = localResults + numLocalAtoms;
04035     } else {
04036       gridResults = localResults;
04037     }
04038 
04039     Vector pairForce = 0.;
04040     Lattice &lattice = patch->flags.lattice;
04041     int g = 0;
04042     if(!simParams->commOnly) {
04043     for ( g=0; g<numGrids; ++g ) {
04044 #ifdef NETWORK_PROGRESS
04045       CmiNetworkProgress();
04046 #endif
04047 
04048 #ifdef NAMD_CUDA
04049       if ( offload ) {
04050         int errfound = 0;
04051         for ( int n=numGridAtoms[g], i=0; i<n; ++i ) {
04052           // Neither isnan() nor x != x worked when testing on Cray; this does.
04053           if ( ((int*)f_data_host)[3*i] == 0x7fffffff ) { errfound = 1; }  // CUDA NaN
04054           gridResults[i].x = f_data_host[3*i];
04055           gridResults[i].y = f_data_host[3*i+1];
04056           gridResults[i].z = f_data_host[3*i+2];
04057         }
04058         if ( errfound ) {
04059           int errcount = 0;
04060           for ( int n=numGridAtoms[g], i=0; i<n; ++i ) {
04061             float f = f_data_host[3*i];
04062             if ( ((int*)f_data_host)[3*i] == 0x7fffffff ) {  // CUDA NaN
04063               ++errcount;
04064               gridResults[i] = 0.;
04065             }
04066           }
04067           iout << iERROR << "Stray PME grid charges detected: "
04068                 << errcount << " atoms on pe " << CkMyPe() << "\n" << endi;
04069         }
04070       } else
04071 #endif // NAMD_CUDA
04072         {
04073           myRealSpace[g]->compute_forces(myMgr->q_arr+g*myMgr->fsize, localGridData[g], gridResults);
04074         }
04075       scale_forces(gridResults, numGridAtoms[g], lattice);
04076       
04077       if (alchOn) {
04078         float scale = 1.;
04079         BigReal elecLambdaUp, elecLambdaDown;
04080         if ( simParams->alchFepWhamOn ) {
04081           if ( simParams->alchFepElecOn ) {
04082             elecLambdaUp = simParams->alchElecLambda;
04083             elecLambdaDown = 1.0 - simParams->alchElecLambda;
04084           }
04085           else {
04086             elecLambdaUp = 0.0;
04087             elecLambdaDown = 1.0;
04088           }
04089         }
04090         else {
04091           BigReal alchLambda = simParams->getCurrentLambda(patch->flags.step);
04092           myMgr->alchLambda = alchLambda;
04093           elecLambdaUp = simParams->getElecLambda(alchLambda);
04094           elecLambdaDown = simParams->getElecLambda(1. - alchLambda);
04095         }
04096         
04097         if ( g == 0 ) scale = elecLambdaUp;
04098         else if ( g == 1 ) scale = elecLambdaDown;
04099         else if ( g == 2 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
04100 
04101         if (alchDecouple) {
04102           if ( g == 2 ) scale = 1 - elecLambdaUp;
04103           else if ( g == 3 ) scale = 1 - elecLambdaDown;
04104           else if ( g == 4 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
04105         }
04106         int nga = 0;
04107         if (!alchDecouple) {
04108           for(int i=0; i<numLocalAtoms; ++i) {
04109             if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
04110               // (g=2: only partition 0)
04111               localResults[i] += gridResults[nga++] * scale;
04112             }
04113           }
04114         }
04115         else {  // alchDecouple
04116           if ( g < 2 ) {
04117             for(int i=0; i<numLocalAtoms; ++i) {
04118               if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
04119                 // g = 0: partition 0 or partition 1
04120                 // g = 1: partition 0 or partition 2
04121                 localResults[i] += gridResults[nga++] * scale;
04122               }
04123             }
04124           }
04125           else {
04126             for(int i=0; i<numLocalAtoms; ++i) {
04127               if ( localPartition[i] == (g-1) || localPartition[i] == (g-4)) {
04128                 // g = 2: partition 1 only
04129                 // g = 3: partition 2 only
04130                 // g = 4: partition 0 only
04131                 localResults[i] += gridResults[nga++] * scale;
04132               }
04133             }
04134           }
04135         }
04136       } else if ( lesOn ) {
04137         float scale = 1.;
04138         if ( alchFepOn ) {
04139           if(simParams->alchFepWhamOn) {
04140             if(simParams->alchFepElecOn) {
04141               if ( g == 0 ) scale = simParams->alchElecLambda;
04142               else if ( g == 1 ) scale = 1. - simParams->alchElecLambda;
04143             }
04144             else {
04145               if ( g == 0 ) scale = 0.0;
04146               else if ( g == 1 ) scale = 1.0;
04147             }
04148           }
04149           else {
04150             BigReal alchLambda = simParams->getCurrentLambda(patch->flags.step);
04151             myMgr->alchLambda = alchLambda;
04152             if ( g == 0 ) scale = alchLambda;
04153             else if ( g == 1 ) scale = 1. - alchLambda;
04154           }
04155         } else if ( lesOn ) {
04156           scale = 1.0 / (float)lesFactor;
04157         }
04158         int nga = 0;
04159         for(int i=0; i<numLocalAtoms; ++i) {
04160           if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
04161             localResults[i] += gridResults[nga++] * scale;
04162           }
04163         }
04164       } else if ( selfOn ) {
04165         PmeParticle *lgd = localGridData[g];
04166         int nga = 0;
04167         for(int i=0; i<numLocalAtoms; ++i) {
04168           if ( localPartition[i] == 1 ) {
04169             pairForce += gridResults[nga];  // should add up to almost zero
04170             localResults[i] += gridResults[nga++];
04171           }
04172         }
04173       } else if ( pairOn ) {
04174         if ( g == 0 ) {
04175           int nga = 0;
04176           for(int i=0; i<numLocalAtoms; ++i) {
04177             if ( localPartition[i] == 1 ) {
04178               pairForce += gridResults[nga];
04179             }
04180             if ( localPartition[i] == 1 || localPartition[i] == 2 ) {
04181               localResults[i] += gridResults[nga++];
04182             }
04183           }
04184         } else if ( g == 1 ) {
04185           int nga = 0;
04186           for(int i=0; i<numLocalAtoms; ++i) {
04187             if ( localPartition[i] == g ) {
04188               pairForce -= gridResults[nga];  // should add up to almost zero
04189               localResults[i] -= gridResults[nga++];
04190             }
04191           }
04192         } else {
04193           int nga = 0;
04194           for(int i=0; i<numLocalAtoms; ++i) {
04195             if ( localPartition[i] == g ) {
04196               localResults[i] -= gridResults[nga++];
04197             }
04198          }
04199         }
04200       }
04201     }
04202     }
04203 
04204     Vector *results_ptr = localResults;
04205     
04206     // add in forces
04207     {
04208       Results *r = forceBox->open();
04209       Force *f = r->f[Results::slow];
04210       int numAtoms = patch->getNumAtoms();
04211 
04212       if ( ! myMgr->strayChargeErrors && ! simParams->commOnly ) {
04213         for(int i=0; i<numAtoms; ++i) {
04214           f[i].x += results_ptr->x;
04215           f[i].y += results_ptr->y;
04216           f[i].z += results_ptr->z;
04217           ++results_ptr;
04218         }
04219       }
04220       forceBox->close(&r);
04221     }
04222 
04223     if ( pairOn || selfOn ) {
04224         ADD_VECTOR_OBJECT(myMgr->reduction,REDUCTION_PAIR_ELECT_FORCE,pairForce);
04225     }
04226 
04227 }


Friends And Related Function Documentation

friend class ComputePmeMgr [friend]

Definition at line 36 of file ComputePme.h.


The documentation for this class was generated from the following files:
Generated on Tue Nov 21 01:17:18 2017 for NAMD by  doxygen 1.4.7