CudaGlobalMasterServer.C

Go to the documentation of this file.

#include "AtomMap.h"
#include "CudaGlobalMasterClient.h"
#include "CudaGlobalMasterServer.h"
#include "CudaGlobalMasterServerKernel.h"
#include "CudaUtils.h"
#include "DeviceCUDA.h"
#include "InfoStream.h"
#include "Lattice.h"
#include "NamdEventsProfiling.h"
#include "NamdTypes.h"
#include "ProcessorPrivate.h"
#include "SequencerCUDA.h"
#include "SimParameters.h"
#include <type_traits>
#include <string>

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

std::string ptr_to_str(const void *ptr) {
  std::ostringstream oss;
  oss << ptr;
  return oss.str();
}
#endif

struct localIDResult {
  bool found;
  // int pe;
  int soaIndex;
};

#if (defined(NAMD_CUDA) || defined(NAMD_HIP)) && defined(NODEGROUP_FORCE_REGISTER)

using CopyListTupleT = CudaGlobalMasterServer::CopyListTuple;
using ClientBufferT = CudaGlobalMasterServer::ClientBuffer;

extern __thread DeviceCUDA *deviceCUDA;

localIDResult queryLocalID(const AtomID globalID,
                           const std::vector<AtomMap *> &atomMapsList,
                           const std::vector<CudaLocalRecord> &localRecords,
                           const int *h_globalToLocalID) {
  // The logic here is similar to ComputeRestraintsCUDA::updateRestrainedAtoms
  localIDResult result{false, -1};
  for (int i = 0; i < atomMapsList.size(); ++i) {
    // Lookup a LocalID from the global atom ID
    const LocalID lid = atomMapsList[i]->localID(globalID);
    if (lid.pid != notUsed) {
      // Atom found
      result.found = true;
      // Mapping from global patch ID to local patch ID
      const int soaPid = h_globalToLocalID[lid.pid];
      // Memory location of the atom = patch start offset + atom local index in the patch
      result.soaIndex = localRecords[soaPid].bufferOffset + lid.index;
#ifdef DEBUGM
      DebugM(1, "Atom " + std::to_string(globalID) + " found in local patch " +
                    std::to_string(lid.pid) + " with local index " +
                    std::to_string(lid.index) + ", SOA patch ID " +
                    std::to_string(soaPid) + ", SOA global index " +
                    std::to_string(result.soaIndex) + "\n");
#endif
      break;
    }
  }
#ifdef DEBUGM
  if (result.found == false) {
    DebugM(3, "Atom " + std::to_string(globalID) + " not found.\n");
  }
#endif
  return result;
}

template <typename F1, typename F2>
void buildCopyList(
    F1 getAtomID, F2 getPosOrForceBuffer,
    std::vector<std::shared_ptr<CudaGlobalMaster::CudaGlobalMasterClient>> &clients,
    const std::vector<std::vector<AtomMap *>> &atomMapsLists,
    const std::vector<std::vector<CudaLocalRecord>> &localRecords,
    const std::vector<int *> h_globalToLocalID,
    const std::vector<int> &sourceDevicesList,
    const std::unordered_map<int, int> &deviceToIndex,
    ClientBufferT *&d_clientBuffers, std::vector<CopyListTupleT, CudaGlobalMasterServer::CudaHostAllocator<CopyListTupleT>> &hostCopyList,
    CopyListTupleT *&d_copyList, cudaStream_t stream, bool* checkUniqueList = nullptr) {
  // Number of devices == number of master PEs == number of
  // "std::vector<AtomMap*>"s
  std::vector<ClientBufferT> clientBuffers;
  hostCopyList.clear();
  const size_t numClients = clients.size();
  const int numDevices = atomMapsLists.size();
  for (size_t i = 0; i < numClients; ++i) {
    const auto &client = clients[i];
    // getPosOrForceBuffer should return the device memory address of the client
    // for copying to/from
    clientBuffers.push_back(ClientBufferT{
        ((*client).*getPosOrForceBuffer)(), client->getMasses(),
        client->getCharges(), client->getTransforms(), client->getVelocities(),
        ((*client).*getAtomID)().size()});
    const auto &requested_atoms = ((*client).*getAtomID)();
    // iout << iINFO << "PE: " << CkMyPe() << ", client atoms address " << &requested_atoms << "\n" << endi;
    // For each client, we try to map the global IDs of its
    // requested atoms to local SOA array indicies.
    // Also, taking the multi-GPU cases into account, we need
    // to bookkeep the device index of the peer arrays.
    for (size_t j = 0; j < requested_atoms.size(); ++j) {
      // The global ID of an atom
      const AtomID gid = requested_atoms[j];
      localIDResult result{false, -1};
      // We find the global ID in all atomMaps associated to all master PEs
      for (int k = 0; k < numDevices; ++k) {
        result = queryLocalID(gid, atomMapsLists[k], localRecords[k],
                              h_globalToLocalID[k]);
        if (result.found) {
          hostCopyList.push_back(CopyListTupleT{
              deviceToIndex.at(sourceDevicesList[k]), result.soaIndex, i, j});
          break;
        }
      }
      // The global ID should be found somewhere. Otherwise there is a bug.
      if (!result.found) {
        const std::string error =
            "Cannot find the local ID in SOA arrays of atom " +
            std::to_string(gid) + " requested by client[" + std::to_string(i) +
            "] (" + client->name() + ")\n";
        NAMD_bug(error.c_str());
      }
    }
  }
  // Sorting the copy list might improve the performance
  std::sort(hostCopyList.begin(), hostCopyList.end(),
            [](const CopyListTupleT &a, const CopyListTupleT &b) {
              return a.m_soa_index < b.m_soa_index;
            });
  if (checkUniqueList) {
    // Check if the list is unique
    auto tmp_list = hostCopyList;
    // std::sort(tmp_list.begin(), tmp_list.end(),
    //           [](const CopyListTupleT &a, const CopyListTupleT &b) {
    //             return a.m_soa_index < b.m_soa_index;
    //           });
    auto last = std::unique(tmp_list.begin(), tmp_list.end(),
                            [](const CopyListTupleT &a, const CopyListTupleT &b){
                              return a.m_soa_index == b.m_soa_index;
                            });
    if (last == tmp_list.end()) {
      *checkUniqueList = true;
    } else {
      *checkUniqueList = false;
    }
  }
  // Copy the copy list structure to the device memory
  size_t copySize = sizeof(CopyListTupleT) * hostCopyList.size();
#ifdef DEBUGM
  DebugM(3, "Will copy " + std::to_string(hostCopyList.size()) + " items.\n");
#endif
  // std::cout << "Trying to assign a copyList with size of " <<
  // std::to_string(copySize) << std::endl;
  if (d_copyList != nullptr) {
    cudaCheck(cudaFree(d_copyList));
  }
  cudaCheck(cudaMalloc(&d_copyList, copySize));
  cudaCheck(cudaMemcpyAsync(d_copyList, hostCopyList.data(), copySize,
                            cudaMemcpyHostToDevice, stream));
  copySize = sizeof(ClientBufferT) * clientBuffers.size();
  if (d_clientBuffers != nullptr) {
    cudaCheck(cudaFree(d_clientBuffers));
  }
  cudaCheck(cudaMalloc(&d_clientBuffers, copySize));
  cudaCheck(cudaMemcpyAsync(d_clientBuffers, clientBuffers.data(), copySize,
                            cudaMemcpyHostToDevice, stream));
}

void copyLatticeToClient(const std::vector<double, CudaGlobalMasterServer::CudaHostAllocator<double>>& lat, int deviceID,
                         std::shared_ptr<CudaGlobalMaster::CudaGlobalMasterClient> client, cudaStream_t stream) {
  const size_t copySize = sizeof(double) * lat.size();
  int savedDevice;
  cudaCheck(cudaGetDevice(&savedDevice));
  cudaCheck(cudaSetDevice(deviceID));
  double *d_lattice = client->getLattice();
  if (d_lattice) {
    cudaCheck(cudaMemcpyAsync(d_lattice, lat.data(), copySize,
                              cudaMemcpyHostToDevice, stream));
  } else {
    const std::string error = "Failed to copy lattice to client " + client->name() +
                              " (lattice requested but the client provides a "
                              "nullptr to the device buffer)";
    NAMD_die(error.c_str());
  }
  cudaCheck(cudaSetDevice(savedDevice));
}

#ifdef DEBUGM
template <typename T>
void debugCopyList(const std::string &name,
                   T &L) {
  std::cout << "CudaGlobalMasterServer: the copylist is " << name << ", with "
            << L.size() << " items.\n";
  for (size_t i = 0; i < L.size(); ++i) {
    fprintf(stdout,
            "i = %lu, deviceIndex = %d, soaIndex = %d, clientIndex = %lu, "
            "clientArrayIndex = %lu\n",
            i, L[i].m_src_dev_index, L[i].m_soa_index, L[i].m_client_index,
            L[i].m_client_atom_pos);
  }
}
#endif

#if 0
void debugClientBuffer(const std::string &name, ClientBufferT *B,
                       size_t numClients) {
  std::vector<ClientBufferT> hB(numClients);
  cudaPointerAttributes attributes;
  cudaCheck(cudaPointerGetAttributes(&attributes, B));
  std::cout << "CudaGlobalMasterServer: the clientBuffer is " << name
            << ", with " << numClients << " items.\n";
  std::cout << "deviceBuffer pointer = " << static_cast<void *>(B) << '\n';
  std::cout << "deviceBuffer attributes:\n";
  std::cout << "memory type: " << attributes.type << std::endl;
  cudaCheck(cudaMemcpy(hB.data(), B, sizeof(ClientBufferT) * numClients,
                       cudaMemcpyDeviceToHost));
  for (size_t i = 0; i < numClients; ++i) {
    fprintf(stdout,
            "i = %lu, d_data = %p, d_mass = %p, d_charge = %p, size = %lu\n", i,
            hB[i].d_data, hB[i].d_mass, hB[i].d_charge, hB[i].sz);
  }
  std::cout << std::endl;
}
#endif

CudaGlobalMasterServer::CudaGlobalMasterServer(int deviceID, int printProfilingFreq)
    : m_device_id(deviceID), m_step(0),
      m_num_devices(deviceCUDA->getNumDevice()), m_clients_changed(0),
      m_atom_maps_changed(0), m_print_profiling_freq(printProfilingFreq),
      reductionGpuResident(nullptr) {
  iout << iINFO << "CudaGlobalMasterServer: initialized on PE " << CkMyPe()
       << " and GPU device " << m_device_id << "\n"
       << endi;
  int savedDevice;
  cudaCheck(cudaGetDevice(&savedDevice));
  cudaCheck(cudaSetDevice(m_device_id));
  cudaCheck(cudaStreamCreate(&m_stream));
#ifdef NAMD_NVTX_ENABLED
  nvtxNameCuStreamA(m_stream, "CudaGlobalMaster stream");
#endif
  m_atom_map_lists.resize(m_num_devices);
  m_src_devs.resize(m_num_devices);
  m_local_records.resize(m_num_devices);
  m_global_to_local_id.resize(m_num_devices);
  const int *allDevices = deviceCUDA->allDevices();
  for (int i = 0; i < m_num_devices; ++i) {
    const int currentDeviceID = allDevices[i];
    m_device_id_to_index[currentDeviceID] = i;
    m_src_devs[i] = currentDeviceID;
  }
  if (m_num_devices > 1) {
    allocatePeerArrays();
  }
  // Copy lists
  m_d_atom_pos_copy_list = nullptr;
  m_d_atom_total_force_copy_list = nullptr;
  m_d_forced_atom_copy_list = nullptr;
  // Client buffers
  m_d_atom_pos_client_buffers = nullptr;
  m_atom_total_force_client_buffers = nullptr;
  m_d_forced_atom_client_buffers = nullptr;
  m_unique_forced_atoms = false;
  cudaCheck(cudaSetDevice(savedDevice));
  constexpr int m_h_lattice_size =
    (sizeof(std::invoke_result_t<decltype(&Lattice::a), Lattice>) +
     sizeof(std::invoke_result_t<decltype(&Lattice::b), Lattice>) +
     sizeof(std::invoke_result_t<decltype(&Lattice::c), Lattice>) +
     sizeof(std::invoke_result_t<decltype(&Lattice::origin), Lattice>)) /
    sizeof(std::invoke_result_t<decltype(&Vector::operator[]), Vector, int>);
  m_h_lattice.resize(m_h_lattice_size);
  // Setup the reductions
  auto* params = Node::Object()->simParameters;
  if (!params->CUDASOAintegrateMode) {
    NAMD_bug("CudaGlobalMasterServer only supports GPU-resident mode");
  }
  reductionGpuResident = ReductionMgr::Object()->willSubmit(REDUCTIONS_GPURESIDENT);
}

