fixed issues with FEP soft-core and cut-off's

Now the cut-off is applied to the PME mesh correction for perturbed
pairs with soft-core. Now the cut-off is applied to VdW with soft-core.
Fixes possible NaN in free energy kernel with r close to 0.
Replaced a confusing grompp warning with PME and soft-core by a note.
Removed an invalid grompp warning about soft-core and twin-range.
Fixes #1146

Change-Id: I79a06b20158df2bce575808b6e31e660163bd307

diff --git a/src/gmxlib/nonbonded/nb_free_energy.c b/src/gmxlib/nonbonded/nb_free_energy.c
index 8fce26eff4169b6cac6ac07eb7c3495655c6a5f6..61277a40d769ca8ef5865fa117464a5e2440eab1 100644 (file)
--- a/src/gmxlib/nonbonded/nb_free_energy.c
+++ b/src/gmxlib/nonbonded/nb_free_energy.c
@@ -108,7 +108,7 @@ gmx_nb_free_energy_kernel(t_nblist *                nlist,
     real *        Vv;
     real *        Vc;
     gmx_bool      bDoForces;
-    real          rcoulomb, rvdw, factor_coul, factor_vdw, sh_invrc6;
+    real          rcoulomb, rvdw, sh_invrc6;
     gmx_bool      bExactElecCutoff, bExactVdwCutoff;
     real          rcutoff, rcutoff2, rswitch, d, d2, swV3, swV4, swV5, swF2, swF3, swF4, sw, dsw, rinvcorr;
 
@@ -361,9 +361,6 @@ gmx_nb_free_energy_kernel(t_nblist *                nlist,
                     rinvV          = pow(rpinvV, 1.0/sc_r_power);
                     rV             = 1.0/rinvV;
 
-                    factor_coul    = (rC <= rcoulomb) ? 1.0 : 0.0;
-                    factor_vdw     = (rV <= rvdw)     ? 1.0 : 0.0;
-
                     if (do_tab)
                     {
                         rtC        = rC*tabscale;
@@ -379,7 +376,14 @@ gmx_nb_free_energy_kernel(t_nblist *                nlist,
                         n1V        = tab_elemsize*n0;
                     }
 
-                    if (qq[i] != 0)
+                    /* With Ewald and soft-core we should put the cut-off on r,
+                     * not on the soft-cored rC, as the real-space and
+                     * reciprocal space contributions should (almost) cancel.
+                     */
+                    if (qq[i] != 0 &&
+                        !(bExactElecCutoff &&
+                          ((icoul != GMX_NBKERNEL_ELEC_EWALD && rC >= rcoulomb) ||
+                           (icoul == GMX_NBKERNEL_ELEC_EWALD && r >= rcoulomb))))
                     {
                         switch (icoul)
                         {
@@ -427,15 +431,10 @@ gmx_nb_free_energy_kernel(t_nblist *                nlist,
                             Vcoul[i]        *= sw;
                             FscalC[i]        = FscalC[i]*sw + Vcoul[i]*dsw;
                         }
-
-                        if (bExactElecCutoff)
-                        {
-                            FscalC[i]        = (rC < rcoulomb) ? FscalC[i] : 0.0;
-                            Vcoul[i]         = (rC < rcoulomb) ? Vcoul[i] : 0.0;
-                        }
                     }
 
-                    if ((c6[i] != 0) || (c12[i] != 0))
+                    if ((c6[i] != 0 || c12[i] != 0) &&
+                        !(bExactVdwCutoff && rV >= rvdw))
                     {
                         switch (ivdw)
                         {
@@ -519,13 +518,21 @@ gmx_nb_free_energy_kernel(t_nblist *                nlist,
 
             Fscal = 0;
 
-            if (icoul == GMX_NBKERNEL_ELEC_EWALD)
+            if (icoul == GMX_NBKERNEL_ELEC_EWALD &&
+                !(bExactElecCutoff && r >= rcoulomb))
             {
                 /* because we compute the softcore normally,
                    we have to remove the ewald short range portion. Done outside of
                    the states loop because this part doesn't depend on the scaled R */
 
-                if (r != 0)
+#ifdef GMX_DOUBLE
+                /* Relative accuracy at R_ERF_R_INACC of 3e-10 */
+#define         R_ERF_R_INACC 0.006
+#else
+                /* Relative accuracy at R_ERF_R_INACC of 2e-5 */
+#define         R_ERF_R_INACC 0.1
+#endif
+                if (ewc*r > R_ERF_R_INACC)
                 {
                     VV    = gmx_erf(ewc*r)*rinv;
                     FF    = rinv*rinv*(VV - ewc*M_2_SQRTPI*exp(-ewc*ewc*rsq));
@@ -533,7 +540,7 @@ gmx_nb_free_energy_kernel(t_nblist *                nlist,
                 else
                 {
                     VV    = ewc*M_2_SQRTPI;
-                    FF    = 0;
+                    FF    = ewc*ewc*ewc*M_2_SQRTPI*(2.0/3.0 - 0.4*ewc*ewc*rsq);
                 }
 
                 for (i = 0; i < NSTATES; i++)

diff --git a/src/kernel/md.c b/src/kernel/md.c
index 631154fe87b3b2b8a182111a3a513f61976de2a2..886f9cbcc6c06db2450e6133e9ef0740d2f4810d 100644 (file)
--- a/src/kernel/md.c
+++ b/src/kernel/md.c
@@ -535,10 +535,13 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                                         repl_ex_nst, repl_ex_nex, repl_ex_seed);
     }
 
-    /* PME tuning is only supported with GPUs or PME nodes and not with rerun */
+    /* PME tuning is only supported with GPUs or PME nodes and not with rerun.
+     * With perturbed charges with soft-core we should not change the cut-off.
+     */
     if ((Flags & MD_TUNEPME) &&
         EEL_PME(fr->eeltype) &&
         ( (fr->cutoff_scheme == ecutsVERLET && fr->nbv->bUseGPU) || !(cr->duty & DUTY_PME)) &&
+        !(ir->efep != efepNO && mdatoms->nChargePerturbed > 0 && ir->fepvals->bScCoul) &&
         !bRerunMD)
     {
         pme_loadbal_init(&pme_loadbal, ir, state->box, fr->ic, fr->pmedata);

diff --git a/src/kernel/readir.c b/src/kernel/readir.c
index eaaa2583b5145917236c3cd9ff591bff10ffa718..e984b2e58ad683b7e4af6a815bdb8fa581d05e5d 100644 (file)
--- a/src/kernel/readir.c
+++ b/src/kernel/readir.c
@@ -580,14 +580,6 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
                 (int)fep->sc_r_power);
         CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
 
-        /* check validity of options */
-        if (fep->n_lambda > 0 && ir->rlist < max(ir->rvdw, ir->rcoulomb))
-        {
-            sprintf(warn_buf,
-                    "For foreign lambda free energy differences it is assumed that the soft-core interactions have no effect beyond the neighborlist cut-off");
-            warning(wi, warn_buf);
-        }
-
         sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
         CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) ||  (fep->init_lambda > 0)));
 
@@ -668,8 +660,16 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
 
         if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
         {
-            sprintf(warn_buf, "With coulomb soft core, the reciprocal space calculation will not necessarily cancel.  It may be necessary to decrease the reciprocal space energy, and increase the cutoff radius to get sufficiently close matches to energies with free energy turned off.");
-            warning(wi, warn_buf);
+            real sigma, lambda, r_sc;
+
+            sigma  = 0.34;
+            /* Maximum estimate for A and B charges equal with lambda power 1 */
+            lambda = 0.5;
+            r_sc   = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
+            sprintf(warn_buf, "With PME there is a minor soft core effect present at the cut-off, proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on energy conservation, but usually other effects dominate. With a common sigma value of %g nm the fraction of the particle-particle potential at the cut-off at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
+                    fep->sc_r_power,
+                    sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
+            warning_note(wi, warn_buf);
         }
 
         /*  Free Energy Checks -- In an ideal world, slow growth and FEP would