PME CUDA recalculate self-energy term 16/4516/1
authorJim Phillips <jim@ks.uiuc.edu>
Thu, 23 Aug 2018 22:30:56 +0000 (17:30 -0500)
committerJim Phillips <jim@ks.uiuc.edu>
Thu, 23 Aug 2018 22:30:56 +0000 (17:30 -0500)
Modify PME CUDA implementation to always recalculate the self-energy term.
This solves the incompatibility with REST2, in which PME CUDA was getting
the energy wrong after a charge rescaling.

Also greatly simplifies code by submitting self energy to reduction system
where it is calculated rather than separately reducing it to a pencil.

Change-Id: I05d8d959004b5e7ed8a274f41936a9400a4a7e1d

src/ComputePmeCUDA.C
src/ComputePmeCUDA.h
src/ComputePmeCUDAMgr.C
src/ComputePmeCUDAMgr.ci
src/ComputePmeCUDAMgr.h
src/PmeSolver.C
src/PmeSolver.ci
src/PmeSolver.h

index 3095411..351cd7a 100644 (file)
@@ -48,6 +48,7 @@ ComputePmeCUDA::~ComputePmeCUDA() {
        PatchMap::Object()->patch(patches[i].patchID)->unregisterForceDeposit(this, &patches[i].forceBox);
     }
   }
+  delete reduction;
   CmiDestroyLock(lock);
 }
 