CudaGlobalMasterServer::~CudaGlobalMasterServer() {
  if (m_d_atom_pos_copy_list != nullptr) {
    cudaCheck(cudaFree(m_d_atom_pos_copy_list));
    m_d_atom_pos_copy_list = nullptr;
  }
  if (m_d_atom_pos_client_buffers != nullptr) {
    cudaCheck(cudaFree(m_d_atom_pos_client_buffers));
    m_d_atom_pos_client_buffers = nullptr;
  }
  if (m_d_atom_total_force_copy_list != nullptr) {
    cudaCheck(cudaFree(m_d_atom_total_force_copy_list));
    m_d_atom_total_force_copy_list = nullptr;
  }
  if (m_atom_total_force_client_buffers != nullptr) {
    cudaCheck(cudaFree(m_atom_total_force_client_buffers));
    m_atom_total_force_client_buffers = nullptr;
  }
  if (m_d_forced_atom_copy_list != nullptr) {
    cudaCheck(cudaFree(m_d_forced_atom_copy_list));
    m_d_forced_atom_copy_list = nullptr;
  }
  if (m_d_forced_atom_client_buffers != nullptr) {
    cudaCheck(cudaFree(m_d_forced_atom_client_buffers));
    m_d_forced_atom_client_buffers = nullptr;
  }
  if (reductionGpuResident != nullptr) {
    delete reductionGpuResident;
  }
  iout << iINFO << "CudaGlobalMasterServer: destructed on PE " << CkMyPe()
       << "\n"
       << endi;
}

void CudaGlobalMasterServer::addClient(
    std::shared_ptr<CudaGlobalMaster::CudaGlobalMasterClient> client) {
  auto it = std::find(m_clients.begin(), m_clients.end(), client);
  if (it == m_clients.end()) {
    iout << iINFO << "CudaGlobalMasterServer: adding client \""
         << client->name() << "\"\n"
         << endi;
    m_clients.push_back(client);
    m_clients_changed = CudaGlobalMasterServer::numCopyLists;
  } else {
    const std::string error =
        "The client \"" + client->name() + "\" are being added twice.\n";
    NAMD_die(error.c_str());
  }
}

void CudaGlobalMasterServer::removeClient(
    std::shared_ptr<CudaGlobalMaster::CudaGlobalMasterClient> client) {
  auto it = std::find(m_clients.begin(), m_clients.end(), client);
  if (it == m_clients.end()) {
    iout << iWARN << "CudaGlobalMasterServer: the client \"" << client->name()
         << "\" is not registered with CudaGlobalMasterServer\n"
         << endi;
  }
  while (it != m_clients.end()) {
    iout << iINFO << "CudaGlobalMasterServer: removing client \""
         << client->name() << "\"\n"
         << endi;
    m_clients.erase(it);
    m_clients_changed = CudaGlobalMasterServer::numCopyLists;
    it = std::find(m_clients.begin(), m_clients.end(), client);
  }
}

