* 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 "network.h"
-#include "gromacs/fileio/filenm.h"
#include <string.h>
-#include "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 "partdec.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++;
}
}
static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md,
const dvec f_pull, int sign, rvec *f)
{
- int i, ii, m, start, end;
+ int i, ii, m;
double wmass, inv_wm;
- start = md->start;
- end = md->homenr + start;
-
inv_wm = pgrp->wscale*pgrp->invtm;
for (i = 0; i < pgrp->nat_loc; i++)
}
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;
for (g = 0; g < pull->ngroup; g++)
{
make_local_pull_group(ga2la, &pull->group[g],
- md->start, md->start+md->homenr);
+ 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,
- int start, int end,
int g, t_pull_group *pg, ivec pulldims,
gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
{
int i, ii, d, nfrozen, ndim;
real m, w, mbd;
double tmass, wmass, wwmass;
- gmx_bool bDomDec;
- gmx_ga2la_t ga2la = NULL;
gmx_groups_t *groups;
gmx_mtop_atomlookup_t alook;
t_atom *atom;
- bDomDec = (cr && DOMAINDECOMP(cr));
- if (bDomDec)
- {
- ga2la = cr->dd->ga2la;
- }
-
if (EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
{
/* There are no masses in the integrator.
{
ii = pg->ind[i];
gmx_mtop_atomnr_to_atom(alook, ii, &atom);
- if (cr && PAR(cr) && !bDomDec && ii >= start && ii < end)
- {
- pg->ind_loc[pg->nat_loc++] = ii;
- }
if (ir->opts.nFreeze)
{
for (d = 0; d < DIM; d++)
pull->bVirial = FALSE;
}
- if (cr && PARTDECOMP(cr))
- {
- pd_at_range(cr, &start, &end);
- }
pull->rbuf = NULL;
pull->dbuf = NULL;
pull->dbuf_cyl = NULL;
}
}
/* Set the indices */
- init_pull_group_index(fplog, cr, start, end, g, pgrp, pull->dim, mtop, ir, lambda);
+ init_pull_group_index(fplog, cr, g, pgrp, pull->dim, mtop, ir, lambda);
if (PULL_CYL(pull) && pgrp->invtm == 0)
{
gmx_fatal(FARGS, "Can not have frozen atoms in a cylinder pull group");
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;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff