ComputeMgr.C

Go to the documentation of this file.

#include "dlloader.h"
#include "CudaGlobalMasterClient.h"
#include "InfoStream.h"
#include "ProcessorPrivate.h"
#include "middle-conv.h"

//#define DEBUGM
#define MIN_DEBUG_LEVEL 3
#include "Debug.h"

#include "BOCgroup.h"
#include "ComputeMgr.decl.h"
#include "ComputeMgr.h"
#include "ProxyMgr.decl.h"
#include "ProxyMgr.h"

#include "Node.h"
#include "ComputeMap.h"
#include "PatchMap.h"
#include "PatchMap.inl"

#include "Compute.h"
#include "ComputeNonbondedUtil.h"
#include "ComputeNonbondedSelf.h"
#include "ComputeNonbondedPair.h"
#include "ComputeAngles.h"
#include "ComputeDihedrals.h"
#include "ComputeImpropers.h"
#include "ComputeThole.h"
#include "ComputeAniso.h"
#include "ComputeCrossterms.h"
// JLai
#include "ComputeGromacsPair.h"
#include "ComputeBonds.h"
#include "ComputeNonbondedCUDAExcl.h"
#include "ComputeFullDirect.h"
#include "ComputeGlobal.h"
#include "ComputeGlobalMsgs.h"
#include "ComputeExt.h"
#include "ComputeQM.h"
#include "ComputeGBISser.h"
#include "ComputeLCPO.h"
#include "ComputeFmmSerial.h"
#include "ComputeMsmSerial.h"
#include "ComputeLjPmeSerial.h"
#include "ComputeMsmMsa.h"
#include "ComputeMsm.h"
#include "ComputeDPMTA.h"
#include "ComputeDPME.h"
#include "ComputeDPMEMsgs.h"
#include "ComputePme.h"
// #ifdef NAMD_CUDA
#include "ComputePmeCUDA.h"
#include "ComputeCUDAMgr.h"
#include "CudaComputeNonbonded.h"
#include "ComputePmeCUDAMgr.h"
// #endif
#include "ComputeEwald.h"
#include "ComputeEField.h"
/* BEGIN gf */
#include "ComputeGridForce.h"
/* END gf */
#include "ComputeStir.h"
#include "ComputeSphericalBC.h"
#include "ComputeCylindricalBC.h"
#include "ComputeTclBC.h"
#include "ComputeRestraints.h"
#include "ComputeConsForce.h"
#include "ComputeConsForceMsgs.h"
#include "WorkDistrib.h"

#include "LdbCoordinator.h"

/* include all of the specific masters we need here */
#include "FreeEnergyEnums.h"
#include "FreeEnergyAssert.h"
#include "FreeEnergyGroup.h"
#include "FreeEnergyVector.h"
#include "FreeEnergyRestrain.h"
#include "FreeEnergyRMgr.h"
#include "FreeEnergyLambda.h"
#include "FreeEnergyLambdMgr.h"

#include "GlobalMasterTest.h"
#include "GlobalMasterIMD.h"
#include "GlobalMasterTcl.h"
#include "GlobalMasterSMD.h"
#include "GlobalMasterTMD.h"
#include "GlobalMasterSymmetry.h"
#include "GlobalMasterEasy.h"
#include "GlobalMasterMisc.h"
#include "GlobalMasterFreeEnergy.h"
#include "GlobalMasterColvars.h"

#include "PatchData.h"
#include "NamdEventsProfiling.h"
#include "DeviceCUDA.h"

#include "CudaGlobalMasterServer.h"
#include "strlib.h"

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
#ifdef WIN32
#define __thread __declspec(thread)
#endif
extern __thread DeviceCUDA *deviceCUDA;
#endif

ComputeMgr::ComputeMgr()
{
    CkpvAccess(BOCclass_group).computeMgr = thisgroup;
    computeGlobalObject = 0;
    computeGlobalResultsMsgSeq = -1;
    computeGlobalResultsMsgMasterSeq = -1;
    computeDPMEObject = 0;
    computeEwaldObject = 0;
    computeNonbondedWorkArrays = new ComputeNonbondedWorkArrays;
    skipSplitting = 0;
    masterServerObject = NULL;
}

ComputeMgr::~ComputeMgr(void)
{
    delete computeNonbondedWorkArrays;
    if (masterServerObject != NULL) delete masterServerObject;
    for (auto& loader: CudaGlobalMasterClientDlloaders) {
      if (loader) {
        iout << iINFO << "Close library " << loader->LibName() << "\n" << endi;
        loader->DLCloseLib();
      }
    }
}

void ComputeMgr::updateComputes(int ep, CkGroupID chareID)
{
    updateComputesReturnEP = ep;
    updateComputesReturnChareID = chareID;
    updateComputesCount = CkNumPes();

    if (CkMyPe())
    {
        NAMD_bug("updateComputes signaled on wrong Pe!");
    }

    CkStartQD(CkIndex_ComputeMgr::updateComputes2((CkQdMsg*)0),&thishandle);
}

void ComputeMgr::updateComputes2(CkQdMsg *msg)
{
    delete msg;

    CProxy_WorkDistrib wd(CkpvAccess(BOCclass_group).workDistrib);
    WorkDistrib  *workDistrib = wd.ckLocalBranch();
    workDistrib->saveComputeMapChanges(CkIndex_ComputeMgr::updateComputes3(),thisgroup);
}

void ComputeMgr::updateComputes3()
{
    if ( skipSplitting ) {
      CProxy_ComputeMgr(thisgroup).updateLocalComputes();
    } else {
      CProxy_ComputeMgr(thisgroup).splitComputes();
      skipSplitting = 1;
    }
}

void ComputeMgr::splitComputes()
{
  if ( ! CkMyRank() ) {
    ComputeMap *computeMap = ComputeMap::Object();
    const int nc = computeMap->numComputes();

    for (int i=0; i<nc; i++) {
      int nnp = computeMap->newNumPartitions(i);
      if ( nnp > 0 ) {
        if ( computeMap->numPartitions(i) != 1 ) {
          CkPrintf("Warning: unable to partition compute %d\n", i);
          computeMap->setNewNumPartitions(i,0);
          continue;
        }
        //CkPrintf("splitting compute %d by %d\n",i,nnp);
        computeMap->setNumPartitions(i,nnp);
        if (computeMap->newNode(i) == -1) {
          computeMap->setNewNode(i,computeMap->node(i));
        }
        for ( int j=1; j<nnp; ++j ) {
          int newcid = computeMap->cloneCompute(i,j);
          //CkPrintf("compute %d partition %d is %d\n",i,j,newcid);
        }
      }
    }
    computeMap->extendPtrs();
  }

  if (!CkMyPe())
  {
    CkStartQD(CkIndex_ComputeMgr::splitComputes2((CkQdMsg*)0), &thishandle);
  }
}

void ComputeMgr::splitComputes2(CkQdMsg *msg)
{
    delete msg;
    CProxy_ComputeMgr(thisgroup).updateLocalComputes();
}

void ComputeMgr::updateLocalComputes()
{
    ComputeMap *computeMap = ComputeMap::Object();
    CProxy_ProxyMgr pm(CkpvAccess(BOCclass_group).proxyMgr);
    ProxyMgr *proxyMgr = pm.ckLocalBranch();
    LdbCoordinator *ldbCoordinator = LdbCoordinator::Object();

     computeFlag.resize(0);

    const int nc = computeMap->numComputes();
    for (int i=0; i<nc; i++) {

        if ( computeMap->node(i) == CkMyPe() &&
             computeMap->newNumPartitions(i) > 1 ) {
           Compute *c = computeMap->compute(i);
           ldbCoordinator->Migrate(c->ldObjHandle,CkMyPe());
           delete c;
           computeMap->registerCompute(i,NULL);
           if ( computeMap->newNode(i) == CkMyPe() ) computeFlag.add(i); 
        } else
        if (computeMap->newNode(i) == CkMyPe() && computeMap->node(i) != CkMyPe())
        {
            computeFlag.add(i);
            for (int n=0; n < computeMap->numPids(i); n++)
            {
                proxyMgr->createProxy(computeMap->pid(i,n));
            }
        }
        else if (computeMap->node(i) == CkMyPe() &&
                 (computeMap->newNode(i) != -1 && computeMap->newNode(i) != CkMyPe() ))
        {
            // CkPrintf("delete compute %d on pe %d\n",i,CkMyPe());
            delete computeMap->compute(i);
            computeMap->registerCompute(i,NULL);
        }
    }

    if (!CkMyPe())
    {
        CkStartQD(CkIndex_ComputeMgr::updateLocalComputes2((CkQdMsg*)0), &thishandle);
    }
}

void
ComputeMgr::updateLocalComputes2(CkQdMsg *msg)
{
    delete msg;
    CProxy_ComputeMgr(thisgroup).updateLocalComputes3();
}

void
ComputeMgr::updateLocalComputes3()
{
    ComputeMap *computeMap = ComputeMap::Object();
    CProxy_ProxyMgr pm(CkpvAccess(BOCclass_group).proxyMgr);
    ProxyMgr *proxyMgr = pm.ckLocalBranch();

    ProxyMgr::nodecount = 0;

    const int nc = computeMap->numComputes();

    if ( ! CkMyRank() ) {
      for (int i=0; i<nc; i++) {
        computeMap->setNewNumPartitions(i,0);
        if (computeMap->newNode(i) != -1) {
          computeMap->setNode(i,computeMap->newNode(i));
          computeMap->setNewNode(i,-1);
        }
      }
    }
 
    for(int i=0; i<computeFlag.size(); i++) createCompute(computeFlag[i], computeMap);
    computeFlag.clear();

    proxyMgr->removeUnusedProxies();

    if (!CkMyPe())
    {
        CkStartQD(CkIndex_ComputeMgr::updateLocalComputes4((CkQdMsg*)0), &thishandle);
    }
}

void
ComputeMgr::updateLocalComputes4(CkQdMsg *msg)
{
    delete msg;
    CProxy_ComputeMgr(thisgroup).updateLocalComputes5();

    // store the latest compute map
           SimParameters *simParams = Node::Object()->simParameters;
    if (simParams->storeComputeMap) {
      ComputeMap *computeMap = ComputeMap::Object();
      computeMap->saveComputeMap(simParams->computeMapFilename);
    }
}

#if 0
int firstphase = 1;
#endif

void
ComputeMgr::updateLocalComputes5()
{
    if ( ! CkMyRank() ) {
      ComputeMap::Object()->checkMap();
      PatchMap::Object()->checkMap();
    }

    // we always use the centralized building of spanning tree
    // distributed building of ST called in Node.C only
    if (proxySendSpanning || proxyRecvSpanning)
        ProxyMgr::Object()->buildProxySpanningTree2();

    // this code needs to be turned on if we want to
    // shift the creation of ST to the load balancer

#if 0
    if (proxySendSpanning || proxyRecvSpanning)
    {
        if (firstphase)
            ProxyMgr::Object()->buildProxySpanningTree2();
        else
            if (CkMyPe() == 0)
                ProxyMgr::Object()->sendSpanningTrees();

        firstphase = 0;
    }
#endif

    if (!CkMyPe())
        CkStartQD(CkIndex_ComputeMgr::doneUpdateLocalComputes(), &thishandle);
}

void ComputeMgr::doneUpdateLocalComputes()
{

//  if (!--updateComputesCount) {
    DebugM(4, "doneUpdateLocalComputes on Pe("<<CkMyPe()<<")\n");
    void *msg = CkAllocMsg(0,0,0);
    CkSendMsgBranch(updateComputesReturnEP,msg,0,updateComputesReturnChareID);
//  }
}

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
// Helper functions for creating and getting pointers to CUDA computes
CudaComputeNonbonded* getCudaComputeNonbonded() {
  return ComputeCUDAMgr::getComputeCUDAMgr()->getCudaComputeNonbonded();
}

CudaComputeNonbonded* createCudaComputeNonbonded(ComputeID c) {
  return ComputeCUDAMgr::getComputeCUDAMgr()->createCudaComputeNonbonded(c);
}

#ifdef BONDED_CUDA
ComputeBondedCUDA* getComputeBondedCUDA() {
  return ComputeCUDAMgr::getComputeCUDAMgr()->getComputeBondedCUDA();
}

ComputeBondedCUDA* createComputeBondedCUDA(ComputeID c, ComputeMgr* computeMgr) {
  return ComputeCUDAMgr::getComputeCUDAMgr()->createComputeBondedCUDA(c, computeMgr);
}
#endif
#endif

//
void
ComputeMgr::createCompute(ComputeID i, ComputeMap *map)
{
    Compute *c;
    PatchID pid2[2];
    PatchIDList pids;
    int trans2[2];
    SimParameters *simParams = Node::Object()->simParameters;

    PatchID pid8[8];
    int trans8[8];
#ifdef NODEGROUP_FORCE_REGISTER      
    CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
    PatchData *patchData = cpdata.ckLocalBranch();
    suspendCounter=&(patchData->suspendCounter);
#endif    

    switch ( map->type(i) )
    {
    case computeNonbondedSelfType:
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
        getCudaComputeNonbonded()->registerComputeSelf(i, map->computeData[i].pids[0].pid);
#else
        c = new ComputeNonbondedSelf(i,map->computeData[i].pids[0].pid,
                                     computeNonbondedWorkArrays,
                                     map->partition(i),map->partition(i)+1,
                                     map->numPartitions(i)); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
#endif
        break;
    case computeLCPOType:
        for (int j = 0; j < 8; j++) {
          pid8[j] = map->computeData[i].pids[j].pid;
          trans8[j] = map->computeData[i].pids[j].trans;
        }
        c = new ComputeLCPO(i,pid8,trans8,
             computeNonbondedWorkArrays,
             map->partition(i),map->partition(i)+1,
             map->numPartitions(i), 8);
        map->registerCompute(i,c);
        c->initialize();
      
        break;
    case computeNonbondedPairType:
        pid2[0] = map->computeData[i].pids[0].pid;
        trans2[0] = map->computeData[i].pids[0].trans;
        pid2[1] = map->computeData[i].pids[1].pid;
        trans2[1] = map->computeData[i].pids[1].trans;
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
        getCudaComputeNonbonded()->registerComputePair(i, pid2, trans2);
#else
        c = new ComputeNonbondedPair(i,pid2,trans2,
                                     computeNonbondedWorkArrays,
                                     map->partition(i),map->partition(i)+1,
                                     map->numPartitions(i)); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
#endif
        break;
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
    case computeNonbondedCUDA2Type:
      c = createCudaComputeNonbonded(i);
      map->registerCompute(i,c);
      // NOTE: initialize() is called at the end of createComputes(),
      //       after all computes have been created
      //c->initialize();
      break;
#ifdef BONDED_CUDA
    case computeBondedCUDAType:
      c = createComputeBondedCUDA(i, this);
      map->registerCompute(i,c);
      break;
#endif
#endif
    case computeExclsType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined(NAMD_HIP))
        if (simParams->bondedCUDA & 16)
        {
          PatchMap::Object()->basePatchIDList(map->computeData[i].node, pids);
          getComputeBondedCUDA()->registerCompute(map->computeData[i].node, map->type(i), pids);
        } else
#endif
        {
        PatchMap::Object()->basePatchIDList(CkMyPe(),pids);
        c = new ComputeExcls(i,pids); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
      }
      break;
    case computeBondsType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined (NAMD_HIP))
        if (simParams->bondedCUDA & 1)
        {
          PatchMap::Object()->basePatchIDList(map->computeData[i].node, pids);
          getComputeBondedCUDA()->registerCompute(map->computeData[i].node, map->type(i), pids);
        } else
#endif
        {
          PatchMap::Object()->basePatchIDList(CkMyPe(),pids);
          c = new ComputeBonds(i,pids); // unknown delete
          map->registerCompute(i,c);
          c->initialize();
        }
        break;
    case computeAnglesType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined (NAMD_HIP))
        if (simParams->bondedCUDA & 2)
        {
          PatchMap::Object()->basePatchIDList(map->computeData[i].node, pids);
          getComputeBondedCUDA()->registerCompute(map->computeData[i].node, map->type(i), pids);
        } else
#endif
        {
          PatchMap::Object()->basePatchIDList(CkMyPe(),pids);
          c = new ComputeAngles(i,pids); // unknown delete
          map->registerCompute(i,c);
          c->initialize();
        }
        break;
    case computeDihedralsType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined (NAMD_HIP))
        if (simParams->bondedCUDA & 4)
        {
          PatchMap::Object()->basePatchIDList(map->computeData[i].node, pids);
          getComputeBondedCUDA()->registerCompute(map->computeData[i].node, map->type(i), pids);
        } else
#endif
        {
          PatchMap::Object()->basePatchIDList(CkMyPe(),pids);
          c = new ComputeDihedrals(i,pids); // unknown delete
          map->registerCompute(i,c);
          c->initialize();
        }
        break;
    case computeImpropersType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined (NAMD_HIP))
        if (simParams->bondedCUDA & 8)
        {
          PatchMap::Object()->basePatchIDList(map->computeData[i].node, pids);
          getComputeBondedCUDA()->registerCompute(map->computeData[i].node, map->type(i), pids);
        } else
#endif
        {
          PatchMap::Object()->basePatchIDList(CkMyPe(),pids);
          c = new ComputeImpropers(i,pids); // unknown delete
          map->registerCompute(i,c);
          c->initialize();
        }
        break;
    case computeTholeType:
        PatchMap::Object()->basePatchIDList(CkMyPe(),pids);
        c = new ComputeThole(i,pids); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeAnisoType:
        PatchMap::Object()->basePatchIDList(CkMyPe(),pids);
        c = new ComputeAniso(i,pids); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeCrosstermsType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined (NAMD_HIP))
        if (simParams->bondedCUDA & 32)
        {
          PatchMap::Object()->basePatchIDList(map->computeData[i].node, pids);
          getComputeBondedCUDA()->registerCompute(map->computeData[i].node, map->type(i), pids);
        } else