void CudaGlobalMasterServer::updateAtomMaps() {
  const int numPes = CkNumPes();
  // Iterate over all PEs
#ifdef DEBUGM
  DebugM(3, "updateAtomMaps: number of PEs = " + std::to_string(numPes) + "\n");
#endif
  for (int i = 0; i < numPes; ++i) {
    // Find the device ID of the i-th PE
    const int peDevice = deviceCUDA->getDeviceIDforPe(i);
    const int j = m_device_id_to_index.at(peDevice);
    // Get the atom maps
    AtomMap *amap = AtomMap::ObjectOnPe(i);
#ifdef DEBUGM
    DebugM(3, "updateAtomMaps: PE " + std::to_string(i) + " atomMap " +
                  ptr_to_str(amap) + " on device " + std::to_string(peDevice) +
                  "\n");
#endif
    m_atom_map_lists[j].push_back(amap);
  }
  const bool multi_gpu = m_num_devices > 1;
  // Iterate over all devices to get the map of global patch ID to local patch
  // ID
  // TODO: I assume all devices are used. Is this correct?
  for (int i = 0; i < m_num_devices; ++i) {
    const int deviceID = m_src_devs[i];
    // Get the master PE
    const int masterPe = deviceCUDA->getMasterPeForDeviceID(deviceID);
    // Get the corresponding SequencerCUDA instance
    const SequencerCUDA *sequencer = SequencerCUDA::ObjectOnPe(masterPe);
    m_global_to_local_id[i] = sequencer->globalToLocalID;
    m_local_records[i] = sequencer->patchData->devData[deviceID].h_localPatches;
    if (multi_gpu) {
      m_h_peer_atom_data.d_pos_x[i] = sequencer->d_pos_x;
      m_h_peer_atom_data.d_pos_y[i] = sequencer->d_pos_y;
      m_h_peer_atom_data.d_pos_z[i] = sequencer->d_pos_z;
      m_h_peer_atom_data.d_vel_x[i] = sequencer->d_vel_x;
      m_h_peer_atom_data.d_vel_y[i] = sequencer->d_vel_y;
      m_h_peer_atom_data.d_vel_z[i] = sequencer->d_vel_z;
      m_h_peer_atom_data.d_mass[i] = sequencer->d_mass;
      m_h_peer_atom_data.d_charge[i] = sequencer->d_charge;
      m_h_peer_atom_data.d_transform[i] = sequencer->d_transform;

      m_h_peer_tf_array.d_f_normal_x[i] = sequencer->d_f_normal_x;
      m_h_peer_tf_array.d_f_normal_y[i] = sequencer->d_f_normal_y;
      m_h_peer_tf_array.d_f_normal_z[i] = sequencer->d_f_normal_z;
      m_h_peer_tf_array.d_f_saved_nbond_x[i] = sequencer->d_f_saved_nbond_x;
      m_h_peer_tf_array.d_f_saved_nbond_y[i] = sequencer->d_f_saved_nbond_y;
      m_h_peer_tf_array.d_f_saved_nbond_z[i] = sequencer->d_f_saved_nbond_z;
      m_h_peer_tf_array.d_f_saved_slow_x[i] = sequencer->d_f_saved_slow_x;
      m_h_peer_tf_array.d_f_saved_slow_y[i] = sequencer->d_f_saved_slow_y;
      m_h_peer_tf_array.d_f_saved_slow_z[i] = sequencer->d_f_saved_slow_z;
      m_h_peer_tf_array.d_atomFixed[i] = sequencer->d_atomFixed;

      m_h_peer_af_array.d_f_applied_x[i] = sequencer->d_f_global_x;
      m_h_peer_af_array.d_f_applied_y[i] = sequencer->d_f_global_y;
      m_h_peer_af_array.d_f_applied_z[i] = sequencer->d_f_global_z;
      m_h_peer_af_array.d_atomFixed[i] = sequencer->d_atomFixed;
    }
#ifdef DEBUGM
    DebugM(3, "updateAtomMaps: device " + std::to_string(deviceID) +
                  ", sequencer " + ptr_to_str(sequencer) + "\n");
#endif
  }
  if (multi_gpu) {
    copyPeerArraysToDevice();
  }
  // Setup the flag to rebuild all copy lists
  m_atom_maps_changed = CudaGlobalMasterServer::numCopyLists;
}

void CudaGlobalMasterServer::setStep(int64_t step) {
  m_step = step;
  // Update the step number
  for (auto& client: m_clients) {
    client->setStep(step);
  }
}

void CudaGlobalMasterServer::communicateToClients(const Lattice* lat) {
  // iout << iINFO << "PE: " << CkMyPe() << ", communicateToClients this address " << this << "\n" << endi;
  const Vector lat_a = lat->a();
  const Vector lat_b = lat->b();
  const Vector lat_c = lat->c();
  const Vector lat_o = lat->origin();
  m_h_lattice[0] = lat_a.x;
  m_h_lattice[1] = lat_a.y;
  m_h_lattice[2] = lat_a.z;
  m_h_lattice[3] = lat_b.x;
  m_h_lattice[4] = lat_b.y;
  m_h_lattice[5] = lat_b.z;
  m_h_lattice[6] = lat_c.x;
  m_h_lattice[7] = lat_c.y;
  m_h_lattice[8] = lat_c.z;
  m_h_lattice[9]  = lat_o.x;
  m_h_lattice[10] = lat_o.y;
  m_h_lattice[11] = lat_o.z;
  bool b_buildAtomsPositionCopyList = false;
  bool b_buildAtomsTotalForcesCopyList = false;
  for (auto&& client: m_clients) {
    // Check all clients if their requested atoms are changed
    if (client->requestedAtomsChanged()) {
      b_buildAtomsPositionCopyList = true;
    }
    // Check all clients if their requested total force atoms are changed
    if (client->requestedTotalForcesAtomsChanged()) {
      b_buildAtomsTotalForcesCopyList = true;
    }
    // Copy lattice to client buffers if necessary
    if (client->requestUpdateLattice()) {
      copyLatticeToClient(m_h_lattice, m_device_id, client, m_stream);
    }
  }
  // Rebuild the list of atoms to copy if necessary
  if (b_buildAtomsPositionCopyList || m_atom_maps_changed > 0 ||
      m_clients_changed > 0) {
    buildAtomsCopyList();
    if (m_atom_maps_changed > 0)
      m_atom_maps_changed--;
    if (m_clients_changed > 0)
      m_clients_changed--;
  }
  if (b_buildAtomsTotalForcesCopyList || m_atom_maps_changed > 0 ||
      m_clients_changed > 0) {
    buildAtomsTotalForcesCopyList();
    if (m_atom_maps_changed > 0)
      m_atom_maps_changed--;
    if (m_clients_changed > 0)
      m_clients_changed--;
  }
  // Check all clients if they request to update the atomic positions
  bool b_copyPositions = false;
  bool b_copyTotalForces = false;
  bool b_copyMasses = false;
  bool b_copyCharges = false;
  bool b_copyTransforms = false;
  bool b_copyVelocities = false;
  for (auto&& client: m_clients) {
    b_copyPositions |= client->requestUpdateAtomPositions();
    b_copyTotalForces |= client->requestUpdateAtomTotalForces();
    b_copyMasses |= client->requestUpdateMasses();
    b_copyCharges |= client->requestUpdateCharges();
    b_copyTransforms |= client->requestUpdateTransforms();
    b_copyVelocities |= client->requestUpdateVelocities();
  }
  // Update the atomic positions if necessary
  if (b_buildAtomsPositionCopyList || b_copyPositions || b_copyMasses ||
      b_copyCharges || b_copyTransforms || b_copyVelocities) {
    copyAtomsToClients(b_copyPositions, b_copyMasses, b_copyCharges,
                       b_copyTransforms, b_copyVelocities);
  }
  if (b_copyTotalForces || b_buildAtomsTotalForcesCopyList) {
    copyTotalForcesToClients();
  }
  for (auto&& client: m_clients) {
    client->onBuffersUpdated();
  }
}

