ComputeLonepairsCUDA.C

Go to the documentation of this file.

#include "ComputeLonepairsCUDA.h"
#include "CudaUtils.h"
#include "ProcessorPrivate.h"
#include "AtomMap.h"
#include "HomePatch.h"
#include "Node.h"
#include "Molecule.h"
#include "DeviceCUDA.h"
#include "ComputeLonepairsCUDAKernel.h"

#include <string>
#include <cmath>

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

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

ComputeLonepairsCUDA::ComputeLonepairsCUDA():
  d_lpcolinear_list(nullptr), lpcolinear_list_size(0),
  d_lprelative_list(nullptr), lprelative_list_size(0),
  d_lpbisector_list(nullptr), lpbisector_list_size(0) {}

ComputeLonepairsCUDA::~ComputeLonepairsCUDA() {
  deallocate_device(&d_lpcolinear_list);
  lpcolinear_list_size = 0;
  deallocate_device(&d_lprelative_list);
  lprelative_list_size = 0;
  deallocate_device(&d_lpbisector_list);
  lpbisector_list_size = 0;
}

int globalAtomIDToSOAID(
  std::vector<AtomMap*> &atomMapsList,
  const std::vector<CudaLocalRecord> &localRecords,
  const int* globalToLocalID, AtomID aid, bool errorOnNotFound) {
  LocalID lid{-1, -1};
  for (int j = 0; j < atomMapsList.size(); ++j) {
    if (atomMapsList[j] == NULL) continue;
    lid = atomMapsList[j]->localID(aid);
    if (lid.pid != -1) break;
  }
  if (errorOnNotFound && (lid.pid == -1)) {
    const std::string error = "Atom " + std::to_string(aid) + " not found in master PE " + std::to_string(deviceCUDA->getMasterPe());
    NAMD_bug(error.c_str());
    return -1;
  }
  const int soaPid = globalToLocalID[lid.pid];
  const int soaIndex = localRecords[soaPid].bufferOffset + lid.index;
  return soaIndex;
}

void ComputeLonepairsCUDA::updateAtoms(
  std::vector<HomePatch*> patchList,
  std::vector<AtomMap*> &atomMapsList,
  const std::vector<CudaLocalRecord> &localRecords,
  const int* h_globalToLocalID, cudaStream_t stream) {
  std::vector<LonepairColinear> lpcolinear_list;
  std::vector<LonepairRelative> lprelative_list;
  std::vector<LonepairBisector> lpbisector_list;
  for (size_t i = 0; i < patchList.size(); ++i) {
    HomePatch *p = patchList[i];
    const int numAtoms = p->getNumAtoms();
    const auto& atomList = p->getAtomList();
    for (int j = 0; j < numAtoms; ++j) {
      const auto mass = atomList[j].mass;
      // Iterate over all atoms with mass < 0.01 to find lone pair IDs
      if (mass < 0.01) {
        const AtomID aid = atomList[j].id;
        // Find the LP host
        Lphost *lph_ptr = Node::Object()->molecule->get_lphost(aid);
        if (lph_ptr == NULL) {
          const std::string error = "rebuildLonePairHostList: no Lphost exists for LP " + std::to_string(aid);
          NAMD_die(error.c_str());
        }
        Lphost lph = *(lph_ptr);
        if (lph.numhosts == 2) { // colinear
          lpcolinear_list.push_back(LonepairColinear{
            globalAtomIDToSOAID(atomMapsList, localRecords, h_globalToLocalID, lph.atom1, true),
            globalAtomIDToSOAID(atomMapsList, localRecords, h_globalToLocalID, lph.atom2, true),
            globalAtomIDToSOAID(atomMapsList, localRecords, h_globalToLocalID, lph.atom3, true),
            lph.distance, lph.angle});
        } else {
          const double distance = lph.distance;
          // Pre-calculate the sine and cosine
          if (distance < 0) { // bisector
            const double dcosa = -distance * std::cos(lph.angle);
            const double dsina = -distance * std::sin(lph.angle);
            const double dsinasint = dsina * std::sin(lph.dihedral);
            const double dsinacost = dsina * std::cos(lph.dihedral);
            lpbisector_list.push_back(LonepairBisector{
              globalAtomIDToSOAID(atomMapsList, localRecords, h_globalToLocalID, lph.atom1, true),
              globalAtomIDToSOAID(atomMapsList, localRecords, h_globalToLocalID, lph.atom2, true),
              globalAtomIDToSOAID(atomMapsList, localRecords, h_globalToLocalID, lph.atom3, true),
              globalAtomIDToSOAID(atomMapsList, localRecords, h_globalToLocalID, lph.atom4, true),
              dcosa, dsinacost, dsinasint});
          } else { // relative
            const double dcosa = distance * std::cos(lph.angle);
            const double dsina = distance * std::sin(lph.angle);
            const double dsinasint = dsina * std::sin(lph.dihedral);
            const double dsinacost = dsina * std::cos(lph.dihedral);
            lprelative_list.push_back(LonepairRelative{
              globalAtomIDToSOAID(atomMapsList, localRecords, h_globalToLocalID, lph.atom1, true),
              globalAtomIDToSOAID(atomMapsList, localRecords, h_globalToLocalID, lph.atom2, true),
              globalAtomIDToSOAID(atomMapsList, localRecords, h_globalToLocalID, lph.atom3, true),
              globalAtomIDToSOAID(atomMapsList, localRecords, h_globalToLocalID, lph.atom4, true),
              dcosa, dsinacost, dsinasint});
          }
        }
      }
    }
  }
  lpcolinear_list_size = lpcolinear_list.size();
  if (!lpcolinear_list.empty()) {
    deallocate_device_async(&d_lpcolinear_list, stream);
    allocate_device_async(&d_lpcolinear_list, lpcolinear_list_size, stream);
    copy_HtoD(lpcolinear_list.data(), d_lpcolinear_list, lpcolinear_list_size, stream);
  }
  lprelative_list_size = lprelative_list.size();
  if (!lprelative_list.empty()) {
    deallocate_device_async(&d_lprelative_list, stream);
    allocate_device_async(&d_lprelative_list, lprelative_list_size, stream);
    copy_HtoD(lprelative_list.data(), d_lprelative_list, lprelative_list_size, stream);
  }
  lpbisector_list_size = lpbisector_list.size();
  if (!lpbisector_list.empty()) {
    deallocate_device_async(&d_lpbisector_list, stream);
    allocate_device_async(&d_lpbisector_list, lpbisector_list_size, stream);
    copy_HtoD(lpbisector_list.data(), d_lpbisector_list, lpbisector_list_size, stream);
  }
  cudaCheck(cudaStreamSynchronize(stream));
}