#endif
        {
          PatchMap::Object()->basePatchIDList(CkMyPe(),pids);
          c = new ComputeCrossterms(i,pids); // unknown delete
          map->registerCompute(i,c);
          c->initialize();
        }
        break;
        // JLai
    case computeGromacsPairType:
        PatchMap::Object()->basePatchIDList(CkMyPe(),pids);
              c = new ComputeGromacsPair(i,pids); // unknown delete
              map->registerCompute(i,c);
              c->initialize();
              break;
  case computeSelfGromacsPairType:
        c = new ComputeSelfGromacsPair(i,map->computeData[i].pids[0].pid); // unknown delete
              map->registerCompute(i,c);
              c->initialize();
              break;
        // End of JLai
    case computeSelfExclsType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined (NAMD_HIP))
        if (simParams->bondedCUDA & 16)
        {
          getComputeBondedCUDA()->registerSelfCompute(map->computeData[i].node, map->type(i), map->computeData[i].pids[0].pid);
        } else
#endif
        {
          c = new ComputeSelfExcls(i,map->computeData[i].pids[0].pid);
          map->registerCompute(i,c);
          c->initialize();
        }
        break;
    case computeSelfBondsType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined (NAMD_HIP))
        if (simParams->bondedCUDA & 1)
        {
          getComputeBondedCUDA()->registerSelfCompute(map->computeData[i].node, map->type(i), map->computeData[i].pids[0].pid);
        } else
#endif
        {
          c = new ComputeSelfBonds(i,map->computeData[i].pids[0].pid);
          map->registerCompute(i,c);
          c->initialize();
        }
        break;
    case computeSelfAnglesType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined (NAMD_HIP))
        if (simParams->bondedCUDA & 2)
        {
          getComputeBondedCUDA()->registerSelfCompute(map->computeData[i].node, map->type(i), map->computeData[i].pids[0].pid);
        } else
#endif
        {
          c = new ComputeSelfAngles(i,map->computeData[i].pids[0].pid);
          map->registerCompute(i,c);
          c->initialize();
        }
        break;
    case computeSelfDihedralsType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined (NAMD_HIP))
        if (simParams->bondedCUDA & 4)
        {
          getComputeBondedCUDA()->registerSelfCompute(map->computeData[i].node, map->type(i), map->computeData[i].pids[0].pid);
        } else
#endif
        {
          c = new ComputeSelfDihedrals(i,map->computeData[i].pids[0].pid);
          map->registerCompute(i,c);
          c->initialize();
        }
        break;
    case computeSelfImpropersType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined (NAMD_HIP))
        if (simParams->bondedCUDA & 8)
        {
          getComputeBondedCUDA()->registerSelfCompute(map->computeData[i].node, map->type(i), map->computeData[i].pids[0].pid);
        } else
#endif
        {
          c = new ComputeSelfImpropers(i,map->computeData[i].pids[0].pid);
          map->registerCompute(i,c);
          c->initialize();
        }
        break;
    case computeSelfTholeType:
        c = new ComputeSelfThole(i,map->computeData[i].pids[0].pid);
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeSelfAnisoType:
        c = new ComputeSelfAniso(i,map->computeData[i].pids[0].pid);
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeSelfCrosstermsType:
#if defined(BONDED_CUDA) && (defined(NAMD_CUDA) || defined (NAMD_HIP))
        if (simParams->bondedCUDA & 32)
        {
          getComputeBondedCUDA()->registerSelfCompute(map->computeData[i].node, map->type(i), map->computeData[i].pids[0].pid);
        } else
#endif
        {
          c = new ComputeSelfCrossterms(i,map->computeData[i].pids[0].pid);
          map->registerCompute(i,c);
          c->initialize();
        }
        break;
#ifdef DPMTA
    case computeDPMTAType:
        c = new ComputeDPMTA(i); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
#endif
#ifdef DPME
    case computeDPMEType:
        c = computeDPMEObject = new ComputeDPME(i,this); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
#endif
    case computePmeType:
        c = new ComputePme(i,map->computeData[i].pids[0].pid); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
    case computePmeCUDAType:
        // PatchMap::Object()->basePatchIDList(CkMyPe(),pids);
        // c = new ComputePmeCUDA(i, pids);
        c = new ComputePmeCUDA(i, map->computeData[i].pids[0].pid);
        map->registerCompute(i,c);
        c->initialize();
        break;
#endif
    case computeEwaldType:
        c = computeEwaldObject = new ComputeEwald(i,this); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeFullDirectType:
        c = new ComputeFullDirect(i); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeGlobalType:
        c = computeGlobalObject = new ComputeGlobal(i,this); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeStirType:
        c = new ComputeStir(i,map->computeData[i].pids[0].pid); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeExtType:
        c = new ComputeExt(i); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeQMType:
        c = new ComputeQM(i);
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeGBISserType: //gbis serial
        c = new ComputeGBISser(i);
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeFmmType: // FMM serial
        c = new ComputeFmmSerial(i);
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeMsmSerialType: // MSM serial
        c = new ComputeMsmSerial(i);
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeLjPmeSerialType: // LJ-PME serial
        c = new ComputeLjPmeSerial(i);
        map->registerCompute(i,c);
        c->initialize();
        break;
#ifdef CHARM_HAS_MSA
    case computeMsmMsaType: // MSM parallel long-range part using MSA
        c = new ComputeMsmMsa(i);
        map->registerCompute(i,c);
        c->initialize();
        break;
#endif
    case computeMsmType: // MSM parallel
        c = new ComputeMsm(i);
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeEFieldType:
        c = new ComputeEField(i,map->computeData[i].pids[0].pid); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
        /* BEGIN gf */
    case computeGridForceType:
        c = new ComputeGridForce(i,map->computeData[i].pids[0].pid);
        map->registerCompute(i,c);
        c->initialize();
        break;
        /* END gf */
    case computeSphericalBCType:
        c = new ComputeSphericalBC(i,map->computeData[i].pids[0].pid); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeCylindricalBCType:
        c = new ComputeCylindricalBC(i,map->computeData[i].pids[0].pid); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeTclBCType:
        c = new ComputeTclBC(i); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeRestraintsType:
        c = new ComputeRestraints(i,map->computeData[i].pids[0].pid); // unknown delete
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeConsForceType:
        c = new ComputeConsForce(i,map->computeData[i].pids[0].pid);
        map->registerCompute(i,c);
        c->initialize();
        break;
    case computeConsTorqueType:
        c = new ComputeConsTorque(i,map->computeData[i].pids[0].pid);
        map->registerCompute(i,c);
        c->initialize();
        break;
    default:
        NAMD_bug("Unknown compute type in ComputeMgr::createCompute().");
        break;
    }
}

void registerUserEventsForAllComputeObjs()
{
#ifdef TRACE_COMPUTE_OBJECTS
    ComputeMap *map = ComputeMap::Object();
    PatchMap *pmap = PatchMap::Object();     
    char user_des[50];
    int p1, p2;
    int adim, bdim, cdim;
    int t1, t2;
    int x1, y1, z1, x2, y2, z2;
    int dx, dy, dz;
    for (int i=0; i<map->numComputes(); i++)
    {
        memset(user_des, 0, 50);
        switch ( map->type(i) )
        {
        case computeNonbondedSelfType:
            sprintf(user_des, "computeNonBondedSelfType_%d_pid_%d", i, map->pid(i,0));
            break;
        case computeLCPOType:
            sprintf(user_des, "computeLCPOType_%d_pid_%d", i, map->pid(i,0));
            break;
        case computeNonbondedPairType:
            adim = pmap->gridsize_a();
            bdim = pmap->gridsize_b();
            cdim = pmap->gridsize_c();
            p1 = map->pid(i, 0);
            t1 = map->trans(i, 0);
            x1 = pmap->index_a(p1) + adim * Lattice::offset_a(t1);
            y1 = pmap->index_b(p1) + bdim * Lattice::offset_b(t1);
            z1 = pmap->index_c(p1) + cdim * Lattice::offset_c(t1);
            p2 = map->pid(i, 1);
            t2 = map->trans(i, 1);
            x2 = pmap->index_a(p2) + adim * Lattice::offset_a(t2);
            y2 = pmap->index_b(p2) + bdim * Lattice::offset_b(t2);
            z2 = pmap->index_c(p2) + cdim * Lattice::offset_c(t2);
            dx = abs(x1-x2);
            dy = abs(y1-y2);
            dz = abs(z1-z2);
            sprintf(user_des, "computeNonBondedPairType_%d(%d,%d,%d)", i, dx,dy,dz);
            break;
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
#ifdef BONDED_CUDA
        case computeBondedCUDAType:
            sprintf(user_des, "computeBondedCUDAType_%d", i);
            break;
#endif
#endif
        case computeExclsType:
            sprintf(user_des, "computeExclsType_%d", i);
            break;
        case computeBondsType:
            sprintf(user_des, "computeBondsType_%d", i);
            break;
        case computeAnglesType:
            sprintf(user_des, "computeAnglesType_%d", i);
            break;
        case computeDihedralsType:
            sprintf(user_des, "computeDihedralsType_%d", i);
            break;
        case computeImpropersType:
            sprintf(user_des, "computeImpropersType_%d", i);
            break;
        case computeTholeType:
            sprintf(user_des, "computeTholeType_%d", i);
            break;
        case computeAnisoType:
            sprintf(user_des, "computeAnisoType_%d", i);
            break;
        case computeCrosstermsType:
            sprintf(user_des, "computeCrosstermsType_%d", i);
            break;
        case computeSelfExclsType:
            sprintf(user_des, "computeSelfExclsType_%d", i);
            break;
        case computeSelfBondsType:
            sprintf(user_des, "computeSelfBondsType_%d", i);
            break;
        case computeSelfAnglesType:
            sprintf(user_des, "computeSelfAnglesType_%d", i);
            break;
        case computeSelfDihedralsType:
            sprintf(user_des, "computeSelfDihedralsType_%d", i);
            break;
        case computeSelfImpropersType:
            sprintf(user_des, "computeSelfImpropersType_%d", i);
            break;
        case computeSelfTholeType:
            sprintf(user_des, "computeSelfTholeType_%d", i);
            break;
        case computeSelfAnisoType:
            sprintf(user_des, "computeSelfAnisoType_%d", i);
            break;
        case computeSelfCrosstermsType:
            sprintf(user_des, "computeSelfCrosstermsType_%d", i);
            break;
#ifdef DPMTA
        case computeDPMTAType:
            sprintf(user_des, "computeDPMTAType_%d", i);
            break;
#endif
#ifdef DPME
        case computeDPMEType:
            sprintf(user_des, "computeDPMEType_%d", i);
            break;
#endif
        case computePmeType:
            sprintf(user_des, "computePMEType_%d", i);
            break;
#if defined(NAMD_CUDA) || defined(NAMD_HIP)
        case computePmeCUDAType:
            sprintf(user_des, "computePMECUDAType_%d", i);
            break;
#endif
        case computeEwaldType:
            sprintf(user_des, "computeEwaldType_%d", i);
            break;
        case computeFullDirectType:
            sprintf(user_des, "computeFullDirectType_%d", i);
            break;
        case computeGlobalType:
            sprintf(user_des, "computeGlobalType_%d", i);
            break;
        case computeStirType:
            sprintf(user_des, "computeStirType_%d", i);
            break;
        case computeExtType:
            sprintf(user_des, "computeExtType_%d", i);
            break;
        case computeQMType:
            sprintf(user_des, "computeQMType_%d", i);
            break;
        case computeEFieldType:
            sprintf(user_des, "computeEFieldType_%d", i);
            break;
            /* BEGIN gf */
        case computeGridForceType:
            sprintf(user_des, "computeGridForceType_%d", i);
            break;
            /* END gf */
        case computeSphericalBCType:
            sprintf(user_des, "computeSphericalBCType_%d", i);
            break;
        case computeCylindricalBCType:
            sprintf(user_des, "computeCylindricalBCType_%d", i);
            break;
        case computeTclBCType:
            sprintf(user_des, "computeTclBCType_%d", i);
            break;
        case computeRestraintsType:
            sprintf(user_des, "computeRestraintsType_%d", i);
            break;
        case computeConsForceType:
            sprintf(user_des, "computeConsForceType_%d", i);
            break;
        case computeConsTorqueType:
            sprintf(user_des, "computeConsTorqueType_%d", i);
            break;
        default:
            NAMD_bug("Unknown compute type in ComputeMgr::registerUserEventForAllComputeObjs().");
            break;
        }
        int user_des_len = strlen(user_des);
        char *user_des_cst = new char[user_des_len+1];
        memcpy(user_des_cst, user_des, user_des_len);
        user_des_cst[user_des_len] = 0;
        //Since the argument in traceRegisterUserEvent is supposed
        //to be a const string which will not be copied inside the
        //function when a new user event is created, user_des_cst 
        //has to be allocated in heap.
        int reEvenId = traceRegisterUserEvent(user_des_cst, TRACE_COMPOBJ_IDOFFSET+i);
        //printf("Register user event (%s) with id (%d)\n", user_des, reEvenId);
    }
#else
    return;
#endif
}

void
ComputeMgr::createComputes(ComputeMap *map)
{
// #ifdef NAMD_CUDA
//     int ComputePmeCUDACounter = 0;
// #endif
    Node *node = Node::Object();
    SimParameters *simParams = node->simParameters;
    int myNode = node->myid();

    if ( simParams->globalForcesOn && !myNode )
    {
        DebugM(4,"Mgr running on Node "<<CkMyPe()<<"\n");
        /* create a master server to allow multiple masters */
        masterServerObject = new GlobalMasterServer(this,
                PatchMap::Object()->numNodesWithPatches());

        #ifdef NODEGROUP_FORCE_REGISTER
        CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
        PatchData *patchData = cpdata.ckLocalBranch();
        patchData->master_mgr = this;
        #endif

        /* create the individual global masters */
        // masterServerObject->addClient(new GlobalMasterTest());
        if (simParams->tclForcesOn)
            masterServerObject->addClient(new GlobalMasterTcl());
        if (simParams->IMDon && ! (simParams->IMDignore || simParams->IMDignoreForces) )
            masterServerObject->addClient(new GlobalMasterIMD());
        // SMD is implemented on GPU resident version of NAMD (NAMD3)
        if (simParams->SMDOn && !simParams->CUDASOAintegrateMode)
            masterServerObject->addClient(
                new GlobalMasterSMD(simParams->SMDk, simParams->SMDk2,
                                    simParams->SMDVel,
                                    simParams->SMDDir, simParams->SMDOutputFreq,
                                    simParams->firstTimestep, simParams->SMDFile,
                                    node->molecule->numAtoms)
            );
            
        if (simParams->symmetryOn && 
          (simParams->firstTimestep < simParams->symmetryLastStep || 
          simParams->symmetryLastStep == -1))
            masterServerObject->addClient(new GlobalMasterSymmetry());    
        if (simParams->TMDOn)
            masterServerObject->addClient(new GlobalMasterTMD());
        if (simParams->miscForcesOn)
            masterServerObject->addClient(new GlobalMasterMisc());
        if ( simParams->freeEnergyOn )
            masterServerObject->addClient(new GlobalMasterFreeEnergy());
        if ( simParams->colvarsOn )
            masterServerObject->addClient(new GlobalMasterColvars());

    }

    if ( !myNode && simParams->IMDon && (simParams->IMDignore || simParams->IMDignoreForces) ) {
      // GlobalMasterIMD constructor saves pointer to node->IMDOutput object
      new GlobalMasterIMD();
    }

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
    bool deviceIsMine = ( deviceCUDA->getMasterPe() == CkMyPe() );
#ifdef BONDED_CUDA
    // Place bonded forces on Pe different from non-bonded forces
    int bondedMasterPe = deviceCUDA->getMasterPe();
    // for (int i=0;i < deviceCUDA->getNumPesSharingDevice();i++) {
    //   int pe = deviceCUDA->getPesSharingDevice(i);
    //   if (pe != deviceCUDA->getMasterPe()) {
    //     bondedMasterPe = pe;
    //   }
    // }
    bool deviceIsMineBonded = (CkMyPe() == bondedMasterPe);
#endif
#endif

    for (int i=0; i < map->nComputes; i++)
    {
        if ( ! ( i % 100 ) )
        {
        }

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
        switch ( map->type(i) )
        {
          // case computePmeCUDAType:
          //   // Only create single ComputePmeCUDA object per Pe
          //  if ( map->computeData[i].node != myNode ) continue;
          //  if (ComputePmeCUDACounter > 0) continue;
          //  ComputePmeCUDACounter++;
          //  break;
          case computeNonbondedSelfType:
            if ( ! deviceIsMine ) continue;
            if ( ! deviceCUDA->device_shared_with_pe(map->computeData[i].node) ) continue;
          break;

          case computeNonbondedPairType:
            if ( ! deviceIsMine ) continue;
            if ( ! deviceCUDA->device_shared_with_pe(map->computeData[i].node) ) continue;
          break;

#ifdef BONDED_CUDA
          case computeSelfBondsType:
          case computeBondsType:
            if (simParams->bondedCUDA & 1) {
              if ( ! deviceIsMineBonded ) continue;
              if ( ! deviceCUDA->device_shared_with_pe(map->computeData[i].node) ) continue;
            } else {
              if ( map->computeData[i].node != myNode ) continue;
            }
          break;

          case computeSelfAnglesType:
          case computeAnglesType:
            if (simParams->bondedCUDA & 2) {
              if ( ! deviceIsMineBonded ) continue;
              if ( ! deviceCUDA->device_shared_with_pe(map->computeData[i].node) ) continue;
            } else {
              if ( map->computeData[i].node != myNode ) continue;
            }
          break;

          case computeSelfDihedralsType:
          case computeDihedralsType:
            if (simParams->bondedCUDA & 4) {
              if ( ! deviceIsMineBonded ) continue;
              if ( ! deviceCUDA->device_shared_with_pe(map->computeData[i].node) ) continue;
            } else {
              if ( map->computeData[i].node != myNode ) continue;
            }
          break;

          case computeSelfImpropersType:
          case computeImpropersType:
            if (simParams->bondedCUDA & 8) {
              if ( ! deviceIsMineBonded ) continue;
              if ( ! deviceCUDA->device_shared_with_pe(map->computeData[i].node) ) continue;
            } else {
              if ( map->computeData[i].node != myNode ) continue;
            }
          break;

          case computeSelfExclsType:
          case computeExclsType:
            if (simParams->bondedCUDA & 16) {
              if ( ! deviceIsMineBonded ) continue;
              if ( ! deviceCUDA->device_shared_with_pe(map->computeData[i].node) ) continue;
            } else {
              if ( map->computeData[i].node != myNode ) continue;
            }
          break;

          case computeSelfCrosstermsType:
          case computeCrosstermsType:
            if (simParams->bondedCUDA & 32) {
              if ( ! deviceIsMineBonded ) continue;
              if ( ! deviceCUDA->device_shared_with_pe(map->computeData[i].node) ) continue;
            } else {
              if ( map->computeData[i].node != myNode ) continue;
            }
          break;

          case computeBondedCUDAType:
            if ( ! deviceIsMineBonded ) continue;
            if ( map->computeData[i].node != myNode ) continue;
          break;
#endif // BONDED_CUDA

          case computeNonbondedCUDA2Type:
            if ( ! deviceIsMine ) continue;
// #ifdef BONDED_CUDA
//           case computeBondedCUDAType:
// #endif
          default:
            if ( map->computeData[i].node != myNode ) continue;
        }
#else // defined(NAMD_CUDA) || defined(NAMD_HIP)
        if ( map->computeData[i].node != myNode ) continue;
#endif
        DebugM(1,"Compute " << i << '\n');
        DebugM(1,"  node = " << map->computeData[i].node << '\n');
        DebugM(1,"  type = " << map->computeData[i].type << '\n');
        DebugM(1,"  numPids = " << map->computeData[i].numPids << '\n');
//         DebugM(1,"  numPidsAllocated = " << map->computeData[i].numPidsAllocated << '\n');
        for (int j=0; j < map->computeData[i].numPids; j++)
        {
            DebugM(1,"  pid " << map->computeData[i].pids[j].pid << '\n');
            if (!((j+1) % 6))
                DebugM(1,'\n');
        }
        DebugM(1,"\n---------------------------------------");
        DebugM(1,"---------------------------------------\n");

        createCompute(i, map);

    }

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
    if (deviceIsMine) {
      getCudaComputeNonbonded()->assignPatches(this);
      getCudaComputeNonbonded()->initialize();
    }
#ifdef BONDED_CUDA
    if (simParams->bondedCUDA) {
      if (deviceIsMineBonded) {
        getComputeBondedCUDA()->initialize();
      }
    }
#endif
#endif
}