bool CudaGlobalMasterServer::requestedTotalForces() const {
  bool result = false;
  for (auto&& client: m_clients) {
    result |= (!(client->getRequestedForcedAtoms().empty()) &&
               client->requestUpdateAtomTotalForces());
  }
  return result;
}

void CudaGlobalMasterServer::buildAtomsCopyList() {
#ifdef DEBUGM
  DebugM(3, "buildAtomsCopyList is called\n");
#endif
  // Save the current device ID
  int savedDevice;
  cudaCheck(cudaGetDevice(&savedDevice));
  cudaCheck(cudaSetDevice(m_device_id));
  buildCopyList(&CudaGlobalMaster::CudaGlobalMasterClient::getRequestedAtoms,
                &CudaGlobalMaster::CudaGlobalMasterClient::getPositions, m_clients,
                m_atom_map_lists, m_local_records, m_global_to_local_id,
                m_src_devs, m_device_id_to_index, m_d_atom_pos_client_buffers,
                m_atom_pos_copy_list, m_d_atom_pos_copy_list, m_stream);
  // cudaCheck(cudaStreamSynchronize(m_stream));
  cudaCheck(cudaSetDevice(savedDevice));
}

void CudaGlobalMasterServer::copyAtomsToClients(bool copyPositions,
                                                bool copyMasses,
                                                bool copyCharges,
                                                bool copyTransforms,
                                                bool copyVelocities) {
#ifdef DEBUGM
  DebugM(1, "copyAtomsToClients is called\n");
#endif
  NAMD_EVENT_START(1, NamdProfileEvent::CUDAGM_ATOMTOCLIENTS);
  // Save the current device ID
  int savedDevice;
  cudaCheck(cudaGetDevice(&savedDevice));
  cudaCheck(cudaSetDevice(m_device_id));
  if (m_num_devices == 1) {
#ifdef DEBUGM
#ifdef MIN_DEBUG_LEVEL
#if MIN_DEBUG_LEVEL <= 1
    debugCopyList("CudaGlobalMasterServer::copyAtomsToClients",
                  m_atom_pos_copy_list);
#endif
#endif
#endif
    const int masterPe = deviceCUDA->getMasterPe();
    const SequencerCUDA *sequencer = SequencerCUDA::ObjectOnPe(masterPe);
    copyAtomsToClientsCUDA(
        copyPositions, copyMasses, copyCharges, copyTransforms, copyVelocities,
        sequencer->d_pos_x, sequencer->d_pos_y, sequencer->d_pos_z,
        sequencer->d_vel_x, sequencer->d_vel_y, sequencer->d_vel_z,
        sequencer->d_transform, sequencer->d_mass, sequencer->d_charge,
        sequencer->myLattice, this->m_d_atom_pos_copy_list,
        this->m_atom_pos_copy_list.size(), this->m_d_atom_pos_client_buffers,
        m_clients.size(), m_stream);
  } else {
    const int masterPe = deviceCUDA->getMasterPe();
    const SequencerCUDA *sequencer = SequencerCUDA::ObjectOnPe(masterPe);
    copyAtomsToClientsCUDAMGPU(
        copyPositions, copyMasses, copyCharges, copyTransforms, copyVelocities,
        (const double **)m_h_peer_atom_data.d_pos_x,
        (const double **)m_h_peer_atom_data.d_pos_y,
        (const double **)m_h_peer_atom_data.d_pos_z,
        (const double **)m_h_peer_atom_data.d_vel_x,
        (const double **)m_h_peer_atom_data.d_vel_y,
        (const double **)m_h_peer_atom_data.d_vel_z,
        (const char3 **)m_h_peer_atom_data.d_transform,
        (const float **)m_h_peer_atom_data.d_mass,
        (const float **)m_h_peer_atom_data.d_charge, sequencer->myLattice,
        this->m_d_atom_pos_copy_list, this->m_atom_pos_copy_list.size(),
        this->m_d_atom_pos_client_buffers, m_clients.size(), m_stream);
  }
  // cudaCheck(cudaStreamSynchronize(m_stream));
  cudaCheck(cudaSetDevice(savedDevice));
  NAMD_EVENT_STOP(1, NamdProfileEvent::CUDAGM_ATOMTOCLIENTS);
}

