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       if (g < 2) {
03178       int nga = 0;
03179       for(int i=0; i<numLocalAtoms; ++i) {
03180         if ( localPartition[i] == 0 || localPartition[i] == (g+1) || localPartition[i] == (g+3)) {
03181           // for FEP/TI: grid 0 gets non-alch + partition 1 + partition 3;
03182           // grid 1 gets non-alch + partition 2 + + partition 4;
03183           lgd[nga++] = localData[i];
03184         }
03185        }
03186        numGridAtoms[g] = nga;
03187        } else {
03188         int nga = 0;
03189         for(int i=0; i<numLocalAtoms; ++i) {
03190         if ( localPartition[i] == 0 ) {
03191           // grid 2 (only if called for with numGrids=3) gets only non-alch
03192           lgd[nga++] = localData[i];
03193         }
03194        }
03195        numGridAtoms[g] = nga;
03196       }
03197     }
03198   } else if ( alchOn && alchDecouple) {
03199     // alchemical decoupling: four grids
03200     // g=0: partition 0 and partition 1
03201     // g=1: partition 0 and partition 2
03202     // g=2: only partition 1 atoms
03203     // g=3: only partition 2 atoms
03204     // plus one grid g=4, only partition 0, if numGrids=5
03205     for ( g=0; g<2; ++g ) {  // same as before for first 2
03206       PmeParticle *lgd = localGridData[g];
03207       int nga = 0;
03208       for(int i=0; i<numLocalAtoms; ++i) {
03209         if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
03210           lgd[nga++] = localData[i];
03211         }
03212       }
03213       numGridAtoms[g] = nga;
03214     }
03215     for (g=2 ; g<4 ; ++g ) {  // only alchemical atoms for these 2
03216       PmeParticle *lgd = localGridData[g];
03217       int nga = 0;
03218       for(int i=0; i<numLocalAtoms; ++i) {
03219         if ( localPartition[i] == (g-1) ) {
03220           lgd[nga++] = localData[i];
03221         }
03222       }
03223       numGridAtoms[g] = nga;
03224     }
03225     for (g=4 ; g<numGrids ; ++g ) {  // only non-alchemical atoms 
03226       // numGrids=5 only if alchElecLambdaStart > 0
03227       PmeParticle *lgd = localGridData[g];
03228       int nga = 0;
03229       for(int i=0; i<numLocalAtoms; ++i) {
03230         if ( localPartition[i] == 0 ) {
03231           lgd[nga++] = localData[i];
03232         }
03233       }
03234       numGridAtoms[g] = nga;
03235     }
03236   } else if ( selfOn ) {
03237     if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 1 failed");
03238     g = 0;
03239     PmeParticle *lgd = localGridData[g];
03240     int nga = 0;
03241     for(int i=0; i<numLocalAtoms; ++i) {
03242       if ( localPartition[i] == 1 ) {
03243         lgd[nga++] = localData[i];
03244       }
03245     }
03246     numGridAtoms[g] = nga;
03247   } else if ( pairOn ) {
03248     if ( numGrids != 3 ) NAMD_bug("ComputePme::doWork assertion 2 failed");
03249     g = 0;
03250     PmeParticle *lgd = localGridData[g];
03251     int nga = 0;
03252     for(int i=0; i<numLocalAtoms; ++i) {
03253       if ( localPartition[i] == 1 || localPartition[i] == 2 ) {
03254         lgd[nga++] = localData[i];
03255       }
03256     }
03257     numGridAtoms[g] = nga;
03258     for ( g=1; g<3; ++g ) {
03259       PmeParticle *lgd = localGridData[g];
03260       int nga = 0;
03261       for(int i=0; i<numLocalAtoms; ++i) {
03262         if ( localPartition[i] == g ) {
03263           lgd[nga++] = localData[i];
03264         }
03265       }
03266       numGridAtoms[g] = nga;
03267     }
03268   } else {
03269     if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 3 failed");
03270     localGridData[0] = localData;
03271     numGridAtoms[0] = numLocalAtoms;
03272   }
03273 
03274  if ( ! myMgr->doWorkCount ) {
03275   myMgr->doWorkCount = myMgr->pmeComputes.size();
03276 
03277 #ifdef NAMD_CUDA
03278  if ( !  offload )
03279 #endif // NAMD_CUDA
03280  {
03281   memset( (void*) myMgr->fz_arr, 0, (myGrid.K3+myGrid.order-1) * sizeof(char) );
03282 
03283   for (int i=0; i<myMgr->q_count; ++i) {
03284     memset( (void*) (myMgr->q_list[i]), 0, (myGrid.K3+myGrid.order-1) * sizeof(float) );
03285   }
03286  }
03287 
03288   for ( g=0; g<numGrids; ++g ) {
03289     myMgr->evir[g] = 0;
03290   }
03291 
03292   myMgr->strayChargeErrors = 0;
03293 
03294   myMgr->compute_sequence = sequence();
03295  }
03296 
03297   if ( sequence() != myMgr->compute_sequence ) NAMD_bug("ComputePme sequence mismatch in doWork()");
03298 
03299   int strayChargeErrors = 0;
03300 
03301   // calculate self energy
03302   BigReal ewaldcof = ComputeNonbondedUtil::ewaldcof;
03303   for ( g=0; g<numGrids; ++g ) {
03304     BigReal selfEnergy = 0;
03305     data_ptr = localGridData[g];
03306     int i;
03307     for(i=0; i<numGridAtoms[g]; ++i)
03308     {
03309       selfEnergy += data_ptr->cg * data_ptr->cg;
03310       ++data_ptr;
03311     }
03312     selfEnergy *= -1. * ewaldcof / SQRT_PI;
03313     myMgr->evir[g][0] += selfEnergy;
03314 
03315     float **q = myMgr->q_arr + g*myMgr->fsize;
03316     char *f = myMgr->f_arr + g*myMgr->fsize;
03317 
03318     scale_coordinates(localGridData[g], numGridAtoms[g], lattice, myGrid);
03319 #ifdef NAMD_CUDA
03320    if ( offload ) {
03321     if ( myMgr->cuda_atoms_alloc == 0 ) {  // first call
03322       int na = myMgr->cuda_atoms_alloc = 1.2 * (myMgr->cuda_atoms_count + 1000);
03323       cuda_errcheck("before malloc atom data for pme");
03324       cudaMallocHost((void**) &(myMgr->a_data_host), 7*na*sizeof(float));
03325       cudaMalloc((void**) &(myMgr->a_data_dev), 7*na*sizeof(float));
03326       cuda_errcheck("malloc atom data for pme");
03327       myMgr->cuda_atoms_count = 0;
03328     }
03329     cuda_atoms_offset = myMgr->cuda_atoms_count;
03330     int n = numGridAtoms[g];
03331     myMgr->cuda_atoms_count += n;
03332     if ( myMgr->cuda_atoms_count > myMgr->cuda_atoms_alloc ) {
03333       CkPrintf("Pe %d expanding CUDA PME atoms allocation because %d > %d\n",
03334                         CkMyPe(), myMgr->cuda_atoms_count, myMgr->cuda_atoms_alloc);
03335       cuda_errcheck("before malloc expanded atom data for pme");
03336       int na = myMgr->cuda_atoms_alloc = 1.2 * (myMgr->cuda_atoms_count + 1000);
03337       const float *a_data_host_old = myMgr->a_data_host;
03338       cudaMallocHost((void**) &(myMgr->a_data_host), 7*na*sizeof(float));
03339       cuda_errcheck("malloc expanded host atom data for pme");
03340       memcpy(myMgr->a_data_host, a_data_host_old, 7*cuda_atoms_offset*sizeof(float));
03341       cudaFreeHost((void*) a_data_host_old);
03342       cuda_errcheck("free expanded host atom data for pme");
03343       cudaFree(myMgr->a_data_dev);
03344       cuda_errcheck("free expanded dev atom data for pme");
03345       cudaMalloc((void**) &(myMgr->a_data_dev), 7*na*sizeof(float));
03346       cuda_errcheck("malloc expanded dev atom data for pme");
03347     }
03348     float *a_data_host = myMgr->a_data_host + 7 * cuda_atoms_offset;
03349     data_ptr = localGridData[g];
03350     double order_1 = myGrid.order - 1;
03351     double K1 = myGrid.K1;
03352     double K2 = myGrid.K2;
03353     double K3 = myGrid.K3;
03354     int found_negative = 0;
03355     for ( int i=0; i<n; ++i ) {
03356       if ( data_ptr[i].x < 0 || data_ptr[i].y < 0 || data_ptr[i].z < 0 ) {
03357         found_negative = 1;
03358         // CkPrintf("low coord: %f %f %f\n", data_ptr[i].x, data_ptr[i].y, data_ptr[i].z);
03359       }
03360       double x_int = (int) data_ptr[i].x;
03361       double y_int = (int) data_ptr[i].y;
03362       double z_int = (int) data_ptr[i].z;
03363       a_data_host[7*i  ] = data_ptr[i].x - x_int;  // subtract in double precision
03364       a_data_host[7*i+1] = data_ptr[i].y - y_int;
03365       a_data_host[7*i+2] = data_ptr[i].z - z_int;
03366       a_data_host[7*i+3] = data_ptr[i].cg;
03367       x_int -= order_1;  if ( x_int < 0 ) x_int += K1;
03368       y_int -= order_1;  if ( y_int < 0 ) y_int += K2;
03369       z_int -= order_1;  if ( z_int < 0 ) z_int += K3;
03370       a_data_host[7*i+4] = x_int;
03371       a_data_host[7*i+5] = y_int;
03372       a_data_host[7*i+6] = z_int;
03373     }
03374     if ( found_negative ) NAMD_bug("found negative atom coordinate in ComputePme::doWork");
03375    } else
03376 #endif // NAMD_CUDA
03377    {
03378     myRealSpace[g]->set_num_atoms(numGridAtoms[g]);
03379     myRealSpace[g]->fill_charges(q, myMgr->q_list, myMgr->q_count, strayChargeErrors, f, myMgr->fz_arr, localGridData[g]);
03380    }
03381   }
03382   myMgr->strayChargeErrors += strayChargeErrors;
03383 
03384 #ifdef TRACE_COMPUTE_OBJECTS
03385     traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+this->cid, traceObjStartTime, CmiWallTimer());
03386 #endif
03387 
03388  if ( --(myMgr->doWorkCount) == 0 ) {
03389 // cudaDeviceSynchronize();  // XXXX
03390 #ifdef NAMD_CUDA
03391   if ( offload ) {
03392     ComputePmeMgr::cuda_submit_charges_args args;
03393     args.mgr = myMgr;
03394     args.lattice = &lattice;
03395     args.sequence = sequence();
03396     CmiLock(ComputePmeMgr::cuda_lock);
03397     if ( ComputePmeMgr::cuda_busy ) {
03398       ComputePmeMgr::cuda_submit_charges_deque.push_back(args);
03399     } else if ( CkMyPe() == deviceCUDA->getMasterPe() ) {
03400       // avoid adding work to nonbonded data preparation pe
03401       args.mgr->cuda_submit_charges(*args.lattice, args.sequence);
03402     } else {
03403       ComputePmeMgr::cuda_busy = true;
03404       while ( 1 ) {
03405         CmiUnlock(ComputePmeMgr::cuda_lock);
03406         args.mgr->cuda_submit_charges(*args.lattice, args.sequence);
03407         CmiLock(ComputePmeMgr::cuda_lock);
03408         if ( ComputePmeMgr::cuda_submit_charges_deque.size() ) {
03409           args = ComputePmeMgr::cuda_submit_charges_deque.front();
03410           ComputePmeMgr::cuda_submit_charges_deque.pop_front();
03411         } else {
03412           ComputePmeMgr::cuda_busy = false;
03413           break;
03414         }
03415       }
03416     }
03417     CmiUnlock(ComputePmeMgr::cuda_lock);
03418   } else
03419 #endif // NAMD_CUDA
03420   {
03421     myMgr->chargeGridReady(lattice,sequence());
03422   }
03423  }
03424  atomsChanged = 0;
03425 }

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 4012 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().

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


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 21 Sep 2020 for NAMD by  doxygen 1.6.1