#if 0
void ComputeMgr:: sendComputeGlobalConfig(ComputeGlobalConfigMsg *msg)
{
    (CProxy_ComputeMgr(CkpvAccess(BOCclass_group).computeMgr)).recvComputeGlobalConfig(msg);
}

void ComputeMgr:: recvComputeGlobalConfig(ComputeGlobalConfigMsg *msg)
{
    if ( computeGlobalObject )
    {
        computeGlobalObject->recvConfig(msg);
    }
    else if ( ! (PatchMap::Object())->numHomePatches() ) delete msg;
    else NAMD_die("ComputeMgr::computeGlobalObject is NULL!");
}
#endif
#ifdef NODEGROUP_FORCE_REGISTER
#endif
void ComputeMgr:: sendComputeGlobalData(ComputeGlobalDataMsg *msg)
{
  NAMD_EVENT_START(1, NamdProfileEvent::GM_SEND_COMP_DATA);        
  //  CkPrintf("*** [%d] Calling sendComputeGlobalData\n", CkMyPe());
  #ifdef NODEGROUP_FORCE_REGISTER
  SimParameters *sp = Node::Object()->simParameters;
  if (sp->CUDASOAintegrate) {
    NAMD_EVENT_START(1, NamdProfileEvent::GM_NODELOCK);        
    CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
    PatchData *patchData = cpdata.ckLocalBranch();
    CmiNodeLock &nl = patchData->nodeLock;
    // atomic access to GlobalMasterServer to simulate queueing
    if (CkMyPe() != 0)
    {
      CmiLock(nl);
      //CkPrintf("*** [%d] Acquired nodelock!\n", CkMyPe());
      patchData->master_mgr->recvComputeGlobalData(msg);
      CmiUnlock(nl);
    }
    NAMD_EVENT_STOP(1, NamdProfileEvent::GM_NODELOCK);
    NAMD_EVENT_START(1, NamdProfileEvent::GM_BARRIER);        
    // Barrier to make sure 0 goes last, since invocation of the clients and
    // message coordination has to happen on PE 0 and the last PE to call
    // recvComputeGlobalData will trigger all of that on itself
    //    CmiNodeBarrier();
    //    CkPrintf("*** sendComputeGlobalData entering barrier 1 on PE %d \n", CkMyPe());
    stowSuspendULT();
    
    NAMD_EVENT_STOP(1, NamdProfileEvent::GM_BARRIER);        
    if (CkMyPe() == 0)
      {
        CmiLock(nl);
        patchData->master_mgr->recvComputeGlobalData(msg);
        CmiUnlock(nl);
      }
    else
    {
      // All PEs other than 0 wait here while the clients run and the global
      // results messages are prepared and copied into their slots (happens from
      // sendComputeGlobalResults on PE0)
      //  CmiNodeBarrier();
      //      CkPrintf("before call to stow %d\n",CkMyPe());
      //      CkPrintf("*** sendComputeGlobalData barrier 3 on PE %d \n", CkMyPe());
      stowSuspendULT();
      //      CkPrintf("*** sendComputeGlobalData out barrier 3 on PE %d \n", CkMyPe());
      //      CkPrintf("returned from call to stow %d\n",CkMyPe());
    }
    // Get the message from the slot for this PE and resume execution
    ComputeGlobalResultsMsg* resultsMsg = CkpvAccess(ComputeGlobalResultsMsg_instance);
    DebugM(3,"["<<CkMyPe()<<"] calling recvComputeGlobalResults\n");
    recvComputeGlobalResults(resultsMsg);
  } else {
  #endif
    CProxy_ComputeMgr cm(CkpvAccess(BOCclass_group).computeMgr);
    DebugM(3,"["<<CkMyPe()<<"] msg to recvComputeGlobalData\n");
    cm[0].recvComputeGlobalData(msg);
  #ifdef NODEGROUP_FORCE_REGISTER
  }
  #endif
  NAMD_EVENT_STOP(1, NamdProfileEvent::GM_SEND_COMP_DATA);
  DebugM(3,"["<<CkMyPe()<<"] done sendComputeGlobalData\n");
}

void ComputeMgr:: recvComputeGlobalData(ComputeGlobalDataMsg *msg)
{
  NAMD_EVENT_START(1, NamdProfileEvent::GM_RECV_COMP_DATA);        
    if (masterServerObject)  // make sure it has been initialized
    {
      DebugM(3, "["<<CkMyPe()<<"] recvComputeGlobalData calling recvData\n");
        masterServerObject->recvData(msg);
    }
    else NAMD_die("ComputeMgr::masterServerObject is NULL!");
    NAMD_EVENT_STOP(1, NamdProfileEvent::GM_RECV_COMP_DATA);        
}