void CudaGlobalMasterServer::copyTotalForcesToClients() {
#ifdef DEBUGM
  DebugM(1, "copyTotalForcesToClients is called\n");
#endif
  const SimParameters *simParams = Node::Object()->simParameters;
  NAMD_EVENT_START(1, NamdProfileEvent::CUDAGM_TOTALFORCETOCLIENTS);
  // Save the current device ID
  int savedDevice;
  cudaCheck(cudaGetDevice(&savedDevice));
  cudaCheck(cudaSetDevice(m_device_id));
  if (m_num_devices == 1) {
#ifdef DEBUGM
#ifdef MIN_DEBUG_LEVEL
#if MIN_DEBUG_LEVEL <= 1
    debugCopyList("CudaGlobalMasterServer::copyTotalForcesToClients",
                  m_atom_total_force_copy_list);
#endif
#endif
#endif
    const int masterPe = deviceCUDA->getMasterPe();
    const SequencerCUDA *sequencer = SequencerCUDA::ObjectOnPe(masterPe);
    copyTotalForcesToClientsCUDA(
        simParams->fixedAtomsOn, sequencer->d_f_normal_x,
        sequencer->d_f_normal_y, sequencer->d_f_normal_z,
        sequencer->d_f_saved_nbond_x, sequencer->d_f_saved_nbond_y,
        sequencer->d_f_saved_nbond_z, sequencer->d_f_saved_slow_x,
        sequencer->d_f_saved_slow_y, sequencer->d_f_saved_slow_z,
        sequencer->d_atomFixed, this->m_d_atom_total_force_copy_list,
        this->m_atom_total_force_copy_list.size(),
        this->m_atom_total_force_client_buffers, m_clients.size(), m_stream);
  } else {
    copyTotalForcesToClientsCUDAMGPU(
        simParams->fixedAtomsOn,
        (const double **)m_d_peer_tf_array.d_f_normal_x,
        (const double **)m_d_peer_tf_array.d_f_normal_y,
        (const double **)m_d_peer_tf_array.d_f_normal_z,
        (const double **)m_d_peer_tf_array.d_f_saved_nbond_x,
        (const double **)m_d_peer_tf_array.d_f_saved_nbond_y,
        (const double **)m_d_peer_tf_array.d_f_saved_nbond_z,
        (const double **)m_d_peer_tf_array.d_f_saved_slow_x,
        (const double **)m_d_peer_tf_array.d_f_saved_slow_y,
        (const double **)m_d_peer_tf_array.d_f_saved_slow_z,
        (const int **)m_d_peer_tf_array.d_atomFixed,
        this->m_d_atom_total_force_copy_list,
        this->m_atom_total_force_copy_list.size(),
        this->m_atom_total_force_client_buffers, m_clients.size(), m_stream);
  }
  // cudaCheck(cudaStreamSynchronize(m_stream));
  cudaCheck(cudaSetDevice(savedDevice));
  NAMD_EVENT_STOP(1, NamdProfileEvent::CUDAGM_TOTALFORCETOCLIENTS);
}

void CudaGlobalMasterServer::buildAtomsTotalForcesCopyList() {
#ifdef DEBUGM
  DebugM(3, "buildAtomsTotalForcesCopyList is called\n");
#endif
  // Save the current device ID
  int savedDevice;
  cudaCheck(cudaGetDevice(&savedDevice));
  cudaCheck(cudaSetDevice(m_device_id));
  buildCopyList(&CudaGlobalMaster::CudaGlobalMasterClient::getRequestedTotalForcesAtoms,
                &CudaGlobalMaster::CudaGlobalMasterClient::getTotalForces, m_clients,
                m_atom_map_lists, m_local_records, m_global_to_local_id,
                m_src_devs, m_device_id_to_index,
                m_atom_total_force_client_buffers, m_atom_total_force_copy_list,
                m_d_atom_total_force_copy_list, m_stream);
  cudaCheck(cudaSetDevice(savedDevice));
}

void CudaGlobalMasterServer::buildForcedAtomsCopyList() {
#ifdef DEBUGM
  DebugM(3, "buildForcedAtomsCopyList is called\n");
#endif
  // Save the current device ID
  int savedDevice;
  cudaCheck(cudaGetDevice(&savedDevice));
  cudaCheck(cudaSetDevice(m_device_id));
  buildCopyList(&CudaGlobalMaster::CudaGlobalMasterClient::getRequestedForcedAtoms,
                &CudaGlobalMaster::CudaGlobalMasterClient::getAppliedForces, m_clients,
                m_atom_map_lists, m_local_records, m_global_to_local_id,
                m_src_devs, m_device_id_to_index,
                m_d_forced_atom_client_buffers, m_forced_atom_copy_list,
                m_d_forced_atom_copy_list, m_stream, &m_unique_forced_atoms);
  cudaCheck(cudaSetDevice(savedDevice));
}

void CudaGlobalMasterServer::calculate() {
  NAMD_EVENT_START(1, NamdProfileEvent::CUDAGM_CALCULATECLIENTS);
  for (auto&& client: m_clients) {
    client->calculate();
  }
  NAMD_EVENT_STOP(1, NamdProfileEvent::CUDAGM_CALCULATECLIENTS);
}

void CudaGlobalMasterServer::communicateToMD(bool doEnergy, bool doVirial) {
  // iout << iINFO << "PE: " << CkMyPe() << ", communicateToMD this address " << this << "\n" << endi;
#ifdef DEBUGM
  DebugM(1, "Calling communicateToMD at step " + std::to_string(m_step));
#endif
  // Call the calculate function of every client
  bool b_buildForcedAtomsCopyList = false;
  for (auto&& client: m_clients) {
    if (client->requestedForcedAtomsChanged()) {
      b_buildForcedAtomsCopyList = true;
    }
  }
  if (b_buildForcedAtomsCopyList || m_atom_maps_changed > 0 ||
      m_clients_changed > 0) {
    buildForcedAtomsCopyList();
    if (m_atom_maps_changed > 0)
      m_atom_maps_changed--;
    if (m_clients_changed > 0)
      m_clients_changed--;
  }
  bool b_copyForcedAtoms = false;
  for (auto&& client: m_clients) {
    if (client->requestUpdateForcedAtoms()) {
      b_copyForcedAtoms = true;
    }
  }
  NAMD_EVENT_START(1, NamdProfileEvent::CUDAGM_ADDGLOBALFORCES);
  if (b_copyForcedAtoms) {
    addGlobalForces();
  }
  auto* currentReduction = getCurrentReduction();
  for (auto& client: m_clients) {
    client->finishReductions(doEnergy, doVirial, currentReduction);
  }
  cudaCheck(cudaStreamSynchronize(m_stream));
  NAMD_EVENT_STOP(1, NamdProfileEvent::CUDAGM_ADDGLOBALFORCES);
}

