diff --git a/src/moves/IntraMoleculeExchange1.h b/src/moves/IntraMoleculeExchange1.h index a67add807..3370dae5f 100644 --- 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 index 123728bb6..d0394ff2b 100644 --- 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 index 42e1873b5..3eb6e3525 100644 --- 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); }