void ComputeMgr:: sendComputeGlobalResults(ComputeGlobalResultsMsg *msg)
{
  NAMD_EVENT_START(1, NamdProfileEvent::GM_SEND_COMP_RESULTS);        
    msg->seq = ++computeGlobalResultsMsgMasterSeq;
    DebugM(3,"["<< CkMyPe()<< "] sendComputeGlobalResults seq "<<msg->seq<<"\n");

  #ifdef NODEGROUP_FORCE_REGISTER
  SimParameters *sp = Node::Object()->simParameters;
  if (sp->CUDASOAintegrate) {
    // Only PE 0 runs this code
    // Copy the message into each PE's slot (Assumes single-node with multicore build)
    for (int pe = 0; pe < CkMyNodeSize(); pe++) {
      if(CkpvAccessOther(ComputeGlobalResultsMsg_instance, pe)!=nullptr)
        {
          // make sure msg delete happens on the same PE as made the msg to
          // avoid unbounded memory pool growth for these unsent messages
          delete CkpvAccessOther(ComputeGlobalResultsMsg_instance, pe);
        }
      CkpvAccessOther(ComputeGlobalResultsMsg_instance, pe) = (ComputeGlobalResultsMsg*)CkCopyMsg((void**)&msg);
    }
    delete msg;
    // Now that copies are done, trigger the barrier to resume the other PEs
    // (most other PEs call this barrier from sendComputeGlobalData)
    //    CkPrintf("this is where we would call awaken\n",CkMyPe());
    //CmiNodeBarrier();
    //    CkPrintf("*** sendComputeGlobalResults entering barrier 2 on PE %d \n", CkMyPe());
    stowSuspendULT();
    //thisProxy.recvComputeGlobalResults(msg);
  } else {
  #endif
    DebugM(3,"["<< CkMyPe() << "] ComputeMgr::sendComputeGlobalResults invoking bcast recvComputeGlobalResults\n");    
    thisProxy.recvComputeGlobalResults(msg);
  #ifdef NODEGROUP_FORCE_REGISTER
  }
  #endif
  NAMD_EVENT_STOP(1, NamdProfileEvent::GM_SEND_COMP_RESULTS);        
}

void ComputeMgr:: enableComputeGlobalResults()
{
  NAMD_EVENT_START(1, NamdProfileEvent::GM_ENABLE_COMP_RESULTS);          
  ++computeGlobalResultsMsgSeq;
  DebugM(3,"["<<CkMyPe() <<"] enableComputeGlobalResults for "<< computeGlobalResultsMsgs.size() <<" messages seq "<< computeGlobalResultsMsgSeq <<"\n");
  for ( int i=0; i<computeGlobalResultsMsgs.size(); ++i ) {
    if ( computeGlobalResultsMsgs[i]->seq == computeGlobalResultsMsgSeq ) {
      ComputeGlobalResultsMsg *msg = computeGlobalResultsMsgs[i];
      computeGlobalResultsMsgs.del(i);
      recvComputeGlobalResults(msg);
      break;
    }
  }
  NAMD_EVENT_STOP(1, NamdProfileEvent::GM_ENABLE_COMP_RESULTS);
  DebugM(3,"["<<CkMyPe() <<"] exiting enableComputeGlobalResults for "<< computeGlobalResultsMsgs.size() <<" messages seq "<< computeGlobalResultsMsgSeq <<"\n");
}

void ComputeMgr:: recvComputeGlobalResults(ComputeGlobalResultsMsg *msg)
{
  NAMD_EVENT_START(1, NamdProfileEvent::GM_RCV_COMP_RESULTS);
  DebugM(3,"[" << CkMyPe() << "] recvComputeGlobalResults msg->seq "<< msg->seq << " computeGlobalResultsMsgSeq " << computeGlobalResultsMsgSeq << "\n");
    if ( computeGlobalObject )
    {
      if ( msg->seq == computeGlobalResultsMsgSeq ) {
        CmiEnableUrgentSend(1);

        computeGlobalObject->recvResults(msg);
        //      CkPrintf("*** past recvResults on PE %d \n", CkMyPe());
        CmiEnableUrgentSend(0);
      } else {
        //      CkPrintf("*** Adding recvComputeGlobalResults on PE %d \n", CkMyPe());
        computeGlobalResultsMsgs.add(msg);
      }
    }
    else if ( ! (PatchMap::Object())->numHomePatches() )
      {
        //      CkPrintf("*** ignoring recvComputeGlobalResults on PE %d due to no home patch\n", CkMyPe());
        delete msg;
      }
    else NAMD_die("ComputeMgr::computeGlobalObject is NULL!");
    NAMD_EVENT_STOP(1, NamdProfileEvent::GM_RCV_COMP_RESULTS);
    //    CkPrintf("*** exiting recvComputeGlobalResults on PE %d \n", CkMyPe());
}

/*
 * Begin Ewald messages
 */
void ComputeMgr:: sendComputeEwaldData(ComputeEwaldMsg *msg)
{
    if (computeEwaldObject)
    {
        int node = computeEwaldObject->getMasterNode();
        CProxy_ComputeMgr cm(CkpvAccess(BOCclass_group).computeMgr);
        cm[node].recvComputeEwaldData(msg);
    }
    else if (!PatchMap::Object()->numHomePatches())
    {
      // CkPrintf("skipping message on Pe(%d)\n", CkMyPe());
        delete msg;
    }
    else NAMD_die("ComputeMgr::computeEwaldObject is NULL!");
}

void ComputeMgr:: recvComputeEwaldData(ComputeEwaldMsg *msg)
{
    if (computeEwaldObject)
        computeEwaldObject->recvData(msg);
    else NAMD_die("ComputeMgr::computeEwaldObject in recvData is NULL!");
}

void ComputeMgr:: sendComputeEwaldResults(ComputeEwaldMsg *msg)
{
    (CProxy_ComputeMgr(CkpvAccess(BOCclass_group).computeMgr)).recvComputeEwaldResults(msg);
}

void ComputeMgr::recvComputeEwaldResults(ComputeEwaldMsg *msg)
{
    if (computeEwaldObject) {
        CmiEnableUrgentSend(1);
        computeEwaldObject->recvResults(msg);
        CmiEnableUrgentSend(0);
    }
    else if ( ! (PatchMap::Object())->numHomePatches() ) delete msg;
    else NAMD_die("ComputeMgr::computeEwaldObject in recvResults is NULL!");
}

void ComputeMgr:: sendComputeDPMEData(ComputeDPMEDataMsg *msg)
{
    if ( computeDPMEObject )
    {
#ifdef DPME
        int node = computeDPMEObject->getMasterNode();
        CProxy_ComputeMgr cm(CkpvAccess(BOCclass_group).computeMgr);
        cm.recvComputeDPMEData(msg,node);
#endif
    }
    else if ( ! (PatchMap::Object())->numHomePatches() ) delete msg;
    else NAMD_die("ComputeMgr::computeDPMEObject is NULL!");
}

void ComputeMgr:: recvComputeDPMEData(ComputeDPMEDataMsg *msg)
{
    if ( computeDPMEObject )
    {
#ifdef DPME
        computeDPMEObject->recvData(msg);
#endif
    }
    else if ( ! (PatchMap::Object())->numHomePatches() ) delete msg;
    else NAMD_die("ComputeMgr::computeDPMEObject is NULL!");
}

void ComputeMgr:: sendComputeDPMEResults(ComputeDPMEResultsMsg *msg, int node)
{
    CProxy_ComputeMgr cm(CkpvAccess(BOCclass_group).computeMgr);
    cm[node].recvComputeDPMEResults(msg);
}

void ComputeMgr:: recvComputeDPMEResults(ComputeDPMEResultsMsg *msg)
{
    if ( computeDPMEObject )
    {
#ifdef DPME
        computeDPMEObject->recvResults(msg);
#endif
    }
    else if ( ! (PatchMap::Object())->numHomePatches() ) delete msg;
    else NAMD_die("ComputeMgr::computeDPMEObject is NULL!");
}

/*
 * Molecule now has only one instance per process, so this must only
 * be done once per process.  

 * TODO: A molecule manager nodegroup would be the natural place
 * for entry methods that alter the molecule like this.
 */
void ComputeMgr::recvComputeConsForceMsg(ComputeConsForceMsg *msg)
{
    Molecule *m = Node::Object()->molecule;
    if(CkMyRank()==0){ // there is only one molecule per process
      delete [] m->consForceIndexes;
      delete [] m->consForce;
      int n = msg->aid.size();
      if (n > 0)
        {
          m->consForceIndexes = new int32[m->numAtoms];
          m->consForce = new Vector[n];
          int i;
          for (i=0; i<m->numAtoms; i++) m->consForceIndexes[i] = -1;
          for (i=0; i<msg->aid.size(); i++)
            {
              m->consForceIndexes[msg->aid[i]] = i;
              m->consForce[i] = msg->f[i];
            }
        }
      else
        {
          m->consForceIndexes = NULL;
          m->consForce = NULL;
        }
    }
    delete msg;
#ifdef NODEGROUP_FORCE_REGISTER
    if(CkMyPe()==0)
      {
        CProxy_PatchData cpdata(CkpvAccess(BOCclass_group).patchData);
        cpdata.setDeviceKernelUpdateCounter();
      }
#endif
}

