ComputeGroupRestraintsCUDA.C

Go to the documentation of this file.

#include "ComputeGroupRestraintsCUDA.h"
#include "ComputeGroupRes1GroupCUDAKernel.h"
#include "ComputeGroupRes2GroupCUDAKernel.h"
#include "ComputeCUDAMgr.h"
#include "SimParameters.h"
#include "Node.h"
#include "Molecule.h"
#include "InfoStream.h"
#define MIN_DEBUG_LEVEL 3
//#define DEBUGM
#include "Debug.h"

#ifdef NODEGROUP_FORCE_REGISTER
GroupRestraintsCUDA::GroupRestraintsCUDA(const GroupRestraintParam *param, bool _mGpuOn, int _numDevices, int _deviceIndex, int grp) {
    resParam = param;
    groupName = param->GetGroupName();
    mGpuOn = _mGpuOn;
    numDevices = _numDevices;
    deviceIndex = _deviceIndex;
    groupRestID = grp;
    splitOverDevicesGrp1=false;
    splitOverDevicesGrp2=false;
    restraintExp = param->GetExponent();
    restraintK = param->GetForce();
    inv_group1_mass = 0.0;
    inv_group2_mass = 0.0;
    const std::vector<int> &group1Index = param->GetGroup1AtomIndex();
    const std::vector<int> &group2Index = param->GetGroup2AtomIndex();
    numRestrainedGroup1 = group1Index.size();
    numRestrainedGroup2 = group2Index.size();
    totalNumRestrained = numRestrainedGroup1 + numRestrainedGroup2; 
    // If we defined a list of atoms for group 1, then we have to
    // calculate the COM for group 1 at every steps
    calcGroup1COM = (numRestrainedGroup1 ? true : false);
    // groupAtomsSOAIndex stores SOAindex of group 1, followed by 
    // SOAindex of group 2 
    groupAtomsSOAIndex.resize(totalNumRestrained);
    useDistMagnitude = param->GetUseDistMagnitude();
    Vector center = param->GetResCenter();
    Vector dir = param->GetResDirection();
    resDirection = make_double3(dir.x, dir.y, dir.z);
    resCenterVec = make_double3(center.x, center.y, center.z);

    // Allocate host and device memory
    allocate_host<double>(&h_resEnergy, 1);
    allocate_host<double3>(&h_diffCOM, 1);
    allocate_host<double3>(&h_group1COM, 1);
    allocate_host<double3>(&h_group2COM, 1);

    allocate_host<double3>(&h_resForce, 1);
    allocate_device<double3>(&d_group1COM, 1);
    allocate_device<double3>(&d_group2COM, 1);
    allocate_device<int>(&d_groupAtomsSOAIndex, totalNumRestrained);
    allocate_device<unsigned int>(&d_tbcatomic, 1);
    // Set the counter to zero
    cudaCheck(cudaMemset(d_tbcatomic, 0, sizeof(unsigned int)));
  if(mGpuOn)
    {// each needs its own peer array
      allocate_device<double3>(&d_peer1COM, sizeof(double3));
      allocate_device<double3>(&d_peer2COM, sizeof(double3));
#ifdef DEBUGM
      allocate_host<double3>(&h_peer1COM, sizeof(double3));
      allocate_host<double3>(&h_peer2COM, sizeof(double3)*numDevices);
#endif
    }

    // Check the atom index and calculate he inverse mass
    Molecule *mol = Node::Object()->molecule;
    int totalAtoms = mol->numAtoms;
    double total_mass = 0.0;

    for (int i = 0; i < numRestrainedGroup2; ++i) {
        int index = group2Index[i]; 
        if (index > -1 && index < totalAtoms) {
            total_mass += mol->atommass(index);    
        } else {
            char err_msg[512];
            sprintf(err_msg, "Group restraints: Bad atom index for %s!"
                " Atom indices must be within [%d, %d].\n", groupName, 0, totalAtoms - 1);
            NAMD_die(err_msg);
        }
    }
    inv_group2_mass = 1.0 / total_mass; 

    // Do we need to calculate COM of group 1, or we have a
    // reference position for it?
    if (calcGroup1COM) {
        total_mass = 0.0;
        for (int i = 0; i < numRestrainedGroup1; ++i) {
            int index = group1Index[i]; 
            if (index > -1 && index < totalAtoms) {
                total_mass += mol->atommass(index);    
            } else {
                char err_msg[512];
                sprintf(err_msg, "Group restraints: Bad atom index for %s!"
                    " Atom indices must be within [%d, %d].\n", groupName, 0, totalAtoms - 1);
                NAMD_die(err_msg);
            }
        }
        inv_group1_mass = 1.0 / total_mass; 
    } else {
        // We defined the reference point for COM of group 1, so no need
        // to calculate it, just copy it.
        // Set the h_group1COM to reference COM position of group 1
        Vector ref = param->GetGroupRes1Position();
        h_group1COM->x = ref.x;
        h_group1COM->y = ref.y;
        h_group1COM->z = ref.z;
    }
}

