Merge release-5-0 into master

[alexxy/gromacs.git] / src / gromacs / pulling / pull.c

diff --git a/src/gromacs/pulling/pull.c b/src/gromacs/pulling/pull.c
index bfadaa93e0b9c7e20cabd3ee82513b97f3b59e5a..1e084858bbb7a117025d75f25ccf1eefe39c07dc 100644 (file)
--- a/src/gromacs/pulling/pull.c
+++ b/src/gromacs/pulling/pull.c
@@ -501,6 +501,11 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
                 r_ij[c][m] = a*pull->coord[c].vec[m];
             }
         }
+
+        if (dnorm2(r_ij[c]) == 0)
+        {
+            gmx_fatal(FARGS, "Distance for pull coordinate %d is zero with constraint pulling, which is not allowed.", c + 1);
+        }
     }
 
     bConverged_all = FALSE;
@@ -740,6 +745,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
             double f_invr;
 
             /* Add the pull contribution to the virial */
+            /* We have already checked above that r_ij[c] != 0 */
             f_invr = pcrd->f_scal/dnorm(r_ij[c]);
 
             for (j = 0; j < DIM; j++)
@@ -784,7 +790,19 @@ static void do_pull_pot(int ePull,
         {
             case epullgDIST:
                 ndr   = dnorm(pcrd->dr);
-                invdr = 1/ndr;
+                if (ndr > 0)
+                {
+                    invdr = 1/ndr;
+                }
+                else
+                {
+                    /* With an harmonic umbrella, the force is 0 at r=0,
+                     * so we can set invdr to any value.
+                     * With a constant force, the force at r=0 is not defined,
+                     * so we zero it (this is anyhow a very rare event).
+                     */
+                    invdr = 0;
+                }
                 if (ePull == epullUMBRELLA)
                 {
                     pcrd->f_scal  =       -k*dev;
@@ -931,6 +949,9 @@ void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md)
         make_local_pull_group(ga2la, &pull->group[g],
                               0, md->homenr);
     }
+
+    /* Since the PBC of atoms might have changed, we need to update the PBC */
+    pull->bSetPBCatoms = TRUE;
 }
 
 static void init_pull_group_index(FILE *fplog, t_commrec *cr,
@@ -1242,6 +1263,9 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[],
         snew(pull->dyna, pull->ncoord);
     }
 
+    /* We still need to initialize the PBC reference coordinates */
+    pull->bSetPBCatoms = TRUE;
+
     /* Only do I/O when we are doing dynamics and if we are the MASTER */
     pull->out_x = NULL;
     pull->out_f = NULL;