void ComputeMgr::recvCudaGlobalMasterCreateMsg(std::vector<std::string> args) {
#ifdef NAMD_CUDA
  Node *node = Node::Object();
  SimParameters *simParams = node->simParameters;
  if (simParams->CUDASOAintegrate && simParams->useCudaGlobal) {
#ifdef NODEGROUP_FORCE_REGISTER
    if (deviceCUDA->getMasterPe() == CkMyPe()) {
      if (deviceCUDA->getIsGlobalDevice()) {
        DebugM(3, "Call recvCudaGlobalMasterCreateMsg on master PE " << CkMyPe() << ".\n");
        ComputeCUDAMgr* cudaMgr = ComputeCUDAMgr::getComputeCUDAMgr();
        cudaMgr->createCudaGlobalMaster();
        std::shared_ptr<CudaGlobalMasterClient> client = nullptr;
        const std::string library_name = args[0];
        // Find to see if library_name has been loaded
        std::shared_ptr<dlloader::DLLoader<CudaGlobalMasterClient>> loader = nullptr;
        for (auto it = CudaGlobalMasterClientDlloaders.begin();
             it != CudaGlobalMasterClientDlloaders.end(); ++it) {
          if ((*it)->LibName() == library_name) {
            loader = (*it);
            break;
          }
        }
        // Create a new loader if not found
        if (loader == nullptr) {
          loader = std::shared_ptr<dlloader::DLLoader<CudaGlobalMasterClient>>(new dlloader::DLLoader<CudaGlobalMasterClient>(library_name));
        }
        try {
          iout << iINFO << "Loading library " << library_name
               << " on PE: " << CkMyPe() << "\n" << endi;
          loader->DLOpenLib();
          client = loader->DLGetInstance();
        } catch (std::exception& e) {
          iout << iERROR << "Cannot load the shared library " << library_name << "\n" << endi;
          NAMD_die(e.what());
        }
        // Try to initialize the client
        try {
          client->initialize(args,
            deviceCUDA->getGlobalDevice(),
            cudaMgr->getCudaGlobalMaster()->getStream());
          client->subscribe(cudaMgr->getCudaGlobalMaster());
          iout << iINFO << "CudaGlobalMaster client \"" << client->name()
                << "\"" << " initialized\n" << endi;
        } catch (std::exception& e) {
          iout << iERROR << "Cannot initialize the CudaGlobalMaster client from "
                << library_name << "\n" << endi;
          NAMD_die(e.what());
        }
        CudaGlobalMasterClientDlloaders.push_back(loader);
      } else {
        DebugM(3, "Skip recvCudaGlobalMasterCreateMsg on master PE " <<
                  CkMyPe() << " that is not scheduled for GPU-resident global master.\n");
      }
    } else {
      DebugM(3, "Skip recvCudaGlobalMasterCreateMsg on non-master PE " << CkMyPe() << ".\n");
    }
#endif // NODEGROUP_FORCE_REGISTER
  } else {
    if (!(simParams->CUDASOAintegrate)) {
      NAMD_die("GPU-resident mode is not enabled.\n");
    }
    if (!(simParams->useCudaGlobal)) {
      NAMD_die("GPU-resident external forces are not enabled.\n");
    }
  }
  // CmiNodeBarrier();
#endif
}

void ComputeMgr::recvCudaGlobalMasterRemoveMsg(std::vector<std::string> args) {
#ifdef NAMD_CUDA
  Node *node = Node::Object();
  SimParameters *simParams = node->simParameters;
  const std::string client_name_to_remove = args[0];
  if (simParams->CUDASOAintegrate && simParams->useCudaGlobal) {
#ifdef NODEGROUP_FORCE_REGISTER
    if (deviceCUDA->getMasterPe() == CkMyPe()) {
      if (deviceCUDA->getIsGlobalDevice()) {
        ComputeCUDAMgr* cudaMgr = ComputeCUDAMgr::getComputeCUDAMgr();
        std::shared_ptr<CudaGlobalMasterServer> gm = cudaMgr->getCudaGlobalMaster();
        if (gm) {
          std::shared_ptr<CudaGlobalMasterClient> c = nullptr;
          const std::vector<std::shared_ptr<CudaGlobalMasterClient>>& clients = gm->getClients();
          for (size_t i = 0; i < clients.size(); ++i) {
            if (client_name_to_remove == clients[i]->name()) {
              c = clients[i];
              break;
            }
          }
          if (c) {
            gm->removeClient(c);
            iout << iINFO << "CudaGlobalMasterClient \""
                 << client_name_to_remove << "\" removed\n" << endi;
          } else {
            const std::string error = "CudaGlobalMasterClient \""
              + client_name_to_remove + "\" not found";
            NAMD_die(error.c_str());
          }
        }
      }
    }
#endif // NODEGROUP_FORCE_REGISTER
  } else {
    if (!(simParams->CUDASOAintegrate)) {
      NAMD_die("GPU-resident mode is not enabled.\n");
    }
    if (!(simParams->useCudaGlobal)) {
      NAMD_die("GPU-resident external forces are not enabled.\n");
    }
  }
#endif
}

void ComputeMgr::recvCudaGlobalMasterUpdateMsg(std::vector<std::string> args) {
  // XXX Should this also be for NAMD_HIP ?
#ifdef NAMD_CUDA
  std::vector<std::string> result_args;
  Node *node = Node::Object();
  SimParameters *simParams = node->simParameters;
  const std::string client_name_to_update = args[0];
  if (simParams->CUDASOAintegrate && simParams->useCudaGlobal) {
#ifdef NODEGROUP_FORCE_REGISTER
    if (deviceCUDA->getMasterPe() == CkMyPe()) {
      if (deviceCUDA->getIsGlobalDevice()) {
        ComputeCUDAMgr* cudaMgr = ComputeCUDAMgr::getComputeCUDAMgr();
        std::shared_ptr<CudaGlobalMasterServer> gm = cudaMgr->getCudaGlobalMaster();
        if (gm) {
          std::shared_ptr<CudaGlobalMasterClient> c = nullptr;
          const std::vector<std::shared_ptr<CudaGlobalMasterClient>>& clients = gm->getClients();
          for (size_t i = 0; i < clients.size(); ++i) {
            if (client_name_to_update == clients[i]->name()) {
              c = clients[i];
              break;
            }
          }
          if (c) {
            result_args.push_back(client_name_to_update);
            result_args.push_back(c->updateFromTCLCommand(args));
            iout << iINFO << "CudaGlobalMasterClient \""
                 << client_name_to_update << "\" updated\n" << endi;
          } else {
            const std::string error = "CudaGlobalMasterClient \""
              + client_name_to_update + "\" not found";
            NAMD_die(error.c_str());
          }
        }
      }
    }
#endif // NODEGROUP_FORCE_REGISTER
  } else {
    if (!(simParams->CUDASOAintegrate)) {
      NAMD_die("GPU-resident mode is not enabled.\n");
    }
    if (!(simParams->useCudaGlobal)) {
      NAMD_die("GPU-resident external forces are not enabled.\n");
    }
  }
  CProxy_ComputeMgr cm(CkpvAccess(BOCclass_group).computeMgr);
  cm[0].recvCudaGlobalMasterUpdateResultMsg(result_args);
#endif
}

void ComputeMgr::recvCudaGlobalMasterUpdateResultMsg(std::vector<std::string> args) {
  if (CkMyPe() == 0) {
    if (!args.empty()) {
      CudaGlobalMasterClientUpdateResults[args[0]] = args[1];
    }
  } else {
    const std::string error =
      "recvCudaGlobalMasterUpdateResultMsg is called on " +
      std::to_string(CkMyPe()) + " but expected on PE 0!\n";
    NAMD_bug(error.c_str());
  }
}

std::string ComputeMgr::getCudaGlobalMasterUpdateResult(const std::string& client_name) const {
  return CudaGlobalMasterClientUpdateResults.at(client_name);
}

void ComputeMgr::sendYieldDevice(int pe) {
    CProxy_ComputeMgr cm(CkpvAccess(BOCclass_group).computeMgr);
    cm[pe].recvYieldDevice(CkMyPe());
}

void ComputeMgr::recvYieldDevice(int pe) {
  // XXX MIC support was only code using YieldDevice functionality
  // computeNonbondedMICObject->recvYieldDevice(pe);
}

#if defined(NAMD_CUDA) || defined(NAMD_HIP)
class CudaComputeNonbondedMsg : public CMessage_CudaComputeNonbondedMsg {
public:
  CudaComputeNonbonded* c;
  int i;
};

void ComputeMgr::sendAssignPatchesOnPe(std::vector<int>& pes, CudaComputeNonbonded* c) {
  for (int i=0;i < pes.size();i++) {
    CudaComputeNonbondedMsg *msg = new CudaComputeNonbondedMsg;
    msg->c = c;
    thisProxy[pes[i]].recvAssignPatchesOnPe(msg);
  }
}

void ComputeMgr::recvAssignPatchesOnPe(CudaComputeNonbondedMsg *msg) {
  msg->c->assignPatchesOnPe();
  delete msg;
}

