A constant force is applied between the centers of mass of two groups.
Thus, the potential is linear.
In this case there is no reference distance of pull rate.
+\item{\textbf{Flat bottom pulling}}
+Like umbrella pulling, but the potential and force are zero for
+negative deviations. This is useful for restraining the distance
+between e.g. two molecules to a certain maximum distance.
\end{enumerate}
\subsubsection{Definition of the center of mass}
from the whole bilayer. The whole bilayer can be taken as reference
group in that case, but it might also be of interest to define the
reaction coordinate for the PMF more locally. The {\tt .mdp} option
-{\tt pull-geometry = cylinder} does not
+{\tt pull-coord?-geometry = cylinder} does not
use all the atoms of the reference group, but instead dynamically only those
-within a cylinder with radius {\tt pull-r1} around the pull vector going
+within a cylinder with radius {\tt pull-cylinder-r} around the pull vector going
through the pull group. This only
works for distances defined in one dimension, and the cylinder is
-oriented with its long axis along this one dimension. A second cylinder
-can be defined with {\tt pull-r0}, with a linear switch function that weighs
-the contribution of atoms between {\tt pull-r0} and {\tt pull-r1} with
-distance. This smooths the effects of atoms moving in and out of the
-cylinder (which causes jumps in the pull forces).
+oriented with its long axis along this one dimension. To avoid jumps in
+the pull force, contributions of atoms are weighted as a function of distance
+(in addition to the mass weighting):
+\bea
+w(r < r_\mathrm{cyl}) & = &
+1-2 \left(\frac{r}{r_\mathrm{cyl}}\right)^2 + \left(\frac{r}{r_\mathrm{cyl}}\right)^4 \\
+w(r \geq r_\mathrm{cyl}) & = & 0
+\eea
+Note that the radial dependence on the weight causes a radial force on
+both cylinder group and the other pull group. This is an undesirable,
+but unavoidable effect. To minimize this effect, the cylinder radius should
+be chosen sufficiently large. The effective mass is 0.47 times that of
+a cylinder with uniform weights and equal to the mass of uniform cylinder
+of 0.79 times the radius.
\begin{figure}
\centerline{\includegraphics[width=6cm]{plots/pullref}}
be ignored (and if present in the :ref:`mdp` file, they unfortunately
generate warnings)
- .. mdp-value:: umbrella
-
- Center of mass pulling using an umbrella potential between the
- reference group and one or more groups.
-
- .. mdp-value:: constraint
-
- Center of mass pulling using a constraint between the reference
- group and one or more groups. The setup is identical to the
- option umbrella, except for the fact that a rigid constraint is
- applied instead of a harmonic potential.
-
- .. mdp-value:: constant-force
+ .. mdp-value:: yes
- Center of mass pulling using a linear potential and therefore a
- constant force. For this option there is no reference position
- and therefore the parameters :mdp:`pull-coord1-init` and
- :mdp:`pull-coord1-rate` are not used.
+ Center of mass pulling will be applied on 1 or more groups using
+ 1 or more pull coordinates.
-.. mdp:: pull-geometry
+.. mdp:: pull-cylinder-r
- .. mdp-value:: distance
+ (1.5) \[nm\]
+ the radius of the cylinder for
+ :mdp:`pull-coord1-geometry` = :mdp-value:`cylinder`
- Pull along the vector connecting the two groups. Components can
- be selected with :mdp:`pull-dim`.
-
- .. mdp-value:: direction
-
- Pull in the direction of :mdp:`pull-coord1-vec`.
+.. mdp:: pull-constr-tol
- .. mdp-value:: direction-periodic
+ (1e-6)
+ the relative constraint tolerance for constraint pulling
- As :mdp-value:`pull-geometry=direction`, but allows the distance
- to be larger than half the box size. With this geometry the box
- should not be dynamic (*e.g.* no pressure scaling) in the pull
- dimensions and the pull force is not added to virial.
+.. mdp:: pull-print-com1
- .. mdp-value:: cylinder
+ .. mdp-value:: no
- Designed for pulling with respect to a layer where the reference
- COM is given by a local cylindrical part of the reference
- group. The pulling is in the direction of
- :mdp:`pull-coord1-vec`. From the reference group a cylinder is
- selected around the axis going through the pull group with
- direction :mdp:`pull-coord1-vec` using two radii. The radius
- :mdp:`pull-r1` gives the radius within which all the relative
- weights are one, between :mdp:`pull-r1` and :mdp:`pull-r0` the
- weights are switched to zero. Mass weighting is also used. 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.
-
-.. mdp:: pull-dim
+ do not print the COM of the first group in each pull coordinate
- (Y Y Y)
- the distance components to be used with :mdp:`pull-geometry`
- distance, and also sets which components are printed to the output
- files
+ .. mdp-value:: yes
-.. mdp:: pull-r1
+ print the COM of the first group in each pull coordinate
- (1) \[nm\]
- the inner radius of the cylinder for :mdp:`pull-geometry` cylinder
+.. mdp:: pull-print-com2
-.. mdp:: pull-r0
+ .. mdp-value:: no
- (1) \[nm\]
- the outer radius of the cylinder for :mdp:`pull-geometry` cylinder
+ do not print the COM of the second group in each pull coordinate
-.. mdp:: pull-constr-tol
+ .. mdp-value:: yes
- (1e-6)
- the relative constraint tolerance for constraint pulling
+ print the COM of the second group in each pull coordinate
-.. mdp:: pull-start
+.. mdp:: pull-print-ref-value
.. mdp-value:: no
- do not modify :mdp:`pull-coord1-init`
+ do not print the reference value for each pull coordinate
.. mdp-value:: yes
- add the COM distance of the starting conformation to
- :mdp:`pull-coord1-init`
+ print the reference value for each pull coordinate
-.. mdp:: pull-print-reference
+.. mdp:: pull-print-components
.. mdp-value:: no
- do not print the COM of the first group in each pull coordinate
-
+ only print the distance for each pull coordinate
+
.. mdp-value:: yes
- print the COM of the first group in each pull coordinate
+ print the distance and Cartesian components selected in
+ :mdp:`pull-coord1-dim`
.. mdp:: pull-nstxout
- (10)
- frequency for writing out the COMs of all the pull group
+ (50)
+ frequency for writing out the COMs of all the pull group (0 is
+ never)
.. mdp:: pull-nstfout
- (1)
+ (50)
frequency for writing out the force of all the pulled group
+ (0 is never)
+
.. mdp:: pull-ngroups
and one should think about what to do with the center of mass
motion.
+.. mdp:: pull-coord1-type:
+
+ .. mdp-value:: umbrella
+
+ Center of mass pulling using an umbrella potential between the
+ reference group and one or more groups.
+
+ .. mdp-value:: constraint
+
+ Center of mass pulling using a constraint between the reference
+ group and one or more groups. The setup is identical to the
+ option umbrella, except for the fact that a rigid constraint is
+ applied instead of a harmonic potential.
+
+ .. mdp-value:: constant-force
+
+ Center of mass pulling using a linear potential and therefore a
+ constant force. For this option there is no reference position
+ and therefore the parameters :mdp:`pull-coord1-init` and
+ :mdp:`pull-coord1-rate` are not used.
+
+ .. mdp-value:: flat-bottom
+
+ At distances beyond :mdp:`pull-coord1-init` a harmonic potential
+ is applied, otherwise no potential is applied.
+
+.. mdp:: pull-coord1-geometry
+
+ .. mdp-value:: distance
+
+ Pull along the vector connecting the two groups. Components can
+ be selected with :mdp:`pull-coord1-dim`.
+
+ .. mdp-value:: direction
+
+ Pull in the direction of :mdp:`pull-coord1-vec`.
+
+ .. mdp-value:: direction-periodic
+
+ As :mdp-value:`direction`, but allows the distance to be larger
+ than half the box size. With this geometry the box should not be
+ dynamic (*e.g.* no pressure scaling) in the pull dimensions and
+ the pull force is not added to virial.
+
+ .. mdp-value:: cylinder
+
+ Designed for pulling with respect to a layer where the reference
+ COM is given by a local cylindrical part of the reference group.
+ The pulling is in the direction of :mdp:`pull-coord1-vec`. From
+ the first of the two groups in :mdp:`pull-coord1-groups` a
+ cylinder is selected around the axis going through the COM of
+ the second group with direction :mdp:`pull-coord1-vec` with
+ radius :mdp:`pull-cylinder-r`. Weights of the atoms decrease
+ continously to zero as the radial distance goes from 0 to
+ :mdp:`pull-cylinder-r` (mass weighting is also used). The radial
+ dependence gives rise to radial forces on both pull groups.
+ Note that the radius 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. This geometry is not supported with constraint
+ pulling.
+
+.. mdp:: pull-coord1-dim
+
+ (Y Y Y)
+ Selects the dimensions that this pull coordinate acts on and that
+ are printed to the output files when
+ :mdp:`pull-print-components` = :mdp-value:`yes`. With
+ :mdp:`pull-coord1-geometry` = :mdp-value:`distance`, only Cartesian
+ components set to Y contribute to the distance. Thus setting this
+ to Y Y N results in a distance in the x/y plane. With other
+ geometries all dimensions with non-zero entries in
+ :mdp:`pull-coord1-vec` should be set to Y, the values for other
+ dimensions only affect the output.
+
.. mdp:: pull-coord1-origin
(0.0 0.0 0.0)
(0.0 0.0 0.0)
The pull direction. :ref:`gmx grompp` normalizes the vector.
+.. mdp:: pull-coord1-start
+
+ .. mdp-value:: no
+
+ do not modify :mdp:`pull-coord1-init`
+
+ .. mdp-value:: yes
+
+ add the COM distance of the starting conformation to
+ :mdp:`pull-coord1-init`
+
.. mdp:: pull-coord1-init
(0.0) \[nm\]
set_constraints(constr, top_local, ir, mdatoms, cr);
}
- if (ir->ePull != epullNO)
+ if (ir->bPull)
{
/* Update the local pull groups */
dd_make_local_pull_groups(dd, ir->pull, mdatoms);
tpxv_Use64BitRandomSeed, /**< change ld_seed from int to gmx_int64_t */
tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, /**< potentials for supporting coarse-grained force fields */
tpxv_InteractiveMolecularDynamics, /**< interactive molecular dynamics (IMD) */
- tpxv_RemoveObsoleteParameters1 /**< remove optimize_fft, dihre_fc, nstcheckpoint */
+ tpxv_RemoveObsoleteParameters1, /**< remove optimize_fft, dihre_fc, nstcheckpoint */
+ tpxv_PullCoordTypeGeom /**< add pull type and geometry per group and flat-bottom */
};
/*! \brief Version number of the file format written to run input
*
* When developing a feature branch that needs to change the run input
* file format, change tpx_tag instead. */
-static const int tpx_version = tpxv_RemoveObsoleteParameters1;
+static const int tpx_version = tpxv_PullCoordTypeGeom;
/* This number should only be increased when you edit the TOPOLOGY section
gmx_fio_do_int(fio, pgrp->pbcatom);
}
-static void do_pull_coord(t_fileio *fio, t_pull_coord *pcrd)
+static void do_pull_coord(t_fileio *fio, t_pull_coord *pcrd, int file_version,
+ int ePullOld, int eGeomOld, ivec dimOld)
{
int i;
gmx_fio_do_int(fio, pcrd->group[0]);
gmx_fio_do_int(fio, pcrd->group[1]);
+ if (file_version >= tpxv_PullCoordTypeGeom)
+ {
+ gmx_fio_do_int(fio, pcrd->eType);
+ gmx_fio_do_int(fio, pcrd->eGeom);
+ gmx_fio_do_ivec(fio, pcrd->dim);
+ }
+ else
+ {
+ pcrd->eType = ePullOld;
+ pcrd->eGeom = eGeomOld;
+ copy_ivec(dimOld, pcrd->dim);
+ }
gmx_fio_do_rvec(fio, pcrd->origin);
gmx_fio_do_rvec(fio, pcrd->vec);
+ if (file_version >= tpxv_PullCoordTypeGeom)
+ {
+ gmx_fio_do_gmx_bool(fio, pcrd->bStart);
+ }
+ else
+ {
+ /* This parameter is only printed, but not actually used by mdrun */
+ pcrd->bStart = FALSE;
+ }
gmx_fio_do_real(fio, pcrd->init);
gmx_fio_do_real(fio, pcrd->rate);
gmx_fio_do_real(fio, pcrd->k);
}
}
-static void do_pull(t_fileio *fio, t_pull *pull, gmx_bool bRead, int file_version)
+static void do_pull(t_fileio *fio, t_pull *pull, gmx_bool bRead,
+ int file_version, int ePullOld)
{
- int g;
+ int eGeomOld;
+ ivec dimOld;
+ int g;
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);
+ if (file_version < tpxv_PullCoordTypeGeom)
+ {
+ real dum;
+
+ gmx_fio_do_int(fio, eGeomOld);
+ gmx_fio_do_ivec(fio, dimOld);
+ /* The inner cylinder radius, now removed */
+ gmx_fio_do_real(fio, dum);
+ }
+ gmx_fio_do_real(fio, pull->cylinder_r);
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->bPrintCOM1);
+ /* With file_version < 95 this value is set below */
+ }
+ if (file_version >= tpxv_PullCoordTypeGeom)
+ {
+ gmx_fio_do_int(fio, pull->bPrintCOM2);
+ gmx_fio_do_int(fio, pull->bPrintRefValue);
+ gmx_fio_do_int(fio, pull->bPrintComp);
+ }
+ else
+ {
+ pull->bPrintCOM2 = FALSE;
+ pull->bPrintRefValue = FALSE;
+ pull->bPrintComp = TRUE;
}
gmx_fio_do_int(fio, pull->nstxout);
gmx_fio_do_int(fio, pull->nstfout);
if (file_version < 95)
{
/* epullgPOS for position pulling, before epullgDIRPBC was removed */
- if (pull->eGeom == epullgDIRPBC)
+ if (eGeomOld == epullgDIRPBC)
{
gmx_fatal(FARGS, "pull-geometry=position is no longer supported");
}
- if (pull->eGeom > epullgDIRPBC)
+ if (eGeomOld > epullgDIRPBC)
{
- pull->eGeom -= 1;
+ eGeomOld -= 1;
}
for (g = 0; g < pull->ngroup; g++)
}
}
- pull->bPrintRef = (pull->group[0].nat > 0);
+ pull->bPrintCOM1 = (pull->group[0].nat > 0);
}
else
{
}
for (g = 0; g < pull->ncoord; g++)
{
- do_pull_coord(fio, &pull->coord[g]);
+ do_pull_coord(fio, &pull->coord[g],
+ file_version, ePullOld, eGeomOld, dimOld);
}
}
}
/* pull stuff */
if (file_version >= 48)
{
- gmx_fio_do_int(fio, ir->ePull);
- if (ir->ePull != epullNO)
+ int ePullOld = 0;
+
+ if (file_version >= tpxv_PullCoordTypeGeom)
+ {
+ gmx_fio_do_gmx_bool(fio, ir->bPull);
+ }
+ else
+ {
+ gmx_fio_do_int(fio, ePullOld);
+ ir->bPull = (ePullOld > 0);
+ /* We removed the first ePull=ePullNo for the enum */
+ ePullOld -= 1;
+ }
+ if (ir->bPull)
{
if (bRead)
{
snew(ir->pull, 1);
}
- do_pull(fio, ir->pull, bRead, file_version);
+ do_pull(fio, ir->pull, bRead, file_version, ePullOld);
}
}
else
{
- ir->ePull = epullNO;
+ ir->bPull = FALSE;
}
/* Enforced rotation */
* 3. returns rcoulomb from the tpr */
static float inspect_tpr(int nfile, t_filenm fnm[], real *rcoulomb)
{
- gmx_bool bPull; /* Is pulling requested in .tpr file? */
gmx_bool bTpi; /* Is test particle insertion requested? */
gmx_bool bFree; /* Is a free energy simulation requested? */
gmx_bool bNM; /* Is a normal mode analysis requested? */
/* Check tpr file for options that trigger extra output files */
read_tpx_state(opt2fn("-s", nfile, fnm), &ir, &state, NULL, &mtop);
- bPull = (epullNO != ir.ePull);
bFree = (efepNO != ir.efep );
bNM = (eiNM == ir.eI );
bSwap = (eswapNO != ir.eSwapCoords);
bTpi = EI_TPI(ir.eI);
/* Set these output files on the tuning command-line */
- if (bPull)
+ if (ir.bPull)
{
setopt("-pf", nfile, fnm);
setopt("-px", nfile, fnm);
* \name Using umbrella pull code since gromacs 4.x
*/
/*!\{*/
- int npullcrds; //!< nr of pull coordinates in tpr file
+ int npullcrds_tot; //!< nr of pull coordinates in tpr file
+ int npullcrds; //!< nr of umbrella pull coordinates for reading
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 ?
+ gmx_bool bPrintComp; //!< Components of pull distance written to pullx.xvg ?
real *k; //!< force constants in tpr file
real *init_dist; //!< reference displacements
real *umbInitDist; //!< reference displacement in umbrella direction
/* printf("Reading %s \n",fn); */
read_tpx_state(fn, &ir, &state, NULL, NULL);
- if (ir.ePull != epullUMBRELLA)
+ if (!ir.bPull)
{
- gmx_fatal(FARGS, "This is not a tpr of an umbrella simulation. Found pull type \"%s\" "
- " (ir.ePull = %d)\n", epull_names[ir.ePull], ir.ePull);
+ gmx_fatal(FARGS, "This is not a tpr with COM pulling");
}
/* nr of pull groups */
- ncrd = ir.pull->ncoord;
+ ncrd = 0;
+ for (i = 0; i < ir.pull->ncoord; i++)
+ {
+ if (ir.pull->coord[i].eType == epullUMBRELLA)
+ {
+ int m;
+
+ if (ncrd == 0)
+ {
+ header->pull_geometry = ir.pull->coord[i].eGeom;
+ copy_ivec(ir.pull->coord[i].dim, header->pull_dim);
+ }
+
+ if (ncrd != i)
+ {
+ /* TODO: remove this restriction */
+ gmx_fatal(FARGS, "Currently tpr files where a non-umbrella pull coordinate is followed by an umbrella pull coordinate are not supported");
+ }
+
+ for (m = 0; m < DIM; m++)
+ {
+ if (ir.pull->coord[i].eGeom != header->pull_geometry)
+ {
+ /* TODO: remove the restriction */
+ gmx_fatal(FARGS, "Currently all umbrella pull coordinates should use the same geometry");
+ }
+
+ if (ir.pull->coord[i].dim[m] != header->pull_dim[m])
+ {
+ /* TODO: remove the restriction */
+ gmx_fatal(FARGS, "Currently all umbrella pull coordinates should operate on the same dimensions");
+ }
+ }
+
+ ncrd++;
+ }
+ }
+
if (ncrd < 1)
{
- gmx_fatal(FARGS, "This is not a tpr of umbrella simulation. Found only %d pull coordinates\n", ncrd);
+ gmx_fatal(FARGS, "This is not a tpr of umbrella simulation. Found only %d umbrella pull coordinates\n", ncrd);
}
- 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];
+ header->npullcrds_tot = ir.pull->ncoord;
+ header->npullcrds = ncrd;
+ header->bPrintRef = ir.pull->bPrintCOM1;
+ if (ir.pull->bPrintCOM2)
+ {
+ gmx_fatal(FARGS, "Reading pull output with printing of the COM of group 2 is currently not supported");
+ }
+ if (ir.pull->bPrintRefValue)
+ {
+ gmx_fatal(FARGS, "Reading pull output with printing of the reference value of the coordinates is currently not supported");
+ }
+ header->bPrintComp = ir.pull->bPrintComp;
+ header->pull_ndim = header->pull_dim[0]+header->pull_dim[1]+header->pull_dim[2];
+ /* We should add a struct for each pull coord with all pull coord data
+ instead of allocating many arrays for each property */
snew(header->k, ncrd);
snew(header->init_dist, ncrd);
snew(header->umbInitDist, ncrd);
* bPrintRef == FALSE: for each pull coordinate: no reference columns, but ndim dx columns
*/
- nColPerCrd = opt->bPullx ? header->pull_ndim : 1;
+ nColPerCrd = 1;
+ if (opt->bPullx && header->bPrintComp)
+ {
+ nColPerCrd += header->pull_ndim;
+ }
quantity = opt->bPullx ? "position" : "force";
if (opt->bPullx && header->bPrintRef)
printf("Reading pull %s file with pull geometry %s and %d pull dimensions\n",
bHaveForce ? "force" : "position", epullg_names[header->pull_geometry],
header->pull_ndim);
+ /* Since this tool code has not updated yet to allow difference pull
+ * dimensions per pull coordinate, we can't easily check the exact
+ * number of expected columns, so we only check what we expect for
+ * the pull coordinates selected for the WHAM analysis.
+ */
printf("Expecting these columns in pull file:\n"
"\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);
+ printf("With %d pull groups, expect %s%d columns (including the time column)\n",
+ header->npullcrds,
+ header->npullcrds < header->npullcrds_tot ? "at least " : "",
+ nColExpect);
bFirst = FALSE;
}
- if (ny != nColExpect)
+ if (ny < nColExpect ||
+ (header->npullcrds == header->npullcrds_tot && ny > nColExpect))
{
- gmx_fatal(FARGS, "Found %d pull coodinates in %s,\n but %d data columns in %s (expected %d)\n"
+ gmx_fatal(FARGS, "Using %d pull coodinates from %s,\n but found %d data columns in %s (expected %s%d)\n"
"\nMaybe you confused options -ix and -if ?\n",
- header->npullcrds, fntpr, ny-1, fn, nColExpect-1);
+ header->npullcrds, fntpr, ny-1, fn,
+ header->npullcrds < header->npullcrds_tot ? "at least " : "",
+ nColExpect-1);
}
if (opt->verbose)
else
{
/* pick the right column from:
- * time ref1[ndim] coord1[ndim] ref2[ndim] coord2[ndim] ...
+ * time ref1[ndim] coord1[1(ndim)] ref2[ndim] coord2[1(ndim)] ...
*/
column = 1 + nColRefEachCrd + g*(nColRefEachCrd+nColPerCrd);
- switch (header->pull_geometry)
- {
- case epullgDIST:
- pos = dist_ndim(y + column, header->pull_ndim, i);
- break;
- case epullgCYL:
- pos = y[column][i];
- break;
- default:
- gmx_fatal(FARGS, "Bad error, this error should have been catched before. Ups.\n");
- }
+ pos = y[column][i];
}
/* printf("crd %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */
"[THISMODULE] is an analysis program that implements the Weighted",
"Histogram Analysis Method (WHAM). It is intended to analyze",
"output files generated by umbrella sampling simulations to ",
- "compute a potential of mean force (PMF).",
+ "compute a potential of mean force (PMF).[PAR]",
+ "",
+ "[THISMODULE] is currently not fully up to date. It only supports pull setups",
+ "where the first pull coordinate(s) is/are umbrella pull coordinates",
+ "and, if multiple coordinates need to be analyzed, all used the same",
+ "geometry and dimensions. In most cases this is not an issue.[PAR]",
"",
"At present, three input modes are supported.",
"",
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
};
const char *epull_names[epullNR+1] = {
- "no", "umbrella", "constraint", "constant-force", NULL
+ "umbrella", "constraint", "constant-force", "flat-bottom", NULL
};
const char *epullg_names[epullgNR+1] = {
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
fprintf(fp, "pull-coord %d:\n", c);
PI("group[0]", pcrd->group[0]);
PI("group[1]", pcrd->group[1]);
+ PS("type", EPULLTYPE(pcrd->eType));
+ PS("geometry", EPULLGEOM(pcrd->eGeom));
+ pr_ivec(fp, indent, "dim", pcrd->dim, DIM, TRUE);
pr_rvec(fp, indent, "origin", pcrd->origin, DIM, TRUE);
pr_rvec(fp, indent, "vec", pcrd->vec, DIM, TRUE);
+ PS("start", EBOOL(pcrd->bStart));
PR("init", pcrd->init);
PR("rate", pcrd->rate);
PR("k", pcrd->k);
{
int g;
- PS("pull-geometry", EPULLGEOM(pull->eGeom));
- pr_ivec(fp, indent, "pull-dim", pull->dim, DIM, TRUE);
- PR("pull-r1", pull->cyl_r1);
- PR("pull-r0", pull->cyl_r0);
+ PR("pull-cylinder-r", pull->cylinder_r);
PR("pull-constr-tol", pull->constr_tol);
- PS("pull-print-reference", EBOOL(pull->bPrintRef));
+ PS("pull-print-COM1", EBOOL(pull->bPrintCOM1));
+ PS("pull-print-COM2", EBOOL(pull->bPrintCOM2));
+ PS("pull-print-ref-value", EBOOL(pull->bPrintRefValue));
+ PS("pull-print-components", EBOOL(pull->bPrintComp));
PI("pull-nstxout", pull->nstxout);
PI("pull-nstfout", pull->nstfout);
PI("pull-ngroups", pull->ngroup);
PR("wall-ewald-zfac", ir->wall_ewald_zfac);
/* COM PULLING */
- PS("pull", EPULLTYPE(ir->ePull));
- if (ir->ePull != epullNO)
+ PS("pull", EBOOL(ir->bPull));
+ if (ir->bPull)
{
pr_pull(fp, indent, ir->pull);
}
}
}
- if (ir->ePull != epullNO)
+ if (ir->bPull)
{
- set_pull_init(ir, sys, state.x, state.box, state.lambda[efptMASS], oenv, opts->pull_start);
+ set_pull_init(ir, sys, state.x, state.box, state.lambda[efptMASS], oenv);
}
if (ir->bRot)
/* COM pulling */
CCTYPE("COM PULLING");
- CTYPE("Pull type: no, umbrella, constraint or constant-force");
- EETYPE("pull", ir->ePull, epull_names);
- if (ir->ePull != epullNO)
+ EETYPE("pull", ir->bPull, yesno_names);
+ if (ir->bPull)
{
snew(ir->pull, 1);
- is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
+ is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, wi);
}
/* Enforced rotation */
}
}
- if (ir->ePull == epullCONSTRAINT)
+ if (ir->bPull)
{
/* Correct nrdf for the COM constraints.
* We correct using the TC and VCM group of the first atom
for (i = 0; i < pull->ncoord; i++)
{
+ if (pull->coord[i].eType != epullCONSTRAINT)
+ {
+ continue;
+ }
+
imin = 1;
for (j = 0; j < 2; j++)
}
}
- if (ir->ePull != epullNO)
+ if (ir->bPull)
{
make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
}
- if (ir->ePull != epullNO)
+ if (ir->bPull)
{
- gmx_bool bPullAbsoluteRef;
+ gmx_bool bWarned;
- bPullAbsoluteRef = FALSE;
- for (i = 0; i < ir->pull->ncoord; i++)
+ bWarned = FALSE;
+ for (i = 0; i < ir->pull->ncoord && !bWarned; 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->pull->coord[i].group[0] == 0 ||
+ ir->pull->coord[i].group[1] == 0)
{
- if (ir->pull->dim[m] && !AbsRef[m])
+ absolute_reference(ir, sys, FALSE, AbsRef);
+ for (m = 0; m < DIM; m++)
{
- warning(wi, "You are using an absolute reference for pulling, but the rest of the system does not have an absolute reference. This will lead to artifacts.");
- break;
+ if (ir->pull->coord[i].dim[m] && !AbsRef[m])
+ {
+ warning(wi, "You are using an absolute reference for pulling, but the rest of the system does not have an absolute reference. This will lead to artifacts.");
+ bWarned = TRUE;
+ break;
+ }
}
}
}
- if (ir->pull->eGeom == epullgDIRPBC)
+ for (i = 0; i < 3; i++)
{
- for (i = 0; i < 3; i++)
+ for (m = 0; m <= i; m++)
{
- for (m = 0; m <= i; m++)
+ if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
+ ir->deform[i][m] != 0)
{
- if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
- ir->deform[i][m] != 0)
+ for (c = 0; c < ir->pull->ncoord; c++)
{
- for (c = 0; c < ir->pull->ncoord; c++)
+ if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
+ ir->pull->coord[c].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);
- }
+ gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
}
}
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
gmx_bool bOrire;
gmx_bool bMorse;
char *wall_atomtype[2];
- gmx_bool pull_start;
char *couple_moltype;
int couple_lam0;
int couple_lam1;
/* Routines In readpull.c */
char **read_pullparams(int *ninp_p, t_inpfile **inp,
- t_pull *pull, gmx_bool *bStart,
+ t_pull *pull,
warninp_t wi);
/* Reads the pull parameters, returns a list of the pull group names */
/* Process the pull coordinates after reading the pull groups */
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);
+ const output_env_t oenv);
/* Prints the initial pull group distances in x.
- * If bStart adds the distance to the initial reference location.
+ * If requested, adds the current distance to the initial reference location.
*/
int str_nelem(const char *str, int maxptr, char *ptr[]);
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/utility/smalloc.h"
-static char pulldim[STRLEN];
-
static void string2dvec(const char buf[], dvec nums)
{
double dum;
}
}
-static void init_pull_coord(t_pull_coord *pcrd, int eGeom,
+static void process_pull_dim(char *dim_buf, ivec dim)
+{
+ int ndim, d, nchar, c;
+ char *ptr, pulldim1[STRLEN];
+ t_pull_coord *pcrd;
+
+ ptr = dim_buf;
+ ndim = 0;
+ for (d = 0; d < DIM; d++)
+ {
+ if (sscanf(ptr, "%s%n", pulldim1, &nchar) != 1)
+ {
+ gmx_fatal(FARGS, "Less than 3 pull dimensions given in pull_dim: '%s'",
+ dim_buf);
+ }
+
+ if (gmx_strncasecmp(pulldim1, "N", 1) == 0)
+ {
+ dim[d] = 0;
+ }
+ else if (gmx_strncasecmp(pulldim1, "Y", 1) == 0)
+ {
+ dim[d] = 1;
+ ndim++;
+ }
+ else
+ {
+ gmx_fatal(FARGS, "Please use Y(ES) or N(O) for pull_dim only (not %s)",
+ pulldim1);
+ }
+ ptr += nchar;
+ }
+ if (ndim == 0)
+ {
+ gmx_fatal(FARGS, "All entries in pull dim are N");
+ }
+}
+
+static void init_pull_coord(t_pull_coord *pcrd,
+ char *dim_buf,
const char *origin_buf, const char *vec_buf,
warninp_t wi)
{
dvec origin, vec;
char buf[STRLEN];
+ if (pcrd->eType == epullCONSTRAINT && pcrd->eGeom == epullgCYL)
+ {
+ gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
+ epull_names[pcrd->eType],
+ epullg_names[pcrd->eGeom],
+ epull_names[epullUMBRELLA]);
+ }
+
+ process_pull_dim(dim_buf, pcrd->dim);
+
string2dvec(origin_buf, origin);
if (pcrd->group[0] != 0 && dnorm(origin) > 0)
{
gmx_fatal(FARGS, "The pull origin can only be set with an absolute reference");
}
- if (eGeom == epullgDIST)
+ if (pcrd->eGeom == epullgDIST)
{
if (pcrd->init < 0)
{
sprintf(buf, "The initial pull distance is negative with geometry %s, while a distance can not be negative. Use geometry %s instead.",
- EPULLGEOM(eGeom), EPULLGEOM(epullgDIR));
+ EPULLGEOM(pcrd->eGeom), EPULLGEOM(epullgDIR));
warning_error(wi, buf);
}
/* TODO: With a positive init but a negative rate things could still
else
{
string2dvec(vec_buf, vec);
- if (eGeom == epullgDIR || eGeom == epullgCYL)
+ if (dnorm2(vec) == 0)
+ {
+ gmx_fatal(FARGS, "With pull geometry %s the pull vector can not be 0,0,0",
+ epullg_names[pcrd->eGeom]);
+ }
+ if (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgCYL)
{
/* Normalize the direction vector */
dsvmul(1/dnorm(vec), vec, vec);
}
char **read_pullparams(int *ninp_p, t_inpfile **inp_p,
- t_pull *pull, gmx_bool *bStart,
+ t_pull *pull,
warninp_t wi)
{
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];
+ char dummy[STRLEN], buf[STRLEN], groups[STRLEN], dim_buf[STRLEN];
+ char 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];
inp = *inp_p;
/* read pull parameters */
- 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");
CTYPE("Cylinder radius for dynamic reaction force groups (nm)");
- 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-cylinder-r", pull->cylinder_r, 1.5);
RTYPE("pull-constr-tol", pull->constr_tol, 1E-6);
- 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);
+ EETYPE("pull-print-com1", pull->bPrintCOM1, yesno_names);
+ EETYPE("pull-print-com2", pull->bPrintCOM2, yesno_names);
+ EETYPE("pull-print-ref-value", pull->bPrintRefValue, yesno_names);
+ EETYPE("pull-print-components", pull->bPrintComp, yesno_names);
+ ITYPE("pull-nstxout", pull->nstxout, 50);
+ ITYPE("pull-nstfout", pull->nstfout, 50);
CTYPE("Number of pull groups");
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");
- }
-
if (pull->ngroup < 1)
{
gmx_fatal(FARGS, "pull-ngroups should be >= 1");
buf, 2);
nerror++;
}
+ sprintf(buf, "pull-coord%d-type", i);
+ EETYPE(buf, pcrd->eType, epull_names);
+ sprintf(buf, "pull-coord%d-geometry", i);
+ EETYPE(buf, pcrd->eGeom, epullg_names);
+ sprintf(buf, "pull-coord%d-dim", i);
+ STYPE(buf, dim_buf, "Y Y Y");
sprintf(buf, "pull-coord%d-origin", i);
- STYPE(buf, origin_buf, "0.0 0.0 0.0");
+ 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");
+ STYPE(buf, vec_buf, "0.0 0.0 0.0");
+ sprintf(buf, "pull-coord%d-start", i);
+ EETYPE(buf, pcrd->bStart, yesno_names);
sprintf(buf, "pull-coord%d-init", i);
- RTYPE(buf, pcrd->init, 0.0);
+ RTYPE(buf, pcrd->init, 0.0);
sprintf(buf, "pull-coord%d-rate", i);
- RTYPE(buf, pcrd->rate, 0.0);
+ RTYPE(buf, pcrd->rate, 0.0);
sprintf(buf, "pull-coord%d-k", i);
- RTYPE(buf, pcrd->k, 0.0);
+ RTYPE(buf, pcrd->k, 0.0);
sprintf(buf, "pull-coord%d-kB", i);
- RTYPE(buf, pcrd->kB, pcrd->k);
+ RTYPE(buf, pcrd->kB, pcrd->k);
/* Initialize the pull coordinate */
- init_pull_coord(pcrd, pull->eGeom, origin_buf, vec_buf, wi);
+ init_pull_coord(pcrd, dim_buf, origin_buf, vec_buf, wi);
}
*ninp_p = ninp;
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)",
void make_pull_coords(t_pull *pull)
{
- int ndim, d, nchar, c;
- char *ptr, pulldim1[STRLEN];
+ int c, d;
t_pull_coord *pcrd;
- ptr = pulldim;
- ndim = 0;
- for (d = 0; d < DIM; d++)
- {
- if (sscanf(ptr, "%s%n", pulldim1, &nchar) != 1)
- {
- gmx_fatal(FARGS, "Less than 3 pull dimensions given in pull_dim: '%s'",
- pulldim);
- }
-
- if (gmx_strncasecmp(pulldim1, "N", 1) == 0)
- {
- pull->dim[d] = 0;
- }
- else if (gmx_strncasecmp(pulldim1, "Y", 1) == 0)
- {
- pull->dim[d] = 1;
- ndim++;
- }
- else
- {
- gmx_fatal(FARGS, "Please use Y(ES) or N(O) for pull_dim only (not %s)",
- pulldim1);
- }
- ptr += nchar;
- }
- if (ndim == 0)
- {
- gmx_fatal(FARGS, "All entries in pull_dim are N");
- }
-
for (c = 0; c < pull->ncoord; c++)
{
pcrd = &pull->coord[c];
gmx_fatal(FARGS, "Identical pull group indices in pull-coord%d-groups", c+1);
}
- if (pull->eGeom == epullgCYL && pcrd->group[0] != 1)
+ if (pcrd->eGeom == epullgCYL)
{
- gmx_fatal(FARGS, "With pull geometry %s, the first pull group should always be 1", EPULLGEOM(pull->eGeom));
+ if (pull->group[pcrd->group[0]].nweight > 0)
+ {
+ gmx_fatal(FARGS, "Weights are not supported for the reference group with cylinder pulling");
+ }
}
- if (pull->eGeom != epullgDIST)
+ if (pcrd->eGeom != epullgDIST)
{
for (d = 0; d < DIM; d++)
{
- if (pcrd->vec[d] != 0 && pull->dim[d] == 0)
+ if (pcrd->vec[d] != 0 && pcrd->dim[d] == 0)
{
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) &&
+ if ((pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgCYL) &&
norm2(pcrd->vec) == 0)
{
gmx_fatal(FARGS, "pull-group%d-vec can not be zero with geometry %s",
- c+1, EPULLGEOM(pull->eGeom));
+ c+1, EPULLGEOM(pcrd->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)
+ const output_env_t oenv)
{
t_mdatoms *md;
t_pull *pull;
fprintf(stderr, "%8d %8d %8d ",
pcrd->group[1], pgrp1->nat, pgrp1->pbcatom+1);
- if (bStart)
+ if (pcrd->bStart)
{
init = pcrd->init;
pcrd->init = 0;
value = pcrd->init + t_start*pcrd->rate + dev;
fprintf(stderr, " %10.3f nm", value);
- if (bStart)
+ if (pcrd->bStart)
{
pcrd->init = value + init;
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
/* Pull stuff */
enum {
- epullNO, epullUMBRELLA, epullCONSTRAINT, epullCONST_F, epullNR
+ epullUMBRELLA, epullCONSTRAINT, epullCONST_F, epullFLATBOTTOM, epullNR
};
enum {
epullgDIST, epullgDIR, epullgCYL, epullgDIRPBC, epullgNR
};
-#define PULL_CYL(pull) ((pull)->eGeom == epullgCYL)
-
/* Enforced rotation groups */
enum {
erotgISO, erotgISOPF,
atom_id pbcatom; /* The reference atom for pbc (global number) */
/* Variables not present in mdp, but used at run time */
+ gmx_bool bCalcCOM; /* Calculate COM? Not if only used as cylinder group */
+ real mwscale; /* mass*weight scaling factor 1/sum w m */
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 *mdw; /* mass*gradient(weight) for atoms */
+ double *dv; /* distance to the other group along vec */
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 */
+ int eType; /* The pull type: umbrella, constraint, ... */
+ int eGeom; /* The pull geometry */
+ ivec dim; /* Used to select components for constraint */
rvec origin; /* The origin for the absolute reference */
rvec vec; /* The pull vector, direction or position */
+ gmx_bool bStart; /* Set init based on the initial structure */
real init; /* Initial reference displacement */
real rate; /* Rate of motion (nm/ps) */
real k; /* force constant */
/* Variables not present in mdp, but used at run time */
dvec dr; /* The distance from the reference group */
+ dvec ffrad; /* conversion factor from vec to radial force */
+ double cyl_dev; /* The deviation from the reference position */
double f_scal; /* Scalar force for directional pulling */
dvec f; /* force due to the pulling/constraining */
} t_pull_coord;
} t_expanded;
typedef struct {
- 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 */
- int npbcdim; /* do pbc in dims 0 <= dim < npbcdim */
- 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_pull_group *group; /* groups to pull/restrain/etc/ */
- t_pull_coord *coord; /* the pull coordinates */
+ int ngroup; /* number of pull groups */
+ int ncoord; /* number of pull coordinates */
+ real cylinder_r; /* radius of cylinder for dynamic COM */
+ real constr_tol; /* absolute tolerance for constraints in (nm) */
+ gmx_bool bPrintCOM1; /* Print coordinates of COM 1 for each coord */
+ gmx_bool bPrintCOM2; /* Print coordinates of COM 2 for each coord */
+ gmx_bool bPrintRefValue; /* Print the reference value for each coord */
+ gmx_bool bPrintComp; /* Print cartesian components for each coord with geometry=distance */
+ int nstxout; /* Output frequency for pull x */
+ int nstfout; /* Output frequency for pull f */
+ int ePBC; /* the boundary conditions */
+ int npbcdim; /* do pbc in dims 0 <= dim < npbcdim */
+ gmx_bool bRefAt; /* do we need reference atoms for a group COM ? */
+ int cosdim; /* dimension for cosine weighting, -1 if none */
+ 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 */
+ gmx_bool bPotential; /* Are there coordinates with potential? */
+ gmx_bool bConstraint; /* Are there constrained coordinates? */
+ gmx_bool bCylinder; /* Is group 0 a cylinder group? */
t_pull_group *dyna; /* dynamic groups for use with local constraints */
gmx_bool bSetPBCatoms; /* Do we need to set x_pbc for the groups? */
int wall_atomtype[2]; /* The atom type for walls */
real wall_density[2]; /* Number density for walls */
real wall_ewald_zfac; /* Scaling factor for the box for Ewald */
- int ePull; /* Type of pulling: no, umbrella or constraint */
+ gmx_bool bPull; /* Do we do COM pulling? */
t_pull *pull; /* The data for center of mass pulling */
gmx_bool bRot; /* Calculate enforced rotation potential(s)? */
t_rot *rot; /* The data for enforced rotation potentials */
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
{
bc_simtempvals(cr, inputrec->simtempvals, inputrec->fepvals->n_lambda);
}
- if (inputrec->ePull != epullNO)
+ if (inputrec->bPull)
{
snew_bc(cr, inputrec->pull, 1);
bc_pull(cr, inputrec->pull);
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
if (econq == econqCoord)
{
- if (ir->ePull == epullCONSTRAINT)
+ if (ir->bPull && ir->pull->bConstraint)
{
if (EI_DYNAMICS(ir->eI))
{
gmx_mtop_ftype_count(mtop, F_CONSTRNC);
nset = gmx_mtop_ftype_count(mtop, F_SETTLE);
- if (ncon+nset == 0 && ir->ePull != epullCONSTRAINT && ed == NULL)
+ if (ncon+nset == 0 &&
+ !(ir->bPull && ir->pull->bConstraint) &&
+ ed == NULL)
{
return NULL;
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
}
else if (i == F_COM_PULL)
{
- md->bEner[i] = (ir->ePull == epullUMBRELLA || ir->ePull == epullCONST_F || ir->bRot);
+ md->bEner[i] = (ir->bPull && ir->pull->bPotential);
}
else if (i == F_ECONSERVED)
{
set_pbc(&pbc, ir->ePBC, box);
dvdl = 0;
enerd->term[F_COM_PULL] +=
- pull_potential(ir->ePull, ir->pull, mdatoms, &pbc,
+ pull_potential(ir->pull, mdatoms, &pbc,
cr, t, lambda[efptRESTRAINT], x, f, vir_force, &dvdl);
enerd->dvdl_lin[efptRESTRAINT] += dvdl;
wallcycle_stop(wcycle, ewcPULLPOT);
clear_rvec(fr->vir_diag_posres);
}
- if (inputrec->ePull == epullCONSTRAINT)
+ if (inputrec->bPull && inputrec->pull->bConstraint)
{
clear_pull_forces(inputrec->pull);
}
}
}
- if (inputrec->ePull == epullUMBRELLA || inputrec->ePull == epullCONST_F)
+ if (inputrec->bPull && inputrec->pull->bPotential)
{
/* Since the COM pulling is always done mass-weighted, no forces are
* applied to vsites and this call can be done after vsite spreading.
clear_rvec(fr->vir_diag_posres);
}
- if (inputrec->ePull == epullCONSTRAINT)
+ if (inputrec->bPull && inputrec->pull->bConstraint)
{
clear_pull_forces(inputrec->pull);
}
}
}
- if (inputrec->ePull == epullUMBRELLA || inputrec->ePull == epullCONST_F)
+ if (inputrec->bPull && inputrec->pull->bPotential)
{
pull_potential_wrapper(cr, inputrec, box, x,
f, vir_force, mdatoms, enerd, lambda, t,
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
}
}
-static void pull_print_coord_dr(FILE *out, ivec dim, const t_pull_coord *pcrd)
+static void pull_print_coord_dr(FILE *out, const t_pull_coord *pcrd,
+ gmx_bool bPrintRefValue, double t,
+ gmx_bool bPrintComponents)
{
- int m;
+ double distance, distance2;
+ int m;
- for (m = 0; m < DIM; m++)
+ if (pcrd->eGeom == epullgDIST)
{
- if (dim[m])
+ /* Geometry=distance: print the length of the distance vector */
+ distance2 = 0;
+ for (m = 0; m < DIM; m++)
{
- fprintf(out, "\t%g", pcrd->dr[m]);
+ if (pcrd->dim[m])
+ {
+ distance2 += pcrd->dr[m]*pcrd->dr[m];
+ }
+ }
+ distance = sqrt(distance2);
+ }
+ else
+ {
+ /* Geometry=direction-like: print distance along the pull vector */
+ distance = 0;
+ for (m = 0; m < DIM; m++)
+ {
+ if (pcrd->dim[m])
+ {
+ distance += pcrd->dr[m]*pcrd->vec[m];
+ }
+ }
+ }
+ fprintf(out, "\t%g", distance);
+
+ if (bPrintRefValue)
+ {
+ fprintf(out, "\t%g", pcrd->init + pcrd->rate*t);
+ }
+
+ if (bPrintComponents)
+ {
+ for (m = 0; m < DIM; m++)
+ {
+ if (pcrd->dim[m])
+ {
+ fprintf(out, "\t%g", pcrd->dr[m]);
+ }
}
}
}
static void pull_print_x(FILE *out, t_pull *pull, double t)
{
- int c;
- const t_pull_coord *pcrd;
+ int c;
+ t_pull_coord *pcrd;
fprintf(out, "%.4f", t);
{
pcrd = &pull->coord[c];
- if (pull->bPrintRef)
+ if (pull->bPrintCOM1)
{
- if (PULL_CYL(pull))
+ if (pcrd->eGeom == epullgCYL)
{
- pull_print_group_x(out, pull->dim, &pull->dyna[c]);
+ pull_print_group_x(out, pcrd->dim, &pull->dyna[c]);
}
else
{
- pull_print_group_x(out, pull->dim, &pull->group[pcrd->group[0]]);
+ pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[0]]);
}
}
- pull_print_coord_dr(out, pull->dim, pcrd);
+ if (pull->bPrintCOM2)
+ {
+ pull_print_group_x(out, pcrd->dim, &pull->group[pcrd->group[1]]);
+ }
+ pull_print_coord_dr(out, pcrd,
+ pull->bPrintRefValue, t, pull->bPrintComp);
}
fprintf(out, "\n");
}
static void pull_print_f(FILE *out, t_pull *pull, double t)
{
- int c, d;
+ int c;
fprintf(out, "%.4f", t);
exvggtXNY, oenv);
}
- snew(setname, 2*pull->ncoord*DIM);
+ /* With default mdp options only the actual distance is printed,
+ * but optionally 2 COMs, the reference distance and distance
+ * components can also be printed.
+ */
+ snew(setname, pull->ncoord*(DIM + DIM + 1 + 1 + DIM));
nsets = 0;
for (c = 0; c < pull->ncoord; c++)
{
if (bCoord)
{
- if (pull->bPrintRef)
+ /* The order of this legend should match the order of printing
+ * the data in print_pull_x above.
+ */
+
+ if (pull->bPrintCOM1)
{
+ /* Legend for reference group position */
for (m = 0; m < DIM; m++)
{
- if (pull->dim[m])
+ if (pull->coord[c].dim[m])
{
- sprintf(buf, "%d %s%c", c+1, "c", 'X'+m);
+ sprintf(buf, "%d %s%d%c", c+1, "c", 1, 'X'+m);
setname[nsets] = gmx_strdup(buf);
nsets++;
}
}
}
- for (m = 0; m < DIM; m++)
+ if (pull->bPrintCOM2)
+ {
+ /* Legend for reference group position */
+ for (m = 0; m < DIM; m++)
+ {
+ if (pull->coord[c].dim[m])
+ {
+ sprintf(buf, "%d %s%d%c", c+1, "c", 2, 'X'+m);
+ setname[nsets] = gmx_strdup(buf);
+ nsets++;
+ }
+ }
+ }
+ /* The pull coord distance */
+ sprintf(buf, "%d", c+1);
+ setname[nsets] = gmx_strdup(buf);
+ nsets++;
+ if (pull->bPrintRefValue)
{
- if (pull->dim[m])
+ sprintf(buf, "%c%d", 'r', c+1);
+ setname[nsets] = gmx_strdup(buf);
+ nsets++;
+ }
+ if (pull->bPrintComp)
+ {
+ /* The pull coord distance components */
+ for (m = 0; m < DIM; m++)
{
- sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
- setname[nsets] = gmx_strdup(buf);
- nsets++;
+ if (pull->coord[c].dim[m])
+ {
+ sprintf(buf, "%d %s%c", c+1, "d", 'X'+m);
+ setname[nsets] = gmx_strdup(buf);
+ nsets++;
+ }
}
}
}
else
{
+ /* For the pull force we always only use one scalar */
sprintf(buf, "%d", c+1);
setname[nsets] = gmx_strdup(buf);
nsets++;
int i, ii, m;
double wmass, inv_wm;
- inv_wm = pgrp->wscale*pgrp->invtm;
+ inv_wm = pgrp->mwscale;
for (i = 0; i < pgrp->nat_loc; i++)
{
}
}
+/* Apply forces in a mass weighted fashion to a cylinder group */
+static void apply_forces_cyl_grp(const t_pull_group *pgrp, double dv_corr,
+ const t_mdatoms *md,
+ const dvec f_pull, double f_scal,
+ int sign, rvec *f)
+{
+ int i, ii, m;
+ double mass, weight, inv_wm, dv_com;
+
+ inv_wm = pgrp->mwscale;
+
+ for (i = 0; i < pgrp->nat_loc; i++)
+ {
+ ii = pgrp->ind_loc[i];
+ mass = md->massT[ii];
+ weight = pgrp->weight_loc[i];
+ /* The stored axial distance from the cylinder center (dv) needs
+ * to be corrected for an offset (dv_corr), which was unknown when
+ * we calculated dv.
+ */
+ dv_com = pgrp->dv[i] + dv_corr;
+
+ /* Here we not only add the pull force working along vec (f_pull),
+ * but also a radial component, due to the dependence of the weights
+ * on the radial distance.
+ */
+ for (m = 0; m < DIM; m++)
+ {
+ f[ii][m] += sign*inv_wm*(mass*weight*f_pull[m] +
+ pgrp->mdw[i][m]*dv_com*f_scal);
+ }
+ }
+}
+
/* Apply forces in a mass weighted fashion */
static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f)
{
{
pcrd = &pull->coord[c];
- if (PULL_CYL(pull))
+ if (pcrd->eGeom == epullgCYL)
{
- apply_forces_grp(&pull->dyna[c], md, pcrd->f, -1, f);
+ dvec f_tot;
+ int m;
+
+ apply_forces_cyl_grp(&pull->dyna[c], pcrd->cyl_dev, md,
+ pcrd->f, pcrd->f_scal, -1, f);
+
+ /* Sum the force along the vector and the radial force */
+ for (m = 0; m < DIM; m++)
+ {
+ f_tot[m] = pcrd->f[m] + pcrd->f_scal*pcrd->ffrad[m];
+ }
+ apply_forces_grp(&pull->group[pcrd->group[1]], md, f_tot, 1, f);
}
else
{
{
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);
}
- apply_forces_grp(&pull->group[pcrd->group[1]], md, pcrd->f, 1, f);
}
}
-static double max_pull_distance2(const t_pull *pull, const t_pbc *pbc)
+static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc)
{
double max_d2;
int m;
max_d2 = GMX_DOUBLE_MAX;
- if (pull->eGeom != epullgDIRPBC)
+ for (m = 0; m < pbc->ndim_ePBC; m++)
{
- for (m = 0; m < pbc->ndim_ePBC; m++)
+ if (pcrd->dim[m] != 0)
{
- if (pull->dim[m] != 0)
- {
- max_d2 = min(max_d2, norm2(pbc->box[m]));
- }
+ max_d2 = min(max_d2, norm2(pbc->box[m]));
}
}
copy_dvec(xref, xrefr);
- if (pull->eGeom == epullgDIRPBC)
+ if (pcrd->eGeom == epullgDIRPBC)
{
for (m = 0; m < DIM; m++)
{
dr2 = 0;
for (m = 0; m < DIM; m++)
{
- dr[m] *= pull->dim[m];
+ dr[m] *= pcrd->dim[m];
dr2 += dr[m]*dr[m];
}
if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2)
pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2));
}
- if (pull->eGeom == epullgDIRPBC)
+ if (pcrd->eGeom == epullgDIRPBC)
{
dvec_inc(dr, dref);
}
double md2;
const t_pull_coord *pcrd;
- if (pull->eGeom == epullgDIRPBC)
+ pcrd = &pull->coord[coord_ind];
+
+ if (pcrd->eGeom == epullgDIRPBC)
{
md2 = -1;
}
else
{
- md2 = max_pull_distance2(pull, pbc);
+ md2 = max_pull_distance2(pcrd, pbc);
}
- 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,
+ pcrd->eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x,
md2,
dr);
}
ref = pcrd->init + pcrd->rate*t;
- switch (pull->eGeom)
+ switch (pcrd->eGeom)
{
case epullgDIST:
/* Pull along the vector between the com's */
double *dr_tot; /* the total update of the coords */
double ref;
dvec vec;
- double d0, inpr;
- double lambda, rm, mass, invdt = 0;
+ double inpr;
+ double lambda, rm, invdt = 0;
gmx_bool bConverged_all, bConverged = FALSE;
int niter = 0, g, c, ii, j, m, max_iter = 100;
double a;
dvec f; /* the pull force */
dvec tmp, tmp3;
- t_pull_group *pdyna, *pgrp0, *pgrp1;
+ t_pull_group *pgrp0, *pgrp1;
t_pull_coord *pcrd;
snew(r_ij, pull->ncoord);
{
copy_dvec(pull->group[g].xp, rnew[g]);
}
- if (PULL_CYL(pull))
- {
- /* 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 (c = 0; c < pull->ncoord; c++)
{
+ pcrd = &pull->coord[c];
+
+ if (pcrd->eType != epullCONSTRAINT)
+ {
+ continue;
+ }
+
get_pull_coord_dr(pull, c, pbc, t, r_ij[c]);
/* Store the difference vector at time t for printing */
- copy_dvec(r_ij[c], pull->coord[c].dr);
+ copy_dvec(r_ij[c], pcrd->dr);
if (debug)
{
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)
+ if (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgDIRPBC)
{
/* Select the component along vec */
a = 0;
for (m = 0; m < DIM; m++)
{
- a += pull->coord[c].vec[m]*r_ij[c][m];
+ a += pcrd->vec[m]*r_ij[c][m];
}
for (m = 0; m < DIM; m++)
{
- r_ij[c][m] = a*pull->coord[c].vec[m];
+ r_ij[c][m] = a*pcrd->vec[m];
}
}
dvec dr0, dr1;
pcrd = &pull->coord[c];
+
+ if (pcrd->eType != epullCONSTRAINT)
+ {
+ continue;
+ }
+
pgrp0 = &pull->group[pcrd->group[0]];
pgrp1 = &pull->group[pcrd->group[1]];
rm = 1.0/(pgrp0->invtm + pgrp1->invtm);
- switch (pull->eGeom)
+ switch (pcrd->eGeom)
{
case epullgDIST:
if (ref <= 0)
for (c = 0; c < pull->ncoord; c++)
{
pcrd = &pull->coord[c];
+
+ if (pcrd->eType != epullCONSTRAINT)
+ {
+ continue;
+ }
+
ref = pcrd->init + pcrd->rate*t;
low_get_pull_coord_dr(pull, pcrd, pbc, t,
rnew[pcrd->group[0]],
-1, unc_ij);
- switch (pull->eGeom)
+ switch (pcrd->eGeom)
{
case epullgDIST:
bConverged = fabs(dnorm(unc_ij) - ref) < pull->constr_tol;
const t_pull_group *pgrp;
dvec dr;
- if (PULL_CYL(pull) && g == pull->coord[0].group[0])
- {
- pgrp = &pull->dyna[0];
- }
- else
- {
- pgrp = &pull->group[g];
- }
+ pgrp = &pull->group[g];
/* 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++)
+
+ if (dnorm2(dr) == 0)
{
- dr[m] *= pull->dim[m];
+ /* No displacement, no update necessary */
+ continue;
}
/* update the atom positions */
/* calculate the constraint forces, used for output and virial only */
for (c = 0; c < pull->ncoord; c++)
{
- pcrd = &pull->coord[c];
+ pcrd = &pull->coord[c];
+
+ if (pcrd->eType != epullCONSTRAINT)
+ {
+ continue;
+ }
+
pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->group[0]].invtm + pull->group[pcrd->group[1]].invtm)*dt*dt);
- if (vir && bMaster)
+ if (vir != NULL && pcrd->eGeom != epullgDIRPBC && bMaster)
{
double f_invr;
}
/* Pulling with a harmonic umbrella potential or constant force */
-static void do_pull_pot(int ePull,
- t_pull *pull, t_pbc *pbc, double t, real lambda,
+static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda,
real *V, tensor vir, real *dVdl)
{
int c, j, m;
- double dev, ndr, invdr;
+ double dev, ndr, invdr = 0;
real k, dkdl;
t_pull_coord *pcrd;
{
pcrd = &pull->coord[c];
+ if (pcrd->eType == epullCONSTRAINT)
+ {
+ continue;
+ }
+
get_pull_coord_distance(pull, c, pbc, t, pcrd->dr, &dev);
k = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB;
dkdl = pcrd->kB - pcrd->k;
- switch (pull->eGeom)
+ if (pcrd->eGeom == epullgDIST)
{
- case epullgDIST:
- ndr = dnorm(pcrd->dr);
- 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;
- *V += 0.5* k*dsqr(dev);
- *dVdl += 0.5*dkdl*dsqr(dev);
- }
- else
- {
- pcrd->f_scal = -k;
- *V += k*ndr;
- *dVdl += dkdl*ndr;
- }
- for (m = 0; m < DIM; m++)
+ ndr = dnorm(pcrd->dr);
+ 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;
+ }
+ }
+ else
+ {
+ ndr = 0;
+ for (m = 0; m < DIM; m++)
+ {
+ ndr += pcrd->vec[m]*pcrd->dr[m];
+ }
+ }
+
+ switch (pcrd->eType)
+ {
+ case epullUMBRELLA:
+ case epullFLATBOTTOM:
+ /* The only difference between an umbrella and a flat-bottom
+ * potential is that a flat-bottom is zero below zero.
+ */
+ if (pcrd->eType == epullFLATBOTTOM && dev < 0)
{
- pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
+ dev = 0;
}
+
+ pcrd->f_scal = -k*dev;
+ *V += 0.5* k*dsqr(dev);
+ *dVdl += 0.5*dkdl*dsqr(dev);
break;
- case epullgDIR:
- case epullgDIRPBC:
- case epullgCYL:
- if (ePull == epullUMBRELLA)
- {
- 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 += pcrd->vec[m]*pcrd->dr[m];
- }
- pcrd->f_scal = -k;
- *V += k*ndr;
- *dVdl += dkdl*ndr;
- }
- for (m = 0; m < DIM; m++)
- {
- pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
- }
+ case epullCONST_F:
+ pcrd->f_scal = -k;
+ *V += k*ndr;
+ *dVdl += dkdl*ndr;
break;
+ default:
+ gmx_incons("Unsupported pull type in do_pull_pot");
+ }
+
+ if (pcrd->eGeom == epullgDIST)
+ {
+ for (m = 0; m < DIM; m++)
+ {
+ pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr;
+ }
+ }
+ else
+ {
+ for (m = 0; m < DIM; m++)
+ {
+ pcrd->f[m] = pcrd->f_scal*pcrd->vec[m];
+ }
}
- if (vir)
+ if (vir != NULL && pcrd->eGeom != epullgDIRPBC)
{
/* Add the pull contribution to the virial */
for (j = 0; j < DIM; j++)
}
}
-real pull_potential(int ePull, t_pull *pull, t_mdatoms *md, t_pbc *pbc,
+real pull_potential(t_pull *pull, t_mdatoms *md, t_pbc *pbc,
t_commrec *cr, double t, real lambda,
rvec *x, rvec *f, tensor vir, real *dvdlambda)
{
pull_calc_coms(cr, pull, md, pbc, t, x, NULL);
- do_pull_pot(ePull, pull, pbc, t, lambda,
- &V, pull->bVirial && MASTER(cr) ? vir : NULL, &dVdl);
+ do_pull_pot(pull, pbc, t, lambda,
+ &V, MASTER(cr) ? vir : NULL, &dVdl);
/* Distribute forces over pulled groups */
apply_forces(pull, md, f);
{
pull_calc_coms(cr, pull, md, pbc, t, x, xp);
- do_constraint(pull, pbc, xp, v, pull->bVirial && MASTER(cr), vir, dt, t);
+ do_constraint(pull, pbc, xp, v, MASTER(cr), vir, dt, t);
}
static void make_local_pull_group(gmx_ga2la_t ga2la,
}
static void init_pull_group_index(FILE *fplog, t_commrec *cr,
- int g, t_pull_group *pg, ivec pulldims,
+ int g, t_pull_group *pg,
+ gmx_bool bConstraint, ivec pulldim_con,
gmx_mtop_t *mtop, t_inputrec *ir, real lambda)
{
int i, ii, d, nfrozen, ndim;
{
ii = pg->ind[i];
gmx_mtop_atomnr_to_atom(alook, ii, &atom);
- if (ir->opts.nFreeze)
+ if (bConstraint && ir->opts.nFreeze)
{
for (d = 0; d < DIM; d++)
{
- if (pulldims[d] && ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
+ if (pulldim_con[d] == 1 &&
+ ir->opts.nFreeze[ggrpnr(groups, egcFREEZE, ii)][d])
{
nfrozen++;
}
if (nfrozen == 0)
{
- /* A value > 0 signals not frozen, it is updated later */
- pg->invtm = 1.0;
+ /* A value != 0 signals not frozen, it is updated later */
+ pg->invtm = -1.0;
}
else
{
ndim = 0;
for (d = 0; d < DIM; d++)
{
- ndim += pulldims[d]*pg->nat;
+ ndim += pulldim_con[d]*pg->nat;
}
if (fplog && nfrozen > 0 && nfrozen < ndim)
{
fprintf(fplog,
"\nWARNING: In pull group %d some, but not all of the degrees of freedom\n"
" that are subject to pulling are frozen.\n"
- " For pulling the whole group will be frozen.\n\n",
+ " For constraint pulling the whole group will be frozen.\n\n",
g);
}
pg->invtm = 0.0;
{
t_pull *pull;
t_pull_group *pgrp;
- int c, g, start = 0, end = 0, m;
+ int c, g, m;
pull = ir->pull;
+ pull->bPotential = FALSE;
+ pull->bConstraint = FALSE;
+ pull->bCylinder = FALSE;
+ for (c = 0; c < pull->ncoord; c++)
+ {
+ t_pull_coord *pcrd;
+
+ pcrd = &pull->coord[c];
+
+ if (pcrd->eType == epullCONSTRAINT)
+ {
+ pull->bConstraint = TRUE;
+ }
+ else
+ {
+ pull->bPotential = TRUE;
+ }
+
+ if (pcrd->eGeom == epullgCYL)
+ {
+ pull->bCylinder = TRUE;
+
+ if (pcrd->eType == epullCONSTRAINT)
+ {
+ gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.",
+ epull_names[pcrd->eType],
+ epullg_names[pcrd->eGeom],
+ epull_names[epullUMBRELLA]);
+ }
+ }
+ else
+ {
+ /* We only need to calculate the plain COM of a group
+ * when it is not only used as a cylinder group.
+ */
+ if (pull->group[pcrd->group[0]].nat > 0)
+ {
+ pull->group[pcrd->group[0]].bCalcCOM = TRUE;
+ }
+ }
+ if (pull->group[pcrd->group[1]].nat > 0)
+ {
+ pull->group[pcrd->group[1]].bCalcCOM = TRUE;
+ }
+ }
+
pull->ePBC = ir->ePBC;
switch (pull->ePBC)
{
}
}
- fprintf(fplog, "\nWill apply %s COM pulling in geometry '%s'\n",
- EPULLTYPE(ir->ePull), EPULLGEOM(pull->eGeom));
+ fprintf(fplog, "\n");
+ if (pull->bPotential)
+ {
+ fprintf(fplog, "Will apply potential COM pulling\n");
+ }
+ if (pull->bConstraint)
+ {
+ fprintf(fplog, "Will apply constraint COM pulling\n");
+ }
fprintf(fplog, "with %d pull coordinate%s and %d group%s\n",
pull->ncoord, pull->ncoord == 1 ? "" : "s",
pull->ngroup, pull->ngroup == 1 ? "" : "s");
}
}
- /* We always add the virial contribution,
- * except for geometry = direction_periodic where this is impossible.
- */
- pull->bVirial = (pull->eGeom != epullgDIRPBC);
- if (getenv("GMX_NO_PULLVIR") != NULL)
- {
- if (fplog)
- {
- fprintf(fplog, "Found env. var., will not add the virial contribution of the COM pull forces\n");
- }
- pull->bVirial = FALSE;
- }
-
pull->rbuf = NULL;
pull->dbuf = NULL;
pull->dbuf_cyl = NULL;
pgrp->epgrppbc = epgrppbcNONE;
if (pgrp->nat > 0)
{
+ gmx_bool bConstraint;
+ ivec pulldim, pulldim_con;
+
+ /* There is an issue when a group is used in multiple coordinates
+ * and constraints are applied in different dimensions with atoms
+ * frozen in some, but not all dimensions.
+ * Since there is only one mass array per group, we can't have
+ * frozen/non-frozen atoms for different coords at the same time.
+ * But since this is a very exotic case, we don't check for this.
+ * A warning is printed in init_pull_group_index.
+ */
+ bConstraint = FALSE;
+ clear_ivec(pulldim);
+ clear_ivec(pulldim_con);
+ for (c = 0; c < pull->ncoord; c++)
+ {
+ if (pull->coord[c].group[0] == g ||
+ pull->coord[c].group[1] == g)
+ {
+ for (m = 0; m < DIM; m++)
+ {
+ if (pull->coord[c].dim[m] == 1)
+ {
+ pulldim[m] = 1;
+
+ if (pull->coord[c].eType == epullCONSTRAINT)
+ {
+ bConstraint = TRUE;
+ pulldim_con[m] = 1;
+ }
+ }
+ }
+ }
+ }
+
/* Determine if we need to take PBC into account for calculating
* the COM's of the pull groups.
*/
for (m = 0; m < pull->npbcdim; m++)
{
- if (pull->dim[m] && pgrp->nat > 1)
+ if (pulldim[m] == 1 && pgrp->nat > 1)
{
if (pgrp->pbcatom >= 0)
{
pgrp->epgrppbc = epgrppbcCOS;
if (pull->cosdim >= 0 && pull->cosdim != m)
{
- gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull_dim)");
+ gmx_fatal(FARGS, "Can only use cosine weighting with pulling in one dimension (use mdp option pull-coord?-dim)");
}
pull->cosdim = m;
}
}
}
/* Set the indices */
- 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");
- }
+ init_pull_group_index(fplog, cr, g, pgrp,
+ bConstraint, pulldim_con,
+ mtop, ir, lambda);
}
else
{
- /* Absolute reference, set the inverse mass to zero */
+ /* Absolute reference, set the inverse mass to zero.
+ * This is only relevant (and used) with constraint pulling.
+ */
pgrp->invtm = 0;
pgrp->wscale = 1;
}
}
- /* if we use dynamic reference groups, do some initialising for them */
- if (PULL_CYL(pull))
+ /* If we use cylinder coordinates, do some initialising for them */
+ if (pull->bCylinder)
{
- if (ir->ePull == epullCONSTRAINT && pull->ncoord > 1)
- {
- /* 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->ncoord);
for (c = 0; c < pull->ncoord; c++)
{
- if (pull->group[pull->coord[c].group[0]].nat == 0)
+ const t_pull_coord *pcrd;
+
+ pcrd = &pull->coord[c];
+
+ if (pcrd->eGeom == epullgCYL)
{
- gmx_fatal(FARGS, "Dynamic reference groups are not supported when using absolute reference!\n");
+ if (pull->group[pcrd->group[0]].nat == 0)
+ {
+ gmx_fatal(FARGS, "A cylinder pull group is not supported when using absolute reference!\n");
+ }
}
}
-
- snew(pull->dyna, pull->ncoord);
}
/* We still need to initialize the PBC reference coordinates */
{
if (pull->nstxout > 0)
{
- pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv, TRUE, Flags);
+ pull->out_x = open_pull_out(opt2fn("-px", nfile, fnm), pull, oenv,
+ TRUE, Flags);
}
if (pull->nstfout > 0)
{
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
/*! \brief Determine the COM pull forces and add them to f, return the potential
*
- * \param[in] ePull Enum defining the type of pulling: umbrella, const force, ...
* \param[in] pull The pull group.
* \param[in] md All atoms.
* \param[in] pbc Information struct about periodicity.
*
* \returns The pull potential energy.
*/
-real pull_potential(int ePull, t_pull *pull, t_mdatoms *md, struct t_pbc *pbc,
+real pull_potential(t_pull *pull, t_mdatoms *md, struct t_pbc *pbc,
t_commrec *cr, double t, real lambda,
rvec *x, rvec *f, tensor vir, real *dvdlambda);
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
rvec *x,
rvec x_pbc)
{
- int a, m;
+ int a;
if (cr != NULL && DOMAINDECOMP(cr))
{
n = 0;
for (g = 0; g < pull->ngroup; g++)
{
- if ((g == 0 && PULL_CYL(pull)) || pull->group[g].pbcatom == -1)
+ if (!pull->group[g].bCalcCOM || pull->group[g].pbcatom == -1)
{
clear_rvec(x_pbc[g]);
}
else
{
pull_set_pbcatom(cr, &pull->group[g], x, x_pbc[g]);
- for (m = 0; m < DIM; m++)
- {
- if (pull->dim[m] == 0)
- {
- x_pbc[g][m] = 0.0;
- }
- }
n++;
}
}
}
}
-/* switch function, x between r and w */
-static real get_weight(real x, real r1, real r0)
-{
- real weight;
-
- if (x >= r0)
- {
- weight = 0;
- }
- else if (x <= r1)
- {
- weight = 1;
- }
- else
- {
- weight = (r0 - x)/(r0 - r1);
- }
-
- return weight;
-}
-
static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md,
- t_pbc *pbc, double t, rvec *x, rvec *xp)
+ t_pbc *pbc, double t, rvec *x)
{
+ /* The size and stride per coord for the reduction buffer */
+ const int stride = 9;
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;
+ double inv_cyl_r2;
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->ncoord*4);
+ snew(pull->dbuf_cyl, pull->ncoord*stride);
}
if (cr && DOMAINDECOMP(cr))
start = 0;
end = md->homenr;
- r0_2 = dsqr(pull->cyl_r0);
+ inv_cyl_r2 = 1/dsqr(pull->cylinder_r);
/* loop over all groups to make a reference group for each*/
for (c = 0; c < pull->ncoord; c++)
{
- pcrd = &pull->coord[c];
+ double sum_a, wmass, wwmass;
+ dvec radf_fac0, radf_fac1;
- /* 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;
- wwmass = 0;
- pdyna->nat_loc = 0;
-
- for (m = 0; m < DIM; m++)
- {
- g_x[m] = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t);
- }
+ pcrd = &pull->coord[c];
+
+ sum_a = 0;
+ wmass = 0;
+ wwmass = 0;
+ clear_dvec(radf_fac0);
+ clear_dvec(radf_fac1);
- /* loop over all atoms in the main ref group */
- for (i = 0; i < pref->nat; i++)
+ if (pcrd->eGeom == epullgCYL)
{
- ii = pref->ind[i];
- if (ga2la)
+ /* 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);
+ pdyna->nat_loc = 0;
+
+ /* We calculate distances with respect to the reference location
+ * of this cylinder group (g_x), which we already have now since
+ * we reduced the other group COM over the ranks. This resolves
+ * any PBC issues and we don't need to use a PBC-atom here.
+ */
+ for (m = 0; m < DIM; m++)
{
- if (!ga2la_get_home(ga2la, pref->ind[i], &ii))
- {
- ii = -1;
- }
+ g_x[m] = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t);
}
- if (ii >= start && ii < end)
+
+ /* loop over all atoms in the main ref group */
+ for (i = 0; i < pref->nat; i++)
{
- pbc_dx_aiuc(pbc, x[ii], g_x, dx);
- inp = iprod(dir, dx);
- dr2 = 0;
- for (m = 0; m < DIM; m++)
+ ii = pref->ind[i];
+ if (ga2la)
{
- dr2 += dsqr(dx[m] - inp*dir[m]);
+ if (!ga2la_get_home(ga2la, pref->ind[i], &ii))
+ {
+ ii = -1;
+ }
}
-
- if (dr2 < r0_2)
+ if (ii >= start && ii < end)
{
- /* add to index, to sum of COM, to weight array */
- if (pdyna->nat_loc >= pdyna->nalloc_loc)
+ double dr2, dr2_rel, inp;
+ dvec dr;
+
+ pbc_dx_aiuc(pbc, x[ii], g_x, dx);
+ inp = iprod(dir, dx);
+ dr2 = 0;
+ for (m = 0; m < DIM; m++)
{
- pdyna->nalloc_loc = over_alloc_large(pdyna->nat_loc+1);
- srenew(pdyna->ind_loc, pdyna->nalloc_loc);
- srenew(pdyna->weight_loc, pdyna->nalloc_loc);
+ /* Determine the radial components */
+ dr[m] = dx[m] - inp*dir[m];
+ dr2 += dr[m]*dr[m];
}
- pdyna->ind_loc[pdyna->nat_loc] = ii;
- mass = md->massT[ii];
- weight = get_weight(sqrt(dr2), pull->cyl_r1, pull->cyl_r0);
- pdyna->weight_loc[pdyna->nat_loc] = weight;
- sum_a += mass*weight*inp;
- if (xp)
+ dr2_rel = dr2*inv_cyl_r2;
+
+ if (dr2_rel < 1)
{
- pbc_dx_aiuc(pbc, xp[ii], g_x, dx);
- inp = iprod(dir, dx);
- sum_ap += mass*weight*inp;
+ double mass, weight, dweight_r;
+ dvec mdw;
+
+ /* add to index, to sum of COM, to weight array */
+ if (pdyna->nat_loc >= pdyna->nalloc_loc)
+ {
+ pdyna->nalloc_loc = over_alloc_large(pdyna->nat_loc+1);
+ srenew(pdyna->ind_loc, pdyna->nalloc_loc);
+ srenew(pdyna->weight_loc, pdyna->nalloc_loc);
+ srenew(pdyna->mdw, pdyna->nalloc_loc);
+ srenew(pdyna->dv, pdyna->nalloc_loc);
+ }
+ pdyna->ind_loc[pdyna->nat_loc] = ii;
+
+ mass = md->massT[ii];
+ /* The radial weight function is 1-2x^2+x^4,
+ * where x=r/cylinder_r. Since this function depends
+ * on the radial component, we also get radial forces
+ * on both groups.
+ */
+ weight = 1 + (-2 + dr2_rel)*dr2_rel;
+ dweight_r = (-4 + 4*dr2_rel)*inv_cyl_r2;
+ pdyna->weight_loc[pdyna->nat_loc] = weight;
+ sum_a += mass*weight*inp;
+ wmass += mass*weight;
+ wwmass += mass*weight*weight;
+ dsvmul(mass*dweight_r, dr, mdw);
+ copy_dvec(mdw, pdyna->mdw[pdyna->nat_loc]);
+ /* Currently we only have the axial component of the
+ * distance (inp) up to an unkown offset. We add this
+ * offset after the reduction needs to determine the
+ * COM of the cylinder group.
+ */
+ pdyna->dv[pdyna->nat_loc] = inp;
+ for (m = 0; m < DIM; m++)
+ {
+ radf_fac0[m] += mdw[m];
+ radf_fac1[m] += mdw[m]*inp;
+ }
+ pdyna->nat_loc++;
}
- wmass += mass*weight;
- wwmass += mass*sqr(weight);
- pdyna->nat_loc++;
}
}
}
- 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;
+ pull->dbuf_cyl[c*stride+0] = wmass;
+ pull->dbuf_cyl[c*stride+1] = wwmass;
+ pull->dbuf_cyl[c*stride+2] = sum_a;
+ pull->dbuf_cyl[c*stride+3] = radf_fac0[XX];
+ pull->dbuf_cyl[c*stride+4] = radf_fac0[YY];
+ pull->dbuf_cyl[c*stride+5] = radf_fac0[ZZ];
+ pull->dbuf_cyl[c*stride+6] = radf_fac1[XX];
+ pull->dbuf_cyl[c*stride+7] = radf_fac1[YY];
+ pull->dbuf_cyl[c*stride+8] = radf_fac1[ZZ];
}
- if (cr && PAR(cr))
+ if (cr != NULL && PAR(cr))
{
- /* Sum the contributions over the nodes */
- gmx_sumd(pull->ncoord*4, pull->dbuf_cyl, cr);
+ /* Sum the contributions over the ranks */
+ gmx_sumd(pull->ncoord*stride, pull->dbuf_cyl, cr);
}
for (c = 0; c < pull->ncoord; c++)
{
pcrd = &pull->coord[c];
- pdyna = &pull->dyna[c];
- pgrp = &pull->group[pcrd->group[1]];
+ if (pcrd->eGeom == epullgCYL)
+ {
+ double wmass, wwmass, inp, dist;
- 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);
+ pdyna = &pull->dyna[c];
+ pgrp = &pull->group[pcrd->group[1]];
- for (m = 0; m < DIM; m++)
- {
- 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)
+ wmass = pull->dbuf_cyl[c*stride+0];
+ wwmass = pull->dbuf_cyl[c*stride+1];
+ pdyna->mwscale = 1.0/wmass;
+ /* Cylinder pulling can't be used with constraints, but we set
+ * wscale and invtm anyhow, in case someone would like to use them.
+ */
+ pdyna->wscale = wmass/wwmass;
+ pdyna->invtm = wwmass/(wmass*wmass);
+
+ /* We store the deviation of the COM from the reference location
+ * used above, since we need it when we apply the radial forces
+ * to the atoms in the cylinder group.
+ */
+ pcrd->cyl_dev = 0;
+ for (m = 0; m < DIM; m++)
{
- pdyna->xp[m] = g_x[m] + pcrd->vec[m]*pull->dbuf_cyl[c*4+3]/wmass;
+ g_x[m] = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t);
+ dist = -pcrd->vec[m]*pull->dbuf_cyl[c*stride+2]*pdyna->mwscale;
+ pdyna->x[m] = g_x[m] - dist;
+ pcrd->cyl_dev += dist;
+ }
+ /* Now we know the exact COM of the cylinder reference group,
+ * we can determine the radial force factor (ffrad) that when
+ * multiplied with the axial pull force will give the radial
+ * force on the pulled (non-cylinder) group.
+ */
+ for (m = 0; m < DIM; m++)
+ {
+ pcrd->ffrad[m] = (pull->dbuf_cyl[c*stride+6+m] +
+ pull->dbuf_cyl[c*stride+3+m]*pcrd->cyl_dev)/wmass;
}
- }
- if (debug)
- {
- fprintf(debug, "Pull cylinder group %d:%8.3f%8.3f%8.3f m:%8.3f\n",
- c, pdyna->x[0], pdyna->x[1],
- pdyna->x[2], 1.0/pdyna->invtm);
+ if (debug)
+ {
+ fprintf(debug, "Pull cylinder group %d:%8.3f%8.3f%8.3f m:%8.3f\n",
+ c, pdyna->x[0], pdyna->x[1],
+ pdyna->x[2], 1.0/pdyna->invtm);
+ fprintf(debug, "ffrad %8.3f %8.3f %8.3f\n",
+ pcrd->ffrad[XX], pcrd->ffrad[YY], pcrd->ffrad[ZZ]);
+ }
}
}
}
{
int g, i, ii, m;
real mass, w, wm, twopi_box = 0;
- double wmass, wwmass, invwmass;
+ double wmass, wwmass;
dvec com, comp;
double cm, sm, cmp, smp, ccm, csm, ssm, csw, snw;
- rvec *xx[2], x_pbc = {0, 0, 0}, dx;
+ rvec x_pbc = {0, 0, 0}, dx;
t_pull_group *pgrp;
if (pull->rbuf == NULL)
ccm = 0;
csm = 0;
ssm = 0;
- if (!(g == 0 && PULL_CYL(pull)))
+ if (pgrp->bCalcCOM)
{
if (pgrp->epgrppbc == epgrppbcREFAT)
{
for (g = 0; g < pull->ngroup; g++)
{
pgrp = &pull->group[g];
- if (pgrp->nat > 0 && !(g == 0 && PULL_CYL(pull)))
+ if (pgrp->nat > 0 && pgrp->bCalcCOM)
{
if (pgrp->epgrppbc != epgrppbcCOS)
{
/* Determine the inverse mass */
- wmass = pull->dbuf[g*3+2][0];
- wwmass = pull->dbuf[g*3+2][1];
- invwmass = 1/wmass;
+ wmass = pull->dbuf[g*3+2][0];
+ wwmass = pull->dbuf[g*3+2][1];
+ pgrp->mwscale = 1.0/wmass;
/* invtm==0 signals a frozen group, so then we should keep it zero */
- if (pgrp->invtm > 0)
+ if (pgrp->invtm != 0)
{
- pgrp->wscale = wmass/wwmass;
- pgrp->invtm = 1.0/(pgrp->wscale*wmass);
+ pgrp->wscale = wmass/wwmass;
+ pgrp->invtm = wwmass/(wmass*wmass);
}
/* Divide by the total mass */
for (m = 0; m < DIM; m++)
{
- pgrp->x[m] = pull->dbuf[g*3 ][m]*invwmass;
+ pgrp->x[m] = pull->dbuf[g*3 ][m]*pgrp->mwscale;
if (xp)
{
- pgrp->xp[m] = pull->dbuf[g*3+1][m]*invwmass;
+ pgrp->xp[m] = pull->dbuf[g*3+1][m]*pgrp->mwscale;
}
if (pgrp->epgrppbc == epgrppbcREFAT)
{
wwmass = (pull->dbuf[g*3+1][0]*csw*csw +
pull->dbuf[g*3+1][1]*csw*snw +
pull->dbuf[g*3+1][2]*snw*snw)/(wmass*wmass);
- pgrp->wscale = wmass/wwmass;
- pgrp->invtm = 1.0/(pgrp->wscale*wmass);
+
+ pgrp->mwscale = 1.0/wmass;
+ pgrp->wscale = wmass/wwmass;
+ pgrp->invtm = wwmass/(wmass*wmass);
/* Set the weights for the local atoms */
csw *= pgrp->invtm;
snw *= pgrp->invtm;
}
}
- if (PULL_CYL(pull))
+ if (pull->bCylinder)
{
/* Calculate the COMs for the cyclinder reference groups */
- make_cyl_refgrps(cr, pull, md, pbc, t, x, xp);
+ make_cyl_refgrps(cr, pull, md, pbc, t, x);
}
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
cmp_real(fp, "inputrec->wall_density[1]", -1, ir1->wall_density[1], ir2->wall_density[1], ftol, abstol);
cmp_real(fp, "inputrec->wall_ewald_zfac", -1, ir1->wall_ewald_zfac, ir2->wall_ewald_zfac, ftol, abstol);
- cmp_int(fp, "inputrec->ePull", -1, ir1->ePull, ir2->ePull);
- if (ir1->ePull == ir2->ePull && ir1->ePull != epullNO)
+ cmp_bool(fp, "inputrec->bPull", -1, ir1->bPull, ir2->bPull);
+ if (ir1->bPull && ir2->bPull)
{
cmp_pull(fp);
}
}
else
{
- if (ir[0].ePull != epullNO)
+ if (ir[0].bPull)
{
comp_pull_AB(stdout, ir->pull, ftol, abstol);
}
step, t,
eprNORMAL, bCompact, mdebin, fcd, groups, &(ir->opts));
}
- if (ir->ePull != epullNO)
+ if (ir->bPull)
{
pull_print_output(ir->pull, step, t);
}
/* Assumes uniform use of the number of OpenMP threads */
walltime_accounting = walltime_accounting_init(gmx_omp_nthreads_get(emntDefault));
- if (inputrec->ePull != epullNO)
+ if (inputrec->bPull)
{
/* Initialize pull code */
init_pull(fplog, inputrec, nfile, fnm, mtop, cr, oenv, inputrec->fepvals->init_lambda,
Flags,
walltime_accounting);
- if (inputrec->ePull != epullNO)
+ if (inputrec->bPull)
{
finish_pull(inputrec->pull);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff