\label{sec:pull}
The pull code applies forces or constraints between the centers
of mass of one or more pairs of groups of atoms.
-There is one reference group and one or more other pull groups.
-Instead of a reference group, one can also use absolute reference
-point in space.
-The most common situation consists of a reference group
-and one pull group. In this case, the two groups are treated
-equivalently.
+Each pull reaction coordinate is called a ``coordinate'' and it operates
+on two pull groups. A pull group can be part of one or more pull
+coordinates. Furthermore, a coordinate can also operate on a single group
+and an absolute reference position in space.
The distance between a pair of groups can be determined
in 1, 2 or 3 dimensions, or can be along a user-defined vector.
The reference distance can be constant or can change linearly with time.
\eeq
\subsubsection{Limitations}
-There is one important limitation:
+There is one theoretical limitation:
strictly speaking, constraint forces can only be calculated between
groups that are not connected by constraints to the rest of the system.
If a group contains part of a molecule of which the bond lengths
Ewald summation can only be used with <b>nwall=2</b>, where one
should use <b><A HREF="#ewald">ewald-geometry</A><tt>=3dc</tt></b>.
The empty layer in the box serves to decrease the unphysical Coulomb
-interaction between periodic images.
+interaction between periodic images.</dd>
</dl>
<A NAME="pull"><br>
Note that the radii should be smaller than half the box size.
For tilted cylinders they should be even smaller than half the box size
since the distance of an atom in the reference group
-from the COM of the pull group has both a radial and an axial component.
-<dt><b>position</b></dt>
-<dd>Pull to the position of the reference group plus
-<b>pull-init</b> + time*<b>pull-rate</b>*<b>pull-vec</b>.</dd>
+from the COM of the pull group has both a radial and an axial component.</dd>
</dl></dd>
<dt><b>pull-dim: (Y Y Y)</b></dt>
<dd>the distance components to be used with geometry <b>distance</b>
<dt><b>yes</b></dt>
<dd>add the COM distance of the starting conformation to <b>pull-init</b></dd>
</dl>
+<dt><b>pull-print-reference: (10)</b></dt>
+<dd><dl compact>
+<dt><b>no</b></dt>
+<dd>do not print the COM of the first group in each pull coordinate</dd>
+<dt><b>yes</b></dt>
+<dd>print the COM of the first group in each pull coordinate</dd>
+</dl>
<dt><b>pull-nstxout: (10)</b></dt>
<dd>frequency for writing out the COMs of all the pull group</dd>
<dt><b>pull-nstfout: (1)</b></dt>
<dd>frequency for writing out the force of all the pulled group</dd>
<dt><b>pull-ngroups: (1)</b></dt>
-<dd>The number of pull groups, not including the reference group.
-If there is only one group, there is no difference in treatment
-of the reference and pulled group (except with the cylinder geometry).
-Below only the pull options for the reference group (ending on 0)
-and the first group (ending on 1) are given,
-further groups work analogously, but with the number 1 replaced
-by the group number.</dd>
-<dt><b>pull-group0: </b></dt>
-<dd>The name of the reference group. When this is empty an absolute reference
-of (0,0,0) is used. With an absolute reference the system is no longer
-translation invariant and one should think about what to do with
-the <A HREF="#run">center of mass motion</A>.</dd>
-<dt><b>pull-weights0: </b></dt>
-<dd>see <b>pull-weights1</b></dd>
-<dt><b>pull-pbcatom0: (0)</b></dt>
-<dd>see <b>pull-pbcatom1</b></dd>
-<dt><b>pull-group1: </b></dt>
-<dd>The name of the pull group.</dd>
-<dt><b>pull-weights1: </b></dt>
+<dd>The number of pull groups, not including the absolute reference group,
+when used. Pull groups can be reused in multiple pull coordinates.
+Below only the pull options for group 1 are given, further groups simply
+increase the group index number.</dd>
+<dt><b>pull-ncoords: (1)</b></dt>
+<dd>The number of pull coordinates. Below only the pull options for
+coordinate 1 are given, further coordinates simply increase the coordinate
+index number.</dd>
+
+<dt><b>pull-group1-name: </b></dt>
+<dd>The name of the pull group, is looked up in the index file
+or in the default groups to obtain the atoms involved.</dd>
+<dt><b>pull-group1-weights: </b></dt>
<dd>Optional relative weights which are multiplied with the masses of the atoms
to give the total weight for the COM. The number should be 0, meaning all 1,
or the number of atoms in the pull group.</dd>
-<dt><b>pull-pbcatom1: (0)</b></dt>
+<dt><b>pull-group1-pbcatom: (0)</b></dt>
<dd>The reference atom for the treatment of periodic boundary conditions
inside the group
(this has no effect on the treatment of the pbc between groups).
This option is only important when the diameter of the pull group
is larger than half the shortest box vector.
For determining the COM, all atoms in the group are put at their periodic image
-which is closest to <b>pull-pbcatom1</b>.
+which is closest to <b>pull-group1-pbcatom</b>.
A value of 0 means that the middle atom (number wise) is used.
This parameter is not used with geometry <b>cylinder</b>.
A value of -1 turns on cosine weighting, which is useful for a group
of molecules in a periodic system, e.g. a water slab (see Engin et al.
J. Chem. Phys. B 2010).</dd>
-<dt><b>pull-vec1: (0.0 0.0 0.0)</b></dt>
+
+<dt><b>pull-coord1-groups: </b></dt>
+<dd>The two groups indices should be given on which this pull coordinate
+will operate. The first index can be 0, in which case an absolute reference
+of <b>pull-coord1-origin</b> is used. With an absolute reference the system
+is no longer translation invariant and one should think about what to do with
+the <A HREF="#run">center of mass motion</A>.</dd>
+<dt><b>pull-coord1-origin: (0.0 0.0 0.0)</b></dt>
+<dd>The pull reference position for use with an absolute reference.</dd>
+<dt><b>pull-coord1-vec: (0.0 0.0 0.0)</b></dt>
<dd>The pull direction. <tt>grompp</tt> normalizes the vector.</dd>
-<dt><b>pull-init1: (0.0) / (0.0 0.0 0.0) [nm]</b></dt>
-<dd>The reference distance at t=0. This is a single value,
-except for geometry <b>position</b> which uses a vector.</dd>
-<dt><b>pull-rate1: (0) [nm/ps]</b></dt>
+<dt><b>pull-coord1-init: (0.0) [nm]</b></dt>
+<dd>The reference distance at t=0.</dd>
+<dt><b>pull-coord1-rate: (0) [nm/ps]</b></dt>
<dd>The rate of change of the reference position.</dd>
-<dt><b>pull-k1: (0) [kJ mol<sup>-1</sup> nm<sup>-2</sup>] / [kJ mol<sup>-1</sup> nm<sup>-1</sup>]</b></dt>
+<dt><b>pull-coord1-k: (0) [kJ mol<sup>-1</sup> nm<sup>-2</sup>] / [kJ mol<sup>-1</sup> nm<sup>-1</sup>]</b></dt>
<dd>The force constant. For umbrella pulling this is the harmonic force
constant in [kJ mol<sup>-1</sup> nm<sup>-2</sup>]. For constant force pulling
this is the force constant of the linear potential, and thus minus (!)
the constant force in [kJ mol<sup>-1</sup> nm<sup>-1</sup>].</dd>
-<dt><b>pull-kB1: (pull-k1) [kJ mol<sup>-1</sup> nm<sup>-2</sup>] / [kJ mol<sup>-1</sup> nm<sup>-1</sup>]</b></dt>
-<dd>As <b>pull-k1</b>, but for state B. This is only used when
+<dt><b>pull-coord1-kB: (pull-k1) [kJ mol<sup>-1</sup> nm<sup>-2</sup>] / [kJ mol<sup>-1</sup> nm<sup>-1</sup>]</b></dt>
+<dd>As <b>pull-coord1-k</b>, but for state B. This is only used when
<A HREF="#free"><b>free-energy</b></A> is turned on.
-The force constant is then (1 - lambda)*<b>pull-k1</b> + lambda*<b>pull-kB1</b>.
+The force constant is then (1 - lambda)*<b>pull-coord1-k</b> + lambda*<b>pull-coord1-kB</b></dt>.
+
</dl>
<A NAME="nmr"><br>
static const char *tpx_tag = TPX_TAG_RELEASE;
/* This number should be increased whenever the file format changes! */
-static const int tpx_version = 94;
+static const int tpx_version = 95;
/* This number should only be increased when you edit the TOPOLOGY section
* or the HEADER of the tpx format.
* Now the higer level routines that do io of the structures and arrays
*
**************************************************************/
-static void do_pullgrp(t_fileio *fio, t_pullgrp *pgrp, gmx_bool bRead,
- int file_version)
+static void do_pullgrp_tpx_pre95(t_fileio *fio,
+ t_pull_group *pgrp,
+ t_pull_coord *pcrd,
+ gmx_bool bRead,
+ int file_version)
{
- int i;
+ int i;
+ rvec tmp;
gmx_fio_do_int(fio, pgrp->nat);
if (bRead)
}
gmx_fio_ndo_real(fio, pgrp->weight, pgrp->nweight);
gmx_fio_do_int(fio, pgrp->pbcatom);
- gmx_fio_do_rvec(fio, pgrp->vec);
- gmx_fio_do_rvec(fio, pgrp->init);
- gmx_fio_do_real(fio, pgrp->rate);
- gmx_fio_do_real(fio, pgrp->k);
+ gmx_fio_do_rvec(fio, pcrd->vec);
+ clear_rvec(pcrd->origin);
+ gmx_fio_do_rvec(fio, tmp);
+ pcrd->init = tmp[0];
+ gmx_fio_do_real(fio, pcrd->rate);
+ gmx_fio_do_real(fio, pcrd->k);
if (file_version >= 56)
{
- gmx_fio_do_real(fio, pgrp->kB);
+ gmx_fio_do_real(fio, pcrd->kB);
}
else
{
- pgrp->kB = pgrp->k;
+ pcrd->kB = pcrd->k;
+ }
+}
+
+static void do_pull_group(t_fileio *fio, t_pull_group *pgrp, gmx_bool bRead)
+{
+ int i;
+
+ gmx_fio_do_int(fio, pgrp->nat);
+ if (bRead)
+ {
+ snew(pgrp->ind, pgrp->nat);
+ }
+ gmx_fio_ndo_int(fio, pgrp->ind, pgrp->nat);
+ gmx_fio_do_int(fio, pgrp->nweight);
+ if (bRead)
+ {
+ snew(pgrp->weight, pgrp->nweight);
}
+ gmx_fio_ndo_real(fio, pgrp->weight, pgrp->nweight);
+ gmx_fio_do_int(fio, pgrp->pbcatom);
+}
+
+static void do_pull_coord(t_fileio *fio, t_pull_coord *pcrd)
+{
+ int i;
+
+ gmx_fio_do_int(fio, pcrd->group[0]);
+ gmx_fio_do_int(fio, pcrd->group[1]);
+ gmx_fio_do_rvec(fio, pcrd->origin);
+ gmx_fio_do_rvec(fio, pcrd->vec);
+ gmx_fio_do_real(fio, pcrd->init);
+ gmx_fio_do_real(fio, pcrd->rate);
+ gmx_fio_do_real(fio, pcrd->k);
+ gmx_fio_do_real(fio, pcrd->kB);
}
static void do_expandedvals(t_fileio *fio, t_expanded *expand, t_lambda *fepvals, gmx_bool bRead, int file_version)
{
int g;
- gmx_fio_do_int(fio, pull->ngrp);
+ if (file_version >= 95)
+ {
+ gmx_fio_do_int(fio, pull->ngroup);
+ }
+ gmx_fio_do_int(fio, pull->ncoord);
+ if (file_version < 95)
+ {
+ pull->ngroup = pull->ncoord + 1;
+ }
gmx_fio_do_int(fio, pull->eGeom);
gmx_fio_do_ivec(fio, pull->dim);
gmx_fio_do_real(fio, pull->cyl_r1);
gmx_fio_do_real(fio, pull->cyl_r0);
gmx_fio_do_real(fio, pull->constr_tol);
+ if (file_version >= 95)
+ {
+ gmx_fio_do_int(fio, pull->bPrintRef);
+ }
gmx_fio_do_int(fio, pull->nstxout);
gmx_fio_do_int(fio, pull->nstfout);
if (bRead)
{
- snew(pull->grp, pull->ngrp+1);
+ snew(pull->group, pull->ngroup);
+ snew(pull->coord, pull->ncoord);
}
- for (g = 0; g < pull->ngrp+1; g++)
+ if (file_version < 95)
{
- do_pullgrp(fio, &pull->grp[g], bRead, file_version);
+ /* epullgPOS for position pulling, before epullgDIRPBC was removed */
+ if (pull->eGeom == epullgDIRPBC)
+ {
+ gmx_fatal(FARGS, "pull-geometry=position is no longer supported");
+ }
+ if (pull->eGeom > epullgDIRPBC)
+ {
+ pull->eGeom -= 1;
+ }
+
+ for (g = 0; g < pull->ngroup; g++)
+ {
+ /* We read and ignore a pull coordinate for group 0 */
+ do_pullgrp_tpx_pre95(fio, &pull->group[g], &pull->coord[max(g-1,0)],
+ bRead, file_version);
+ if (g > 0)
+ {
+ pull->coord[g-1].group[0] = 0;
+ pull->coord[g-1].group[1] = g;
+ }
+ }
+
+ pull->bPrintRef = (pull->group[0].nat > 0);
+ }
+ else
+ {
+ for (g = 0; g < pull->ngroup; g++)
+ {
+ do_pull_group(fio, &pull->group[g], bRead);
+ }
+ for (g = 0; g < pull->ncoord; g++)
+ {
+ do_pull_coord(fio, &pull->coord[g]);
+ }
}
}
* \name Using umbrella pull code since gromacs 4.x
*/
/*!\{*/
- int npullgrps; //!< nr of pull groups in tpr file
- int pull_geometry; //!< such as distance, position
+ int npullcrds; //!< nr of pull coordinates in tpr file
+ int pull_geometry; //!< such as distance, direction
ivec pull_dim; //!< pull dimension with geometry distance
int pull_ndim; //!< nr of pull_dim != 0
+ gmx_bool bPrintRef; //!< Coordinates of reference groups written to pullx.xvg ?
real *k; //!< force constants in tpr file
- rvec *init_dist; //!< reference displacements
+ real *init_dist; //!< reference displacements
real *umbInitDist; //!< reference displacement in umbrella direction
/*!\}*/
/*!
void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions *opt)
{
t_inputrec ir;
- int i, ngrp, d;
+ int i, ncrd;
t_state state;
static int first = 1;
}
/* nr of pull groups */
- ngrp = ir.pull->ngrp;
- if (ngrp < 1)
+ ncrd = ir.pull->ncoord;
+ if (ncrd < 1)
{
- gmx_fatal(FARGS, "This is not a tpr of umbrella simulation. Found only %d pull groups\n", ngrp);
+ gmx_fatal(FARGS, "This is not a tpr of umbrella simulation. Found only %d pull coordinates\n", ncrd);
}
- header->npullgrps = ir.pull->ngrp;
+ header->npullcrds = ir.pull->ncoord;
header->pull_geometry = ir.pull->eGeom;
+ header->bPrintRef = ir.pull->bPrintRef;
copy_ivec(ir.pull->dim, header->pull_dim);
header->pull_ndim = header->pull_dim[0]+header->pull_dim[1]+header->pull_dim[2];
- if (header->pull_geometry == epullgPOS && header->pull_ndim > 1)
- {
- gmx_fatal(FARGS, "Found pull geometry 'position' and more than 1 pull dimension (%d).\n"
- "Hence, the pull potential does not correspond to a one-dimensional umbrella potential.\n"
- "If you have some special umbrella setup you may want to write your own pdo files\n"
- "and feed them into g_wham. Check g_wham -h !\n", header->pull_ndim);
- }
- snew(header->k, ngrp);
- snew(header->init_dist, ngrp);
- snew(header->umbInitDist, ngrp);
+ snew(header->k, ncrd);
+ snew(header->init_dist, ncrd);
+ snew(header->umbInitDist, ncrd);
/* only z-direction with epullgCYL? */
if (header->pull_geometry == epullgCYL)
}
}
- for (i = 0; i < ngrp; i++)
+ for (i = 0; i < ncrd; i++)
{
- header->k[i] = ir.pull->grp[i+1].k;
+ header->k[i] = ir.pull->coord[i].k;
if (header->k[i] == 0.0)
{
- gmx_fatal(FARGS, "Pull group %d has force constant of of 0.0 in %s.\n"
+ gmx_fatal(FARGS, "Pull coordinate %d has force constant of of 0.0 in %s.\n"
"That doesn't seem to be an Umbrella tpr.\n",
i, fn);
}
- copy_rvec(ir.pull->grp[i+1].init, header->init_dist[i]);
+ header->init_dist[i] = ir.pull->coord[i].init;
/* initial distance to reference */
switch (header->pull_geometry)
{
- case epullgPOS:
- for (d = 0; d < DIM; d++)
- {
- if (header->pull_dim[d])
- {
- header->umbInitDist[i] = header->init_dist[i][d];
- }
- }
- break;
case epullgCYL:
- /* umbrella distance stored in init_dist[i][0] for geometry cylinder (not in ...[i][ZZ]) */
+ /* umbrella distance stored in init_dist[i] for geometry cylinder (not in ...[i][ZZ]) */
case epullgDIST:
case epullgDIR:
case epullgDIRPBC:
- header->umbInitDist[i] = header->init_dist[i][0];
+ header->umbInitDist[i] = header->init_dist[i];
break;
default:
gmx_fatal(FARGS, "Pull geometry %s not supported\n", epullg_names[header->pull_geometry]);
if (opt->verbose || first)
{
- printf("File %s, %d groups, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n",
- fn, header->npullgrps, epullg_names[header->pull_geometry],
+ printf("File %s, %d coordinates, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n"
+ "\tPull group coordinates%s expected in pullx files.\n",
+ fn, header->npullcrds, epullg_names[header->pull_geometry],
int2YN(header->pull_dim[0]), int2YN(header->pull_dim[1]), int2YN(header->pull_dim[2]),
- header->pull_ndim);
- for (i = 0; i < ngrp; i++)
+ header->pull_ndim, (header->bPrintRef ? "" : " not"));
+ for (i = 0; i < ncrd; i++)
{
- printf("\tgrp %d) k = %-5g position = %g\n", i, header->k[i], header->umbInitDist[i]);
+ printf("\tcrd %d) k = %-5g position = %g\n", i, header->k[i], header->umbInitDist[i]);
}
}
if (!opt->verbose && first)
{
- printf("\tUse option -v to see this output for all input tpr files\n");
+ printf("\tUse option -v to see this output for all input tpr files\n\n");
}
first = 0;
t_groupselection *groupsel)
{
double **y = 0, pos = 0., t, force, time0 = 0., dt;
- int ny, nt, bins, ibin, i, g, gUsed, dstep = 1, nColPerGrp, nColRefOnce, nColRefEachGrp, nColExpect, ntot;
+ int ny, nt, bins, ibin, i, g, gUsed, dstep = 1, nColPerCrd, nColRefEachCrd, nColExpect, ntot, column;
real min, max, minfound = 1e20, maxfound = -1e20;
gmx_bool dt_ok, timeok, bHaveForce;
const char *quantity;
static gmx_bool bFirst = TRUE;
/*
- in force output pullf.xvg:
- No reference, one column per pull group
- in position output pullx.xvg (not cylinder)
- ndim reference, ndim columns per pull group
- in position output pullx.xvg (in geometry cylinder):
- ndim*2 columns per pull group (ndim for ref, ndim for group)
+ * Data columns in pull output:
+ * - in force output pullf.xvg:
+ * No reference columns, one column per pull coordinate
+ *
+ * - in position output pullx.xvg
+ * bPrintRef == TRUE: for each pull coordinate: ndim reference columns, and ndim dx columns
+ * bPrintRef == FALSE: for each pull coordinate: no reference columns, but ndim dx columns
*/
- nColPerGrp = opt->bPullx ? header->pull_ndim : 1;
+ nColPerCrd = opt->bPullx ? header->pull_ndim : 1;
quantity = opt->bPullx ? "position" : "force";
- if (opt->bPullx)
+ if (opt->bPullx && header->bPrintRef)
{
- if (header->pull_geometry == epullgCYL)
- {
- /* Geometry cylinder -> reference group before each pull group */
- nColRefEachGrp = header->pull_ndim;
- nColRefOnce = 0;
- }
- else
- {
- /* Geometry NOT cylinder -> reference group only once after time column */
- nColRefEachGrp = 0;
- nColRefOnce = header->pull_ndim;
- }
+ nColRefEachCrd = header->pull_ndim;
}
- else /* read forces, no reference groups */
+ else
{
- nColRefEachGrp = 0;
- nColRefOnce = 0;
+ nColRefEachCrd = 0;
}
- nColExpect = 1 + nColRefOnce + header->npullgrps*(nColRefEachGrp+nColPerGrp);
+ nColExpect = 1 + header->npullcrds*(nColRefEachCrd+nColPerCrd);
bHaveForce = opt->bPullf;
/* With geometry "distance" or "distance_periodic", only force reading is supported so far.
- That avoids the somewhat tedious extraction of the right columns from the pullx files
- to compute the distances projection on the vector. Sorry for the laziness. */
+ That avoids the somewhat tedious extraction of the reaction coordinate from the pullx files.
+ Sorry for the laziness, this is a To-do. */
if ( (header->pull_geometry == epullgDIR || header->pull_geometry == epullgDIRPBC)
&& opt->bPullx)
{
bHaveForce ? "force" : "position", epullg_names[header->pull_geometry],
header->pull_ndim);
printf("Expecting these columns in pull file:\n"
- "\t%d reference columns for all pull groups together\n"
- "\t%d reference columns for each individual pull group\n"
- "\t%d data columns for each pull group\n", nColRefOnce, nColRefEachGrp, nColPerGrp);
- printf("With %d pull groups, expect %d columns (including the time column)\n", header->npullgrps, nColExpect);
+ "\t%d reference columns for each individual pull coordinate\n"
+ "\t%d data columns for each pull coordinate\n", nColRefEachCrd, nColPerCrd);
+ printf("With %d pull groups, expect %d columns (including the time column)\n", header->npullcrds, nColExpect);
bFirst = FALSE;
}
if (ny != nColExpect)
{
- gmx_fatal(FARGS, "Found %d pull groups in %s,\n but %d data columns in %s (expected %d)\n"
+ gmx_fatal(FARGS, "Found %d pull coodinates in %s,\n but %d data columns in %s (expected %d)\n"
"\nMaybe you confused options -ix and -if ?\n",
- header->npullgrps, fntpr, ny-1, fn, nColExpect-1);
+ header->npullcrds, fntpr, ny-1, fn, nColExpect-1);
}
if (opt->verbose)
{
- printf("Found %d times and %d %s sets %s\n", nt, (ny-1)/nColPerGrp, quantity, fn);
+ printf("Found %d times and %d %s sets %s\n", nt, (ny-1)/nColPerCrd, quantity, fn);
}
if (!bGetMinMax)
if (groupsel)
{
/* Use only groups selected with option -is file */
- if (header->npullgrps != groupsel->n)
+ if (header->npullcrds != groupsel->n)
{
gmx_fatal(FARGS, "tpr file contains %d pull groups, but expected %d from group selection file\n",
- header->npullgrps, groupsel->n);
+ header->npullcrds, groupsel->n);
}
window->nPull = groupsel->nUse;
}
else
{
- window->nPull = header->npullgrps;
+ window->nPull = header->npullcrds;
}
window->nBin = bins;
/* Copying umbrella center and force const is more involved since not
all pull groups from header (tpr file) may be used in window variable */
- for (g = 0, gUsed = 0; g < header->npullgrps; ++g)
+ for (g = 0, gUsed = 0; g < header->npullcrds; ++g)
{
if (groupsel && (groupsel->bUse[g] == FALSE))
{
/*if (opt->verbose)
printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n",
t,opt->tmin, opt->tmax, dt_ok,timeok); */
-
if (timeok)
{
/* Note: if groupsel == NULL:
* only groups with groupsel.bUse[g]==TRUE are stored. gUsed is not always equal g
*/
gUsed = -1;
- for (g = 0; g < header->npullgrps; ++g)
+ for (g = 0; g < header->npullcrds; ++g)
{
/* was this group selected for application in WHAM? */
if (groupsel && (groupsel->bUse[g] == FALSE))
if (bHaveForce)
{
- /* y has 1 time column y[0] and one column per force y[1],...,y[nGrps] */
+ /* y has 1 time column y[0] and one column per force y[1],...,y[nCrds] */
force = y[g+1][i];
pos = -force/header->k[g] + header->umbInitDist[g];
}
else
{
+ /* pick the right column from:
+ * time ref1[ndim] coord1[ndim] ref2[ndim] coord2[ndim] ...
+ */
+ column = 1 + nColRefEachCrd + g*(nColRefEachCrd+nColPerCrd);
switch (header->pull_geometry)
{
case epullgDIST:
- /* y has 1 time column y[0] and nColPerGrps columns per pull group;
- Distance to reference: */
- /* pos=dist_ndim(y+1+nColRef+g*nColPerGrp,header->pull_ndim,i); gmx 4.0 */
- pos = dist_ndim(y + 1 + nColRefOnce + g*nColPerGrp, header->pull_ndim, i);
+ pos = dist_ndim(y + column, header->pull_ndim, i);
break;
- case epullgPOS:
- /* Columns
- Time ref[ndim] group1[ndim] group2[ndim] ... */
case epullgCYL:
- /* Columns
- Time ref1[ndim] group1[ndim] ref2[ndim] group2[ndim] ... */
-
- /* * with geometry==position, we have the reference once (nColRefOnce==ndim), but
- no extra reference group columns before each group (nColRefEachGrp==0)
-
- * with geometry==cylinder, we have no initial ref group column (nColRefOnce==0),
- but ndim ref group colums before every group (nColRefEachGrp==ndim)
- Distance to reference: */
- pos = y[1 + nColRefOnce + nColRefEachGrp + g*(nColRefEachGrp+nColPerGrp)][i];
+ pos = y[column][i];
break;
default:
gmx_fatal(FARGS, "Bad error, this error should have been catched before. Ups.\n");
}
}
- /* printf("grp %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
+ /* printf("crd %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
if (bGetMinMax)
{
if (pos < minfound)
}
}
-static void bc_pullgrp(const t_commrec *cr, t_pullgrp *pgrp)
+static void bc_pull_group(const t_commrec *cr, t_pull_group *pgrp)
{
block_bc(cr, *pgrp);
if (pgrp->nat > 0)
int g;
block_bc(cr, *pull);
- snew_bc(cr, pull->grp, pull->ngrp+1);
- for (g = 0; g < pull->ngrp+1; g++)
+ snew_bc(cr, pull->group, pull->ngroup);
+ for (g = 0; g < pull->ngroup; g++)
{
- bc_pullgrp(cr, &pull->grp[g]);
+ bc_pull_group(cr, &pull->group[g]);
}
+ snew_bc(cr, pull->coord, pull->ncoord);
+ nblock_bc(cr, pull->ncoord, pull->coord);
}
static void bc_rotgrp(const t_commrec *cr, t_rotgrp *rotg)
};
const char *epullg_names[epullgNR+1] = {
- "distance", "direction", "cylinder", "position", "direction-periodic", NULL
+ "distance", "direction", "cylinder", "direction-periodic", NULL
};
const char *erotg_names[erotgNR+1] = {
#define PR(t, s) pr_real(fp, indent, t, s)
#define PD(t, s) pr_double(fp, indent, t, s)
-static void pr_pullgrp(FILE *fp, int indent, int g, t_pullgrp *pg)
+static void pr_pull_group(FILE *fp, int indent, int g, t_pull_group *pgrp)
{
pr_indent(fp, indent);
fprintf(fp, "pull-group %d:\n", g);
indent += 2;
- pr_ivec_block(fp, indent, "atom", pg->ind, pg->nat, TRUE);
- pr_rvec(fp, indent, "weight", pg->weight, pg->nweight, TRUE);
- PI("pbcatom", pg->pbcatom);
- pr_rvec(fp, indent, "vec", pg->vec, DIM, TRUE);
- pr_rvec(fp, indent, "init", pg->init, DIM, TRUE);
- PR("rate", pg->rate);
- PR("k", pg->k);
- PR("kB", pg->kB);
+ pr_ivec_block(fp, indent, "atom", pgrp->ind, pgrp->nat, TRUE);
+ pr_rvec(fp, indent, "weight", pgrp->weight, pgrp->nweight, TRUE);
+ PI("pbcatom", pgrp->pbcatom);
+}
+
+static void pr_pull_coord(FILE *fp, int indent, int c, t_pull_coord *pcrd)
+{
+ pr_indent(fp, indent);
+ fprintf(fp, "pull-coord %d:\n", c);
+ PI("group[0]", pcrd->group[0]);
+ PI("group[1]", pcrd->group[1]);
+ pr_rvec(fp, indent, "origin", pcrd->origin, DIM, TRUE);
+ pr_rvec(fp, indent, "vec", pcrd->vec, DIM, TRUE);
+ PR("init", pcrd->init);
+ PR("rate", pcrd->rate);
+ PR("k", pcrd->k);
+ PR("kB", pcrd->kB);
}
static void pr_simtempvals(FILE *fp, int indent, t_simtemp *simtemp, int n_lambda)
PR("pull-r1", pull->cyl_r1);
PR("pull-r0", pull->cyl_r0);
PR("pull-constr-tol", pull->constr_tol);
+ PS("pull-bPrintRef", EBOOL(pull->bPrintRef));
PI("pull-nstxout", pull->nstxout);
PI("pull-nstfout", pull->nstfout);
- PI("pull-ngrp", pull->ngrp);
- for (g = 0; g < pull->ngrp+1; g++)
+ PI("pull-ngroup", pull->ngroup);
+ for (g = 0; g < pull->ngroup; g++)
+ {
+ pr_pull_group(fp, indent, g, &pull->group[g]);
+ }
+ PI("pull-ncoord", pull->ncoord);
+ for (g = 0; g < pull->ncoord; g++)
{
- pr_pullgrp(fp, indent, g, &pull->grp[g]);
+ pr_pull_coord(fp, indent, g, &pull->coord[g]);
}
}
done_blocka(&(top->excls));
}
-static void done_pullgrp(t_pullgrp *pgrp)
+static void done_pull_group(t_pull_group *pgrp)
{
- sfree(pgrp->ind);
- sfree(pgrp->ind_loc);
- sfree(pgrp->weight);
- sfree(pgrp->weight_loc);
+ if (pgrp->nat > 0)
+ {
+ sfree(pgrp->ind);
+ sfree(pgrp->ind_loc);
+ sfree(pgrp->weight);
+ sfree(pgrp->weight_loc);
+ }
}
static void done_pull(t_pull *pull)
{
int i;
- for (i = 0; i < pull->ngrp+1; i++)
+ for (i = 0; i < pull->ngroup+1; i++)
{
- done_pullgrp(pull->grp);
- done_pullgrp(pull->dyna);
+ done_pull_group(pull->group);
+ done_pull_group(pull->dyna);
}
}
sfree(dumstr[1]);
}
-static int search_QMstring(char *s, int ng, const char *gn[])
+static int search_QMstring(const char *s, int ng, const char *gn[])
{
/* same as normal search_string, but this one searches QM strings */
int i;
} /* search_QMstring */
-
-int search_string(char *s, int ng, char *gn[])
+/* We would like gn to be const as well, but C doesn't allow this */
+int search_string(const char *s, int ng, char *gn[])
{
int i;
t_grpopts *opts;
gmx_groups_t *groups;
t_pull *pull;
- int natoms, ai, aj, i, j, d, g, imin, jmin, nc;
+ int natoms, ai, aj, i, j, d, g, imin, jmin;
t_iatom *ia;
int *nrdf2, *na_vcm, na_tot;
double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
* to determine the optimal nrdf assignment.
*/
pull = ir->pull;
- if (pull->eGeom == epullgPOS)
+
+ for (i = 0; i < pull->ncoord; i++)
{
- nc = 0;
- for (i = 0; i < DIM; i++)
+ imin = 1;
+
+ for (j = 0; j < 2; j++)
{
- if (pull->dim[i])
+ const t_pull_group *pgrp;
+
+ pgrp = &pull->group[pull->coord[i].group[j]];
+
+ if (pgrp->nat > 0)
{
- nc++;
+ /* Subtract 1/2 dof from each group */
+ ai = pgrp->ind[0];
+ nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
+ nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
+ if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
+ {
+ gmx_fatal(FARGS, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups->grps[egcTC].nm_ind[ggrpnr(groups, egcTC, ai)]]);
+ }
}
- }
- }
- else
- {
- nc = 1;
- }
- for (i = 0; i < pull->ngrp; i++)
- {
- imin = 2*nc;
- if (pull->grp[0].nat > 0)
- {
- /* Subtract 1/2 dof from the reference group */
- ai = pull->grp[0].ind[0];
- if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1)
+ else
{
- nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5;
- nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5;
- imin--;
+ /* We need to subtract the whole DOF from group j=1 */
+ imin += 1;
}
}
- /* Subtract 1/2 dof from the pulled group */
- ai = pull->grp[1+i].ind[0];
- nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
- nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
- if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
- {
- gmx_fatal(FARGS, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups->grps[egcTC].nm_ind[ggrpnr(groups, egcTC, ai)]]);
- }
}
}
if (ir->ePull != epullNO)
{
make_pull_groups(ir->pull, pull_grp, grps, gnames);
+
+ make_pull_coords(ir->pull);
}
if (ir->bRot)
warninp_t wi)
{
char err_buf[256];
- int i, m, g, nmol, npct;
+ int i, m, c, nmol, npct;
gmx_bool bCharge, bAcc;
real gdt_max, *mgrp, mt;
rvec acc;
if (ir->ePull != epullNO)
{
- if (ir->pull->grp[0].nat == 0)
+ gmx_bool bPullAbsoluteRef;
+
+ bPullAbsoluteRef = FALSE;
+ for (i = 0; i < ir->pull->ncoord; i++)
+ {
+ bPullAbsoluteRef = bPullAbsoluteRef ||
+ ir->pull->coord[i].group[0] == 0 ||
+ ir->pull->coord[i].group[1] == 0;
+ }
+ if (bPullAbsoluteRef)
{
absolute_reference(ir, sys, FALSE, AbsRef);
for (m = 0; m < DIM; m++)
if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
ir->deform[i][m] != 0)
{
- for (g = 1; g < ir->pull->ngrp; g++)
+ for (c = 0; c < ir->pull->ncoord; c++)
{
- if (ir->pull->grp[g].vec[m] != 0)
+ if (ir->pull->coord[c].vec[m] != 0)
{
gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
}
/* Validate inputrec data.
* Fatal errors will be added to nerror.
*/
-extern int search_string(char *s, int ng, char *gn[]);
+extern int search_string(const char *s, int ng, char *gn[]);
/* Returns the index of string s in the index groups */
extern void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr,
warninp_t wi);
/* Reads the pull parameters, returns a list of the pull group names */
-extern void make_pull_groups(t_pull *pull, char **pgnames,
- t_blocka *grps, char **gnames);
-/* Process the pull parameters after reading the index groups */
+extern void make_pull_groups(t_pull *pull,
+ char **pgnames,
+ const t_blocka *grps, char **gnames);
+/* Process the pull group parameters after reading the index groups */
+
+extern void make_pull_coords(t_pull *pull);
+/* Process the pull coordinates after reading the pull groups */
extern void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real lambda,
const output_env_t oenv, gmx_bool bStart);
static char pulldim[STRLEN];
-static void string2dvec(char buf[], dvec nums)
+static void string2dvec(const char buf[], dvec nums)
{
- if (sscanf(buf, "%lf%lf%lf", &nums[0], &nums[1], &nums[2]) != 3)
+ double dum;
+
+ if (sscanf(buf, "%lf%lf%lf%lf", &nums[0], &nums[1], &nums[2], &dum) != 3)
{
gmx_fatal(FARGS, "Expected three numbers at input line %s", buf);
}
}
-static void init_pullgrp(t_pullgrp *pg, char *wbuf,
- gmx_bool bRef, int eGeom, char *s_vec)
+static void init_pull_group(t_pull_group *pg,
+ const char *wbuf)
{
double d;
int n, m;
- dvec vec;
pg->nweight = 0;
while (sscanf(wbuf, "%lf %n", &d, &n) == 1)
pg->weight[pg->nweight++] = d;
wbuf += n;
}
- if (!bRef)
+}
+
+static void init_pull_coord(t_pull_coord *pcrd, int eGeom,
+ const char *origin_buf, const char *vec_buf)
+{
+ int m;
+ dvec origin, vec;
+
+ string2dvec(origin_buf, origin);
+ if (pcrd->group[0] != 0 && dnorm(origin) > 0)
{
- if (eGeom == epullgDIST)
- {
- clear_dvec(vec);
- }
- else
- {
- string2dvec(s_vec, vec);
- if (eGeom == epullgDIR || eGeom == epullgCYL ||
- (eGeom == epullgPOS && dnorm(vec) != 0))
- {
- /* Normalize the direction vector */
- dsvmul(1/dnorm(vec), vec, vec);
- }
- }
- for (m = 0; m < DIM; m++)
+ gmx_fatal(FARGS,"The pull origin can only be set with an absolute reference");
+ }
+
+ if (eGeom == epullgDIST)
+ {
+ clear_dvec(vec);
+ }
+ else
+ {
+ string2dvec(vec_buf, vec);
+ if (eGeom == epullgDIR || eGeom == epullgCYL)
{
- pg->vec[m] = vec[m];
+ /* Normalize the direction vector */
+ dsvmul(1/dnorm(vec), vec, vec);
}
}
+ for (m = 0; m < DIM; m++)
+ {
+ pcrd->origin[m] = origin[m];
+ pcrd->vec[m] = vec[m];
+ }
}
char **read_pullparams(int *ninp_p, t_inpfile **inp_p,
t_pull *pull, gmx_bool *bStart,
warninp_t wi)
{
- int ninp, nerror = 0, i, nchar, ndim, nscan, m;
- t_inpfile *inp;
- const char *tmp;
- char **grpbuf;
- char dummy[STRLEN], buf[STRLEN], init[STRLEN];
- const char *init_def1 = "0.0", *init_def3 = "0.0 0.0 0.0";
- char wbuf[STRLEN], VecTemp[STRLEN];
- dvec vec;
+ int ninp, nerror = 0, i, nchar, nscan, m, idum;
+ t_inpfile *inp;
+ const char *tmp;
+ char **grpbuf;
+ char dummy[STRLEN], buf[STRLEN], groups[STRLEN], init[STRLEN];
+ const char *init_def1 = "0.0", *init_def3 = "0.0 0.0 0.0";
+ char wbuf[STRLEN], origin_buf[STRLEN], vec_buf[STRLEN];
- t_pullgrp *pgrp;
+ t_pull_group *pgrp;
+ t_pull_coord *pcrd;
ninp = *ninp_p;
inp = *inp_p;
/* read pull parameters */
- CTYPE("Pull geometry: distance, direction, cylinder or position");
- EETYPE("pull_geometry", pull->eGeom, epullg_names);
+ CTYPE("Pull geometry: distance, direction, direction-periodic or cylinder");
+ EETYPE("pull-geometry", pull->eGeom, epullg_names);
CTYPE("Select components for the pull vector. default: Y Y Y");
- STYPE("pull_dim", pulldim, "Y Y Y");
+ STYPE("pull-dim", pulldim, "Y Y Y");
CTYPE("Cylinder radius for dynamic reaction force groups (nm)");
- RTYPE("pull_r1", pull->cyl_r1, 1.0);
+ RTYPE("pull-r1", pull->cyl_r1, 1.0);
CTYPE("Switch from r1 to r0 in case of dynamic reaction force");
- RTYPE("pull_r0", pull->cyl_r0, 1.5);
+ RTYPE("pull-r0", pull->cyl_r0, 1.5);
RTYPE("pull_constr_tol", pull->constr_tol, 1E-6);
- EETYPE("pull_start", *bStart, yesno_names);
- ITYPE("pull_nstxout", pull->nstxout, 10);
- ITYPE("pull_nstfout", pull->nstfout, 1);
+ EETYPE("pull-start", *bStart, yesno_names);
+ EETYPE("pull-print-reference", pull->bPrintRef, yesno_names);
+ ITYPE("pull-nstxout", pull->nstxout, 10);
+ ITYPE("pull-nstfout", pull->nstfout, 1);
CTYPE("Number of pull groups");
- ITYPE("pull_ngroups", pull->ngrp, 1);
+ ITYPE("pull-ngroups", pull->ngroup, 1);
+ CTYPE("Number of pull coordinates");
+ ITYPE("pull-ncoords", pull->ncoord, 1);
if (pull->cyl_r1 > pull->cyl_r0)
{
- warning_error(wi, "pull_r1 > pull_r0");
+ warning_error(wi, "pull-r1 > pull_r0");
}
- if (pull->ngrp < 1)
+ if (pull->ngroup < 1)
{
- gmx_fatal(FARGS, "pull_ngroups should be >= 1");
+ gmx_fatal(FARGS, "pull-ngroups should be >= 1");
}
+ /* We always add an absolute reference group (index 0), even if not used */
+ pull->ngroup += 1;
- snew(pull->grp, pull->ngrp+1);
-
- if (pull->eGeom == epullgPOS)
- {
- ndim = 3;
- }
- else
+ if (pull->ncoord < 1)
{
- ndim = 1;
+ gmx_fatal(FARGS, "pull-ncoords should be >= 1");
}
+ snew(pull->group, pull->ngroup);
+
+ snew(pull->coord, pull->ncoord);
+
/* pull group options */
CTYPE("Group name, weight (default all 1), vector, init, rate (nm/ps), kJ/(mol*nm^2)");
+
/* Read the pull groups */
- snew(grpbuf, pull->ngrp+1);
- for (i = 0; i < pull->ngrp+1; i++)
+ snew(grpbuf, pull->ngroup);
+ /* Group 0 is the absolute reference, we don't read anything for 0 */
+ for (i = 1; i < pull->ngroup; i++)
{
- pgrp = &pull->grp[i];
+ pgrp = &pull->group[i];
snew(grpbuf[i], STRLEN);
- sprintf(buf, "pull_group%d", i);
+ sprintf(buf, "pull-group%d-name", i);
STYPE(buf, grpbuf[i], "");
- sprintf(buf, "pull_weights%d", i);
+ sprintf(buf, "pull-group%d-weights", i);
STYPE(buf, wbuf, "");
- sprintf(buf, "pull_pbcatom%d", i);
+ sprintf(buf, "pull-group%d-pbcatom", i);
ITYPE(buf, pgrp->pbcatom, 0);
- if (i > 0)
- {
- sprintf(buf, "pull_vec%d", i);
- STYPE(buf, VecTemp, "0.0 0.0 0.0");
- sprintf(buf, "pull_init%d", i);
- STYPE(buf, init, ndim == 1 ? init_def1 : init_def3);
- nscan = sscanf(init, "%lf %lf %lf", &vec[0], &vec[1], &vec[2]);
- if (nscan != ndim)
- {
- fprintf(stderr, "ERROR: %s should have %d components\n", buf, ndim);
- nerror++;
- }
- for (m = 0; m < DIM; m++)
- {
- pgrp->init[m] = (m < ndim ? vec[m] : 0.0);
- }
- sprintf(buf, "pull_rate%d", i);
- RTYPE(buf, pgrp->rate, 0.0);
- sprintf(buf, "pull_k%d", i);
- RTYPE(buf, pgrp->k, 0.0);
- sprintf(buf, "pull_kB%d", i);
- RTYPE(buf, pgrp->kB, pgrp->k);
- }
/* Initialize the pull group */
- init_pullgrp(pgrp, wbuf, i == 0, pull->eGeom, VecTemp);
+ init_pull_group(pgrp, wbuf);
+ }
+
+ /* Read the pull coordinates */
+ for (i = 1; i < pull->ncoord + 1; i++)
+ {
+ pcrd = &pull->coord[i-1];
+ sprintf(buf, "pull-coord%d-groups", i);
+ STYPE(buf, groups, "");
+ nscan = sscanf(groups, "%d %d %d", &pcrd->group[0], &pcrd->group[1], &idum);
+ if (nscan != 2)
+ {
+ fprintf(stderr, "ERROR: %s should have %d components\n", buf, 2);
+ nerror++;
+ }
+ sprintf(buf, "pull-coord%d-origin", i);
+ STYPE(buf, origin_buf, "0.0 0.0 0.0");
+ sprintf(buf, "pull-coord%d-vec", i);
+ STYPE(buf, vec_buf, "0.0 0.0 0.0");
+ sprintf(buf, "pull-coord%d-init", i);
+ RTYPE(buf, pcrd->init, 0.0);
+ sprintf(buf, "pull-coord%d-rate", i);
+ RTYPE(buf, pcrd->rate, 0.0);
+ sprintf(buf, "pull-coord%d-k", i);
+ RTYPE(buf, pcrd->k, 0.0);
+ sprintf(buf, "pull-coord%d-kB", i);
+ RTYPE(buf, pcrd->kB, pcrd->k);
+
+ /* Initialize the pull coordinate */
+ init_pull_coord(pcrd, pull->eGeom, origin_buf, vec_buf);
}
*ninp_p = ninp;
return grpbuf;
}
-void make_pull_groups(t_pull *pull, char **pgnames, t_blocka *grps, char **gnames)
+void make_pull_groups(t_pull *pull,
+ char **pgnames,
+ const t_blocka *grps, char **gnames)
{
- int d, nchar, g, ig = -1, i;
- char *ptr, pulldim1[STRLEN];
- t_pullgrp *pgrp;
+ int g, ig = -1, i;
+ t_pull_group *pgrp;
- ptr = pulldim;
- i = 0;
+ /* Absolute reference group (might not be used) is special */
+ pgrp = &pull->group[0];
+ pgrp->nat = 0;
+ pgrp->pbcatom = -1;
+
+ for (g = 1; g < pull->ngroup; g++)
+ {
+ pgrp = &pull->group[g];
+
+ ig = search_string(pgnames[g], grps->nr, gnames);
+ pgrp->nat = grps->index[ig+1] - grps->index[ig];
+
+ fprintf(stderr, "Pull group %d '%s' has %d atoms\n",
+ g, pgnames[g], pgrp->nat);
+
+ if (pgrp->nat == 0)
+ {
+ gmx_fatal(FARGS, "Pull group %d '%s' is empty", g, pgnames[g]);
+ }
+
+ snew(pgrp->ind, pgrp->nat);
+ for (i = 0; i < pgrp->nat; i++)
+ {
+ pgrp->ind[i] = grps->a[grps->index[ig]+i];
+ }
+
+ if (pull->eGeom == epullgCYL && g == 1 && pgrp->nweight > 0)
+ {
+ gmx_fatal(FARGS, "Weights are not supported for the reference group with cylinder pulling");
+ }
+ if (pgrp->nweight > 0 && pgrp->nweight != pgrp->nat)
+ {
+ gmx_fatal(FARGS, "Number of weights (%d) for pull group %d '%s' does not match the number of atoms (%d)",
+ pgrp->nweight, g, pgnames[g], pgrp->nat);
+ }
+
+ if (pgrp->nat == 1)
+ {
+ /* No pbc is required for this group */
+ pgrp->pbcatom = -1;
+ }
+ else
+ {
+ if (pgrp->pbcatom > 0)
+ {
+ pgrp->pbcatom -= 1;
+ }
+ else if (pgrp->pbcatom == 0)
+ {
+ pgrp->pbcatom = pgrp->ind[(pgrp->nat-1)/2];
+ }
+ else
+ {
+ /* Use cosine weighting */
+ pgrp->pbcatom = -1;
+ }
+ }
+ }
+}
+
+void make_pull_coords(t_pull *pull)
+{
+ int ndim, d, nchar, c;
+ char *ptr, pulldim1[STRLEN];
+ t_pull_coord *pcrd;
+
+ ptr = pulldim;
+ ndim = 0;
for (d = 0; d < DIM; d++)
{
if (sscanf(ptr, "%s%n", pulldim1, &nchar) != 1)
else if (gmx_strncasecmp(pulldim1, "Y", 1) == 0)
{
pull->dim[d] = 1;
- i++;
+ ndim++;
}
else
{
}
ptr += nchar;
}
- if (i == 0)
+ if (ndim == 0)
{
gmx_fatal(FARGS, "All entries in pull_dim are N");
}
- for (g = 0; g < pull->ngrp+1; g++)
+ for (c = 0; c < pull->ncoord; c++)
{
- pgrp = &pull->grp[g];
- if (g == 0 && strcmp(pgnames[g], "") == 0)
+ pcrd = &pull->coord[c];
+
+ if (pcrd->group[0] < 0 || pcrd->group[0] >= pull->ngroup ||
+ pcrd->group[1] < 0 || pcrd->group[1] >= pull->ngroup)
{
- pgrp->nat = 0;
+ gmx_fatal(FARGS, "Pull group index in pull-coord%d-groups out of range, should be between %d and %d", c+1, 0, pull->ngroup+1);
}
- else
+
+ if (pcrd->group[0] == pcrd->group[1])
{
- ig = search_string(pgnames[g], grps->nr, gnames);
- pgrp->nat = grps->index[ig+1] - grps->index[ig];
+ gmx_fatal(FARGS, "Identical pull group indices in pull-coord%d-groups", c+1);
}
- if (pgrp->nat > 0)
- {
- fprintf(stderr, "Pull group %d '%s' has %d atoms\n",
- g, pgnames[g], pgrp->nat);
- snew(pgrp->ind, pgrp->nat);
- for (i = 0; i < pgrp->nat; i++)
- {
- pgrp->ind[i] = grps->a[grps->index[ig]+i];
- }
-
- if (pull->eGeom == epullgCYL && g == 0 && pgrp->nweight > 0)
- {
- gmx_fatal(FARGS, "Weights are not supported for the reference group with cylinder pulling");
- }
- if (pgrp->nweight > 0 && pgrp->nweight != pgrp->nat)
- {
- gmx_fatal(FARGS, "Number of weights (%d) for pull group %d '%s' does not match the number of atoms (%d)",
- pgrp->nweight, g, pgnames[g], pgrp->nat);
- }
- if (pgrp->nat == 1)
- {
- /* No pbc is required for this group */
- pgrp->pbcatom = -1;
- }
- else
- {
- if (pgrp->pbcatom > 0)
- {
- pgrp->pbcatom -= 1;
- }
- else if (pgrp->pbcatom == 0)
- {
- pgrp->pbcatom = pgrp->ind[(pgrp->nat-1)/2];
- }
- else
- {
- /* Use cosine weighting */
- pgrp->pbcatom = -1;
- }
- }
+ if (pull->eGeom == epullgCYL && pcrd->group[0] != 1)
+ {
+ gmx_fatal(FARGS, "With pull geometry %s, the first pull group should always be 1", EPULLGEOM(pull->eGeom));
+ }
- if (g > 0 && pull->eGeom != epullgDIST)
+ if (pull->eGeom != epullgDIST)
+ {
+ for (d = 0; d < DIM; d++)
{
- for (d = 0; d < DIM; d++)
+ if (pcrd->vec[d] != 0 && pull->dim[d] == 0)
{
- if (pgrp->vec[d] != 0 && pull->dim[d] == 0)
- {
- gmx_fatal(FARGS, "ERROR: pull_vec%d has non-zero %c-component while pull_dim in N\n", g, 'x'+d);
- }
+ gmx_fatal(FARGS, "ERROR: pull-group%d-vec has non-zero %c-component while pull_dim for the %c-dimension is N\n", c+1, 'x'+d, 'x'+d);
}
}
-
- if ((pull->eGeom == epullgDIR || pull->eGeom == epullgCYL) &&
- g > 0 && norm2(pgrp->vec) == 0)
- {
- gmx_fatal(FARGS, "pull_vec%d can not be zero with geometry %s",
- g, EPULLGEOM(pull->eGeom));
- }
- if ((pull->eGeom == epullgPOS) && pgrp->rate != 0 &&
- g > 0 && norm2(pgrp->vec) == 0)
- {
- gmx_fatal(FARGS, "pull_vec%d can not be zero with geometry %s and non-zero rate",
- g, EPULLGEOM(pull->eGeom));
- }
}
- else
+
+ if ((pull->eGeom == epullgDIR || pull->eGeom == epullgCYL) &&
+ norm2(pcrd->vec) == 0)
{
- if (g == 0)
- {
- if (pull->eGeom == epullgCYL)
- {
- gmx_fatal(FARGS, "Absolute reference groups are not supported with geometry %s", EPULLGEOM(pull->eGeom));
- }
- }
- else
- {
- gmx_fatal(FARGS, "Pull group %d '%s' is empty", g, pgnames[g]);
- }
- pgrp->pbcatom = -1;
+ gmx_fatal(FARGS, "pull-group%d-vec can not be zero with geometry %s",
+ c+1, EPULLGEOM(pull->eGeom));
}
}
}
void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real lambda,
const output_env_t oenv, gmx_bool bStart)
{
- t_mdatoms *md;
- t_pull *pull;
- t_pullgrp *pgrp;
- t_pbc pbc;
- int ndim, g, m;
- double t_start, tinvrate;
- rvec init;
- dvec dr, dev;
+ t_mdatoms *md;
+ t_pull *pull;
+ t_pull_coord *pcrd;
+ t_pull_group *pgrp0, *pgrp1;
+ t_pbc pbc;
+ int c, m;
+ double t_start, tinvrate;
+ real init;
+ dvec dr;
+ double dev;
init_pull(NULL, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0);
md = init_mdatoms(NULL, mtop, ir->efep);
update_mdatoms(md, lambda);
}
pull = ir->pull;
- if (pull->eGeom == epullgPOS)
- {
- ndim = 3;
- }
- else
- {
- ndim = 1;
- }
set_pbc(&pbc, ir->ePBC, box);
pull_calc_coms(NULL, pull, md, &pbc, t_start, x, NULL);
fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n");
- for (g = 0; g < pull->ngrp+1; g++)
+ for (c = 0; c < pull->ncoord; c++)
{
- pgrp = &pull->grp[g];
- fprintf(stderr, "%8d %8d %8d ", g, pgrp->nat, pgrp->pbcatom+1);
- copy_rvec(pgrp->init, init);
- clear_rvec(pgrp->init);
- if (g > 0)
+ pcrd = &pull->coord[c];
+
+ pgrp0 = &pull->group[pcrd->group[0]];
+ pgrp1 = &pull->group[pcrd->group[1]];
+ fprintf(stderr, "%8d %8d %8d\n",
+ pcrd->group[0], pgrp0->nat, pgrp0->pbcatom+1);
+ fprintf(stderr, "%8d %8d %8d ",
+ pcrd->group[1], pgrp1->nat, pgrp1->pbcatom+1);
+
+ init = pcrd->init;
+ pcrd->init = 0;
+
+ if (pcrd->rate == 0)
{
- if (pgrp->rate == 0)
- {
- tinvrate = 0;
- }
- else
- {
- tinvrate = t_start/pgrp->rate;
- }
- get_pullgrp_distance(pull, &pbc, g, 0, dr, dev);
- for (m = 0; m < DIM; m++)
- {
- if (m < ndim)
- {
- fprintf(stderr, " %6.3f", dev[m]);
- }
- else
- {
- fprintf(stderr, " ");
- }
- }
- fprintf(stderr, " ");
- for (m = 0; m < DIM; m++)
- {
- if (bStart)
- {
- pgrp->init[m] = init[m] + dev[m]
- - tinvrate*(pull->eGeom == epullgPOS ? pgrp->vec[m] : 1);
- }
- else
- {
- pgrp->init[m] = init[m];
- }
- if (m < ndim)
- {
- fprintf(stderr, " %6.3f", pgrp->init[m]);
- }
- else
- {
- fprintf(stderr, " ");
- }
- }
+ tinvrate = 0;
+ }
+ else
+ {
+ tinvrate = t_start/pcrd->rate;
+ }
+ get_pull_coord_distance(pull, c, &pbc, 0, dr, &dev);
+ fprintf(stderr, " %6.3f ", dev);
+
+ if (bStart)
+ {
+ pcrd->init = init + dev - tinvrate;
+ }
+ else
+ {
+ pcrd->init = init;
}
- fprintf(stderr, "\n");
+ fprintf(stderr, " %6.3f\n", pcrd->init);
}
}
/* This file contains datatypes and function declarations necessary
for mdrun to interface with the pull code */
-/* Get the distance to the reference and deviation for pull group g */
-void get_pullgrp_distance(t_pull *pull, t_pbc *pbc, int g, double t,
- dvec dr, dvec dev);
+/* Get the distance to the reference and deviation for pull coord coord_ind */
+void get_pull_coord_distance(const t_pull *pull,
+ int coord_ind,
+ const t_pbc *pbc, double t,
+ dvec dr, double *dev);
/* Set the all the pull forces to zero */
void clear_pull_forces(t_pull *pull);
};
enum {
- epullgDIST, epullgDIR, epullgCYL, epullgPOS, epullgDIRPBC, epullgNR
+ epullgDIST, epullgDIR, epullgCYL, epullgDIRPBC, epullgNR
};
#define PULL_CYL(pull) ((pull)->eGeom == epullgCYL)
real *weight_loc; /* Weights for the local indices */
int epgrppbc; /* The type of pbc for this pull group, see enum above */
atom_id pbcatom; /* The reference atom for pbc (global number) */
- rvec vec; /* The pull vector, direction or position */
- rvec init; /* Initial reference displacement */
- real rate; /* Rate of motion (nm/ps) */
- real k; /* force constant */
- real kB; /* force constant for state B */
+
+ /* Variables not present in mdp, but used at run time */
real wscale; /* scaling factor for the weights: sum w m/sum w w m */
real invtm; /* inverse total mass of the group: 1/wscale sum w m */
dvec x; /* center of mass before update */
dvec xp; /* center of mass after update before constraining */
+} t_pull_group;
+
+typedef struct {
+ int group[2]; /* The pull groups, index in group in t_pull */
+ rvec origin; /* The origin for the absolute reference */
+ rvec vec; /* The pull vector, direction or position */
+ real init; /* Initial reference displacement */
+ real rate; /* Rate of motion (nm/ps) */
+ real k; /* force constant */
+ real kB; /* force constant for state B */
+
+ /* Variables not present in mdp, but used at run time */
dvec dr; /* The distance from the reference group */
double f_scal; /* Scalar force for directional pulling */
dvec f; /* force due to the pulling/constraining */
-} t_pullgrp;
+} t_pull_coord;
typedef struct {
int eSimTempScale; /* simulated temperature scaling; linear or exponential */
} t_expanded;
typedef struct {
- int ngrp; /* number of groups */
+ int ngroup; /* number of pull groups */
+ int ncoord; /* number of pull coordinates */
int eGeom; /* pull geometry */
ivec dim; /* used to select components for constraint */
real cyl_r1; /* radius of cylinder for dynamic COM */
real cyl_r0; /* radius of cylinder including switch length */
real constr_tol; /* absolute tolerance for constraints in (nm) */
+ gmx_bool bPrintRef; /* Print coordinates of the first group in each coord */
int nstxout; /* Output frequency for pull x */
int nstfout; /* Output frequency for pull f */
int ePBC; /* the boundary conditions */
gmx_bool bRefAt; /* do we need reference atoms for a group COM ? */
int cosdim; /* dimension for cosine weighting, -1 if none */
gmx_bool bVirial; /* do we need to add the pull virial? */
- t_pullgrp *grp; /* groups to pull/restrain/etc/ */
- t_pullgrp *dyna; /* dynamic groups for use with local constraints */
+ t_pull_group *group; /* groups to pull/restrain/etc/ */
+ t_pull_coord *coord; /* the pull coordinates */
+
+ /* Variables not present in mdp, but used at run time */
+ t_pull_group *dyna; /* dynamic groups for use with local constraints */
rvec *rbuf; /* COM calculation buffer */
dvec *dbuf; /* COM calculation buffer */
double *dbuf_cyl; /* cylinder ref. groups COM calculation buffer */
#include "copyrite.h"
#include "macros.h"
-static void pull_print_x_grp(FILE *out, gmx_bool bRef, ivec dim, t_pullgrp *pgrp)
+static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp)
{
int m;
{
if (dim[m])
{
- fprintf(out, "\t%g", bRef ? pgrp->x[m] : pgrp->dr[m]);
+ fprintf(out, "\t%g", pgrp->x[m]);
}
}
}
-static void pull_print_x(FILE *out, t_pull *pull, double t)
+static void pull_print_coord_dr(FILE *out, ivec dim, const t_pull_coord *pcrd)
{
- int g;
-
- fprintf(out, "%.4f", t);
+ int m;
- if (PULL_CYL(pull))
+ for (m = 0; m < DIM; m++)
{
- for (g = 1; g < 1+pull->ngrp; g++)
+ if (dim[m])
{
- pull_print_x_grp(out, TRUE, pull->dim, &pull->dyna[g]);
- pull_print_x_grp(out, FALSE, pull->dim, &pull->grp[g]);
+ fprintf(out, "\t%g", pcrd->dr[m]);
}
}
- else
+}
+
+static void pull_print_x(FILE *out, t_pull *pull, double t)
+{
+ int c;
+ const t_pull_coord *pcrd;
+
+ fprintf(out, "%.4f", t);
+
+ for (c = 0; c < pull->ncoord; c++)
{
- for (g = 0; g < 1+pull->ngrp; g++)
+ pcrd = &pull->coord[c];
+
+ if (pull->bPrintRef)
{
- if (pull->grp[g].nat > 0)
+ if (PULL_CYL(pull))
{
- pull_print_x_grp(out, g == 0, pull->dim, &pull->grp[g]);
+ pull_print_group_x(out, pull->dim, &pull->dyna[c]);
+ }
+ else
+ {
+ pull_print_group_x(out, pull->dim, &pull->group[pcrd->group[0]]);
}
}
+ pull_print_coord_dr(out, pull->dim, pcrd);
}
fprintf(out, "\n");
}
static void pull_print_f(FILE *out, t_pull *pull, double t)
{
- int g, d;
+ int c, d;
fprintf(out, "%.4f", t);
- for (g = 1; g < 1+pull->ngrp; g++)
+ for (c = 0; c < pull->ncoord; c++)
{
- if (pull->eGeom == epullgPOS)
- {
- for (d = 0; d < DIM; d++)
- {
- if (pull->dim[d])
- {
- fprintf(out, "\t%g", pull->grp[g].f[d]);
- }
- }
- }
- else
- {
- fprintf(out, "\t%g", pull->grp[g].f_scal);
- }
+ fprintf(out, "\t%g", pull->coord[c].f_scal);
}
fprintf(out, "\n");
}
gmx_bool bCoord, unsigned long Flags)
{
FILE *fp;
- int nsets, g, m;
+ int nsets, c, m;
char **setname, buf[10];
if (Flags & MD_APPENDFILES)
exvggtXNY, oenv);
}
- snew(setname, (1+pull->ngrp)*DIM);
+ snew(setname, 2*pull->ncoord*DIM);
nsets = 0;
- for (g = 0; g < 1+pull->ngrp; g++)
+ for (c = 0; c < pull->ncoord; c++)
{
- if (pull->grp[g].nat > 0 &&
- (g > 0 || (bCoord && !PULL_CYL(pull))))
+ if (bCoord)
{
- if (bCoord || pull->eGeom == epullgPOS)
+ if (pull->bPrintRef)
{
- if (PULL_CYL(pull))
- {
- for (m = 0; m < DIM; m++)
- {
- if (pull->dim[m])
- {
- sprintf(buf, "%d %s%c", g, "c", 'X'+m);
- setname[nsets] = strdup(buf);
- nsets++;
- }
- }
- }
for (m = 0; m < DIM; m++)
{
if (pull->dim[m])
{
- sprintf(buf, "%d %s%c",
- g, (bCoord && g > 0) ? "d" : "", 'X'+m);
+ sprintf(buf, "%d %s%c", c+1, "c", 'X'+m);
setname[nsets] = strdup(buf);
nsets++;
}
}
}
- else
+ for (m = 0; m < DIM; m++)
{
- sprintf(buf, "%d", g);
- setname[nsets] = strdup(buf);
- nsets++;
+ if (pull->dim[m])
+ {
+ sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
+ setname[nsets] = strdup(buf);
+ nsets++;
+ }
}
}
+ else
+ {
+ sprintf(buf, "%d", c+1);
+ setname[nsets] = strdup(buf);
+ nsets++;
+ }
}
- if (bCoord || nsets > 1)
+ if (nsets > 1)
{
xvgr_legend(fp, nsets, (const char**)setname, oenv);
}
- for (g = 0; g < nsets; g++)
+ for (c = 0; c < nsets; c++)
{
- sfree(setname[g]);
+ sfree(setname[c]);
}
sfree(setname);
}
}
/* Apply forces in a mass weighted fashion */
-static void apply_forces_grp(t_pullgrp *pgrp, t_mdatoms * md,
- dvec f_pull, int sign, rvec *f)
+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;
double wmass, inv_wm;
/* Apply forces in a mass weighted fashion */
static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f)
{
- int i;
- t_pullgrp *pgrp;
+ int c;
+ const t_pull_coord *pcrd;
- for (i = 1; i < pull->ngrp+1; i++)
+ for (c = 0; c < pull->ncoord; c++)
{
- pgrp = &(pull->grp[i]);
- apply_forces_grp(pgrp, md, pgrp->f, 1, f);
- if (pull->grp[0].nat)
+ pcrd = &pull->coord[c];
+
+ if (PULL_CYL(pull))
{
- if (PULL_CYL(pull))
- {
- apply_forces_grp(&(pull->dyna[i]), md, pgrp->f, -1, f);
- }
- else
+ apply_forces_grp(&pull->dyna[c], md, pcrd->f, -1, f);
+ }
+ else
+ {
+ if (pull->group[pcrd->group[0]].nat > 0)
{
- apply_forces_grp(&(pull->grp[0]), md, pgrp->f, -1, f);
+ apply_forces_grp(&pull->group[pcrd->group[0]], md, pcrd->f, -1, f);
}
}
+ apply_forces_grp(&pull->group[pcrd->group[1]], md, pcrd->f, 1, f);
}
}
return 0.25*max_d2;
}
-static void get_pullgrps_dr(const t_pull *pull, const t_pbc *pbc, int g, double t,
- dvec xg, dvec xref, double max_dist2,
- dvec dr)
+static void low_get_pull_coord_dr(const t_pull *pull,
+ const t_pull_coord *pcrd,
+ const t_pbc *pbc, double t,
+ dvec xg, dvec xref, double max_dist2,
+ dvec dr)
{
- t_pullgrp *pref, *pgrp;
- int m;
- dvec xrefr, dref = {0, 0, 0};
- double dr2;
+ const t_pull_group *pgrp0, *pgrp1;
+ int m;
+ dvec xrefr, dref = {0, 0, 0};
+ double dr2;
+
+ pgrp0 = &pull->group[pcrd->group[0]];
+ pgrp1 = &pull->group[pcrd->group[1]];
- pgrp = &pull->grp[g];
+ /* Only the first group can be an absolute reference, in that case nat=0 */
+ if (pgrp0->nat == 0)
+ {
+ for (m = 0; m < DIM; m++)
+ {
+ xref[m] = pcrd->origin[m];
+ }
+ }
copy_dvec(xref, xrefr);
{
for (m = 0; m < DIM; m++)
{
- dref[m] = (pgrp->init[0] + pgrp->rate*t)*pull->grp[g].vec[m];
+ dref[m] = (pcrd->init + pcrd->rate*t)*pcrd->vec[m];
}
/* Add the reference position, so we use the correct periodic image */
dvec_inc(xrefr, dref);
}
if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
{
- gmx_fatal(FARGS, "Distance of pull group %d (%f nm) is larger than 0.49 times the box size (%f)", g, sqrt(dr2), sqrt(max_dist2));
+ gmx_fatal(FARGS, "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the box size (%f)",
+ pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2));
}
if (pull->eGeom == epullgDIRPBC)
}
}
-static void get_pullgrp_dr(const t_pull *pull, const t_pbc *pbc, int g, double t,
- dvec dr)
+static void get_pull_coord_dr(const t_pull *pull,
+ int coord_ind,
+ const t_pbc *pbc, double t,
+ dvec dr)
{
- double md2;
+ double md2;
+ const t_pull_coord *pcrd;
if (pull->eGeom == epullgDIRPBC)
{
md2 = max_pull_distance2(pull, pbc);
}
- get_pullgrps_dr(pull, pbc, g, t,
- pull->grp[g].x,
- PULL_CYL(pull) ? pull->dyna[g].x : pull->grp[0].x,
- md2,
- dr);
+ pcrd = &pull->coord[coord_ind];
+
+ low_get_pull_coord_dr(pull, pcrd, pbc, t,
+ pull->group[pcrd->group[1]].x,
+ PULL_CYL(pull) ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x,
+ md2,
+ dr);
}
-void get_pullgrp_distance(t_pull *pull, t_pbc *pbc, int g, double t,
- dvec dr, dvec dev)
+void get_pull_coord_distance(const t_pull *pull,
+ int coord_ind,
+ const t_pbc *pbc, double t,
+ dvec dr, double *dev)
{
static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety,
but is fairly benign */
- t_pullgrp *pgrp;
- int m;
- dvec ref;
- double drs, inpr;
+ const t_pull_coord *pcrd;
+ int m;
+ double ref, drs, inpr;
- pgrp = &pull->grp[g];
+ pcrd = &pull->coord[coord_ind];
- get_pullgrp_dr(pull, pbc, g, t, dr);
+ get_pull_coord_dr(pull, coord_ind, pbc, t, dr);
- if (pull->eGeom == epullgPOS)
- {
- for (m = 0; m < DIM; m++)
- {
- ref[m] = pgrp->init[m] + pgrp->rate*t*pgrp->vec[m];
- }
- }
- else
- {
- ref[0] = pgrp->init[0] + pgrp->rate*t;
- }
+ ref = pcrd->init + pcrd->rate*t;
switch (pull->eGeom)
{
case epullgDIST:
/* Pull along the vector between the com's */
- if (ref[0] < 0 && !bWarned)
+ if (ref < 0 && !bWarned)
{
- fprintf(stderr, "\nPull reference distance for group %d is negative (%f)\n", g, ref[0]);
+ fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, ref);
bWarned = TRUE;
}
drs = dnorm(dr);
/* With no vector we can not determine the direction for the force,
* so we set the force to zero.
*/
- dev[0] = 0;
+ *dev = 0;
}
else
{
/* Determine the deviation */
- dev[0] = drs - ref[0];
+ *dev = drs - ref;
}
break;
case epullgDIR:
inpr = 0;
for (m = 0; m < DIM; m++)
{
- inpr += pgrp->vec[m]*dr[m];
- }
- dev[0] = inpr - ref[0];
- break;
- case epullgPOS:
- /* Determine the difference of dr and ref along each dimension */
- for (m = 0; m < DIM; m++)
- {
- dev[m] = (dr[m] - ref[m])*pull->dim[m];
+ inpr += pcrd->vec[m]*dr[m];
}
+ *dev = inpr - ref;
break;
}
}
* It can happen that multiple constraint steps need to be applied
* and therefore the constraint forces need to be accumulated.
*/
- for (i = 0; i < 1+pull->ngrp; i++)
+ for (i = 0; i < pull->ncoord; i++)
{
- clear_dvec(pull->grp[i].f);
- pull->grp[i].f_scal = 0;
+ clear_dvec(pull->coord[i].f);
+ pull->coord[i].f_scal = 0;
}
}
dvec *r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */
dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */
-
- dvec *rinew; /* current 'new' position of group i */
- dvec *rjnew; /* current 'new' position of group j */
- dvec ref, vec;
+ dvec *rnew; /* current 'new' positions of the groups */
+ double *dr_tot; /* the total update of the coords */
+ double ref;
+ dvec vec;
double d0, inpr;
double lambda, rm, mass, invdt = 0;
gmx_bool bConverged_all, bConverged = FALSE;
- int niter = 0, g, ii, j, m, max_iter = 100;
- double q, a, b, c; /* for solving the quadratic equation,
- see Num. Recipes in C ed 2 p. 184 */
- dvec *dr; /* correction for group i */
- dvec ref_dr; /* correction for group j */
+ int niter = 0, g, c, ii, j, m, max_iter = 100;
+ double a;
dvec f; /* the pull force */
dvec tmp, tmp3;
- t_pullgrp *pdyna, *pgrp, *pref;
+ t_pull_group *pdyna, *pgrp0, *pgrp1;
+ t_pull_coord *pcrd;
- snew(r_ij, pull->ngrp+1);
- if (PULL_CYL(pull))
- {
- snew(rjnew, pull->ngrp+1);
- }
- else
- {
- snew(rjnew, 1);
- }
- snew(dr, pull->ngrp+1);
- snew(rinew, pull->ngrp+1);
+ snew(r_ij, pull->ncoord);
+ snew(dr_tot, pull->ncoord);
+
+ snew(rnew, pull->ngroup);
/* copy the current unconstrained positions for use in iterations. We
iterate until rinew[i] and rjnew[j] obey the constraints. Then
rinew - pull.x_unc[i] is the correction dr to group i */
- for (g = 1; g < 1+pull->ngrp; g++)
+ for (g = 0; g < pull->ngroup; g++)
{
- copy_dvec(pull->grp[g].xp, rinew[g]);
+ copy_dvec(pull->group[g].xp, rnew[g]);
}
if (PULL_CYL(pull))
{
- for (g = 1; g < 1+pull->ngrp; g++)
- {
- copy_dvec(pull->dyna[g].xp, rjnew[g]);
- }
- }
- else
- {
- copy_dvec(pull->grp[0].xp, rjnew[0]);
+ /* There is only one pull coordinate and reference group */
+ copy_dvec(pull->dyna[0].xp, rnew[pull->coord[0].group[0]]);
}
/* Determine the constraint directions from the old positions */
- for (g = 1; g < 1+pull->ngrp; g++)
+ for (c = 0; c < pull->ncoord; c++)
{
- get_pullgrp_dr(pull, pbc, g, t, r_ij[g]);
+ get_pull_coord_dr(pull, c, pbc, t, r_ij[c]);
/* Store the difference vector at time t for printing */
- copy_dvec(r_ij[g], pull->grp[g].dr);
+ copy_dvec(r_ij[c], pull->coord[c].dr);
if (debug)
{
- fprintf(debug, "Pull group %d dr %f %f %f\n",
- g, r_ij[g][XX], r_ij[g][YY], r_ij[g][ZZ]);
+ fprintf(debug, "Pull coord %d dr %f %f %f\n",
+ c, r_ij[c][XX], r_ij[c][YY], r_ij[c][ZZ]);
}
if (pull->eGeom == epullgDIR || pull->eGeom == epullgDIRPBC)
a = 0;
for (m = 0; m < DIM; m++)
{
- a += pull->grp[g].vec[m]*r_ij[g][m];
+ a += pull->coord[c].vec[m]*r_ij[c][m];
}
for (m = 0; m < DIM; m++)
{
- r_ij[g][m] = a*pull->grp[g].vec[m];
+ r_ij[c][m] = a*pull->coord[c].vec[m];
}
}
}
bConverged_all = TRUE;
/* loop over all constraints */
- for (g = 1; g < 1+pull->ngrp; g++)
+ for (c = 0; c < pull->ncoord; c++)
{
- pgrp = &pull->grp[g];
- if (PULL_CYL(pull))
- {
- pref = &pull->dyna[g];
- }
- else
- {
- pref = &pull->grp[0];
- }
+ dvec dr0, dr1;
+
+ pcrd = &pull->coord[c];
+ pgrp0 = &pull->group[pcrd->group[0]];
+ pgrp1 = &pull->group[pcrd->group[1]];
/* Get the current difference vector */
- get_pullgrps_dr(pull, pbc, g, t, rinew[g], rjnew[PULL_CYL(pull) ? g : 0],
- -1, unc_ij);
+ low_get_pull_coord_dr(pull, pcrd, pbc, t,
+ rnew[pcrd->group[1]],
+ rnew[pcrd->group[0]],
+ -1, unc_ij);
- if (pull->eGeom == epullgPOS)
- {
- for (m = 0; m < DIM; m++)
- {
- ref[m] = pgrp->init[m] + pgrp->rate*t*pgrp->vec[m];
- }
- }
- else
- {
- ref[0] = pgrp->init[0] + pgrp->rate*t;
- /* Keep the compiler happy */
- ref[1] = 0;
- ref[2] = 0;
- }
+ ref = pcrd->init + pcrd->rate*t;
if (debug)
{
- fprintf(debug, "Pull group %d, iteration %d\n", g, niter);
+ fprintf(debug, "Pull coord %d, iteration %d\n", c, niter);
}
- rm = 1.0/(pull->grp[g].invtm + pref->invtm);
+ rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
switch (pull->eGeom)
{
case epullgDIST:
- if (ref[0] <= 0)
+ if (ref <= 0)
{
- gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", g, ref[0]);
+ gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, ref);
}
- a = diprod(r_ij[g], r_ij[g]);
- b = diprod(unc_ij, r_ij[g])*2;
- c = diprod(unc_ij, unc_ij) - dsqr(ref[0]);
-
- if (b < 0)
- {
- q = -0.5*(b - sqrt(b*b - 4*a*c));
- lambda = -q/a;
- }
- else
{
- q = -0.5*(b + sqrt(b*b - 4*a*c));
- lambda = -c/q;
- }
+ double q, c_a, c_b, c_c;
- if (debug)
- {
- fprintf(debug,
- "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
- a, b, c, lambda);
+ c_a = diprod(r_ij[c], r_ij[c]);
+ c_b = diprod(unc_ij, r_ij[c])*2;
+ c_c = diprod(unc_ij, unc_ij) - dsqr(ref);
+
+ if (c_b < 0)
+ {
+ q = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c));
+ lambda = -q/c_a;
+ }
+ else
+ {
+ q = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c));
+ lambda = -c_c/q;
+ }
+
+ if (debug)
+ {
+ fprintf(debug,
+ "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n",
+ c_a, c_b, c_c, lambda);
+ }
}
/* The position corrections dr due to the constraints */
- dsvmul(-lambda*rm*pgrp->invtm, r_ij[g], dr[g]);
- dsvmul( lambda*rm*pref->invtm, r_ij[g], ref_dr);
+ dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1);
+ dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0);
+ dr_tot[c] += -lambda*dnorm(r_ij[c]);
break;
case epullgDIR:
case epullgDIRPBC:
a = 0;
for (m = 0; m < DIM; m++)
{
- vec[m] = pgrp->vec[m];
+ vec[m] = pcrd->vec[m];
a += unc_ij[m]*vec[m];
}
/* Select only the component along the vector */
dsvmul(a, vec, unc_ij);
- lambda = a - ref[0];
+ lambda = a - ref;
if (debug)
{
fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda);
}
/* The position corrections dr due to the constraints */
- dsvmul(-lambda*rm*pull->grp[g].invtm, vec, dr[g]);
- dsvmul( lambda*rm* pref->invtm, vec, ref_dr);
- break;
- case epullgPOS:
- for (m = 0; m < DIM; m++)
- {
- if (pull->dim[m])
- {
- lambda = r_ij[g][m] - ref[m];
- /* The position corrections dr due to the constraints */
- dr[g][m] = -lambda*rm*pull->grp[g].invtm;
- ref_dr[m] = lambda*rm*pref->invtm;
- }
- else
- {
- dr[g][m] = 0;
- ref_dr[m] = 0;
- }
- }
+ dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1);
+ dsvmul( lambda*rm*pgrp0->invtm, vec, dr0);
+ dr_tot[c] += -lambda;
break;
}
/* DEBUG */
if (debug)
{
- j = (PULL_CYL(pull) ? g : 0);
- get_pullgrps_dr(pull, pbc, g, t, rinew[g], rjnew[j], -1, tmp);
- get_pullgrps_dr(pull, pbc, g, t, dr[g], ref_dr, -1, tmp3);
+ int g0, g1;
+
+ g0 = pcrd->group[0];
+ g1 = pcrd->group[1];
+ low_get_pull_coord_dr(pull, pcrd, pbc, t, rnew[g1], rnew[g0], -1, tmp);
+ low_get_pull_coord_dr(pull, pcrd, pbc, t, dr1, dr0, -1, tmp3);
fprintf(debug,
"Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
- rinew[g][0], rinew[g][1], rinew[g][2],
- rjnew[j][0], rjnew[j][1], rjnew[j][2], dnorm(tmp));
- if (pull->eGeom == epullgPOS)
- {
- fprintf(debug,
- "Pull ref %8.5f %8.5f %8.5f\n",
- pgrp->vec[0], pgrp->vec[1], pgrp->vec[2]);
- }
- else
- {
- fprintf(debug,
- "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f %8.5f %8.5f\n",
- "", "", "", "", "", "", ref[0], ref[1], ref[2]);
- }
+ rnew[g0][0], rnew[g0][1], rnew[g0][2],
+ rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp));
+ fprintf(debug,
+ "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n",
+ "", "", "", "", "", "", ref);
fprintf(debug,
"Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n",
- dr[g][0], dr[g][1], dr[g][2],
- ref_dr[0], ref_dr[1], ref_dr[2],
+ dr0[0], dr0[1], dr0[2],
+ dr1[0], dr1[1], dr1[2],
dnorm(tmp3));
- fprintf(debug,
- "Pull cor %10.7f %10.7f %10.7f\n",
- dr[g][0], dr[g][1], dr[g][2]);
} /* END DEBUG */
/* Update the COMs with dr */
- dvec_inc(rinew[g], dr[g]);
- dvec_inc(rjnew[PULL_CYL(pull) ? g : 0], ref_dr);
+ dvec_inc(rnew[pcrd->group[1]], dr1);
+ dvec_inc(rnew[pcrd->group[0]], dr0);
}
/* Check if all constraints are fullfilled now */
- for (g = 1; g < 1+pull->ngrp; g++)
+ for (c = 0; c < pull->ncoord; c++)
{
- pgrp = &pull->grp[g];
+ pcrd = &pull->coord[c];
- get_pullgrps_dr(pull, pbc, g, t, rinew[g], rjnew[PULL_CYL(pull) ? g : 0],
- -1, unc_ij);
+ low_get_pull_coord_dr(pull, pcrd, pbc, t,
+ rnew[pcrd->group[1]],
+ rnew[pcrd->group[0]],
+ -1, unc_ij);
switch (pull->eGeom)
{
case epullgDIST:
- bConverged = fabs(dnorm(unc_ij) - ref[0]) < pull->constr_tol;
+ bConverged = fabs(dnorm(unc_ij) - ref) < pull->constr_tol;
break;
case epullgDIR:
case epullgDIRPBC:
case epullgCYL:
for (m = 0; m < DIM; m++)
{
- vec[m] = pgrp->vec[m];
+ vec[m] = pcrd->vec[m];
}
inpr = diprod(unc_ij, vec);
dsvmul(inpr, vec, unc_ij);
bConverged =
- fabs(diprod(unc_ij, vec) - ref[0]) < pull->constr_tol;
- break;
- case epullgPOS:
- bConverged = TRUE;
- for (m = 0; m < DIM; m++)
- {
- if (pull->dim[m] &&
- fabs(unc_ij[m] - ref[m]) >= pull->constr_tol)
- {
- bConverged = FALSE;
- }
- }
+ fabs(diprod(unc_ij, vec) - ref) < pull->constr_tol;
break;
}
if (debug)
{
fprintf(debug, "NOT CONVERGED YET: Group %d:"
- "d_ref = %f %f %f, current d = %f\n",
- g, ref[0], ref[1], ref[2], dnorm(unc_ij));
+ "d_ref = %f, current d = %f\n",
+ g, ref, dnorm(unc_ij));
}
bConverged_all = FALSE;
invdt = 1/dt;
}
- /* update the normal groups */
- for (g = 1; g < 1+pull->ngrp; g++)
+ /* update atoms in the groups */
+ for (g = 0; g < pull->ngroup; g++)
{
- pgrp = &pull->grp[g];
- /* get the final dr and constraint force for group i */
- dvec_sub(rinew[g], pgrp->xp, dr[g]);
- /* select components of dr */
- for (m = 0; m < DIM; m++)
+ const t_pull_group *pgrp;
+ dvec dr;
+
+ if (PULL_CYL(pull) && g == pull->coord[0].group[0])
{
- dr[g][m] *= pull->dim[m];
+ pgrp = &pull->dyna[0];
}
- dsvmul(1.0/(pgrp->invtm*dt*dt), dr[g], f);
- dvec_inc(pgrp->f, f);
- switch (pull->eGeom)
+ else
{
- case epullgDIST:
- for (m = 0; m < DIM; m++)
- {
- pgrp->f_scal += r_ij[g][m]*f[m]/dnorm(r_ij[g]);
- }
- break;
- case epullgDIR:
- case epullgDIRPBC:
- case epullgCYL:
- for (m = 0; m < DIM; m++)
- {
- pgrp->f_scal += pgrp->vec[m]*f[m];
- }
- break;
- case epullgPOS:
- break;
+ pgrp = &pull->group[g];
}
- if (vir && bMaster)
+ /* get the final constraint displacement dr for group g */
+ dvec_sub(rnew[g], pgrp->xp, dr);
+ /* select components of dr */
+ for (m = 0; m < DIM; m++)
{
- /* Add the pull contribution to the virial */
- for (j = 0; j < DIM; j++)
- {
- for (m = 0; m < DIM; m++)
- {
- vir[j][m] -= 0.5*f[j]*r_ij[g][m];
- }
- }
+ dr[m] *= pull->dim[m];
}
/* update the atom positions */
- copy_dvec(dr[g], tmp);
+ copy_dvec(dr, tmp);
for (j = 0; j < pgrp->nat_loc; j++)
{
ii = pgrp->ind_loc[j];
if (pgrp->weight_loc)
{
- dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr[g], tmp);
+ dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp);
}
for (m = 0; m < DIM; m++)
{
}
}
- /* update the reference groups */
- if (PULL_CYL(pull))
+ /* calculate the constraint forces, used for output and virial only */
+ for (c = 0; c < pull->ncoord; c++)
{
- /* update the dynamic reference groups */
- for (g = 1; g < 1+pull->ngrp; g++)
- {
- pdyna = &pull->dyna[g];
- dvec_sub(rjnew[g], pdyna->xp, ref_dr);
- /* select components of ref_dr */
- for (m = 0; m < DIM; m++)
- {
- ref_dr[m] *= pull->dim[m];
- }
+ pcrd = &pull->coord[c];
+ pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->group[0]].invtm + pull->group[pcrd->group[1]].invtm)*dt*dt);
- for (j = 0; j < pdyna->nat_loc; j++)
- {
- /* reset the atoms with dr, weighted by w_i */
- dsvmul(pdyna->wscale*pdyna->weight_loc[j], ref_dr, tmp);
- ii = pdyna->ind_loc[j];
- for (m = 0; m < DIM; m++)
- {
- x[ii][m] += tmp[m];
- }
- if (v)
- {
- for (m = 0; m < DIM; m++)
- {
- v[ii][m] += invdt*tmp[m];
- }
- }
- }
- }
- }
- else
- {
- pgrp = &pull->grp[0];
- /* update the reference group */
- dvec_sub(rjnew[0], pgrp->xp, ref_dr);
- /* select components of ref_dr */
- for (m = 0; m < DIM; m++)
+ if (vir && bMaster)
{
- ref_dr[m] *= pull->dim[m];
- }
+ double f_invr;
- copy_dvec(ref_dr, tmp);
- for (j = 0; j < pgrp->nat_loc; j++)
- {
- ii = pgrp->ind_loc[j];
- if (pgrp->weight_loc)
- {
- dsvmul(pgrp->wscale*pgrp->weight_loc[j], ref_dr, tmp);
- }
- for (m = 0; m < DIM; m++)
- {
- x[ii][m] += tmp[m];
- }
- if (v)
+ /* Add the pull contribution to the virial */
+ f_invr = pcrd->f_scal/dnorm(r_ij[c]);
+
+ for (j = 0; j < DIM; j++)
{
for (m = 0; m < DIM; m++)
{
- v[ii][m] += invdt*tmp[m];
+ vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m];
}
}
}
/* finished! I hope. Give back some memory */
sfree(r_ij);
- sfree(rinew);
- sfree(rjnew);
- sfree(dr);
+ sfree(dr_tot);
+ sfree(rnew);
}
/* Pulling with a harmonic umbrella potential or constant force */
t_pull *pull, t_pbc *pbc, double t, real lambda,
real *V, tensor vir, real *dVdl)
{
- int g, j, m;
- dvec dev;
- double ndr, invdr;
- real k, dkdl;
- t_pullgrp *pgrp;
+ int c, j, m;
+ double dev, ndr, invdr;
+ real k, dkdl;
+ t_pull_coord *pcrd;
- /* loop over the groups that are being pulled */
+ /* loop over the pull coordinates */
*V = 0;
*dVdl = 0;
- for (g = 1; g < 1+pull->ngrp; g++)
+ for (c = 0; c < pull->ncoord; c++)
{
- pgrp = &pull->grp[g];
- get_pullgrp_distance(pull, pbc, g, t, pgrp->dr, dev);
+ pcrd = &pull->coord[c];
+
+ get_pull_coord_distance(pull, c, pbc, t, pcrd->dr, &dev);
- k = (1.0 - lambda)*pgrp->k + lambda*pgrp->kB;
- dkdl = pgrp->kB - pgrp->k;
+ k = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB;
+ dkdl = pcrd->kB - pcrd->k;
switch (pull->eGeom)
{
case epullgDIST:
- ndr = dnorm(pgrp->dr);
+ ndr = dnorm(pcrd->dr);
invdr = 1/ndr;
if (ePull == epullUMBRELLA)
{
- pgrp->f_scal = -k*dev[0];
- *V += 0.5* k*dsqr(dev[0]);
- *dVdl += 0.5*dkdl*dsqr(dev[0]);
+ pcrd->f_scal = -k*dev;
+ *V += 0.5* k*dsqr(dev);
+ *dVdl += 0.5*dkdl*dsqr(dev);
}
else
{
- pgrp->f_scal = -k;
+ pcrd->f_scal = -k;
*V += k*ndr;
*dVdl += dkdl*ndr;
}
for (m = 0; m < DIM; m++)
{
- pgrp->f[m] = pgrp->f_scal*pgrp->dr[m]*invdr;
+ pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
}
break;
case epullgDIR:
case epullgCYL:
if (ePull == epullUMBRELLA)
{
- pgrp->f_scal = -k*dev[0];
- *V += 0.5* k*dsqr(dev[0]);
- *dVdl += 0.5*dkdl*dsqr(dev[0]);
+ pcrd->f_scal = -k*dev;
+ *V += 0.5* k*dsqr(dev);
+ *dVdl += 0.5*dkdl*dsqr(dev);
}
else
{
ndr = 0;
for (m = 0; m < DIM; m++)
{
- ndr += pgrp->vec[m]*pgrp->dr[m];
+ ndr += pcrd->vec[m]*pcrd->dr[m];
}
- pgrp->f_scal = -k;
+ pcrd->f_scal = -k;
*V += k*ndr;
*dVdl += dkdl*ndr;
}
for (m = 0; m < DIM; m++)
{
- pgrp->f[m] = pgrp->f_scal*pgrp->vec[m];
- }
- break;
- case epullgPOS:
- for (m = 0; m < DIM; m++)
- {
- if (ePull == epullUMBRELLA)
- {
- pgrp->f[m] = -k*dev[m];
- *V += 0.5* k*dsqr(dev[m]);
- *dVdl += 0.5*dkdl*dsqr(dev[m]);
- }
- else
- {
- pgrp->f[m] = -k*pull->dim[m];
- *V += k*pgrp->dr[m]*pull->dim[m];
- *dVdl += dkdl*pgrp->dr[m]*pull->dim[m];
- }
+ pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
}
break;
}
{
for (m = 0; m < DIM; m++)
{
- vir[j][m] -= 0.5*pgrp->f[j]*pgrp->dr[m];
+ vir[j][m] -= 0.5*pcrd->f[j]*pcrd->dr[m];
}
}
}
}
static void make_local_pull_group(gmx_ga2la_t ga2la,
- t_pullgrp *pg, int start, int end)
+ t_pull_group *pg, int start, int end)
{
int i, ii;
ga2la = NULL;
}
- if (pull->grp[0].nat > 0)
+ for (g = 0; g < pull->ngroup; g++)
{
- make_local_pull_group(ga2la, &pull->grp[0], md->start, md->start+md->homenr);
- }
- for (g = 1; g < 1+pull->ngrp; g++)
- {
- make_local_pull_group(ga2la, &pull->grp[g], md->start, md->start+md->homenr);
+ make_local_pull_group(ga2la, &pull->group[g],
+ md->start, md->start+md->homenr);
}
}
static void init_pull_group_index(FILE *fplog, t_commrec *cr,
int start, int end,
- int g, t_pullgrp *pg, ivec pulldims,
+ int g, t_pull_group *pg, ivec pulldims,
gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
{
int i, ii, d, nfrozen, ndim;
gmx_bool bOutFile, unsigned long Flags)
{
t_pull *pull;
- t_pullgrp *pgrp;
- int g, start = 0, end = 0, m;
- gmx_bool bCite;
+ t_pull_group *pgrp;
+ int c, g, start = 0, end = 0, m;
pull = ir->pull;
if (fplog)
{
- fprintf(fplog, "\nWill apply %s COM pulling in geometry '%s'\n",
- EPULLTYPE(ir->ePull), EPULLGEOM(pull->eGeom));
- if (pull->grp[0].nat > 0)
+ gmx_bool bAbs, bCos;
+
+ bAbs = FALSE;
+ for (c = 0; c < pull->ncoord; c++)
{
- fprintf(fplog, "between a reference group and %d group%s\n",
- pull->ngrp, pull->ngrp == 1 ? "" : "s");
+ if (pull->group[pull->coord[c].group[0]].nat == 0 ||
+ pull->group[pull->coord[c].group[1]].nat == 0)
+ {
+ bAbs = TRUE;
+ }
}
- else
+
+ fprintf(fplog, "\nWill apply %s COM pulling in geometry '%s'\n",
+ EPULLTYPE(ir->ePull), EPULLGEOM(pull->eGeom));
+ fprintf(fplog, "with %d pull coordinate%s and %d group%s\n",
+ pull->ncoord, pull->ncoord == 1 ? "" : "s",
+ pull->ngroup, pull->ngroup == 1 ? "" : "s");
+ if (bAbs)
{
- fprintf(fplog, "with an absolute reference on %d group%s\n",
- pull->ngrp, pull->ngrp == 1 ? "" : "s");
+ fprintf(fplog, "with an absolute reference\n");
}
- bCite = FALSE;
- for (g = 0; g < pull->ngrp+1; g++)
+ bCos = FALSE;
+ for (g = 0; g < pull->ngroup; g++)
{
- if (pull->grp[g].nat > 1 &&
- pull->grp[g].pbcatom < 0)
+ if (pull->group[g].nat > 1 &&
+ pull->group[g].pbcatom < 0)
{
/* We are using cosine weighting */
fprintf(fplog, "Cosine weighting is used for group %d\n", g);
- bCite = TRUE;
+ bCos = TRUE;
}
}
- if (bCite)
+ if (bCos)
{
please_cite(fplog, "Engin2010");
}
pull->dbuf_cyl = NULL;
pull->bRefAt = FALSE;
pull->cosdim = -1;
- for (g = 0; g < pull->ngrp+1; g++)
+ for (g = 0; g < pull->ngroup; g++)
{
- pgrp = &pull->grp[g];
+ pgrp = &pull->group[g];
pgrp->epgrppbc = epgrppbcNONE;
if (pgrp->nat > 0)
{
/* if we use dynamic reference groups, do some initialising for them */
if (PULL_CYL(pull))
{
- if (pull->grp[0].nat == 0)
+ if (ir->ePull == epullCONSTRAINT && pull->ncoord > 1)
{
- gmx_fatal(FARGS, "Dynamic reference groups are not supported when using absolute reference!\n");
+ /* We can't easily update the single reference group with multiple
+ * constraints. This would require recalculating COMs.
+ */
+ gmx_fatal(FARGS, "Constraint COM pulling supports only one coordinate with geometry=cylinder, you can use umbrella pulling with multiple coordinates");
}
- snew(pull->dyna, pull->ngrp+1);
+
+ for (c = 0; c < pull->ncoord; c++)
+ {
+ if (pull->group[pull->coord[c].group[0]].nat == 0)
+ {
+ gmx_fatal(FARGS, "Dynamic reference groups are not supported when using absolute reference!\n");
+ }
+ }
+
+ snew(pull->dyna, pull->ncoord);
}
/* Only do I/O when we are doing dynamics and if we are the MASTER */
#include "pull.h"
#include "gmx_ga2la.h"
-static void pull_set_pbcatom(t_commrec *cr, t_pullgrp *pg,
+static void pull_set_pbcatom(t_commrec *cr, t_pull_group *pgrp,
t_mdatoms *md, rvec *x,
rvec x_pbc)
{
{
if (DOMAINDECOMP(cr))
{
- if (!ga2la_get_home(cr->dd->ga2la, pg->pbcatom, &a))
+ if (!ga2la_get_home(cr->dd->ga2la, pgrp->pbcatom, &a))
{
a = -1;
}
}
else
{
- a = pg->pbcatom;
+ a = pgrp->pbcatom;
}
if (a >= md->start && a < md->start+md->homenr)
}
else
{
- copy_rvec(x[pg->pbcatom], x_pbc);
+ copy_rvec(x[pgrp->pbcatom], x_pbc);
}
}
int g, n, m;
n = 0;
- for (g = 0; g < 1+pull->ngrp; g++)
+ for (g = 0; g < pull->ngroup; g++)
{
- if ((g == 0 && PULL_CYL(pull)) || pull->grp[g].pbcatom == -1)
+ if ((g == 0 && PULL_CYL(pull)) || pull->group[g].pbcatom == -1)
{
clear_rvec(x_pbc[g]);
}
else
{
- pull_set_pbcatom(cr, &pull->grp[g], md, x, x_pbc[g]);
+ pull_set_pbcatom(cr, &pull->group[g], md, x, x_pbc[g]);
for (m = 0; m < DIM; m++)
{
if (pull->dim[m] == 0)
if (cr && PAR(cr) && n > 0)
{
/* Sum over the nodes to get x_pbc from the home node of pbcatom */
- gmx_sum((1+pull->ngrp)*DIM, x_pbc[0], cr);
+ gmx_sum(pull->ngroup*DIM, x_pbc[0], cr);
}
}
static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md,
t_pbc *pbc, double t, rvec *x, rvec *xp)
{
- int g, i, ii, m, start, end;
- rvec g_x, dx, dir;
- double r0_2, sum_a, sum_ap, dr2, mass, weight, wmass, wwmass, inp;
- t_pullgrp *pref, *pgrp, *pdyna;
- gmx_ga2la_t ga2la = NULL;
+ int c, i, ii, m, start, end;
+ rvec g_x, dx, dir;
+ double r0_2, sum_a, sum_ap, dr2, mass, weight, wmass, wwmass, inp;
+ t_pull_coord *pcrd;
+ t_pull_group *pref, *pgrp, *pdyna;
+ gmx_ga2la_t ga2la = NULL;
if (pull->dbuf_cyl == NULL)
{
- snew(pull->dbuf_cyl, pull->ngrp*4);
+ snew(pull->dbuf_cyl, pull->ncoord*4);
}
if (cr && DOMAINDECOMP(cr))
r0_2 = dsqr(pull->cyl_r0);
/* loop over all groups to make a reference group for each*/
- pref = &pull->grp[0];
- for (g = 1; g < 1+pull->ngrp; g++)
+ for (c = 0; c < pull->ncoord; c++)
{
- pgrp = &pull->grp[g];
- pdyna = &pull->dyna[g];
- copy_rvec(pgrp->vec, dir);
+ pcrd = &pull->coord[c];
+
+ /* pref will be the same group for all pull coordinates */
+ pref = &pull->group[pcrd->group[0]];
+ pgrp = &pull->group[pcrd->group[1]];
+ pdyna = &pull->dyna[c];
+ copy_rvec(pcrd->vec, dir);
sum_a = 0;
sum_ap = 0;
wmass = 0;
for (m = 0; m < DIM; m++)
{
- g_x[m] = pgrp->x[m] - pgrp->vec[m]*(pgrp->init[0] + pgrp->rate*t);
+ g_x[m] = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t);
}
/* loop over all atoms in the main ref group */
for (i = 0; i < pref->nat; i++)
{
- ii = pull->grp[0].ind[i];
+ ii = pref->ind[i];
if (ga2la)
{
if (!ga2la_get_home(ga2la, pref->ind[i], &ii))
}
}
}
- pull->dbuf_cyl[(g-1)*4+0] = wmass;
- pull->dbuf_cyl[(g-1)*4+1] = wwmass;
- pull->dbuf_cyl[(g-1)*4+2] = sum_a;
- pull->dbuf_cyl[(g-1)*4+3] = sum_ap;
+ pull->dbuf_cyl[c*4+0] = wmass;
+ pull->dbuf_cyl[c*4+1] = wwmass;
+ pull->dbuf_cyl[c*4+2] = sum_a;
+ pull->dbuf_cyl[c*4+3] = sum_ap;
}
if (cr && PAR(cr))
{
/* Sum the contributions over the nodes */
- gmx_sumd(pull->ngrp*4, pull->dbuf_cyl, cr);
+ gmx_sumd(pull->ncoord*4, pull->dbuf_cyl, cr);
}
- for (g = 1; g < 1+pull->ngrp; g++)
+ for (c = 0; c < pull->ncoord; c++)
{
- pgrp = &pull->grp[g];
- pdyna = &pull->dyna[g];
+ pcrd = &pull->coord[c];
+
+ pdyna = &pull->dyna[c];
+ pgrp = &pull->group[pcrd->group[1]];
- wmass = pull->dbuf_cyl[(g-1)*4+0];
- wwmass = pull->dbuf_cyl[(g-1)*4+1];
+ wmass = pull->dbuf_cyl[c*4+0];
+ wwmass = pull->dbuf_cyl[c*4+1];
pdyna->wscale = wmass/wwmass;
pdyna->invtm = 1.0/(pdyna->wscale*wmass);
for (m = 0; m < DIM; m++)
{
- g_x[m] = pgrp->x[m] - pgrp->vec[m]*(pgrp->init[0] + pgrp->rate*t);
- pdyna->x[m] = g_x[m] + pgrp->vec[m]*pull->dbuf_cyl[(g-1)*4+2]/wmass;
+ g_x[m] = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t);
+ pdyna->x[m] = g_x[m] + pcrd->vec[m]*pull->dbuf_cyl[c*4+2]/wmass;
if (xp)
{
- pdyna->xp[m] = g_x[m] + pgrp->vec[m]*pull->dbuf_cyl[(g-1)*4+3]/wmass;
+ pdyna->xp[m] = g_x[m] + pcrd->vec[m]*pull->dbuf_cyl[c*4+3]/wmass;
}
}
if (debug)
{
fprintf(debug, "Pull cylinder group %d:%8.3f%8.3f%8.3f m:%8.3f\n",
- g, pdyna->x[0], pdyna->x[1],
+ c, pdyna->x[0], pdyna->x[1],
pdyna->x[2], 1.0/pdyna->invtm);
}
}
t_pull *pull, t_mdatoms *md, t_pbc *pbc, double t,
rvec x[], rvec *xp)
{
- int g, i, ii, m;
- real mass, w, wm, twopi_box = 0;
- double wmass, wwmass, invwmass;
- dvec com, comp;
- double cm, sm, cmp, smp, ccm, csm, ssm, csw, snw;
- rvec *xx[2], x_pbc = {0, 0, 0}, dx;
- t_pullgrp *pgrp;
+ int g, i, ii, m;
+ real mass, w, wm, twopi_box = 0;
+ double wmass, wwmass, invwmass;
+ dvec com, comp;
+ double cm, sm, cmp, smp, ccm, csm, ssm, csw, snw;
+ rvec *xx[2], x_pbc = {0, 0, 0}, dx;
+ t_pull_group *pgrp;
if (pull->rbuf == NULL)
{
- snew(pull->rbuf, 1+pull->ngrp);
+ snew(pull->rbuf, pull->ngroup);
}
if (pull->dbuf == NULL)
{
- snew(pull->dbuf, 3*(1+pull->ngrp));
+ snew(pull->dbuf, 3*pull->ngroup);
}
if (pull->bRefAt)
twopi_box = 2.0*M_PI/pbc->box[pull->cosdim][pull->cosdim];
}
- for (g = 0; g < 1+pull->ngrp; g++)
+ for (g = 0; g < pull->ngroup; g++)
{
- pgrp = &pull->grp[g];
+ pgrp = &pull->group[g];
clear_dvec(com);
clear_dvec(comp);
wmass = 0;
if (cr && PAR(cr))
{
/* Sum the contributions over the nodes */
- gmx_sumd((1+pull->ngrp)*3*DIM, pull->dbuf[0], cr);
+ gmx_sumd(pull->ngroup*3*DIM, pull->dbuf[0], cr);
}
- for (g = 0; g < 1+pull->ngrp; g++)
+ for (g = 0; g < pull->ngroup; g++)
{
- pgrp = &pull->grp[g];
+ pgrp = &pull->group[g];
if (pgrp->nat > 0 && !(g == 0 && PULL_CYL(pull)))
{
if (pgrp->epgrppbc != epgrppbcCOS)
{
int i;
- for (i = 0; i < pull->ngrp+1; i++)
+ for (i = 0; i < pull->ncoord; i++)
{
- fprintf(fp, "comparing pull group %d\n", i);
- cmp_real(fp, "pullgrp->k", -1, pull->grp[i].k, pull->grp[i].kB, ftol, abstol);
+ fprintf(fp, "comparing pull coord %d\n", i);
+ cmp_real(fp, "pull-coord->k", -1, pull->coord[i].k, pull->coord[i].kB, ftol, abstol);
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff