From dbe80640f78830b973220930b4b8d8fe064e42f4 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Tue, 26 Aug 2014 14:56:09 +0200 Subject: [PATCH] COM pulling options per coord, improved cylinder The COM pull type (potential/constraint), geometry, dimensions and pull-start can now be set differently for each pull coordinate. Also added a flat-bottom potential. Cylinder geometry now uses a smooth radial weight function and adds radial forces so it now conserves energy. Fixes #1590. There is now also an option for printing the COM of the second group, as well as the reference value for all pull coordinates. Writing the distance components to the pull coordinate output is now optional (off by default). gmx wham currently only supports input with the old restrictions (and it checks for this). These restrctions should be removed. Change-Id: I86e443b225accb85cc7fe0aa6c23d937f15b2efb --- docs/manual/special.tex | 27 +- docs/user-guide/mdp-options.rst | 192 ++++--- src/gromacs/domdec/domdec.cpp | 2 +- src/gromacs/fileio/tpxio.c | 96 +++- src/gromacs/gmxana/gmx_tune_pme.c | 4 +- src/gromacs/gmxana/gmx_wham.cpp | 116 +++-- src/gromacs/gmxlib/names.c | 4 +- src/gromacs/gmxlib/txtdump.c | 20 +- src/gromacs/gmxpreprocess/grompp.c | 4 +- src/gromacs/gmxpreprocess/readir.c | 63 +-- src/gromacs/gmxpreprocess/readir.h | 9 +- src/gromacs/gmxpreprocess/readpull.c | 175 ++++--- src/gromacs/legacyheaders/types/enums.h | 6 +- src/gromacs/legacyheaders/types/inputrec.h | 48 +- src/gromacs/mdlib/broadcaststructs.cpp | 4 +- src/gromacs/mdlib/constr.cpp | 8 +- src/gromacs/mdlib/mdebin.c | 4 +- src/gromacs/mdlib/sim_util.cpp | 10 +- src/gromacs/pulling/pull.c | 576 ++++++++++++++------- src/gromacs/pulling/pull.h | 5 +- src/gromacs/pulling/pullutil.c | 293 ++++++----- src/gromacs/tools/compare.c | 8 +- src/programs/mdrun/md.cpp | 2 +- src/programs/mdrun/runner.cpp | 4 +- 24 files changed, 1078 insertions(+), 602 deletions(-) diff --git a/docs/manual/special.tex b/docs/manual/special.tex index ec8a53142e..d79a680e58 100644 --- a/docs/manual/special.tex +++ b/docs/manual/special.tex @@ -214,6 +214,10 @@ exact if only two groups are constrained. 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} @@ -234,16 +238,25 @@ to calculate the PMF of a lipid as function of its distance 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}} diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst index 86fdfef610..b1585e6bf7 100644 --- a/docs/user-guide/mdp-options.rst +++ b/docs/user-guide/mdp-options.rst @@ -1583,112 +1583,75 @@ applicable pulling coordinate. 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 @@ -1741,6 +1704,82 @@ applicable pulling coordinate. 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) @@ -1751,6 +1790,17 @@ applicable pulling coordinate. (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\] diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp index e593445135..00b643ea05 100644 --- a/src/gromacs/domdec/domdec.cpp +++ b/src/gromacs/domdec/domdec.cpp @@ -9854,7 +9854,7 @@ void dd_partition_system(FILE *fplog, 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); diff --git a/src/gromacs/fileio/tpxio.c b/src/gromacs/fileio/tpxio.c index 0ca368e5c7..84cb9d478e 100644 --- a/src/gromacs/fileio/tpxio.c +++ b/src/gromacs/fileio/tpxio.c @@ -90,7 +90,8 @@ enum tpxv { 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 @@ -104,7 +105,7 @@ enum tpxv { * * 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 @@ -287,14 +288,36 @@ static void do_pull_group(t_fileio *fio, t_pull_group *pgrp, gmx_bool bRead) 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); @@ -603,9 +626,12 @@ static void do_fepvals(t_fileio *fio, t_lambda *fepvals, gmx_bool bRead, int fil } } -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) { @@ -616,14 +642,33 @@ static void do_pull(t_fileio *fio, t_pull *pull, gmx_bool bRead, int file_versio { 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); @@ -635,13 +680,13 @@ static void do_pull(t_fileio *fio, t_pull *pull, gmx_bool bRead, int file_versio 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++) @@ -656,7 +701,7 @@ static void do_pull(t_fileio *fio, t_pull *pull, gmx_bool bRead, int file_versio } } - pull->bPrintRef = (pull->group[0].nat > 0); + pull->bPrintCOM1 = (pull->group[0].nat > 0); } else { @@ -666,7 +711,8 @@ static void do_pull(t_fileio *fio, t_pull *pull, gmx_bool bRead, int file_versio } 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); } } } @@ -1495,19 +1541,31 @@ static void do_inputrec(t_fileio *fio, t_inputrec *ir, gmx_bool bRead, /* 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 */ diff --git a/src/gromacs/gmxana/gmx_tune_pme.c b/src/gromacs/gmxana/gmx_tune_pme.c index d5b4c667ed..4c1f7745b0 100644 --- a/src/gromacs/gmxana/gmx_tune_pme.c +++ b/src/gromacs/gmxana/gmx_tune_pme.c @@ -2043,7 +2043,6 @@ static void setopt(const char *opt, int nfile, t_filenm fnm[]) * 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? */ @@ -2055,14 +2054,13 @@ static float inspect_tpr(int nfile, t_filenm fnm[], real *rcoulomb) /* 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); diff --git a/src/gromacs/gmxana/gmx_wham.cpp b/src/gromacs/gmxana/gmx_wham.cpp index b3ef53bcfd..a1968ef659 100644 --- a/src/gromacs/gmxana/gmx_wham.cpp +++ b/src/gromacs/gmxana/gmx_wham.cpp @@ -120,11 +120,13 @@ typedef struct * \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 @@ -2002,24 +2004,70 @@ void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions /* 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); @@ -2120,7 +2168,11 @@ void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header, * 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) @@ -2160,17 +2212,28 @@ void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header, 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) @@ -2324,20 +2387,10 @@ void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header, 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); */ @@ -3129,7 +3182,12 @@ int gmx_wham(int argc, char *argv[]) "[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.", "", diff --git a/src/gromacs/gmxlib/names.c b/src/gromacs/gmxlib/names.c index 6e63f2fb62..3bae62a817 100644 --- a/src/gromacs/gmxlib/names.c +++ b/src/gromacs/gmxlib/names.c @@ -3,7 +3,7 @@ * * 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. @@ -219,7 +219,7 @@ const char *ewt_names[ewtNR+1] = { }; 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] = { diff --git a/src/gromacs/gmxlib/txtdump.c b/src/gromacs/gmxlib/txtdump.c index c22bdb311d..ff8bbf07fb 100644 --- a/src/gromacs/gmxlib/txtdump.c +++ b/src/gromacs/gmxlib/txtdump.c @@ -3,7 +3,7 @@ * * 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. @@ -634,8 +634,12 @@ static void pr_pull_coord(FILE *fp, int indent, int c, t_pull_coord *pcrd) 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); @@ -755,12 +759,12 @@ static void pr_pull(FILE *fp, int indent, t_pull *pull) { 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); @@ -1024,8 +1028,8 @@ void pr_inputrec(FILE *fp, int indent, const char *title, t_inputrec *ir, 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); } diff --git a/src/gromacs/gmxpreprocess/grompp.c b/src/gromacs/gmxpreprocess/grompp.c index c7cfab8f7f..83b4638ae4 100644 --- a/src/gromacs/gmxpreprocess/grompp.c +++ b/src/gromacs/gmxpreprocess/grompp.c @@ -2034,9 +2034,9 @@ int gmx_grompp(int argc, char *argv[]) } } - 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) diff --git a/src/gromacs/gmxpreprocess/readir.c b/src/gromacs/gmxpreprocess/readir.c index d3c06878d8..4e69abd8b3 100644 --- a/src/gromacs/gmxpreprocess/readir.c +++ b/src/gromacs/gmxpreprocess/readir.c @@ -2082,12 +2082,11 @@ void get_ir(const char *mdparin, const char *mdparout, /* 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 */ @@ -2840,7 +2839,7 @@ static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames) } } - 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 @@ -2852,6 +2851,11 @@ static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames) for (i = 0; i < pull->ncoord; i++) { + if (pull->coord[i].eType != epullCONSTRAINT) + { + continue; + } + imin = 1; for (j = 0; j < 2; j++) @@ -3477,7 +3481,7 @@ void do_index(const char* mdparin, const char *ndx, } } - if (ir->ePull != epullNO) + if (ir->bPull) { make_pull_groups(ir->pull, is->pull_grp, grps, gnames); @@ -4235,45 +4239,42 @@ void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys, 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); } } } diff --git a/src/gromacs/gmxpreprocess/readir.h b/src/gromacs/gmxpreprocess/readir.h index dcfd3fc2fc..49b6de6f73 100644 --- a/src/gromacs/gmxpreprocess/readir.h +++ b/src/gromacs/gmxpreprocess/readir.h @@ -3,7 +3,7 @@ * * 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. @@ -66,7 +66,6 @@ typedef struct { gmx_bool bOrire; gmx_bool bMorse; char *wall_atomtype[2]; - gmx_bool pull_start; char *couple_moltype; int couple_lam0; int couple_lam1; @@ -128,7 +127,7 @@ void do_index(const char* mdparin, /* 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 */ @@ -141,9 +140,9 @@ void make_pull_coords(t_pull *pull); /* 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[]); diff --git a/src/gromacs/gmxpreprocess/readpull.c b/src/gromacs/gmxpreprocess/readpull.c index eeab07e031..dde971b5c5 100644 --- a/src/gromacs/gmxpreprocess/readpull.c +++ b/src/gromacs/gmxpreprocess/readpull.c @@ -3,7 +3,7 @@ * * 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. @@ -54,8 +54,6 @@ #include "gromacs/utility/smalloc.h" -static char pulldim[STRLEN]; - static void string2dvec(const char buf[], dvec nums) { double dum; @@ -84,7 +82,46 @@ static void init_pull_group(t_pull_group *pg, } } -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) { @@ -92,18 +129,28 @@ static void init_pull_coord(t_pull_coord *pcrd, int eGeom, 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 @@ -119,7 +166,12 @@ static void init_pull_coord(t_pull_coord *pcrd, int eGeom, 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); @@ -133,14 +185,15 @@ static void init_pull_coord(t_pull_coord *pcrd, int eGeom, } 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]; @@ -151,29 +204,20 @@ char **read_pullparams(int *ninp_p, t_inpfile **inp_p, 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"); @@ -224,21 +268,29 @@ char **read_pullparams(int *ninp_p, t_inpfile **inp_p, 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; @@ -285,10 +337,6 @@ void make_pull_groups(t_pull *pull, 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)", @@ -321,41 +369,9 @@ void make_pull_groups(t_pull *pull, 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]; @@ -371,33 +387,36 @@ void make_pull_coords(t_pull *pull) 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; @@ -437,7 +456,7 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l fprintf(stderr, "%8d %8d %8d ", pcrd->group[1], pgrp1->nat, pgrp1->pbcatom+1); - if (bStart) + if (pcrd->bStart) { init = pcrd->init; pcrd->init = 0; @@ -448,7 +467,7 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l value = pcrd->init + t_start*pcrd->rate + dev; fprintf(stderr, " %10.3f nm", value); - if (bStart) + if (pcrd->bStart) { pcrd->init = value + init; } diff --git a/src/gromacs/legacyheaders/types/enums.h b/src/gromacs/legacyheaders/types/enums.h index 8756469d3d..92bdcb8237 100644 --- a/src/gromacs/legacyheaders/types/enums.h +++ b/src/gromacs/legacyheaders/types/enums.h @@ -3,7 +3,7 @@ * * 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. @@ -291,15 +291,13 @@ enum { /* 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, diff --git a/src/gromacs/legacyheaders/types/inputrec.h b/src/gromacs/legacyheaders/types/inputrec.h index cbfc3e3110..0bd0c85cd7 100644 --- a/src/gromacs/legacyheaders/types/inputrec.h +++ b/src/gromacs/legacyheaders/types/inputrec.h @@ -115,16 +115,24 @@ typedef struct { 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 */ @@ -132,6 +140,8 @@ typedef struct { /* 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; @@ -206,25 +216,27 @@ typedef struct { } 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? */ @@ -444,7 +456,7 @@ typedef struct { 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 */ diff --git a/src/gromacs/mdlib/broadcaststructs.cpp b/src/gromacs/mdlib/broadcaststructs.cpp index 98b476bc85..c13b7185fe 100644 --- a/src/gromacs/mdlib/broadcaststructs.cpp +++ b/src/gromacs/mdlib/broadcaststructs.cpp @@ -3,7 +3,7 @@ * * 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. @@ -687,7 +687,7 @@ static void bc_inputrec(const t_commrec *cr, t_inputrec *inputrec) { 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); diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index ad36ec2a7f..ef655c062d 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -3,7 +3,7 @@ * * 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. @@ -614,7 +614,7 @@ gmx_bool constrain(FILE *fplog, gmx_bool bLog, gmx_bool bEner, if (econq == econqCoord) { - if (ir->ePull == epullCONSTRAINT) + if (ir->bPull && ir->pull->bConstraint) { if (EI_DYNAMICS(ir->eI)) { @@ -1171,7 +1171,9 @@ gmx_constr_t init_constraints(FILE *fplog, 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; } diff --git a/src/gromacs/mdlib/mdebin.c b/src/gromacs/mdlib/mdebin.c index 84b1db0f3d..d46d50f637 100644 --- a/src/gromacs/mdlib/mdebin.c +++ b/src/gromacs/mdlib/mdebin.c @@ -3,7 +3,7 @@ * * 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. @@ -291,7 +291,7 @@ t_mdebin *init_mdebin(ener_file_t fp_ene, } 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) { diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 2df9b2a73e..38d043e28f 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -331,7 +331,7 @@ static void pull_potential_wrapper(t_commrec *cr, 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); @@ -1154,7 +1154,7 @@ void do_force_cutsVERLET(FILE *fplog, t_commrec *cr, clear_rvec(fr->vir_diag_posres); } - if (inputrec->ePull == epullCONSTRAINT) + if (inputrec->bPull && inputrec->pull->bConstraint) { clear_pull_forces(inputrec->pull); } @@ -1456,7 +1456,7 @@ void do_force_cutsVERLET(FILE *fplog, t_commrec *cr, } } - 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. @@ -1799,7 +1799,7 @@ void do_force_cutsGROUP(FILE *fplog, t_commrec *cr, clear_rvec(fr->vir_diag_posres); } - if (inputrec->ePull == epullCONSTRAINT) + if (inputrec->bPull && inputrec->pull->bConstraint) { clear_pull_forces(inputrec->pull); } @@ -1920,7 +1920,7 @@ void do_force_cutsGROUP(FILE *fplog, t_commrec *cr, } } - 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, diff --git a/src/gromacs/pulling/pull.c b/src/gromacs/pulling/pull.c index 575e07ca87..fca9f171c9 100644 --- a/src/gromacs/pulling/pull.c +++ b/src/gromacs/pulling/pull.c @@ -3,7 +3,7 @@ * * 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. @@ -73,23 +73,61 @@ static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp) } } -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); @@ -97,25 +135,30 @@ static void pull_print_x(FILE *out, t_pull *pull, double 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); @@ -164,36 +207,73 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv 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++; @@ -220,7 +300,7 @@ static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md, 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++) { @@ -238,6 +318,40 @@ static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md, } } +/* 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) { @@ -248,9 +362,20 @@ 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 { @@ -258,26 +383,23 @@ static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *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); } - 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])); } } @@ -309,7 +431,7 @@ static void low_get_pull_coord_dr(const t_pull *pull, copy_dvec(xref, xrefr); - if (pull->eGeom == epullgDIRPBC) + if (pcrd->eGeom == epullgDIRPBC) { for (m = 0; m < DIM; m++) { @@ -323,7 +445,7 @@ static void low_get_pull_coord_dr(const t_pull *pull, 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) @@ -332,7 +454,7 @@ static void low_get_pull_coord_dr(const t_pull *pull, pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2)); } - if (pull->eGeom == epullgDIRPBC) + if (pcrd->eGeom == epullgDIRPBC) { dvec_inc(dr, dref); } @@ -346,20 +468,20 @@ static void get_pull_coord_dr(const t_pull *pull, 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); } @@ -381,7 +503,7 @@ void get_pull_coord_distance(const t_pull *pull, ref = pcrd->init + pcrd->rate*t; - switch (pull->eGeom) + switch (pcrd->eGeom) { case epullgDIST: /* Pull along the vector between the com's */ @@ -446,14 +568,14 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, 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); @@ -468,35 +590,37 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, { 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]; } } @@ -517,6 +641,12 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, dvec dr0, dr1; pcrd = &pull->coord[c]; + + if (pcrd->eType != epullCONSTRAINT) + { + continue; + } + pgrp0 = &pull->group[pcrd->group[0]]; pgrp1 = &pull->group[pcrd->group[1]]; @@ -535,7 +665,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, rm = 1.0/(pgrp0->invtm + pgrp1->invtm); - switch (pull->eGeom) + switch (pcrd->eGeom) { case epullgDIST: if (ref <= 0) @@ -631,6 +761,12 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, 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, @@ -638,7 +774,7 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, rnew[pcrd->group[0]], -1, unc_ij); - switch (pull->eGeom) + switch (pcrd->eGeom) { case epullgDIST: bConverged = fabs(dnorm(unc_ij) - ref) < pull->constr_tol; @@ -692,21 +828,15 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, 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 */ @@ -735,10 +865,16 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, /* 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; @@ -763,12 +899,11 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, } /* 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; @@ -779,73 +914,83 @@ static void do_pull_pot(int ePull, { 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++) @@ -859,7 +1004,7 @@ static void do_pull_pot(int ePull, } } -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) { @@ -867,8 +1012,8 @@ real pull_potential(int ePull, t_pull *pull, t_mdatoms *md, t_pbc *pbc, 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); @@ -887,7 +1032,7 @@ void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc, { 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, @@ -953,7 +1098,8 @@ void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md) } 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; @@ -1009,11 +1155,12 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, { 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++; } @@ -1092,22 +1239,22 @@ static void init_pull_group_index(FILE *fplog, t_commrec *cr, 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; @@ -1121,10 +1268,56 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], { 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) { @@ -1147,8 +1340,15 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], } } - 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"); @@ -1173,19 +1373,6 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], } } - /* 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; @@ -1197,12 +1384,47 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], 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) { @@ -1218,47 +1440,46 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], 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 */ @@ -1271,7 +1492,8 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], { 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) { diff --git a/src/gromacs/pulling/pull.h b/src/gromacs/pulling/pull.h index 3a5f8d6b57..7b23de006a 100644 --- a/src/gromacs/pulling/pull.h +++ b/src/gromacs/pulling/pull.h @@ -3,7 +3,7 @@ * * 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. @@ -83,7 +83,6 @@ void clear_pull_forces(t_pull *pull); /*! \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. @@ -97,7 +96,7 @@ void clear_pull_forces(t_pull *pull); * * \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); diff --git a/src/gromacs/pulling/pullutil.c b/src/gromacs/pulling/pullutil.c index 53072c1ef2..5f231846dd 100644 --- a/src/gromacs/pulling/pullutil.c +++ b/src/gromacs/pulling/pullutil.c @@ -3,7 +3,7 @@ * * 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. @@ -56,7 +56,7 @@ static void pull_set_pbcatom(t_commrec *cr, t_pull_group *pgrp, rvec *x, rvec x_pbc) { - int a, m; + int a; if (cr != NULL && DOMAINDECOMP(cr)) { @@ -84,20 +84,13 @@ static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull, 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++; } } @@ -109,40 +102,21 @@ static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull, } } -/* 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)) @@ -153,115 +127,182 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, 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]); + } } } } @@ -285,10 +326,10 @@ void pull_calc_coms(t_commrec *cr, { 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) @@ -343,7 +384,7 @@ void pull_calc_coms(t_commrec *cr, ccm = 0; csm = 0; ssm = 0; - if (!(g == 0 && PULL_CYL(pull))) + if (pgrp->bCalcCOM) { if (pgrp->epgrppbc == epgrppbcREFAT) { @@ -456,27 +497,27 @@ void pull_calc_coms(t_commrec *cr, 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) { @@ -499,8 +540,10 @@ void pull_calc_coms(t_commrec *cr, 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; @@ -525,9 +568,9 @@ void pull_calc_coms(t_commrec *cr, } } - 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); } } diff --git a/src/gromacs/tools/compare.c b/src/gromacs/tools/compare.c index 89038a4c24..7be156ab53 100644 --- a/src/gromacs/tools/compare.c +++ b/src/gromacs/tools/compare.c @@ -3,7 +3,7 @@ * * 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. @@ -836,8 +836,8 @@ static void cmp_inputrec(FILE *fp, t_inputrec *ir1, t_inputrec *ir2, real ftol, 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); } @@ -1010,7 +1010,7 @@ void comp_tpx(const char *fn1, const char *fn2, } else { - if (ir[0].ePull != epullNO) + if (ir[0].bPull) { comp_pull_AB(stdout, ir->pull, ftol, abstol); } diff --git a/src/programs/mdrun/md.cpp b/src/programs/mdrun/md.cpp index f5168fc2f0..6a92f66f7d 100644 --- a/src/programs/mdrun/md.cpp +++ b/src/programs/mdrun/md.cpp @@ -1625,7 +1625,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[], step, t, eprNORMAL, bCompact, mdebin, fcd, groups, &(ir->opts)); } - if (ir->ePull != epullNO) + if (ir->bPull) { pull_print_output(ir->pull, step, t); } diff --git a/src/programs/mdrun/runner.cpp b/src/programs/mdrun/runner.cpp index 8004e6a7ed..10c6cc47fa 100644 --- a/src/programs/mdrun/runner.cpp +++ b/src/programs/mdrun/runner.cpp @@ -1589,7 +1589,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt, /* 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, @@ -1639,7 +1639,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt, Flags, walltime_accounting); - if (inputrec->ePull != epullNO) + if (inputrec->bPull) { finish_pull(inputrec->pull); } -- 2.22.0