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 2650 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.

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

ComputePme::~ComputePme (  )  [virtual]

Definition at line 2923 of file ComputePme.C.

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


Member Function Documentation

void ComputePme::atomUpdate (  )  [virtual]

Reimplemented from Compute.

Definition at line 2648 of file ComputePme.C.

02648 { atomsChanged = 1; }

void ComputePme::doQMWork (  ) 

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

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

void ComputePme::doWork (  )  [virtual]

Reimplemented from Compute.

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

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

void ComputePme::initialize (  )  [virtual]

Reimplemented from Compute.

Definition at line 2689 of file ComputePme.C.

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

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

int ComputePme::noWork (  )  [virtual]

Reimplemented from Compute.

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

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

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

Definition at line 56 of file ComputePme.h.

00056 { myMgr = mgr; }

void ComputePme::ungridForces (  ) 

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

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


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 Mon Feb 19 01:17:18 2018 for NAMD by  doxygen 1.4.7