ComputePme Class Reference

#include <ComputePme.h>

Inheritance diagram for ComputePme:

Compute ComputePmeUtil 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 46 of file ComputePme.h.


Constructor & Destructor Documentation

ComputePme::ComputePme ( ComputeID  c,
PatchID  pid 
)

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

02651                                                : Compute(c), patchID(pid)
02652 {
02653   DebugM(4,"ComputePme created.\n");
02654   basePriority = PME_PRIORITY;
02655   setNumPatches(1);
02656 
02657   CProxy_ComputePmeMgr::ckLocalBranch(
02658         CkpvAccess(BOCclass_group).computePmeMgr)->addCompute(this);
02659 
02660   SimParameters *simParams = Node::Object()->simParameters;
02661 
02662   qmForcesOn =  simParams->qmForcesOn;
02663   offload = simParams->PMEOffload;
02664 
02665   numGridsMax = numGrids;
02666 
02667   myGrid.K1 = simParams->PMEGridSizeX;
02668   myGrid.K2 = simParams->PMEGridSizeY;
02669   myGrid.K3 = simParams->PMEGridSizeZ;
02670   myGrid.order = simParams->PMEInterpOrder;
02671   myGrid.dim2 = myGrid.K2;
02672   myGrid.dim3 = 2 * (myGrid.K3/2 + 1);
02673 
02674 #ifdef NAMD_CUDA
02675   cuda_atoms_offset = 0;
02676   f_data_host = 0;
02677   f_data_dev = 0;
02678  if ( ! offload )
02679 #endif
02680  {
02681   for ( int g=0; g<numGrids; ++g ) myRealSpace[g] = new PmeRealSpace(myGrid);
02682  }
02683 
02684   atomsChanged = 0;
02685   
02686   qmLoclIndx = 0;
02687   qmLocalCharges = 0;
02688 }

ComputePme::~ComputePme (  )  [virtual]

Definition at line 2924 of file ComputePme.C.

02925 {
02926 #ifdef NAMD_CUDA
02927   if ( ! offload )
02928 #endif
02929   {
02930     for ( int g=0; g<numGridsMax; ++g ) delete myRealSpace[g];
02931   }
02932 }


Member Function Documentation

void ComputePme::atomUpdate (  )  [virtual]

Reimplemented from Compute.

Definition at line 2649 of file ComputePme.C.

02649 { atomsChanged = 1; }

void ComputePme::doQMWork (  ) 

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

03017                           {
03018     
03019 //     iout << CkMyPe() << ") ----> PME doQMWork.\n" << endi ;
03020     
03021     
03022     int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
03023     const Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
03024     const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
03025     const Real *qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
03026     
03027     const CompAtomExt *xExt = patch->getCompAtomExtInfo();
03028     
03029     // Determine number of qm atoms in this patch for the current step.
03030     numLocalQMAtoms = 0;
03031     for (int paIter=0; paIter<patch->getNumAtoms(); paIter++) {
03032         if ( qmAtomGroup[xExt[paIter].id] != 0 ) {
03033             numLocalQMAtoms++;
03034         }
03035     }
03036     
03037     // We prepare a charge vector with QM charges for use in the PME calculation.
03038     
03039     // Clears data from last step, if there is any.
03040     if (qmLoclIndx != 0)
03041         delete [] qmLoclIndx;
03042     if (qmLocalCharges != 0)
03043         delete [] qmLocalCharges;
03044     
03045     qmLoclIndx = new int[numLocalQMAtoms] ;
03046     qmLocalCharges = new Real[numLocalQMAtoms] ;
03047     
03048     // I am assuming there will be (in general) more QM atoms among all QM groups
03049     // than MM atoms in a patch.
03050     int procAtms = 0;
03051     
03052     for (int paIter=0; paIter<patch->getNumAtoms(); paIter++) {
03053         
03054         for (int i=0; i<numQMAtms; i++) {
03055             
03056             if (qmAtmIndx[i] == xExt[paIter].id) {
03057                 
03058                 qmLoclIndx[procAtms] = paIter ;
03059                 qmLocalCharges[procAtms] = qmAtmChrg[i];
03060                 
03061                 procAtms++;
03062                 break;
03063             }
03064             
03065         }
03066         
03067         if (procAtms == numLocalQMAtoms)
03068             break;
03069     }
03070     
03071     doWork();
03072     return ;
03073 }

void ComputePme::doWork (  )  [virtual]

Reimplemented from Compute.

Definition at line 3075 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, ComputePmeUtil::selfOn, ComputePmeMgr::submitReductions(), ungridForces(), ComputePmeMgr::ungridForcesCount, PmeParticle::x, x, PmeParticle::y, and PmeParticle::z.

Referenced by doQMWork().

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

void ComputePme::initialize (  )  [virtual]

Reimplemented from Compute.

Definition at line 2690 of file ComputePme.C.

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

02690                             {
02691   if (!(patch = PatchMap::Object()->patch(patchID))) {
02692     NAMD_bug("ComputePme used with unknown patch.");
02693   }
02694   positionBox = patch->registerPositionPickup(this);
02695   avgPositionBox = patch->registerAvgPositionPickup(this);
02696   forceBox = patch->registerForceDeposit(this);
02697 #ifdef NAMD_CUDA
02698  if ( offload ) {
02699   myMgr->cuda_atoms_count += patch->getNumAtoms();
02700  }
02701 #endif
02702 }

int ComputePme::noWork (  )  [virtual]

Reimplemented from Compute.

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

02978                        {
02979 
02980   if ( patch->flags.doFullElectrostatics ) {
02981     // In QM/MM simulations, atom charges form QM regions need special treatment.
02982     if ( qmForcesOn ) {
02983         return 1;
02984     }
02985     if ( ! myMgr->ungridForcesCount && ! myMgr->recipEvirCount ) return 0;  // work to do, enqueue as usual
02986     myMgr->heldComputes.add(this);
02987     return 1;  // don't enqueue yet
02988   }
02989 
02990   positionBox->skip();
02991   forceBox->skip();
02992 
02993   if ( ++(myMgr->noWorkCount) == myMgr->pmeComputes.size() ) {
02994     myMgr->noWorkCount = 0;
02995     myMgr->reduction->submit();
02996   }
02997 
02998   atomsChanged = 0;
02999 
03000   return 1;  // no work for this step
03001 }

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

Definition at line 56 of file ComputePme.h.

00056 { myMgr = mgr; }

void ComputePme::ungridForces (  ) 

Definition at line 3984 of file ComputePme.C.

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

Referenced by doWork().

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


Friends And Related Function Documentation

friend class ComputePmeMgr [friend]

Definition at line 58 of file ComputePme.h.


The documentation for this class was generated from the following files:
Generated on Wed Jun 20 01:17:18 2018 for NAMD by  doxygen 1.4.7