* 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>
#include <stdlib.h>
-#include "gromacs/utility/futil.h"
-#include "index.h"
-#include "gromacs/fileio/gmxfio.h"
-#include "vec.h"
-#include "typedefs.h"
-#include "types/commrec.h"
-#include "network.h"
-#include "gromacs/fileio/filenm.h"
#include <string.h>
-#include "gromacs/utility/smalloc.h"
+
+#include "gromacs/utility/futil.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/legacyheaders/network.h"
#include "pull.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/pbcutil/pbc.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"
#include "gromacs/fileio/xvgr.h"
-#include "names.h"
-#include "pbc.h"
-#include "mtop_util.h"
-#include "mdrun.h"
-#include "gmx_ga2la.h"
-#include "copyrite.h"
-#include "macros.h"
-#include "vec.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/smalloc.h"
static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp)
{
if (pull->dim[m])
{
sprintf(buf, "%d %s%c", c+1, "c", 'X'+m);
- setname[nsets] = strdup(buf);
+ setname[nsets] = gmx_strdup(buf);
nsets++;
}
}
if (pull->dim[m])
{
sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
- setname[nsets] = strdup(buf);
+ setname[nsets] = gmx_strdup(buf);
nsets++;
}
}
else
{
sprintf(buf, "%d", c+1);
- setname[nsets] = strdup(buf);
+ setname[nsets] = gmx_strdup(buf);
nsets++;
}
}
}
if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
{
- gmx_fatal(FARGS, "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the box size (%f)",
+ gmx_fatal(FARGS, "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the box size (%f).\nYou might want to consider using \"pull-geometry = direction-periodic\" instead.\n",
pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2));
}
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;
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]],
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++)
{
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;
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,
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;