From 32eedeb6b3ad77cb016d0179d818c4d80cbb972c Mon Sep 17 00:00:00 2001
From: Michele Invernizzi <inve.michele@gmail.com>
Date: Wed, 8 Jul 2020 21:57:33 +0200
Subject: [PATCH] fixed exponential overflow of rct in case of non-tempered
 metad

---
 src/bias/MetaD.cpp | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/bias/MetaD.cpp b/src/bias/MetaD.cpp
index e5f24d9fb..ef0262ba3 100644
--- a/src/bias/MetaD.cpp
+++ b/src/bias/MetaD.cpp
@@ -1946,21 +1946,21 @@ void MetaD::computeReweightingFactor()
     minusBetaF=1;
     minusBetaFplusV=0;
   }
-  const double big_number=minusBetaF*BiasGrid_->getMaxValue(); //to avoid exp overflow
+  max_bias_=BiasGrid_->getMaxValue(); //to avoid exp overflow
 
   const unsigned rank=comm.Get_rank();
   const unsigned stride=comm.Get_size();
   for (Grid::index_t t=rank; t<BiasGrid_->getSize(); t+=stride) {
     const double val=BiasGrid_->getValue(t);
-    Z_0+=std::exp(minusBetaF*val-big_number);
-    Z_V+=std::exp(minusBetaFplusV*val-big_number);
+    Z_0+=std::exp(minusBetaF*(val-max_bias_));
+    Z_V+=std::exp(minusBetaFplusV*(val-max_bias_));
   }
   if (stride>1) {
     comm.Sum(Z_0);
     comm.Sum(Z_V);
   }
 
-  reweight_factor_=kbt_*std::log(Z_0/Z_V);
+  reweight_factor_=kbt_*std::log(Z_0/Z_V)+max_bias_;
   getPntrToComponent("rct")->set(reweight_factor_);
 }
 
-- 
GitLab