* 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/fileio/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 "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/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 "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++;
}
}
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;
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;