ComputePme Class Reference

#include <ComputePme.h>

Inheritance diagram for ComputePme:

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)
Public Member Functions inherited from Compute
	Compute (ComputeID)
int	type ()
virtual	~Compute ()
void	setNumPatches (int n)
int	getNumPatches ()
virtual void	patchReady (PatchID, int doneMigration, int seq)
virtual void	finishPatch (int)
int	sequence (void)
int	priority (void)
int	getGBISPhase (void)
virtual void	gbisP2PatchReady (PatchID, int seq)
virtual void	gbisP3PatchReady (PatchID, int seq)
Public Member Functions inherited from ComputePmeUtil
	ComputePmeUtil ()
	~ComputePmeUtil ()

Friends
class	ComputePmeMgr

Additional Inherited Members
Static Public Member Functions inherited from ComputePmeUtil
static void	select (void)
Public Attributes inherited from Compute
const ComputeID	cid
LDObjHandle	ldObjHandle
LocalWorkMsg *const	localWorkMsg
Static Public Attributes inherited from ComputePmeUtil
static int	numGrids
static Bool	alchOn
static Bool	alchFepOn
static Bool	alchThermIntOn
static Bool	alchDecouple
static BigReal	alchElecLambdaStart
static Bool	lesOn
static int	lesFactor
static Bool	pairOn
static Bool	selfOn
static Bool	LJPMEOn
Protected Member Functions inherited from Compute
void	enqueueWork ()
Protected Attributes inherited from Compute
int	computeType
int	basePriority
int	gbisPhase
int	gbisPhasePriority [3]

Detailed Description

Definition at line 48 of file ComputePme.h.

Constructor & Destructor Documentation

ComputePme()

ComputePme::ComputePme	(	ComputeID	c,
		PatchID	pid
	)

Definition at line 2712 of file ComputePme.C.

References Compute::basePriority, DebugM, PmeGrid::dim2, PmeGrid::dim3, PmeGrid::K1, PmeGrid::K2, PmeGrid::K3, ComputePmeUtil::numGrids, Node::Object(), PmeGrid::order, PME_PRIORITY, Compute::setNumPatches(), Node::simParameters, and simParams.

                                               : Compute(c), patchID(pid)
{
  DebugM(4,"ComputePme created.\n");
  basePriority = PME_PRIORITY;
  setNumPatches(1);

  CProxy_ComputePmeMgr::ckLocalBranch(
        CkpvAccess(BOCclass_group).computePmeMgr)->addCompute(this);

  SimParameters *simParams = Node::Object()->simParameters;

  qmForcesOn =  simParams->qmForcesOn;
  offload = simParams->PMEOffload;

  numGridsMax = numGrids;

  myGrid.K1 = simParams->PMEGridSizeX;
  myGrid.K2 = simParams->PMEGridSizeY;
  myGrid.K3 = simParams->PMEGridSizeZ;
  myGrid.order = simParams->PMEInterpOrder;
  myGrid.dim2 = myGrid.K2;
  myGrid.dim3 = 2 * (myGrid.K3/2 + 1);

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
  cuda_atoms_offset = 0;
  f_data_host = 0;
  f_data_dev = 0;
 if ( ! offload )
#endif
 {
  for ( int g=0; g<numGrids; ++g ) myRealSpace[g] = new PmeRealSpace(myGrid);
 }

  atomsChanged = 0;
  
  qmLoclIndx = 0;
  qmLocalCharges = 0;
}

~ComputePme()

ComputePme::~ComputePme ( )

virtual

Definition at line 2985 of file ComputePme.C.

{
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
  if ( ! offload )
#endif
  {
    for ( int g=0; g<numGridsMax; ++g ) delete myRealSpace[g];
  }
}

Member Function Documentation

atomUpdate()

void ComputePme::atomUpdate ( void  )

virtual

Reimplemented from Compute.

Definition at line 2710 of file ComputePme.C.

{ atomsChanged = 1; }

doQMWork()

void ComputePme::doQMWork

(

)

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

                          {
    
//     iout << CkMyPe() << ") ----> PME doQMWork.\n" << endi ;
    
    
    int numQMAtms = Node::Object()->molecule->get_numQMAtoms();
    const Real *qmAtmChrg = Node::Object()->molecule->get_qmAtmChrg() ;
    const int *qmAtmIndx = Node::Object()->molecule->get_qmAtmIndx() ;
    const Real *qmAtomGroup = Node::Object()->molecule->get_qmAtomGroup() ;
    
    const CompAtomExt *xExt = patch->getCompAtomExtInfo();
    
    // Determine number of qm atoms in this patch for the current step.
    numLocalQMAtoms = 0;
    for (int paIter=0; paIter<patch->getNumAtoms(); paIter++) {
        if ( qmAtomGroup[xExt[paIter].id] != 0 ) {
            numLocalQMAtoms++;
        }
    }
    
    // We prepare a charge vector with QM charges for use in the PME calculation.
    
    // Clears data from last step, if there is any.
    if (qmLoclIndx != 0)
        delete [] qmLoclIndx;
    if (qmLocalCharges != 0)
        delete [] qmLocalCharges;
    
    qmLoclIndx = new int[numLocalQMAtoms] ;
    qmLocalCharges = new Real[numLocalQMAtoms] ;
    
    // I am assuming there will be (in general) more QM atoms among all QM groups
    // than MM atoms in a patch.
    int procAtms = 0;
    
    for (int paIter=0; paIter<patch->getNumAtoms(); paIter++) {
        
        for (int i=0; i<numQMAtms; i++) {
            
            if (qmAtmIndx[i] == xExt[paIter].id) {
                
                qmLoclIndx[procAtms] = paIter ;
                qmLocalCharges[procAtms] = qmAtmChrg[i];
                
                procAtms++;
                break;
            }
            
        }
        
        if (procAtms == numLocalQMAtoms)
            break;
    }
    
    doWork();
    return ;
}

doWork()

void ComputePme::doWork ( void  )

virtual

Reimplemented from Compute.

Definition at line 3136 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, ComputePmeMgr::cuda_busy, cuda_errcheck(), ComputePmeMgr::cuda_lock, ComputePmeMgr::cuda_submit_charges(), ComputePmeMgr::cuda_submit_charges_deque, DebugM, deviceCUDA, ComputeNonbondedUtil::dielectric_1, CompAtomExt::dispcoef, Flags::doMolly, ComputeNonbondedUtil::ewaldcof, PmeRealSpace::fill_charges(), Patch::flags, Patch::getCompAtomExtInfo(), DeviceCUDA::getDeviceID(), DeviceCUDA::getMasterPe(), Patch::getNumAtoms(), PmeGrid::K1, PmeGrid::K2, PmeGrid::K3, Flags::lattice, ComputePmeMgr::cuda_submit_charges_args::lattice, ComputePmeUtil::lesOn, ComputeNonbondedUtil::LJewaldcof, ComputePmeUtil::LJPMEOn, ComputePmeMgr::cuda_submit_charges_args::mgr, NAMD_bug(), ComputePmeUtil::numGrids, 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, Compute::sequence(), ComputePmeMgr::cuda_submit_charges_args::sequence, PmeRealSpace::set_num_atoms(), ResizeArray< Elem >::size(), SQRT_PI, ComputePmeMgr::submitReductions(), TRACE_COMPOBJ_IDOFFSET, ungridForces(), PmeParticle::x, Vector::x, PmeParticle::y, Vector::y, PmeParticle::z, and Vector::z.

Referenced by doQMWork().

{
  DebugM(4,"Entering ComputePme::doWork().\n");

  if ( basePriority >= COMPUTE_HOME_PRIORITY ) {
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
    basePriority = ( offload ? PME_OFFLOAD_PRIORITY : PME_PRIORITY );
#else
    basePriority = PME_PRIORITY;
#endif
    ungridForces();
    // CkPrintf("doWork 2 pe %d %d %d\n", CkMyPe(), myMgr->ungridForcesCount, myMgr->recipEvirCount);
    if ( ! --(myMgr->ungridForcesCount) && ! myMgr->recipEvirCount ) myMgr->submitReductions();
    return;
  }
  basePriority = COMPUTE_HOME_PRIORITY + PATCH_PRIORITY(patchID);
  // CkPrintf("doWork 1 pe %d %d %d\n", CkMyPe(), myMgr->ungridForcesCount, myMgr->recipEvirCount);

#ifdef TRACE_COMPUTE_OBJECTS
    double traceObjStartTime = CmiWallTimer();
#endif

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
  if ( offload ) cudaSetDevice(deviceCUDA->getDeviceID());
#endif

  // allocate storage
  numLocalAtoms = patch->getNumAtoms();

  Lattice &lattice = patch->flags.lattice;

  // For more than one grid, allocate base grid plus numGrids extra storage.
  // Each storage segment can hold the max atoms possible, numLocalAtoms.
  // Copy coordinate data from position box into base grid plus create 
  // auxiliary array with partition number. For alchemy, copy coordinate data
  // from base grid into each extra buffer depending on partition values.
  // Storage is all PmeParticle: x,y,z,q double precision.
  int extraGrids = 0;
  if ( ! LJPMEOn && (numGrids > 1 || selfOn) ) {
    extraGrids = 1;
  }

  localData_alloc.resize(numLocalAtoms*(numGrids+extraGrids));
  localData = localData_alloc.begin();
  localPartition_alloc.resize(numLocalAtoms);
  localPartition = localPartition_alloc.begin();

  // We have local buffers: base, 0, 1, ..., numGrids-1 (for numGrids > 1).
  // localGridData points to the "0, 1, ..., numGrids-1" buffers.
  int g;
  for ( g=0; g<numGrids; ++g ) {
    localGridData[g] = localData + numLocalAtoms*(g+extraGrids);
  }

  // get positions and charges
  PmeParticle * data_ptr = localData;
  unsigned char * part_ptr = localPartition;
  const BigReal coulomb_sqrt = sqrt( COULOMB * ComputeNonbondedUtil::scaling
                                * ComputeNonbondedUtil::dielectric_1 );
  {
    CompAtom *x = positionBox->open();
    // CompAtomExt *xExt = patch->getCompAtomExtInfo();
    if ( patch->flags.doMolly ) {
      positionBox->close(&x);
      x = avgPositionBox->open();
    }
    int numAtoms = patch->getNumAtoms();

    for(int i=0; i<numAtoms; ++i)
    {
      data_ptr->x = x[i].position.x;
      data_ptr->y = x[i].position.y;
      data_ptr->z = x[i].position.z;
      data_ptr->cg = coulomb_sqrt * x[i].charge;
      ++data_ptr;
      *part_ptr = x[i].partition;
      ++part_ptr;
    }

    // QM loop to overwrite charges of QM atoms.
    // They are zero for NAMD, but are updated in ComputeQM.
    if ( qmForcesOn ) {
        
        for(int i=0; i<numLocalQMAtoms; ++i)
        {
          localData[qmLoclIndx[i]].cg = coulomb_sqrt * qmLocalCharges[i];
        }
        
    }
    
    if ( patch->flags.doMolly ) { avgPositionBox->close(&x); }
    else { positionBox->close(&x); }
  }

  // copy to other grids if needed
  if ( (alchOn && (!alchDecouple)) || lesOn ) {
    for ( g=0; g<numGrids; ++g ) {
      PmeParticle *lgd = localGridData[g];
      if (g < 2) {
      int nga = 0;
      for(int i=0; i<numLocalAtoms; ++i) {
        if ( localPartition[i] == 0 || localPartition[i] == (g+1) || localPartition[i] == (g+3)) {
          // for FEP/TI: grid 0 gets non-alch + partition 1 + partition 3;
          // grid 1 gets non-alch + partition 2 + + partition 4;
          lgd[nga++] = localData[i];
        }
       }
       numGridAtoms[g] = nga;
       } else {
        int nga = 0;
        for(int i=0; i<numLocalAtoms; ++i) {
        if ( localPartition[i] == 0 ) {
          // grid 2 (only if called for with numGrids=3) gets only non-alch
          lgd[nga++] = localData[i];
        }
       }
       numGridAtoms[g] = nga;
      }
    }
  } else if ( alchOn && alchDecouple) {
    // alchemical decoupling: four grids
    // g=0: partition 0 and partition 1
    // g=1: partition 0 and partition 2
    // g=2: only partition 1 atoms
    // g=3: only partition 2 atoms
    // plus one grid g=4, only partition 0, if numGrids=5
    for ( g=0; g<2; ++g ) {  // same as before for first 2
      PmeParticle *lgd = localGridData[g];
      int nga = 0;
      for(int i=0; i<numLocalAtoms; ++i) {
        if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
          lgd[nga++] = localData[i];
        }
      }
      numGridAtoms[g] = nga;
    }
    for (g=2 ; g<4 ; ++g ) {  // only alchemical atoms for these 2
      PmeParticle *lgd = localGridData[g];
      int nga = 0;
      for(int i=0; i<numLocalAtoms; ++i) {
        if ( localPartition[i] == (g-1) ) {
          lgd[nga++] = localData[i];
        }
      }
      numGridAtoms[g] = nga;
    }
    for (g=4 ; g<numGrids ; ++g ) {  // only non-alchemical atoms 
      // numGrids=5 only if alchElecLambdaStart > 0
      PmeParticle *lgd = localGridData[g];
      int nga = 0;
      for(int i=0; i<numLocalAtoms; ++i) {
        if ( localPartition[i] == 0 ) {
          lgd[nga++] = localData[i];
        }
      }
      numGridAtoms[g] = nga;
    }
  } else if ( selfOn ) {
    if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 1 failed");
    g = 0;
    PmeParticle *lgd = localGridData[g];
    int nga = 0;
    for(int i=0; i<numLocalAtoms; ++i) {
      if ( localPartition[i] == 1 ) {
        lgd[nga++] = localData[i];
      }
    }
    numGridAtoms[g] = nga;
  } else if ( pairOn ) {
    if ( numGrids != 3 ) NAMD_bug("ComputePme::doWork assertion 2 failed");
    g = 0;
    PmeParticle *lgd = localGridData[g];
    int nga = 0;
    for(int i=0; i<numLocalAtoms; ++i) {
      if ( localPartition[i] == 1 || localPartition[i] == 2 ) {
        lgd[nga++] = localData[i];
      }
    }
    numGridAtoms[g] = nga;
    for ( g=1; g<3; ++g ) {
      PmeParticle *lgd = localGridData[g];
      int nga = 0;
      for(int i=0; i<numLocalAtoms; ++i) {
        if ( localPartition[i] == g ) {
          lgd[nga++] = localData[i];
        }
      }
      numGridAtoms[g] = nga;
    }
  } else if ( LJPMEOn ) {
    const CompAtomExt *xExt = patch->getCompAtomExtInfo(); // for dispersion coef
    if ( numGrids != 2 ) NAMD_bug("ComputePme::doWork assertion for LJ-PME failed");
    // Reset localGridData pointers and set atom counts
    // localGridData[1] = localGridData[0];
    // localGridData[0] = localData;
    numGridAtoms[0] = numLocalAtoms;
    numGridAtoms[1] = numLocalAtoms;
    PmeParticle *lgd = localGridData[1];
    for (int i=0;  i < numLocalAtoms;  ++i) {
      lgd[i].x = localData[i].x;
      lgd[i].y = localData[i].y;
      lgd[i].z = localData[i].z;
      lgd[i].cg = xExt[i].dispcoef;  // no scaling needed for dispersion
    }
  } else {
    // This else handles the numGrids==1 case.
    // In this case, localGridData[0] and numGridAtoms[0] aren't set to
    // usable values, so we reset them to point to the base buffer.
    // Expect the calculation to be done on localGridData[0..numGrids],
    // each buffer containing numGridAtoms[0..numGrids].
    if ( numGrids != 1 ) NAMD_bug("ComputePme::doWork assertion 3 failed");
    // localGridData[0] = localData;
    numGridAtoms[0] = numLocalAtoms;
  }

 if ( ! myMgr->doWorkCount ) {
  myMgr->doWorkCount = myMgr->pmeComputes.size();

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
 if ( !  offload )
#endif // NAMD_CUDA
 {
  memset( (void*) myMgr->fz_arr, 0, (myGrid.K3+myGrid.order-1) * sizeof(char) );

  for (int i=0; i<myMgr->q_count; ++i) {
    memset( (void*) (myMgr->q_list[i]), 0, (myGrid.K3+myGrid.order-1) * sizeof(float) );
  }
 }

  for ( g=0; g<numGrids; ++g ) {
    myMgr->evir[g] = 0;
  }

  myMgr->strayChargeErrors = 0;

  myMgr->compute_sequence = sequence();
 }

  if ( sequence() != myMgr->compute_sequence ) NAMD_bug("ComputePme sequence mismatch in doWork()");

  int strayChargeErrors = 0;

  // XXX need self energy for LJ-PME
  // calculate self energy
  const BigReal ewaldcof = ComputeNonbondedUtil::ewaldcof;
  for ( g=0; g<numGrids; ++g ) {
    BigReal selfEnergy = 0;
    data_ptr = localGridData[g];
    for (int i=0; i<numGridAtoms[g]; ++i) {
      selfEnergy += data_ptr->cg * data_ptr->cg;
      ++data_ptr;
    }
    if ( LJPMEOn && 1==g ) {
      const BigReal LJewaldcof = ComputeNonbondedUtil::LJewaldcof;
      double alpha6 = LJewaldcof * LJewaldcof * LJewaldcof;
      alpha6 = alpha6 * alpha6;
      selfEnergy *= (1./12.) * alpha6;
    } else {
      selfEnergy *= -1. * ewaldcof / SQRT_PI;
    }
    myMgr->evir[g][0] += selfEnergy;

    float **q = myMgr->q_arr + g*myMgr->fsize;
    char *f = myMgr->f_arr + g*myMgr->fsize;

    scale_coordinates(localGridData[g], numGridAtoms[g], lattice, myGrid);
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
   if ( offload ) {
    if ( myMgr->cuda_atoms_alloc == 0 ) {  // first call
      int na = myMgr->cuda_atoms_alloc = 1.2 * (myMgr->cuda_atoms_count + 1000);
      cuda_errcheck("before malloc atom data for pme");
      cudaMallocHost((void**) &(myMgr->a_data_host), 7*na*sizeof(float));
      cudaMalloc((void**) &(myMgr->a_data_dev), 7*na*sizeof(float));
      cuda_errcheck("malloc atom data for pme");
      myMgr->cuda_atoms_count = 0;
    }
    cuda_atoms_offset = myMgr->cuda_atoms_count;
    int n = numGridAtoms[g];
    myMgr->cuda_atoms_count += n;
    if ( myMgr->cuda_atoms_count > myMgr->cuda_atoms_alloc ) {
      CkPrintf("Pe %d expanding CUDA PME atoms allocation because %d > %d\n",
                        CkMyPe(), myMgr->cuda_atoms_count, myMgr->cuda_atoms_alloc);
      cuda_errcheck("before malloc expanded atom data for pme");
      int na = myMgr->cuda_atoms_alloc = 1.2 * (myMgr->cuda_atoms_count + 1000);
      const float *a_data_host_old = myMgr->a_data_host;
      cudaMallocHost((void**) &(myMgr->a_data_host), 7*na*sizeof(float));
      cuda_errcheck("malloc expanded host atom data for pme");
      memcpy(myMgr->a_data_host, a_data_host_old, 7*cuda_atoms_offset*sizeof(float));
      cudaFreeHost((void*) a_data_host_old);
      cuda_errcheck("free expanded host atom data for pme");
      cudaFree(myMgr->a_data_dev);
      cuda_errcheck("free expanded dev atom data for pme");
      cudaMalloc((void**) &(myMgr->a_data_dev), 7*na*sizeof(float));
      cuda_errcheck("malloc expanded dev atom data for pme");
    }
    float *a_data_host = myMgr->a_data_host + 7 * cuda_atoms_offset;
    data_ptr = localGridData[g];
    double order_1 = myGrid.order - 1;
    double K1 = myGrid.K1;
    double K2 = myGrid.K2;
    double K3 = myGrid.K3;
    int found_negative = 0;
    for ( int i=0; i<n; ++i ) {
      if ( data_ptr[i].x < 0 || data_ptr[i].y < 0 || data_ptr[i].z < 0 ) {
        found_negative = 1;
        // CkPrintf("low coord: %f %f %f\n", data_ptr[i].x, data_ptr[i].y, data_ptr[i].z);
      }
      double x_int = (int) data_ptr[i].x;
      double y_int = (int) data_ptr[i].y;
      double z_int = (int) data_ptr[i].z;
      a_data_host[7*i  ] = data_ptr[i].x - x_int;  // subtract in double precision
      a_data_host[7*i+1] = data_ptr[i].y - y_int;
      a_data_host[7*i+2] = data_ptr[i].z - z_int;
      a_data_host[7*i+3] = data_ptr[i].cg;
      x_int -= order_1;  if ( x_int < 0 ) x_int += K1;
      y_int -= order_1;  if ( y_int < 0 ) y_int += K2;
      z_int -= order_1;  if ( z_int < 0 ) z_int += K3;
      a_data_host[7*i+4] = x_int;
      a_data_host[7*i+5] = y_int;
      a_data_host[7*i+6] = z_int;
    }
    if ( found_negative ) NAMD_bug("found negative atom coordinate in ComputePme::doWork");
   } else
#endif // NAMD_CUDA
   {
    myRealSpace[g]->set_num_atoms(numGridAtoms[g]);
    myRealSpace[g]->fill_charges(q, myMgr->q_list, myMgr->q_count, strayChargeErrors, f, myMgr->fz_arr, localGridData[g]);
   }
  }
  myMgr->strayChargeErrors += strayChargeErrors;

#ifdef TRACE_COMPUTE_OBJECTS
    traceUserBracketEvent(TRACE_COMPOBJ_IDOFFSET+this->cid, traceObjStartTime, CmiWallTimer());
#endif

 if ( --(myMgr->doWorkCount) == 0 ) {
// cudaDeviceSynchronize();  // XXXX
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
  if ( offload ) {
    ComputePmeMgr::cuda_submit_charges_args args;
    args.mgr = myMgr;
    args.lattice = &lattice;
    args.sequence = sequence();
    CmiLock(ComputePmeMgr::cuda_lock);
    if ( ComputePmeMgr::cuda_busy ) {
      ComputePmeMgr::cuda_submit_charges_deque.push_back(args);
    } else if ( CkMyPe() == deviceCUDA->getMasterPe() ) {
      // avoid adding work to nonbonded data preparation pe
      args.mgr->cuda_submit_charges(*args.lattice, args.sequence);
    } else {
      ComputePmeMgr::cuda_busy = true;
      while ( 1 ) {
        CmiUnlock(ComputePmeMgr::cuda_lock);
        args.mgr->cuda_submit_charges(*args.lattice, args.sequence);
        CmiLock(ComputePmeMgr::cuda_lock);
        if ( ComputePmeMgr::cuda_submit_charges_deque.size() ) {
          args = ComputePmeMgr::cuda_submit_charges_deque.front();
          ComputePmeMgr::cuda_submit_charges_deque.pop_front();
        } else {
          ComputePmeMgr::cuda_busy = false;
          break;
        }
      }
    }
    CmiUnlock(ComputePmeMgr::cuda_lock);
  } else
#endif // NAMD_CUDA
  {
    myMgr->chargeGridReady(lattice,sequence());
  }
 }
 atomsChanged = 0;
}

initialize()

void ComputePme::initialize ( void  )

virtual

Reimplemented from Compute.

Definition at line 2751 of file ComputePme.C.

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

                            {
  if (!(patch = PatchMap::Object()->patch(patchID))) {
    NAMD_bug("ComputePme used with unknown patch.");
  }
  positionBox = patch->registerPositionPickup(this);
  avgPositionBox = patch->registerAvgPositionPickup(this);
  forceBox = patch->registerForceDeposit(this);
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
 if ( offload ) {
  myMgr->cuda_atoms_count += patch->getNumAtoms();
 }
#endif
}

noWork()

int ComputePme::noWork ( )

virtual

Reimplemented from Compute.

Definition at line 3039 of file ComputePme.C.

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

                       {

  if ( patch->flags.doFullElectrostatics ) {
    // In QM/MM simulations, atom charges form QM regions need special treatment.
    if ( qmForcesOn ) {
        return 1;
    }
    if ( ! myMgr->ungridForcesCount && ! myMgr->recipEvirCount ) return 0;  // work to do, enqueue as usual
    myMgr->heldComputes.add(this);
    return 1;  // don't enqueue yet
  }

  positionBox->skip();
  forceBox->skip();

  if ( ++(myMgr->noWorkCount) == myMgr->pmeComputes.size() ) {
    myMgr->noWorkCount = 0;
    myMgr->reduction->submit();
  }

  atomsChanged = 0;

  return 1;  // no work for this step
}

setMgr()

void ComputePme::setMgr ( ComputePmeMgr *  mgr )

inline

Definition at line 58 of file ComputePme.h.

{ myMgr = mgr; }

ungridForces()

void ComputePme::ungridForces

(

)

Definition at line 4090 of file ComputePme.C.

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

Referenced by doWork().

                              {

    if ( sequence() != myMgr->compute_sequence ) NAMD_bug("ComputePme sequence mismatch in ungridForces()");
 
    SimParameters *simParams = Node::Object()->simParameters;

    localResults_alloc.resize(numLocalAtoms* ((numGrids>1 || selfOn)?2:1));
    Vector *localResults = localResults_alloc.begin();
    Vector *gridResults;

    if ( alchOn || lesOn || selfOn || pairOn ) {
      for(int i=0; i<numLocalAtoms; ++i) { localResults[i] = 0.; }
      gridResults = localResults + numLocalAtoms;
    } else {
      gridResults = localResults;
    }

    Vector pairForce = 0.;
    Lattice &lattice = patch->flags.lattice;
    int g = 0;
    if(!simParams->commOnly) {
    for ( g=0; g<numGrids; ++g ) {
#ifdef NETWORK_PROGRESS
      CmiNetworkProgress();
#endif

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
      if ( offload ) {
        int errfound = 0;
        for ( int n=numGridAtoms[g], i=0; i<n; ++i ) {
          // Neither isnan() nor x != x worked when testing on Cray; this does.
          if ( ((int*)f_data_host)[3*i] == 0x7fffffff ) { errfound = 1; }  // CUDA NaN
          gridResults[i].x = f_data_host[3*i];
          gridResults[i].y = f_data_host[3*i+1];
          gridResults[i].z = f_data_host[3*i+2];
        }
        if ( errfound ) {
          int errcount = 0;
          for ( int n=numGridAtoms[g], i=0; i<n; ++i ) {
            float f = f_data_host[3*i];
            if ( ((int*)f_data_host)[3*i] == 0x7fffffff ) {  // CUDA NaN
              ++errcount;
              gridResults[i] = 0.;
            }
          }
          iout << iERROR << "Stray PME grid charges detected: "
                << errcount << " atoms on pe " << CkMyPe() << "\n" << endi;
        }
      } else
#endif // NAMD_CUDA
        {
          myRealSpace[g]->compute_forces(myMgr->q_arr+g*myMgr->fsize, localGridData[g], gridResults);
        }
      scale_forces(gridResults, numGridAtoms[g], lattice);
      
      if (LJPMEOn) {
        if (0==g) {
          // finished loop g==0, next loop gathers
          // LJ-PME force contributions into upper buffer
          gridResults += numLocalAtoms;
        } else {
          // sum LJ-PME forces into electrostatic forces buffer
          for (int i=0;  i < numLocalAtoms;  i++) {
            localResults[i] += gridResults[i];
          }
        }
      } else if (alchOn) {
        float scale = 1.;
        BigReal elecLambdaUp, elecLambdaDown;
        BigReal alchLambda = simParams->getCurrentLambda(patch->flags.step);
        myMgr->alchLambda = alchLambda;
        BigReal alchLambda2 = simParams->getCurrentLambda2(patch->flags.step);
        myMgr->alchLambda2 = alchLambda2;
        elecLambdaUp = simParams->getElecLambda(alchLambda);
        elecLambdaDown = simParams->getElecLambda(1. - alchLambda);
        
        if ( g == 0 ) scale = elecLambdaUp;
        else if ( g == 1 ) scale = elecLambdaDown;
        else if ( g == 2 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);

        if (alchDecouple) {
          if ( g == 2 ) scale = 1 - elecLambdaUp;
          else if ( g == 3 ) scale = 1 - elecLambdaDown;
          else if ( g == 4 ) scale = (elecLambdaUp + elecLambdaDown - 1)*(-1);
        }
        int nga = 0;
        if (!alchDecouple) {
         if (g < 2 ) {
          for(int i=0; i<numLocalAtoms; ++i) {
            if ( localPartition[i] == 0 || localPartition[i] == (g+1) || localPartition[i] == (g+3) ) {
              // (g=0: only partition 0 and partiton 1 and partion 3)
              // (g=1: only partition 0 and partiton 2 and partion 4) 
              localResults[i] += gridResults[nga++] * scale;
            }
           }
         } else {
          for(int i=0; i<numLocalAtoms; ++i) {
            if ( localPartition[i] == 0 ) {
              // (g=2: only partition 0)
              localResults[i] += gridResults[nga++] * scale;
            }
          }
        }              
        } else {  // alchDecouple
          if ( g < 2 ) {
            for(int i=0; i<numLocalAtoms; ++i) {
              if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
                // g = 0: partition 0 or partition 1
                // g = 1: partition 0 or partition 2
                localResults[i] += gridResults[nga++] * scale;
              }
            }
          }
          else {
            for(int i=0; i<numLocalAtoms; ++i) {
              if ( localPartition[i] == (g-1) || localPartition[i] == (g-4)) {
                // g = 2: partition 1 only
                // g = 3: partition 2 only
                // g = 4: partition 0 only
                localResults[i] += gridResults[nga++] * scale;
              }
            }
          }
        }
      } else if ( lesOn ) {
        float scale = 1.;
        if ( alchFepOn ) { 
          BigReal alchLambda = simParams->getCurrentLambda(patch->flags.step);
          myMgr->alchLambda = alchLambda;
          BigReal alchLambda2 = simParams->getCurrentLambda2(patch->flags.step);
          myMgr->alchLambda2 = alchLambda2;
          if ( g == 0 ) scale = alchLambda;
          else if ( g == 1 ) scale = 1. - alchLambda;
        } else if ( lesOn ) {
          scale = 1.0 / (float)lesFactor;
        }
        int nga = 0;
        for(int i=0; i<numLocalAtoms; ++i) {
          if ( localPartition[i] == 0 || localPartition[i] == (g+1) ) {
            localResults[i] += gridResults[nga++] * scale;
          }
        }
      } else if ( selfOn ) {
        PmeParticle *lgd = localGridData[g];
        int nga = 0;
        for(int i=0; i<numLocalAtoms; ++i) {
          if ( localPartition[i] == 1 ) {
            pairForce += gridResults[nga];  // should add up to almost zero
            localResults[i] += gridResults[nga++];
          }
        }
      } else if ( pairOn ) {
        if ( g == 0 ) {
          int nga = 0;
          for(int i=0; i<numLocalAtoms; ++i) {
            if ( localPartition[i] == 1 ) {
              pairForce += gridResults[nga];
            }
            if ( localPartition[i] == 1 || localPartition[i] == 2 ) {
              localResults[i] += gridResults[nga++];
            }
          }
        } else if ( g == 1 ) {
          int nga = 0;
          for(int i=0; i<numLocalAtoms; ++i) {
            if ( localPartition[i] == g ) {
              pairForce -= gridResults[nga];  // should add up to almost zero
              localResults[i] -= gridResults[nga++];
            }
          }
        } else {
          int nga = 0;
          for(int i=0; i<numLocalAtoms; ++i) {
            if ( localPartition[i] == g ) {
              localResults[i] -= gridResults[nga++];
            }
         }
        }
      }
    }
    }

    Vector *results_ptr = localResults;
    
    // add in forces
    {
      Results *r = forceBox->open();
      Force *f = r->f[Results::slow];
      int numAtoms = patch->getNumAtoms();

      if ( ! myMgr->strayChargeErrors && ! simParams->commOnly ) {
        for(int i=0; i<numAtoms; ++i) {
          f[i].x += results_ptr->x;
          f[i].y += results_ptr->y;
          f[i].z += results_ptr->z;
          ++results_ptr;
        }
      }
      forceBox->close(&r);
    }

    if ( pairOn || selfOn ) {
        ADD_VECTOR_OBJECT(myMgr->reduction,REDUCTION_PAIR_ELECT_FORCE,pairForce);
    }

}

Friends And Related Function Documentation

ComputePmeMgr

friend class ComputePmeMgr

friend

Definition at line 60 of file ComputePme.h.

---

The documentation for this class was generated from the following files:

• ComputePme.h
• ComputePme.C