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 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, SimParameters::PMEGridSizeX, SimParameters::PMEGridSizeY, SimParameters::PMEGridSizeZ, SimParameters::PMEInterpOrder, SimParameters::PMEOffload, SimParameters::qmForcesOn, Compute::setNumPatches(), and Node::simParameters.

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   numGridsMax = numGrids;
02682 
02683   myGrid.K1 = simParams->PMEGridSizeX;
02684   myGrid.K2 = simParams->PMEGridSizeY;
02685   myGrid.K3 = simParams->PMEGridSizeZ;
02686   myGrid.order = simParams->PMEInterpOrder;
02687   myGrid.dim2 = myGrid.K2;
02688   myGrid.dim3 = 2 * (myGrid.K3/2 + 1);
02689 
02690 #ifdef NAMD_CUDA
02691   cuda_atoms_offset = 0;
02692   f_data_host = 0;
02693   f_data_dev = 0;
02694  if ( ! offload )
02695 #endif
02696  {
02697   for ( int g=0; g<numGrids; ++g ) myRealSpace[g] = new PmeRealSpace(myGrid);
02698  }
02699 
02700   atomsChanged = 0;
02701   
02702   qmLoclIndx = 0;
02703   qmLocalCharges = 0;
02704 }

ComputePme::~ComputePme (  )  [virtual]

Definition at line 2940 of file ComputePme.C.

02941 {
02942 #ifdef NAMD_CUDA
02943   if ( ! offload )
02944 #endif
02945   {
02946     for ( int g=0; g<numGridsMax; ++g ) delete myRealSpace[g];
02947   }
02948 }


Member Function Documentation

void ComputePme::atomUpdate ( void   )  [virtual]

Reimplemented from Compute.

Definition at line 2665 of file ComputePme.C.

02665 { atomsChanged = 1; }

void ComputePme::doQMWork (  ) 

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

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

void ComputePme::doWork ( void   )  [virtual]

Reimplemented from Compute.

Definition at line 3091 of file ComputePme.C.

References ComputePmeMgr::a_data_dev, ComputePmeMgr::a_data_host, ComputePmeUtil::alchDecouple, ComputePmeUtil::alchOn, Compute::basePriority, ResizeArray< Elem >::begin(), PmeParticle::cg, CompAtom::charge, ComputePmeMgr::chargeGridReady(), Compute::cid, Box< Owner, Data >::close(), COMPUTE_HOME_PRIORITY, COULOMB, ComputePmeMgr::cuda_atoms_alloc, ComputePmeMgr::cuda_atoms_count, cuda_errcheck(), ComputePmeMgr::cuda_submit_charges(), DebugM, deviceCUDA, ComputeNonbondedUtil::dielectric_1, Flags::doMolly, ComputeNonbondedUtil::ewaldcof, PmeRealSpace::fill_charges(), Patch::flags, DeviceCUDA::getDeviceID(), DeviceCUDA::getMasterPe(), Patch::getNumAtoms(), PmeGrid::K1, PmeGrid::K2, PmeGrid::K3, ComputePmeMgr::cuda_submit_charges_args::lattice, Flags::lattice, ComputePmeUtil::lesOn, ComputePmeMgr::cuda_submit_charges_args::mgr, NAMD_bug(), Box< Owner, Data >::open(), PmeGrid::order, ComputePmeUtil::pairOn, CompAtom::partition, PATCH_PRIORITY, PME_OFFLOAD_PRIORITY, PME_PRIORITY, ComputePmeMgr::pmeComputes, CompAtom::position, ResizeArray< Elem >::resize(), scale_coordinates(), ComputeNonbondedUtil::scaling, ComputePmeUtil::selfOn, ComputePmeMgr::cuda_submit_charges_args::sequence, Compute::sequence(), PmeRealSpace::set_num_atoms(), ResizeArray< Elem >::size(), SQRT_PI, ComputePmeMgr::submitReductions(), TRACE_COMPOBJ_IDOFFSET, ungridForces(), Vector::x, PmeParticle::x, Vector::y, PmeParticle::y, Vector::z, and PmeParticle::z.

Referenced by doQMWork().

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

void ComputePme::initialize ( void   )  [virtual]

Reimplemented from Compute.

Definition at line 2706 of file ComputePme.C.

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

02706                             {
02707   if (!(patch = PatchMap::Object()->patch(patchID))) {
02708     NAMD_bug("ComputePme used with unknown patch.");
02709   }
02710   positionBox = patch->registerPositionPickup(this);
02711   avgPositionBox = patch->registerAvgPositionPickup(this);
02712   forceBox = patch->registerForceDeposit(this);
02713 #ifdef NAMD_CUDA
02714  if ( offload ) {
02715   myMgr->cuda_atoms_count += patch->getNumAtoms();
02716  }
02717 #endif
02718 }

int ComputePme::noWork (  )  [virtual]

Reimplemented from Compute.

Definition at line 2994 of file ComputePme.C.

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

02994                        {
02995 
02996   if ( patch->flags.doFullElectrostatics ) {
02997     // In QM/MM simulations, atom charges form QM regions need special treatment.
02998     if ( qmForcesOn ) {
02999         return 1;
03000     }
03001     if ( ! myMgr->ungridForcesCount && ! myMgr->recipEvirCount ) return 0;  // work to do, enqueue as usual
03002     myMgr->heldComputes.add(this);
03003     return 1;  // don't enqueue yet
03004   }
03005 
03006   positionBox->skip();
03007   forceBox->skip();
03008 
03009   if ( ++(myMgr->noWorkCount) == myMgr->pmeComputes.size() ) {
03010     myMgr->noWorkCount = 0;
03011     myMgr->reduction->submit();
03012   }
03013 
03014   atomsChanged = 0;
03015 
03016   return 1;  // no work for this step
03017 }

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

Definition at line 56 of file ComputePme.h.

00056 { myMgr = mgr; }

void ComputePme::ungridForces (  ) 

Definition at line 4002 of file ComputePme.C.

References ADD_VECTOR_OBJECT, ComputePmeUtil::alchDecouple, ComputePmeUtil::alchFepOn, ComputePmeUtil::alchOn, ResizeArray< Elem >::begin(), Box< Owner, Data >::close(), SimParameters::commOnly, PmeRealSpace::compute_forces(), endi(), Results::f, Patch::flags, SimParameters::getCurrentLambda(), SimParameters::getCurrentLambda2(), SimParameters::getElecLambda(), Patch::getNumAtoms(), iERROR(), iout, Flags::lattice, ComputePmeUtil::lesFactor, ComputePmeUtil::lesOn, NAMD_bug(), Node::Object(), Box< Owner, Data >::open(), ComputePmeUtil::pairOn, ResizeArray< Elem >::resize(), scale_forces(), ComputePmeUtil::selfOn, Compute::sequence(), Node::simParameters, Results::slow, Flags::step, Vector::x, Vector::y, and Vector::z.

Referenced by doWork().

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


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 11 Nov 2019 for NAMD by  doxygen 1.6.1