GroupRestraintsCUDA::~GroupRestraintsCUDA() {
    deallocate_host<double>(&h_resEnergy);
    deallocate_host<double3>(&h_resForce);
    deallocate_host<double3>(&h_diffCOM);
    deallocate_host<double3>(&h_group1COM);
    deallocate_host<double3>(&h_group2COM);
    deallocate_device<double3>(&d_group1COM); 
    deallocate_device<double3>(&d_group2COM); 
    deallocate_device<int>(&d_groupAtomsSOAIndex);
    deallocate_device<unsigned int>(&d_tbcatomic);
}

void GroupRestraintsCUDA::updateAtoms(
      std::vector<AtomMap*> &atomMapsList,
      std::vector<CudaLocalRecord> &localRecords,
      const int *h_globalToLocalID) {

    // If we need to calculate COM of group 1, we have to store
    // index to SOA data structures for group 1
    numRestrainedGroup1Local = numRestrainedGroup1;
    numRestrainedGroup2Local = numRestrainedGroup2;
    if(mGpuOn)
      groupAtomsSOAIndex.clear();
    if(calcGroup1COM) {
        const std::vector<int> &group1Index = resParam->GetGroup1AtomIndex();
        DebugM(3, "[" << CkMyPe() << "]" <<" updateAtoms grp1 "<< group1Index.size() << "\n" << endi);
        if(numRestrainedGroup1 != group1Index.size()) {
            char err_msg[512];
            sprintf("Number of atoms in group 1 restraint for '%s' is changed!", groupName);
            NAMD_bug(err_msg);
        }
        // Map the global index to local position in SOA data structure 
        for(int i = 0 ; i < numRestrainedGroup1; ++i){
            int gid = group1Index[i];
            LocalID lid;
            DebugM(2, "[" << CkMyPe() << "]" <<" updateAtoms grp1 looking for "<<  gid << "\n" << endi);
            // Search for a valid localID in all atoms
            for(int j = 0 ; j < atomMapsList.size(); ++j){
                lid = atomMapsList[j]->localID(gid);
                if( lid.pid != -1) {
                    break;
                } 
            }
            // Fields of lid need to be != -1, bc the atom needs to be somewhere
            //  otherwise we have a bug
            if(lid.pid == -1){
              if(!mGpuOn)
                NAMD_bug("LocalAtomID not found in patchMap");
              else
                splitOverDevicesGrp1=true;
            }
            else
              {
                // Converts global patch ID to its local position in our SOA data structures
                int soaPid = h_globalToLocalID[lid.pid];
                if(soaPid>=0)
                  {
                    int soaIndex = localRecords[soaPid].bufferOffset + lid.index;
                    DebugM(2, "[" << CkMyPe() << "]" <<" updateAtoms grp1 found gid "<< gid << " as soa "<<  soaIndex << "\n" << endi);
                    if(mGpuOn)
                      groupAtomsSOAIndex.push_back(soaIndex);
                    else
                      groupAtomsSOAIndex[i] = soaIndex;
                  }
              }
        }
        DebugM(3, "[" << CkMyPe() << "]" <<" updateAtoms grp1 soa "<<  groupAtomsSOAIndex.size() << "\n" << endi);
        // Sort vector for better coalesce memory access. Just sort only group 1
        if(mGpuOn)
          {
            numRestrainedGroup1Local = groupAtomsSOAIndex.size();
          }
        if(numRestrainedGroup1Local>0)
          std::sort(groupAtomsSOAIndex.begin(), groupAtomsSOAIndex.begin()+numRestrainedGroup1Local);
    }

    // We always calculate the COM of group 2, so we store
    // SOAIndex of group 2, after group 1 index
    const std::vector<int> &group2Index = resParam->GetGroup2AtomIndex();
    if(numRestrainedGroup2 != group2Index.size()) {
        char err_msg[512];
        sprintf("Number of atoms in group 2 restraint for '%s' is changed!", groupName);
        NAMD_bug(err_msg);
    }
    DebugM(4, "[" << CkMyPe() << "]" <<" updateAtoms grp2 "<< group2Index.size() << "\n" << endi);
    // Map the global index to local position in SOA data structure 
    for(int i = 0 ; i < numRestrainedGroup2; ++i){
        int gid = group2Index[i];
        LocalID lid;
        DebugM(2, "[" << CkMyPe() << "]" <<" updateAtoms grp2 looking for "<<  gid << "\n" << endi);
        // Search for a valid localID in all atoms
        for(int j = 0 ; j < atomMapsList.size(); ++j){
            lid = atomMapsList[j]->localID(gid);
            if( lid.pid != -1) {
                break;
            } 
        }
        // Fields of lid need to be != -1, bc the atom needs to be somewhere
        //  otherwise we have a bug
        if(lid.pid == -1){
          if(!mGpuOn)
            NAMD_bug("LocalAtomID not found in patchMap");
          else
            splitOverDevicesGrp2=true;
        }
        else
          {

        // Converts global patch ID to its local position in our SOA data structures
            int soaPid = h_globalToLocalID[lid.pid];
            if(soaPid>=0)
              {
                int soaIndex = localRecords[soaPid].bufferOffset + lid.index;
                DebugM(2, "[" << CkMyPe() << "]" <<" updateAtoms grp2 found "<<  gid << " as SOA "<< soaIndex << "\n" << endi);
                // store the index for group 2, after group 1 index
                if(mGpuOn)
                  groupAtomsSOAIndex.push_back(soaIndex);
                else
                  groupAtomsSOAIndex[i + numRestrainedGroup1] = soaIndex;
                DebugM(2, "[" << CkMyPe() << "]" <<" updateAtoms grp2 SOA soaidx "<< soaIndex<< " size now "<< groupAtomsSOAIndex.size() <<"\n" << endi);
              }
          }
    }
    numRestrainedGroup2Local = groupAtomsSOAIndex.size() - numRestrainedGroup1Local;
    DebugM(3, "[" << CkMyPe() << "]" <<" updateAtoms grp2 SOA now "<< groupAtomsSOAIndex.size() << " numRestrainedGroup1Local " << numRestrainedGroup1Local << " numRestrainedGroup2Local " << numRestrainedGroup2Local <<"\n" << endi);
    // if there are local atoms, send indices to the device
    if(groupAtomsSOAIndex.size()>0){
      // Sort vector for better coalesce memory access. Sort only for group 2
      if(numRestrainedGroup2Local>0)
        std::sort(groupAtomsSOAIndex.begin() + numRestrainedGroup1Local, groupAtomsSOAIndex.end());
      // Update the SOA index in device
      copy_HtoD<int>(groupAtomsSOAIndex.data(), d_groupAtomsSOAIndex, groupAtomsSOAIndex.size());
    }
}
void GroupRestraintsCUDA::initPeerCOM(double3** d_peerCOM1G,
                                 double3** d_peerCOM2G, cudaStream_t stream){
  DebugM(3, "[" << CkMyPe() << "]" << " initPeerCOM\n" << endi);
  initPeerCOMmgpuG(numDevices, deviceIndex, d_peerCOM1G, d_peer1COM, stream);
  initPeerCOMmgpuG(numDevices, deviceIndex, d_peerCOM2G, d_peer2COM, stream);
}

void GroupRestraintsCUDA::doCOM_mgpu(
        const Lattice &lat,
        const char3*  d_transform,
        const float*  d_mass,
        const double* d_pos_x,
        const double* d_pos_y,
        const double* d_pos_z,
        const int     grp,
        cudaStream_t  stream) {
  groupRestID = grp;
  ComputeCUDAMgr* cudaMgr    = ComputeCUDAMgr::getComputeCUDAMgr();
  DebugM(3, "[" << CkMyPe() << "]" << " doCOM_mgpu g1count "<< numRestrainedGroup1Local<< " g2count "<<numRestrainedGroup2Local  << endi);
  if(numRestrainedGroup1Local>0 && calcGroup1COM)
    {
      // if compute COM1 and have some of the group's atom, compute distCOM grp1
        computeCOMMgpu(numRestrainedGroup1Local,lat,
                    d_mass, d_pos_x, d_pos_y, d_pos_z,
                       d_transform, this->d_groupAtomsSOAIndex,
                       this->d_peer1COM,
                       cudaMgr->curGrp1COM[grp],
                       this->d_tbcatomic,
                       numDevices, deviceIndex,
                    stream);
#ifdef DEBUGM
        cudaCheck(cudaStreamSynchronize(stream));
        copy_DtoH_sync<double3>(d_peer1COM, h_peer1COM, 1);
        DebugM(3, "gid " << groupRestID <<" deviceIndex "<< deviceIndex << " g1 COM " <<h_peer1COM->x*inv_group1_mass<<", "
               <<h_peer1COM->y*inv_group1_mass<<", "
               <<h_peer2COM->z*inv_group1_mass<<"\n"<<endi);
#endif
    }
  if(numRestrainedGroup2Local>0)
    {
      // if we have some of the group's atoms, compute distCOM grp2
        computeCOMMgpu(numRestrainedGroup2Local,lat,
                    d_mass, d_pos_x, d_pos_y, d_pos_z,
                       d_transform,
                       this->d_groupAtomsSOAIndex+numRestrainedGroup1Local,
                       this->d_peer2COM,
                       cudaMgr->curGrp2COM[grp],
                       this->d_tbcatomic,  numDevices, deviceIndex,
                    stream);
#ifdef DEBUGM
        cudaCheck(cudaStreamSynchronize(stream));
        copy_DtoH_sync<double3>(d_peer2COM, h_peer2COM, 1);
        DebugM(3, " gid "<< groupRestID << " deviceIndex "<< deviceIndex << " g2 COM " <<h_peer2COM[0].x*inv_group2_mass<<", "
               <<h_peer2COM[0].y*inv_group2_mass<<", "
               <<h_peer2COM[0].z*inv_group2_mass<<"\n"<<endi);
#endif
    }
}

void GroupRestraintsCUDA::doForce(
    const int timeStep,
    const int doEnergy,
    const int doVirial,
    const int doOutput,
    const int grp,
    const Lattice &lat, 
    const char3*  d_transform, 
    const float*  d_mass, 
    const double* d_pos_x,
    const double* d_pos_y,
    const double* d_pos_z,
    double*       d_f_normal_x,
    double*       d_f_normal_y,
    double*       d_f_normal_z,
    cudaTensor*   d_virial,
    double*       h_extEnergy, 
    double3*      h_extForce, 
    cudaTensor*   h_extVirial,
    double3**      d_peer1COM,
    double3**      d_peer2COM,
    cudaStream_t  stream) {
  DebugM(3, "[" << CkMyPe() << "]" << " doForce g1 "<< numRestrainedGroup1Local <<" g2 "<< numRestrainedGroup2Local <<"\n" << endi);
    // distributed COM for mGpu is not effectively the same as the 1Group
    // case, because the force computation uses both group1 and group2
    ComputeCUDAMgr* cudaMgr    = ComputeCUDAMgr::getComputeCUDAMgr();
    if (calcGroup1COM) {
        computeGroupRestraint_2Group(
        mGpuOn,
        useDistMagnitude,
        doEnergy,
        doVirial,
        numRestrainedGroup1Local,
        numRestrainedGroup1Local + numRestrainedGroup2Local,
        this->restraintExp,
        this->restraintK,
        this->resCenterVec,
        this->resDirection,
        this->inv_group1_mass,
        this->inv_group2_mass,
        this->d_groupAtomsSOAIndex,
        lat,
        d_transform,
        d_mass,
        d_pos_x,
        d_pos_y,
        d_pos_z,
        d_f_normal_x,
        d_f_normal_y,
        d_f_normal_z,
        d_virial,
        h_extVirial,
        this->h_resEnergy,
        this->h_resForce,
        this->h_group1COM,
        this->h_group2COM,
        this->h_diffCOM,
        this->d_group1COM,
        this->d_group2COM,
        d_peer1COM,
        d_peer2COM,
        this->d_tbcatomic,
        numDevices,
        stream);
    } else {
      if(numRestrainedGroup2Local>0)
        {
          computeGroupRestraint_1Group(
                                         mGpuOn,
                                         useDistMagnitude,
                                         doEnergy,
                                         doVirial,
                                         numRestrainedGroup2Local,
                                         this->restraintExp,
                                         this->restraintK,
                                         this->resCenterVec,
                                         this->resDirection,
                                         this->inv_group2_mass,
                                         this->d_groupAtomsSOAIndex,
                                         lat,
                                         d_transform,
                                         d_mass,
                                         d_pos_x,
                                         d_pos_y,
                                         d_pos_z,
                                         d_f_normal_x,
                                         d_f_normal_y,
                                         d_f_normal_z,
                                         d_virial,
                                         h_extVirial,
                                         this->h_resEnergy,
                                         this->h_resForce,
                                         this->h_group1COM,
                                         this->h_group2COM,
                                         this->h_diffCOM,
                                         d_peer2COM,
                                         this->d_group2COM,
                                         this->d_tbcatomic,
                                         numDevices,
                                         stream);
        }
    }
    if(doOutput && (numDevices==1 || deviceIndex == cudaMgr->reducerGroupRestraintDevice)){
        cudaCheck(cudaStreamSynchronize(stream));
        // sum up external energy and virial from this group
        h_extEnergy[0] += h_resEnergy[0];
        // If we have restraint to reference point, then we have net force
        // otherwise the net force is zero for restraining two atom groups
        if (!calcGroup1COM) {
            h_extForce->x += h_resForce->x;
            h_extForce->y += h_resForce->y;
            h_extForce->z += h_resForce->z;
        }
      
        char msg[1024];
        sprintf(msg,"GRES: %9d %14s %14.4f %14.4f %14.4f %19.4f %14.4f %14.4f %14.4f\n",
            timeStep, groupName, h_diffCOM->x, h_diffCOM->y, h_diffCOM->z,
            h_resForce->x, h_resForce->y, h_resForce->z, h_resEnergy[0]);
        iout << msg << endi;

        // {
        //     printf("!!!Accu. exForce: %14.8f %14.8f %14.8f, Vir.x: %14.8f %14.8f %14.8f,"
        //     "Vir.y: %14.8f %14.8f %14.8f, Vir.z: %14.8f %14.8f %14.8f \n",
        //     h_extForce->x, h_extForce->y, h_extForce->z,
        //     h_extVirial->xx, h_extVirial->xy, h_extVirial->xz,
        //     h_extVirial->yx, h_extVirial->yy, h_extVirial->yz,
        //     h_extVirial->zx, h_extVirial->zy, h_extVirial->zz);
        // }
    }
}

// ###########################################################################
// # ComputeGroupRestraintsCUDA functions
// ###########################################################################

ComputeGroupRestraintsCUDA::ComputeGroupRestraintsCUDA(const int outputFreq,
                                                       const GroupRestraintList &resList, bool _mGpuOn, int _numDevices, int _deviceIndex) {
    gResOutputFreq = outputFreq;
    mGpuOn = _mGpuOn;
    numDevices = _numDevices;
    deviceIndex = _deviceIndex;
    const std::map<std::string, GroupRestraintParam*> & groupMap = resList.GetGroupResMap();
    ComputeCUDAMgr* cudaMgr    = ComputeCUDAMgr::getComputeCUDAMgr();
    for (auto it = groupMap.begin(); it != groupMap.end(); ++it) {
      GroupRestraintsCUDA * gResCUDA = new GroupRestraintsCUDA(it->second, mGpuOn, numDevices, deviceIndex, restraintsCUDAList.size());
        restraintsCUDAList.push_back(gResCUDA);
    }
}

ComputeGroupRestraintsCUDA::~ComputeGroupRestraintsCUDA() {
    int numGroup = restraintsCUDAList.size(); 
    for (int i = 0; i < numGroup; ++i) {
        delete restraintsCUDAList[i]; 
    }
    restraintsCUDAList.clear(); 
}

void ComputeGroupRestraintsCUDA::updateAtoms(
        std::vector<AtomMap*> &atomMapsList,
        std::vector<CudaLocalRecord> &localRecords,
        const int *h_globalToLocalID) {

    int numGroup = restraintsCUDAList.size();
    ComputeCUDAMgr* cudaMgr    = ComputeCUDAMgr::getComputeCUDAMgr();
    for (int i = 0; i < numGroup; ++i) {
        restraintsCUDAList[i]->updateAtoms(atomMapsList, localRecords, h_globalToLocalID);
        if(mGpuOn && restraintsCUDAList[i]->numRestrainedGroup1Local + restraintsCUDAList[i]->numRestrainedGroup2Local >0)
          {
            // last one wins
            cudaMgr->reducerGroupRestraintDevice.store(deviceIndex);
          }
    }
}

void ComputeGroupRestraintsCUDA::initPeerCOM(cudaStream_t stream){
  DebugM(3, "[" << CkMyPe() << "]" << " initPeerCOM\n" << endi);
    int numGroup = restraintsCUDAList.size();
    ComputeCUDAMgr* cudaMgr    = ComputeCUDAMgr::getComputeCUDAMgr();
    for (int i = 0; i < numGroup; ++i) {
      restraintsCUDAList[i]->initPeerCOM(cudaMgr->curGrp1COM[i], cudaMgr->curGrp2COM[i], stream);
    }

}


void ComputeGroupRestraintsCUDA::doCOM_mgpu(
        const Lattice &lat,
        const char3*  d_transform,
        const float*  d_mass,
        const double* d_pos_x,
        const double* d_pos_y,
        const double* d_pos_z,
        cudaStream_t  stream) {
  int numGroup = restraintsCUDAList.size();
  DebugM(3, "[" << CkMyPe() << "]" << " doCOM_mgpu "<< numGroup <<"\n" << endi);
  for (int gIdx = 0; gIdx < numGroup; ++gIdx) {
    //each element will determine if it needs a distributed COM
    restraintsCUDAList[gIdx]->doCOM_mgpu(
                                         lat, d_transform,
                                         d_mass, d_pos_x, d_pos_y, d_pos_z, gIdx,
                                         stream);
  }
}

void ComputeGroupRestraintsCUDA::doForce(
        const int timeStep, 
        const int doEnergy,
        const int doVirial,
        const Lattice &lat, 
        const char3*  d_transform, 
        const float*  d_mass, 
        const double* d_pos_x,
        const double* d_pos_y,
        const double* d_pos_z,
        double*       d_f_normal_x,
        double*       d_f_normal_y,
        double*       d_f_normal_z,
        cudaTensor*   d_virial,
        double*       h_extEnergy, 
        double3*      h_extForce, 
        cudaTensor*   h_extVirial, 
        cudaStream_t  stream) {

    const int doOutput = (timeStep % gResOutputFreq) == 0;
    // Since output freq is same as energyOutputFrq, we need to calculate virial
    // for outputting energy data
    int doVirCalc = (doOutput ? 1 : doVirial);
    int numGroup = restraintsCUDAList.size();

    // Reset the values before we add the energy, force, and virial value
    // for each restraint group
    h_extEnergy[0]  = 0.0;
    h_extForce->x   = 0.0;
    h_extForce->y   = 0.0; 
    h_extForce->z   = 0.0;  
    h_extVirial->xx = 0.0;
    h_extVirial->xy = 0.0;
    h_extVirial->xz = 0.0;
    h_extVirial->yx = 0.0;
    h_extVirial->yy = 0.0;
    h_extVirial->yz = 0.0;
    h_extVirial->zx = 0.0;
    h_extVirial->zy = 0.0;
    h_extVirial->zz = 0.0;
    ComputeCUDAMgr* cudaMgr    = ComputeCUDAMgr::getComputeCUDAMgr();
    if(doOutput && (numDevices==1 || deviceIndex == cudaMgr->reducerGroupRestraintDevice)){
      if(timeStep % (100 * gResOutputFreq) == 0) {
        char msg[1024];
        sprintf(msg,"\nGRES_TITLE: %3s %14s %14s %14s %14s %19s %14s %14s %14s\n",
            "TS", "GROUP_NAME", "DISTANCE.X", "DISTANCE.Y", "DISTANCE.Z",
            "FORCE.X", "FORCE.Y", "FORCE.Z", "ENERGY");
        iout << msg << endi;
      }
    }

    for (int gIdx = 0; gIdx < numGroup; ++gIdx) { 
        restraintsCUDAList[gIdx]->doForce( 
            timeStep, doEnergy, doVirCalc, doOutput, gIdx, lat, d_transform,
            d_mass, d_pos_x, d_pos_y, d_pos_z,
            d_f_normal_x, d_f_normal_y, d_f_normal_z, d_virial,
            h_extEnergy, h_extForce, h_extVirial, cudaMgr->curGrp1COM[gIdx],  cudaMgr->curGrp2COM[gIdx], stream);
    }
    if(doOutput && (numDevices==1 || deviceIndex == cudaMgr->reducerGroupRestraintDevice)) {
        iout <<"\n" << endi;;
    }
}
#endif