void ComputeMgr::sendSkipPatchesOnPe(std::vector<int>& pes, CudaComputeNonbonded* c) {
  for (int i=0;i < pes.size();i++) {
    CudaComputeNonbondedMsg *msg = new CudaComputeNonbondedMsg;
    msg->c = c;
    thisProxy[pes[i]].recvSkipPatchesOnPe(msg);
  }
}

void ComputeMgr::recvSkipPatchesOnPe(CudaComputeNonbondedMsg *msg) {
  msg->c->skipPatchesOnPe();
  delete msg;
}

void ComputeMgr::sendFinishPatchesOnPe(std::vector<int>& pes, CudaComputeNonbonded* c) {
  for (int i=0;i < pes.size();i++) {
    CudaComputeNonbondedMsg *msg = new (PRIORITY_SIZE) CudaComputeNonbondedMsg;
    SET_PRIORITY(msg, c->sequence(), COMPUTE_PROXY_PRIORITY);
    msg->c = c;
    thisProxy[pes[i]].recvFinishPatchesOnPe(msg);
  }
}

void ComputeMgr::recvFinishPatchesOnPe(CudaComputeNonbondedMsg *msg) {
  msg->c->finishPatchesOnPe();
  delete msg;
}

void ComputeMgr::sendFinishPatchOnPe(int pe, CudaComputeNonbonded* c, int i, PatchID patchID) {
  CudaComputeNonbondedMsg *msg = new (PRIORITY_SIZE) CudaComputeNonbondedMsg;
  SET_PRIORITY(msg, c->sequence(), COMPUTE_PROXY_PRIORITY + PATCH_PRIORITY(patchID));
  msg->c = c;
  msg->i = i;
  thisProxy[pe].recvFinishPatchOnPe(msg);
}

void ComputeMgr::recvFinishPatchOnPe(CudaComputeNonbondedMsg *msg) {
  msg->c->finishPatchOnPe(msg->i);
  delete msg;
}

void ComputeMgr::sendOpenBoxesOnPe(std::vector<int>& pes, CudaComputeNonbonded* c) {
  for (int i=0;i < pes.size();i++) {
    CudaComputeNonbondedMsg *msg = new (PRIORITY_SIZE) CudaComputeNonbondedMsg;
    SET_PRIORITY(msg, c->sequence(), PROXY_DATA_PRIORITY+1); // after bonded
    msg->c = c;
    thisProxy[pes[i]].recvOpenBoxesOnPe(msg);
  }
}

void ComputeMgr::recvOpenBoxesOnPe(CudaComputeNonbondedMsg *msg) {
  msg->c->openBoxesOnPe();
  delete msg;
}

void ComputeMgr::sendFinishReductions(int pe, CudaComputeNonbonded* c) {
  CudaComputeNonbondedMsg *msg = new CudaComputeNonbondedMsg;
  msg->c = c;
  thisProxy[pe].recvFinishReductions(msg);
}

void ComputeMgr::recvFinishReductions(CudaComputeNonbondedMsg *msg) {
  msg->c->finishReductions();
  delete msg;
}

void ComputeMgr::sendMessageEnqueueWork(int pe, CudaComputeNonbonded* c) {
  CudaComputeNonbondedMsg *msg = new CudaComputeNonbondedMsg;
  msg->c = c;
  thisProxy[pe].recvMessageEnqueueWork(msg);
}

void ComputeMgr::recvMessageEnqueueWork(CudaComputeNonbondedMsg *msg) {
  msg->c->messageEnqueueWork();
  delete msg;
}

void ComputeMgr::sendLaunchWork(int pe, CudaComputeNonbonded* c) {
  CudaComputeNonbondedMsg *msg = new CudaComputeNonbondedMsg;
  msg->c = c;
  thisProxy[pe].recvLaunchWork(msg);
}

void ComputeMgr::recvLaunchWork(CudaComputeNonbondedMsg *msg) {
  msg->c->launchWork();
  delete msg;
}

void ComputeMgr::sendUnregisterBoxesOnPe(std::vector<int>& pes, CudaComputeNonbonded* c) {
  for (int i=0;i < pes.size();i++) {
    CudaComputeNonbondedMsg *msg = new CudaComputeNonbondedMsg;
    msg->c = c;
    thisProxy[pes[i]].recvUnregisterBoxesOnPe(msg);
  }
}

void ComputeMgr::recvUnregisterBoxesOnPe(CudaComputeNonbondedMsg *msg) {
  msg->c->unregisterBoxesOnPe();
  delete msg;
}

#ifdef BONDED_CUDA

class ComputeBondedCUDAMsg : public CMessage_ComputeBondedCUDAMsg {
public:
  ComputeBondedCUDA* c;
  int i;
};

void ComputeMgr::sendAssignPatchesOnPe(std::vector<int>& pes, ComputeBondedCUDA* c) {
  for (int i=0;i < pes.size();i++) {
    ComputeBondedCUDAMsg *msg = new ComputeBondedCUDAMsg;
    msg->c = c;
    thisProxy[pes[i]].recvAssignPatchesOnPe(msg);
  }
}

void ComputeMgr::recvAssignPatchesOnPe(ComputeBondedCUDAMsg *msg) {
  msg->c->assignPatchesOnPe();
  delete msg;
}

void ComputeMgr::sendMessageEnqueueWork(int pe, ComputeBondedCUDA* c) {
  ComputeBondedCUDAMsg *msg = new ComputeBondedCUDAMsg;
  msg->c = c;
  thisProxy[pe].recvMessageEnqueueWork(msg);
}

void ComputeMgr::recvMessageEnqueueWork(ComputeBondedCUDAMsg *msg) {
  msg->c->messageEnqueueWork();
  delete msg;
}

void ComputeMgr::sendOpenBoxesOnPe(std::vector<int>& pes, ComputeBondedCUDA* c) {
  for (int i=0;i < pes.size();i++) {
    ComputeBondedCUDAMsg *msg = new (PRIORITY_SIZE) ComputeBondedCUDAMsg;
    SET_PRIORITY(msg, c->sequence(), PROXY_DATA_PRIORITY);
    msg->c = c;
    thisProxy[pes[i]].recvOpenBoxesOnPe(msg);
  }
}

void ComputeMgr::recvOpenBoxesOnPe(ComputeBondedCUDAMsg *msg) {
  msg->c->openBoxesOnPe();
  delete msg;
}

void ComputeMgr::sendLoadTuplesOnPe(std::vector<int>& pes, ComputeBondedCUDA* c) {
  for (int i=0;i < pes.size();i++) {
    ComputeBondedCUDAMsg *msg = new ComputeBondedCUDAMsg;
    msg->c = c;
    thisProxy[pes[i]].recvLoadTuplesOnPe(msg);
  }
}

void ComputeMgr::recvLoadTuplesOnPe(ComputeBondedCUDAMsg *msg) {
  msg->c->loadTuplesOnPe();
  delete msg;
}

void ComputeMgr::sendLaunchWork(int pe, ComputeBondedCUDA* c) {
  ComputeBondedCUDAMsg *msg = new ComputeBondedCUDAMsg;
  msg->c = c;
  thisProxy[pe].recvLaunchWork(msg);
}

void ComputeMgr::recvLaunchWork(ComputeBondedCUDAMsg *msg) {
  msg->c->launchWork();
  delete msg;
}

void ComputeMgr::sendFinishPatchesOnPe(std::vector<int>& pes, ComputeBondedCUDA* c) {
  for (int i=0;i < pes.size();i++) {
    ComputeBondedCUDAMsg *msg = new (PRIORITY_SIZE) ComputeBondedCUDAMsg;
    SET_PRIORITY(msg, c->sequence(), COMPUTE_PROXY_PRIORITY);
    msg->c = c;
    thisProxy[pes[i]].recvFinishPatchesOnPe(msg);
  }
}

void ComputeMgr::recvFinishPatchesOnPe(ComputeBondedCUDAMsg *msg) {
  msg->c->finishPatchesOnPe();
  delete msg;
}

void ComputeMgr::sendFinishReductions(int pe, ComputeBondedCUDA* c) {
  ComputeBondedCUDAMsg *msg = new ComputeBondedCUDAMsg;
  msg->c = c;
  thisProxy[pe].recvFinishReductions(msg);
}

void ComputeMgr::recvFinishReductions(ComputeBondedCUDAMsg *msg) {
  msg->c->finishReductions();
  delete msg;
}

void ComputeMgr::sendUnregisterBoxesOnPe(std::vector<int>& pes, ComputeBondedCUDA* c) {
  for (int i=0;i < pes.size();i++) {
    ComputeBondedCUDAMsg *msg = new ComputeBondedCUDAMsg;
    msg->c = c;
    thisProxy[pes[i]].recvUnregisterBoxesOnPe(msg);
  }
}

void ComputeMgr::recvUnregisterBoxesOnPe(ComputeBondedCUDAMsg *msg) {
  msg->c->unregisterBoxesOnPe();
  delete msg;
}

#endif // BONDED_CUDA

#endif // NAMD_CUDA

#include "ComputeMgr.def.h"