Remove unnecessary config.h includes

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

diff --git a/src/gromacs/pulling/pull.c b/src/gromacs/pulling/pull.c
index 3ab4fa5d2fbb04cfb34fc20ac41b34e33c2ecaba..e84ad19714c9de45ac48b6138a26823ec83e0744 100644 (file)
--- a/src/gromacs/pulling/pull.c
+++ b/src/gromacs/pulling/pull.c
@@ -34,9 +34,7 @@
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
 
 #include <math.h>
 #include <stdio.h>
@@ -44,18 +42,17 @@
 #include <string.h>
 
 #include "gromacs/utility/futil.h"
-#include "index.h"
-#include "typedefs.h"
-#include "types/commrec.h"
-#include "network.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/legacyheaders/network.h"
 #include "pull.h"
-#include "names.h"
+#include "gromacs/legacyheaders/names.h"
 #include "gromacs/pbcutil/pbc.h"
-#include "mtop_util.h"
-#include "mdrun.h"
-#include "gmx_ga2la.h"
-#include "copyrite.h"
-#include "macros.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/legacyheaders/mdrun.h"
+#include "gromacs/legacyheaders/gmx_ga2la.h"
+#include "gromacs/legacyheaders/copyrite.h"
+#include "gromacs/legacyheaders/macros.h"
 
 #include "gromacs/fileio/filenm.h"
 #include "gromacs/fileio/gmxfio.h"
@@ -180,7 +177,7 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv
                         if (pull->dim[m])
                         {
                             sprintf(buf, "%d %s%c", c+1, "c", 'X'+m);
-                            setname[nsets] = strdup(buf);
+                            setname[nsets] = gmx_strdup(buf);
                             nsets++;
                         }
                     }
@@ -190,7 +187,7 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv
                     if (pull->dim[m])
                     {
                         sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
-                        setname[nsets] = strdup(buf);
+                        setname[nsets] = gmx_strdup(buf);
                         nsets++;
                     }
                 }
@@ -198,7 +195,7 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv
             else
             {
                 sprintf(buf, "%d", c+1);
-                setname[nsets] = strdup(buf);
+                setname[nsets] = gmx_strdup(buf);
                 nsets++;
             }
         }
@@ -502,6 +499,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;
@@ -629,6 +631,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc,
         for (c = 0; c < pull->ncoord; c++)
         {
             pcrd = &pull->coord[c];
+            ref  = pcrd->init + pcrd->rate*t;
 
             low_get_pull_coord_dr(pull, pcrd, pbc, t,
                                   rnew[pcrd->group[1]],
@@ -740,6 +743,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 +788,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 +947,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 +1261,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;