void ComputeLonepairsCUDA::reposition(
  double*  d_pos_x,
  double*  d_pos_y,
  double*  d_pos_z,
  cudaStream_t stream) const {
  if (lpcolinear_list_size > 0) {
    repositionColinear(d_pos_x, d_pos_y, d_pos_z, d_lpcolinear_list, lpcolinear_list_size, stream);
  }
  if (lpbisector_list_size > 0) {
    repositionBisector(d_pos_x, d_pos_y, d_pos_z, d_lpbisector_list, lpbisector_list_size, stream);
  }
  if (lprelative_list_size > 0) {
    repositionRelative(d_pos_x, d_pos_y, d_pos_z, d_lprelative_list, lprelative_list_size, stream);
  }
}

void ComputeLonepairsCUDA::redistributeForce(
  double*        d_f_normal_x,
  double*        d_f_normal_y,
  double*        d_f_normal_z,
  double*        d_f_nbond_x,
  double*        d_f_nbond_y,
  double*        d_f_nbond_z,
  double*        d_f_slow_x,
  double*        d_f_slow_y,
  double*        d_f_slow_z,
  cudaTensor*    d_virial_normal,
  cudaTensor*    d_virial_nbond,
  cudaTensor*    d_virial_slow,
  const double*  d_pos_x,
  const double*  d_pos_y,
  const double*  d_pos_z,
  const int      maxForceNumber,
  const int      doVirial,
  cudaStream_t stream) const {
  if (lpcolinear_list_size > 0) {
    redistributeForceColinear(
      d_f_normal_x, d_f_normal_y, d_f_normal_z,
      d_f_nbond_x, d_f_nbond_y, d_f_nbond_z,
      d_f_slow_x, d_f_slow_y, d_f_slow_z,
      d_virial_normal, d_virial_nbond, d_virial_slow,
      d_pos_x, d_pos_y, d_pos_z, maxForceNumber,
      doVirial, d_lpcolinear_list, lpcolinear_list_size,
      stream);
  }
  if (lpbisector_list_size > 0) {
    redistributeForceBisector(
      d_f_normal_x, d_f_normal_y, d_f_normal_z,
      d_f_nbond_x, d_f_nbond_y, d_f_nbond_z,
      d_f_slow_x, d_f_slow_y, d_f_slow_z,
      d_virial_normal, d_virial_nbond, d_virial_slow,
      d_pos_x, d_pos_y, d_pos_z, maxForceNumber,
      doVirial, d_lpbisector_list, lpbisector_list_size,
      stream);
  }
  if (lprelative_list_size > 0) {
    redistributeForceRelative(
      d_f_normal_x, d_f_normal_y, d_f_normal_z,
      d_f_nbond_x, d_f_nbond_y, d_f_nbond_z,
      d_f_slow_x, d_f_slow_y, d_f_slow_z,
      d_virial_normal, d_virial_nbond, d_virial_slow,
      d_pos_x, d_pos_y, d_pos_z, maxForceNumber,
      doVirial, d_lprelative_list, lprelative_list_size,
      stream);
  }
}

#endif // NAMD_CUDA