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

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

ComputePme::~ComputePme (  )  [virtual]

Definition at line 2960 of file ComputePme.C.

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


Member Function Documentation

void ComputePme::atomUpdate (  )  [virtual]

Reimplemented from Compute.

Definition at line 2663 of file ComputePme.C.

02663 { atomsChanged = 1; }

void ComputePme::doQMWork (  ) 

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

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

void ComputePme::doWork (  )  [virtual]

Reimplemented from Compute.

Definition at line 3111 of file ComputePme.C.

References Compute::basePriority, ResizeArray< Elem >::begin(), PmeParticle::cg, Box< Owner, Data >::close(), COMPUTE_HOME_PRIORITY, COULOMB, DebugM, ComputeNonbondedUtil::dielectric_1, Flags::doMolly, Patch::flags, 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().

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

void ComputePme::initialize (  )  [virtual]

Reimplemented from Compute.

Definition at line 2727 of file ComputePme.C.

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

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

int ComputePme::noWork (  )  [virtual]

Reimplemented from Compute.

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

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

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

Definition at line 34 of file ComputePme.h.

00034 { myMgr = mgr; }

void ComputePme::ungridForces (  ) 

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

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


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 Mon Sep 25 01:17:17 2017 for NAMD by  doxygen 1.4.7