@@ -65,7 +66,7 @@ void ComputePmeCUDA::initialize() {
   if (simParams->pairInteractionOn) NAMD_bug("ComputePmeCUDA::ComputePmeCUDA, pairInteractionOn not yet implemented");
 
   sendAtomsDone = false;
-  selfEnergyDone = false;
+  reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
   // basePriority = PME_PRIORITY;
   patchCounter = getNumPatches();
 
@@ -123,6 +124,8 @@ int ComputePmeCUDA::noWork() {
 
   if (patches[0].patch->flags.doFullElectrostatics) return 0;
 
+  reduction->submit();
+
   for (int i=0;i < getNumPatches();i++) {
     patches[i].positionBox->skip();
     patches[i].forceBox->skip();
@@ -167,6 +170,8 @@ void ComputePmeCUDA::sendAtoms() {
   double r3y = recip3.y;
   double r3z = recip3.z;
 
+  double selfEnergy = 0.;
+
   for (int i=0;i < getNumPatches();i++) {
     if (patches[i].pmeForceMsg != NULL)
       NAMD_bug("ComputePmeCUDA::sendAtoms, pmeForceMsg is not empty");
@@ -192,14 +197,7 @@ void ComputePmeCUDA::sendAtoms() {
 
     int numAtoms = patches[i].patch->getNumAtoms();
 
-    if (!selfEnergyDone) {
-      double selfEnergy = calcSelfEnergy(numAtoms, x);
-      // Send self-energy to one of the pencils, doesn't matter which one, the self energy is always
-      // delivered to (0,0) pencil
-      PmeSelfEnergyMsg *msg = new PmeSelfEnergyMsg();
-      msg->energy = selfEnergy;
-      computePmeCUDAMgrProxy[patches[i].homePencilNode].recvSelfEnergy(msg);
-    }
+    if ( doEnergy ) selfEnergy += calcSelfEnergy(numAtoms, x);
 
     // const Vector ucenter = patches[i].patch->lattice.unscale(patchMap->center(patches[i].patchID));
     // const BigReal recip11 = patches[i].patch->lattice.a_r().x;
@@ -276,9 +274,8 @@ void ComputePmeCUDA::sendAtoms() {
       patches[i].positionBox->close(&x);
   }
 
-  if (!selfEnergyDone) {
-    selfEnergyDone = true;
-  }
+  reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += selfEnergy;
+  reduction->submit();
 
 }
 
index 0d8cdd2..acb715e 100644 (file)
@@ -59,7 +59,7 @@ private:
 
   std::vector<PatchRecord> patches;
 
-  bool selfEnergyDone;
+  SubmitReduction *reduction;
 
   bool sendAtomsDone;
 
@@ -74,4 +74,4 @@ private:
 };
 #endif // NAMD_CUDA
 
-#endif // COMPUTEPMECUDA_H
\ No newline at end of file
+#endif // COMPUTEPMECUDA_H
index 4513f6f..ba8ad0c 100644 (file)
@@ -1702,25 +1702,6 @@ void ComputePmeCUDADevice::sendForcesToPatch(PmeForceMsg *forceMsg) {
     wdProxy[pe].enqueuePme(lmsg);
   }
 }
-
-//
-// Received self energy from ComputePmeCUDA computes
-// NOTE: This should only happen once at the start of the simulation since self energy is constant
-//
-void ComputePmeCUDAMgr::recvSelfEnergy(PmeSelfEnergyMsg *msg) {
-  // NOTE: msg is deleted in PmePencilXYZ / PmePencilX
-  switch(pmePencilType) {
-    case 1:
-    pmePencilZ(0,0,0).recvSelfEnergy(msg);
-    break;
-    case 2:
-    pmePencilZ(0,0,0).recvSelfEnergy(msg);
-    break;
-    case 3:
-    pmePencilXYZ[0].recvSelfEnergy(msg);
-    break;
-  }  
-}
 #endif // NAMD_CUDA
 
 #include "ComputePmeCUDAMgr.def.h"
index 866f457..75032a1 100644 (file)
@@ -58,7 +58,6 @@ module ComputePmeCUDAMgr {
     entry void recvPencils(CProxy_CudaPmePencilXYZ xyz);
     entry void recvPencils(CProxy_CudaPmePencilXY xy, CProxy_CudaPmePencilZ z);
     entry void recvPencils(CProxy_CudaPmePencilX x, CProxy_CudaPmePencilY y, CProxy_CudaPmePencilZ z);
-    entry void recvSelfEnergy(PmeSelfEnergyMsg *msg);
 
     entry void recvDevices(RecvDeviceMsg* msg);
     entry void recvAtomFiler(CProxy_PmeAtomFiler filer);
index e81cb2d..d1f4d90 100644 (file)
@@ -448,7 +448,6 @@ public:
   void recvPencils(CProxy_CudaPmePencilXYZ xyz);
   void recvPencils(CProxy_CudaPmePencilXY xy, CProxy_CudaPmePencilZ z);
   void recvPencils(CProxy_CudaPmePencilX x, CProxy_CudaPmePencilY y, CProxy_CudaPmePencilZ z);
-  void recvSelfEnergy(PmeSelfEnergyMsg *msg);
 
   void createDevicesAndAtomFiler();
   void recvDevices(RecvDeviceMsg* msg);
@@ -536,4 +535,4 @@ class ComputePmeCUDAMgr : public CBase_ComputePmeCUDAMgr {
 };
 #endif // NAMD_CUDA
 
-#endif // COMPUTEPMECUDAMGR_H
\ No newline at end of file
+#endif // COMPUTEPMECUDAMGR_H
index 16ac216..51c451e 100644 (file)
@@ -21,8 +21,6 @@ PmePencilXYZ::PmePencilXYZ() {
   fftCompute = NULL;
   pmeKSpaceCompute = NULL;
   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
-  numSelfEnergyRecv = 0;
-  selfEnergy = 0.0;
   doEnergy = false;
   doVirial = false;
 }
@@ -76,7 +74,7 @@ void PmePencilXYZ::submitReductions() {
   double virial[9];
   double energy = pmeKSpaceCompute->getEnergy();
   pmeKSpaceCompute->getVirial(virial);
-  reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += energy + selfEnergy;
+  reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += energy;
   reduction->item(REDUCTION_VIRIAL_SLOW_XX) += virial[0];
   reduction->item(REDUCTION_VIRIAL_SLOW_XY) += virial[1];
   reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += virial[2];
@@ -372,8 +370,6 @@ PmePencilZ::PmePencilZ() {
   pmeTranspose = NULL;
   pmeKSpaceCompute = NULL;
   reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
-  numSelfEnergyRecv = 0;
-  selfEnergy = 0.0;
   doEnergy = false;
   doVirial = false;
   numStrayAtoms = 0;
@@ -436,7 +432,7 @@ void PmePencilZ::submitReductions() {
   double energy = pmeKSpaceCompute->getEnergy();
   // fprintf(stderr, "PmePencilZ::submitReductions(), numStrayAtoms %d\n", numStrayAtoms);
   pmeKSpaceCompute->getVirial(virial);
-  reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += energy + selfEnergy;
+  reduction->item(REDUCTION_ELECT_ENERGY_SLOW) += energy;
   reduction->item(REDUCTION_VIRIAL_SLOW_XX) += virial[0];
   reduction->item(REDUCTION_VIRIAL_SLOW_XY) += virial[1];
   reduction->item(REDUCTION_VIRIAL_SLOW_XZ) += virial[2];
index fbdd1c3..390dcb8 100644 (file)
@@ -1,7 +1,5 @@
 module PmeSolver {
 
-  message PmeSelfEnergyMsg;
-
   message PmeStartMsg;
 
   message PmeRunMsg;
@@ -22,24 +20,12 @@ module PmeSolver {
 
   array[1D] PmePencilXYZ {
     entry PmePencilXYZ();
-    entry void recvSelfEnergy(PmeSelfEnergyMsg *msg) {
-      serial {
-        selfEnergy += msg->energy;
-        delete msg;
-        numSelfEnergyRecv++;
-        if (numSelfEnergyRecv == PatchMap::Object()->numPatches()) {
-          selfEnergyReady();
-        }
-      }
-    };
-    entry void selfEnergyReady();
     entry void chargeGridReady(PmeRunMsg *msg);
     entry void skip();
     entry void start(PmeStartMsg *pmeStartMsg) {
       serial "initFFT" {
         initFFT(pmeStartMsg); 
         delete pmeStartMsg;
-        firstIter = true;
       }
       while (true) {
         when chargeGridReady(PmeRunMsg *msg) serial "forwardFFT" {
@@ -52,12 +38,7 @@ module PmeSolver {
         }
         serial "forwardDone" { forwardDone(); }
         serial "backwardFFT" { backwardFFT(); }
-        if (firstIter) {
-          // On the first iteration, wait until self energy has been received
-          when selfEnergyReady() serial "backwardDoneFirst" { firstIter = false; backwardDone(); }
-        } else {
-          serial "backwardDone" { backwardDone(); }
-        }
+        serial "backwardDone" { backwardDone(); }
       }
     };
   };
@@ -180,19 +161,6 @@ module PmeSolver {
 
   array[3D] PmePencilZ {
     entry PmePencilZ();
-    entry void recvSelfEnergy(PmeSelfEnergyMsg *msg) {
-      serial {
-        if (thisIndex.x != 0 || thisIndex.y != 0 || thisIndex.z != 0)
-          NAMD_bug("PmePencilZ::recvSelfEnergy can only be called on (0,0,0) element");
-        selfEnergy += msg->energy;
-        delete msg;
-        numSelfEnergyRecv++;
-        if (numSelfEnergyRecv == PatchMap::Object()->numPatches()) {
-          selfEnergyReady();
-        }
-      }
-    };
-    entry void selfEnergyReady();
     entry void recvBlock(PmeBlockMsg *msg);
     entry void skip();
     entry void start(PmeStartMsg *pmeStartMsg) {
@@ -200,10 +168,6 @@ module PmeSolver {
         initFFT(pmeStartMsg); 
         delete pmeStartMsg;
         start();
-        if (thisIndex.x == 0 && thisIndex.y == 0)
-          firstIter = true;
-        else
-          firstIter = false;
       }
       while (true) {
         for (imsg=0;imsg < pmeGrid.zBlocks;++imsg) {
@@ -220,12 +184,7 @@ module PmeSolver {
           forwardDone();
         }
         serial "backwardFFT" { backwardFFT(); }
-        if (firstIter) {
-          // On the first iteration on the (0,0) element, wait until self energy has been received
-          when selfEnergyReady() serial "backwardDoneFirst" { firstIter = false; backwardDone(); }
-        } else {
-          serial "backwardDone" { backwardDone(); }
-        }
+        serial "backwardDone" { backwardDone(); }
       }
     };
   };
index d137cbe..cf0d6e5 100644 (file)
@@ -112,11 +112,6 @@ public:
   Lattice lattice;
 };
 
-class PmeSelfEnergyMsg : public CMessage_PmeSelfEnergyMsg {
-public:
-  double energy;
-};
-
 class PmeDoneMsg : public CMessage_PmeDoneMsg {
 public:
   PmeDoneMsg(int i, int j) : i(i), j(j) {}
@@ -155,9 +150,6 @@ private:
   void forwardDone();
   void initFFT(PmeStartMsg *msg);
 
-  int numSelfEnergyRecv;
-  bool firstIter;
-  double selfEnergy;
   SubmitReduction* reduction;
 
 };
@@ -272,9 +264,6 @@ private:
   virtual void recvDataFromY(PmeBlockMsg *msg);
   virtual void start();
 
-  int numSelfEnergyRecv;
-  bool firstIter;
-  double selfEnergy;
   SubmitReduction* reduction;
 
 };