bool CudaGlobalMasterServer::willAddGlobalForces() const {
  bool result = false;
  for (auto&& client: m_clients) {
    result |= (!(client->getRequestedForcedAtoms().empty()) &&
               client->requestUpdateForcedAtoms());
  }
  return result;
}

void CudaGlobalMasterServer::addGlobalForces() {
#ifdef DEBUGM
  DebugM(1, "Calling addGlobalForces at step " + std::to_string(m_step));
#endif
  const SimParameters *simParams = Node::Object()->simParameters;
  // Save the current device ID
  int savedDevice;
  cudaCheck(cudaGetDevice(&savedDevice));
  cudaCheck(cudaSetDevice(m_device_id));
  if (m_num_devices == 1) {
    const int masterPe = deviceCUDA->getMasterPe();
    const SequencerCUDA *sequencer = SequencerCUDA::ObjectOnPe(masterPe);
    addGlobalForcesFromClients(
        simParams->fixedAtomsOn, m_unique_forced_atoms, sequencer->d_f_global_x,
        sequencer->d_f_global_y, sequencer->d_f_global_z,
        sequencer->d_atomFixed, this->m_d_forced_atom_copy_list,
        this->m_forced_atom_copy_list.size(),
        this->m_d_forced_atom_client_buffers, m_clients.size(), m_stream);
  } else {
    addGlobalForcesFromClientsMGPU(
        simParams->fixedAtomsOn, m_unique_forced_atoms, m_d_peer_af_array.d_f_applied_x,
        m_d_peer_af_array.d_f_applied_y, m_d_peer_af_array.d_f_applied_z,
        (const int **)m_d_peer_af_array.d_atomFixed,
        this->m_d_forced_atom_copy_list, this->m_forced_atom_copy_list.size(),
        this->m_d_forced_atom_client_buffers, m_clients.size(), m_stream);
  }
  // cudaCheck(cudaStreamSynchronize(m_stream));
  cudaCheck(cudaSetDevice(savedDevice));
}

void CudaGlobalMasterServer::allocatePeerArrays() {
#ifdef DEBUGM
  DebugM(3, "CudaGlobalMasterServer::allocatePeerArrays");
#endif
  allocate_host<double *>(&(m_h_peer_atom_data.d_pos_x), m_num_devices);
  allocate_host<double *>(&(m_h_peer_atom_data.d_pos_y), m_num_devices);
  allocate_host<double *>(&(m_h_peer_atom_data.d_pos_z), m_num_devices);
  allocate_host<double *>(&(m_h_peer_atom_data.d_vel_x), m_num_devices);
  allocate_host<double *>(&(m_h_peer_atom_data.d_vel_y), m_num_devices);
  allocate_host<double *>(&(m_h_peer_atom_data.d_vel_z), m_num_devices);
  allocate_host<float *>(&(m_h_peer_atom_data.d_mass), m_num_devices);
  allocate_host<float *>(&(m_h_peer_atom_data.d_charge), m_num_devices);
  allocate_host<char3 *>(&(m_h_peer_atom_data.d_transform), m_num_devices);

  allocate_host<tf_type *>(&(m_h_peer_tf_array.d_f_normal_x), m_num_devices);
  allocate_host<tf_type *>(&(m_h_peer_tf_array.d_f_normal_y), m_num_devices);
  allocate_host<tf_type *>(&(m_h_peer_tf_array.d_f_normal_z), m_num_devices);
  allocate_host<tf_type *>(&(m_h_peer_tf_array.d_f_saved_nbond_x),
                           m_num_devices);
  allocate_host<tf_type *>(&(m_h_peer_tf_array.d_f_saved_nbond_y),
                           m_num_devices);
  allocate_host<tf_type *>(&(m_h_peer_tf_array.d_f_saved_nbond_z),
                           m_num_devices);
  allocate_host<tf_type *>(&(m_h_peer_tf_array.d_f_saved_slow_x),
                           m_num_devices);
  allocate_host<tf_type *>(&(m_h_peer_tf_array.d_f_saved_slow_y),
                           m_num_devices);
  allocate_host<tf_type *>(&(m_h_peer_tf_array.d_f_saved_slow_z),
                           m_num_devices);
  allocate_host<int *>(&(m_h_peer_tf_array.d_atomFixed), m_num_devices);

  allocate_host<tf_type *>(&(m_h_peer_af_array.d_f_applied_x), m_num_devices);
  allocate_host<tf_type *>(&(m_h_peer_af_array.d_f_applied_y), m_num_devices);
  allocate_host<tf_type *>(&(m_h_peer_af_array.d_f_applied_z), m_num_devices);
  allocate_host<int *>(&(m_h_peer_af_array.d_atomFixed), m_num_devices);

  // Save the current device ID
  int savedDevice;
  cudaCheck(cudaGetDevice(&savedDevice));
  cudaCheck(cudaSetDevice(m_device_id));
  allocate_device<double *>(&(m_d_peer_atom_data.d_pos_x), m_num_devices);
  allocate_device<double *>(&(m_d_peer_atom_data.d_pos_y), m_num_devices);
  allocate_device<double *>(&(m_d_peer_atom_data.d_pos_z), m_num_devices);
  allocate_device<double *>(&(m_d_peer_atom_data.d_vel_x), m_num_devices);
  allocate_device<double *>(&(m_d_peer_atom_data.d_vel_y), m_num_devices);
  allocate_device<double *>(&(m_d_peer_atom_data.d_vel_z), m_num_devices);
  allocate_device<float *>(&(m_d_peer_atom_data.d_mass), m_num_devices);
  allocate_device<float *>(&(m_d_peer_atom_data.d_charge), m_num_devices);
  allocate_device<char3 *>(&(m_d_peer_atom_data.d_transform), m_num_devices);

  allocate_device<tf_type *>(&(m_d_peer_tf_array.d_f_normal_x), m_num_devices);
  allocate_device<tf_type *>(&(m_d_peer_tf_array.d_f_normal_y), m_num_devices);
  allocate_device<tf_type *>(&(m_d_peer_tf_array.d_f_normal_z), m_num_devices);
  allocate_device<tf_type *>(&(m_d_peer_tf_array.d_f_saved_nbond_x),
                             m_num_devices);
  allocate_device<tf_type *>(&(m_d_peer_tf_array.d_f_saved_nbond_y),
                             m_num_devices);
  allocate_device<tf_type *>(&(m_d_peer_tf_array.d_f_saved_nbond_z),
                             m_num_devices);
  allocate_device<tf_type *>(&(m_d_peer_tf_array.d_f_saved_slow_x),
                             m_num_devices);
  allocate_device<tf_type *>(&(m_d_peer_tf_array.d_f_saved_slow_y),
                             m_num_devices);
  allocate_device<tf_type *>(&(m_d_peer_tf_array.d_f_saved_slow_z),
                             m_num_devices);
  allocate_device<int *>(&(m_d_peer_tf_array.d_atomFixed), m_num_devices);

  allocate_device<tf_type *>(&(m_d_peer_af_array.d_f_applied_x), m_num_devices);
  allocate_device<tf_type *>(&(m_d_peer_af_array.d_f_applied_y), m_num_devices);
  allocate_device<tf_type *>(&(m_d_peer_af_array.d_f_applied_z), m_num_devices);
  allocate_device<int *>(&(m_d_peer_af_array.d_atomFixed), m_num_devices);
  cudaCheck(cudaSetDevice(savedDevice));
}

void CudaGlobalMasterServer::copyPeerArraysToDevice() {
#ifdef DEBUGM
  DebugM(3, "CudaGlobalMasterServer::copyPeerArraysToDevice");
#endif
  // Save the current device ID
  int savedDevice;
  cudaCheck(cudaGetDevice(&savedDevice));
  cudaCheck(cudaSetDevice(m_device_id));
  copy_HtoD(m_h_peer_atom_data.d_pos_x, m_d_peer_atom_data.d_pos_x,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_atom_data.d_pos_y, m_d_peer_atom_data.d_pos_y,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_atom_data.d_pos_z, m_d_peer_atom_data.d_pos_z,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_atom_data.d_vel_x, m_d_peer_atom_data.d_vel_x,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_atom_data.d_vel_y, m_d_peer_atom_data.d_vel_y,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_atom_data.d_vel_z, m_d_peer_atom_data.d_vel_z,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_atom_data.d_mass, m_d_peer_atom_data.d_mass, m_num_devices,
            m_stream);
  copy_HtoD(m_h_peer_atom_data.d_charge, m_d_peer_atom_data.d_charge,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_atom_data.d_transform, m_d_peer_atom_data.d_transform,
            m_num_devices, m_stream);

  copy_HtoD(m_h_peer_tf_array.d_f_normal_x, m_d_peer_tf_array.d_f_normal_x,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_tf_array.d_f_normal_y, m_d_peer_tf_array.d_f_normal_y,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_tf_array.d_f_normal_z, m_d_peer_tf_array.d_f_normal_z,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_tf_array.d_f_saved_nbond_x,
            m_d_peer_tf_array.d_f_saved_nbond_x, m_num_devices, m_stream);
  copy_HtoD(m_h_peer_tf_array.d_f_saved_nbond_y,
            m_d_peer_tf_array.d_f_saved_nbond_y, m_num_devices, m_stream);
  copy_HtoD(m_h_peer_tf_array.d_f_saved_nbond_z,
            m_d_peer_tf_array.d_f_saved_nbond_z, m_num_devices, m_stream);
  copy_HtoD(m_h_peer_tf_array.d_f_saved_slow_x,
            m_d_peer_tf_array.d_f_saved_slow_x, m_num_devices, m_stream);
  copy_HtoD(m_h_peer_tf_array.d_f_saved_slow_y,
            m_d_peer_tf_array.d_f_saved_slow_y, m_num_devices, m_stream);
  copy_HtoD(m_h_peer_tf_array.d_f_saved_slow_z,
            m_d_peer_tf_array.d_f_saved_slow_z, m_num_devices, m_stream);
  copy_HtoD(m_h_peer_tf_array.d_atomFixed, m_d_peer_tf_array.d_atomFixed,
            m_num_devices, m_stream);

  copy_HtoD(m_h_peer_af_array.d_f_applied_x, m_d_peer_af_array.d_f_applied_x,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_af_array.d_f_applied_y, m_d_peer_af_array.d_f_applied_y,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_af_array.d_f_applied_z, m_d_peer_af_array.d_f_applied_z,
            m_num_devices, m_stream);
  copy_HtoD(m_h_peer_af_array.d_atomFixed, m_d_peer_af_array.d_atomFixed,
            m_num_devices, m_stream);

  // cudaCheck(cudaStreamSynchronize(m_stream));
  cudaCheck(cudaSetDevice(savedDevice));
}

void CudaGlobalMasterServer::finishReductions() {
#ifdef DEBUGM
  DebugM(1, "Calling finishReductions at step " + std::to_string(m_step));
#endif
  // CkPrintf("Calling CudaGlobalMasterServer::finishReductions\n");
  auto* currentReduction = getCurrentReduction();
  currentReduction->submit();
}

#else

CudaGlobalMasterServer::CudaGlobalMasterServer(
  int deviceID, int printProfilingFreq /* = -1 */
) {
  NAMD_die(
      "CudaGlobalMasterServer requires to build NAMD with CUDA support.\n");
}
#endif // (defined(NAMD_CUDA) || defined(NAMD_HIP)) && defined(NODEGROUP_FORCE_REGISTER)