From 4da8ff3d219c6b63f8ad0b34ff97ca3f3553318e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 7 Jun 2020 15:47:46 -0400 Subject: [PATCH] avoid division by zero --- src/USER-QTB/fix_qbmsst.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/USER-QTB/fix_qbmsst.cpp b/src/USER-QTB/fix_qbmsst.cpp index 34350c9b05..bb396201cb 100644 --- a/src/USER-QTB/fix_qbmsst.cpp +++ b/src/USER-QTB/fix_qbmsst.cpp @@ -520,8 +520,8 @@ void FixQBMSST::initial_integrate(int /*vflag*/) // decide if the qtb temperature need to be updated or not if (counter_l == 0) { t_current -= dtv*fric_coef*eta*beta*(old_eavg-e0)/(3*ntotal*boltz); + if (t_current <= 0.0) break; old_eavg = 0;//clear old energy average - if (t_current < 0.0) t_current=0; // load omega_H with calculated spectrum at a specific temperature (corrected spectrum), omega_H is the Fourier transformation of time_H for (int k = 0; k < 2*N_f; k++) { @@ -529,7 +529,7 @@ void FixQBMSST::initial_integrate(int /*vflag*/) if(k == N_f) { omega_H[k]=sqrt(force->boltz * t_current); } else { - double energy_k= force->hplanck * fabs(f_k); + double energy_k= force->hplanck * fabs(f_k); omega_H[k]=sqrt( energy_k * (0.5+1.0/( exp(energy_k/(force->boltz * t_current)) - 1.0 )) ); omega_H[k]*=alpha*sin((k-N_f)*MY_PI/(2*alpha*N_f))/sin((k-N_f)*MY_PI/(2*N_f)); }