Merge pull request #426 from GOMC-WSU/IntraMolExchange

Fix #425 use the correct recip energy for IntraMEMC move acceptance
GOMC-WSU · Feb 22, 2022 · ab60a57 · ab60a57
2 parents 563acda + 07f4f74
commit ab60a57
Show file tree

Hide file tree

Showing 3 changed files with 16 additions and 14 deletions.
diff --git a/src/moves/IntraMoleculeExchange1.h b/src/moves/IntraMoleculeExchange1.h
@@ -46,8 +46,7 @@ class IntraMoleculeExchange1 : public MoveBase
     numInCavB = 0;
     numSCavA = 0;
     numSCavB = 0;
-    recipDiffA = 0.0;
-    recipDiffB = 0.0;
+    recipDiff = 0.0;
     sourceBox = 0;
     volCav = 0.0;
 
@@ -127,7 +126,7 @@ class IntraMoleculeExchange1 : public MoveBase
 
   int exDiff, exchangeRatio;
   double volCav, lastAccept;
-  double W_recip, recipDiffA, recipDiffB, correctDiff;
+  double W_recip, recipDiff, correctDiff;
 
   XYZ centerA, centerB, cavity;
   XYZArray cavA, invCavA, cavB, invCavB;
@@ -492,16 +491,18 @@ inline void IntraMoleculeExchange1::CalcEn()
 {
   GOMC_EVENT_START(1, GomcProfileEvent::CALC_EN_INTRA_MEMC);
   W_recip = 1.0;
-  recipDiffA = 0.0, recipDiffB = 0.0;
+  recipDiff = 0.0;
   correctDiff = 0.0;
 
   if(!overlap) {
-    recipDiffA = calcEwald->MolExchangeReciprocal(newMolA, oldMolA, molIndexA, molIndexA, true);
-    recipDiffB = calcEwald->MolExchangeReciprocal(newMolB, oldMolB, molIndexB, molIndexB, false);
+    //MolExchangeReciprocal returns the total change in recip energy. It accumulates with each
+    //call, so we should use only the last of the two.
+    recipDiff = calcEwald->MolExchangeReciprocal(newMolA, oldMolA, molIndexA, molIndexA, true);
+    recipDiff = calcEwald->MolExchangeReciprocal(newMolB, oldMolB, molIndexB, molIndexB, false);
 
     //No need to contribute the self and correction energy since insertion
     //and deletion are rigid body
-    W_recip = exp(-1.0 * ffRef.beta * (recipDiffA + recipDiffB));
+    W_recip = exp(-ffRef.beta * recipDiff);
   }
   GOMC_EVENT_STOP(1, GomcProfileEvent::CALC_EN_INTRA_MEMC);
 }
@@ -588,7 +589,7 @@ inline void IntraMoleculeExchange1::Accept(const uint rejectState,
       }
 
       // Add Reciprocal energy
-      sysPotRef.boxEnergy[sourceBox].recip += recipDiffB;
+      sysPotRef.boxEnergy[sourceBox].recip += recipDiff;
 
       // Add correction energy
       sysPotRef.boxEnergy[sourceBox].correction += correctDiff;

diff --git a/src/moves/IntraMoleculeExchange2.h b/src/moves/IntraMoleculeExchange2.h
@@ -342,7 +342,7 @@ inline uint IntraMoleculeExchange2::Transform()
 
 inline void IntraMoleculeExchange2::CalcEn()
 {
-  // Updates recipDiffA and recipDiffB and updates sum new arrays at the same time
+  // Updates recipDiff and updates sum new arrays at the same time
   IntraMoleculeExchange1::CalcEn();
 }
 

diff --git a/src/moves/IntraMoleculeExchange3.h b/src/moves/IntraMoleculeExchange3.h
@@ -296,7 +296,7 @@ inline void IntraMoleculeExchange3::CalcEn()
 {
   GOMC_EVENT_START(1, GomcProfileEvent::CALC_EN_INTRA_MEMC);
   W_recip = 1.0;
-  recipDiffA = 0.0, recipDiffB = 0.0;
+  recipDiff = 0.0;
   correctDiff = 0.0;
   //No need to calculate the correction term for kindS since it is
   // inserted rigid body. We just need it for kindL
@@ -305,11 +305,12 @@ inline void IntraMoleculeExchange3::CalcEn()
       correctDiff += calcEwald->SwapCorrection(newMolB[n], molIndexB[n]);
       correctDiff -= calcEwald->SwapCorrection(oldMolB[n], molIndexB[n]);
     }
-    recipDiffA = calcEwald->MolExchangeReciprocal(newMolA, oldMolA, molIndexA, molIndexA, true);
-    recipDiffB = calcEwald->MolExchangeReciprocal(newMolB, oldMolB, molIndexB, molIndexB, false);
+    //MolExchangeReciprocal returns the total change in recip energy. It accumulates with each
+    //call, so we should use only the last of the two.
+    recipDiff = calcEwald->MolExchangeReciprocal(newMolA, oldMolA, molIndexA, molIndexA, true);
+    recipDiff = calcEwald->MolExchangeReciprocal(newMolB, oldMolB, molIndexB, molIndexB, false);
 
-    W_recip = exp(-1.0 * ffRef.beta * (recipDiffA + recipDiffB +
-                                       correctDiff));
+    W_recip = exp(-ffRef.beta * (recipDiff + correctDiff));
   }
   GOMC_EVENT_STOP(1, GomcProfileEvent::CALC_EN_INTRA_MEMC);
 }