Merge branch release-2018

[alexxy/gromacs.git] / src / gromacs / mdlib / tpi.cpp

diff --git a/src/gromacs/mdlib/tpi.cpp b/src/gromacs/mdlib/tpi.cpp
index 5b15d1734e4b98a63db75ee7f9f31a1087dc2272..bf25ece3882365159419203bef1c6697d0f74c5f 100644 (file)
--- a/src/gromacs/mdlib/tpi.cpp
+++ b/src/gromacs/mdlib/tpi.cpp
@@ -127,7 +127,9 @@ namespace gmx
 {
 
 /*! \brief Do test particle insertion.
-    \copydoc integrator_t (FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
+    \copydoc integrator_t (FILE *fplog, t_commrec *cr,
+                           const gmx_multi_sim_t *,
+                           const gmx::MDLogger &mdlog,
                            int nfile, const t_filenm fnm[],
                            const gmx_output_env_t *oenv,
                            const MdrunOptions &mdrunOptions,
@@ -144,7 +146,9 @@ namespace gmx
                            gmx_membed_t gmx_unused *membed,
                            gmx_walltime_accounting_t walltime_accounting)
  */
-double do_tpi(FILE *fplog, t_commrec *cr, const gmx::MDLogger gmx_unused &mdlog,
+double do_tpi(FILE *fplog, t_commrec *cr,
+              const gmx_multisim_t *ms,
+              const gmx::MDLogger gmx_unused &mdlog,
               int nfile, const t_filenm fnm[],
               const gmx_output_env_t *oenv,
               const MdrunOptions &mdrunOptions,
@@ -660,12 +664,12 @@ double do_tpi(FILE *fplog, t_commrec *cr, const gmx::MDLogger gmx_unused &mdlog,
              * out of the box. */
             /* Make do_force do a single node force calculation */
             cr->nnodes = 1;
-            do_force(fplog, cr, inputrec,
+            do_force(fplog, cr, ms, inputrec,
                      step, nrnb, wcycle, top, &top_global->groups,
                      state_global->box, state_global->x, &state_global->hist,
                      f, force_vir, mdatoms, enerd, fcd,
                      state_global->lambda,
-                     nullptr, fr, nullptr, mu_tot, t, nullptr, FALSE,
+                     nullptr, fr, nullptr, mu_tot, t, nullptr,
                      GMX_FORCE_NONBONDED | GMX_FORCE_ENERGY |
                      (bNS ? GMX_FORCE_DYNAMICBOX | GMX_FORCE_NS : 0) |
                      (bStateChanged ? GMX_FORCE_STATECHANGED : 0),
@@ -715,7 +719,8 @@ double do_tpi(FILE *fplog, t_commrec *cr, const gmx::MDLogger gmx_unused &mdlog,
             }
             else
             {
-                embU      = exp(-beta*epot);
+                // Exponent argument is fine in SP range, but output can be in DP range
+                embU      = exp(static_cast<double>(-beta*epot));
                 sum_embU += embU;
                 /* Determine the weighted energy contributions of each energy group */
                 e                = 0;
@@ -846,7 +851,13 @@ double do_tpi(FILE *fplog, t_commrec *cr, const gmx::MDLogger gmx_unused &mdlog,
     {
         fprintf(fplog, "\n");
         fprintf(fplog, "  <V>  = %12.5e nm^3\n", V_all/frame);
-        fprintf(fplog, "  <mu> = %12.5e kJ/mol\n", -log(VembU_all/V_all)/beta);
+        const double mu = -log(VembU_all/V_all)/beta;
+        fprintf(fplog, "  <mu> = %12.5e kJ/mol\n", mu);
+
+        if (!std::isfinite(mu))
+        {
+            fprintf(fplog, "\nThe computed chemical potential is not finite - consider increasing the number of steps and/or the number of frames to insert into.\n");
+        }
     }
 
     /* Write the Boltzmann factor histogram */