From cc819af8c6691d70a84f76ae622f604f83ebf440 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Tue, 26 Nov 2013 09:41:01 +0100 Subject: [PATCH] Removed restriction of one reference pull group The pull code now uses separate pull groups and coordinates, which removes any restriction on group pairing. The pull geometry "position" has been removed, but the origin option replaces this functionaliy in case of an absolute reference. Removed triple use of the init option. Refs #1346 Change-Id: I90736ac6b6f18e3bf26a344072172a69627566e9 --- manual/special.tex | 12 +- share/html/online/mdp_opt.html | 77 ++- src/gromacs/fileio/tpxio.c | 113 ++- src/gromacs/gmxana/gmx_wham.cpp | 159 ++--- src/gromacs/gmxlib/mvdata.c | 10 +- src/gromacs/gmxlib/names.c | 2 +- src/gromacs/gmxlib/txtdump.c | 39 +- src/gromacs/gmxlib/typedefs.c | 19 +- src/gromacs/gmxpreprocess/readir.c | 78 +-- src/gromacs/gmxpreprocess/readir.h | 12 +- src/gromacs/gmxpreprocess/readpull.c | 475 ++++++------- src/gromacs/legacyheaders/pull.h | 8 +- src/gromacs/legacyheaders/types/enums.h | 2 +- src/gromacs/legacyheaders/types/inputrec.h | 32 +- src/gromacs/mdlib/pull.c | 759 +++++++++------------ src/gromacs/mdlib/pullutil.c | 104 +-- src/programs/gmx/tpbcmp.c | 6 +- 17 files changed, 946 insertions(+), 961 deletions(-) diff --git a/manual/special.tex b/manual/special.tex index db543530d9..000a6ce882 100644 --- a/manual/special.tex +++ b/manual/special.tex @@ -178,12 +178,10 @@ a canonical ensemble of the initial state A and $\beta=1/k_B T$. \label{sec:pull} The pull code applies forces or constraints between the centers of mass of one or more pairs of groups of atoms. -There is one reference group and one or more other pull groups. -Instead of a reference group, one can also use absolute reference -point in space. -The most common situation consists of a reference group -and one pull group. In this case, the two groups are treated -equivalently. +Each pull reaction coordinate is called a ``coordinate'' and it operates +on two pull groups. A pull group can be part of one or more pull +coordinates. Furthermore, a coordinate can also operate on a single group +and an absolute reference position in space. The distance between a pair of groups can be determined in 1, 2 or 3 dimensions, or can be along a user-defined vector. The reference distance can be constant or can change linearly with time. @@ -295,7 +293,7 @@ as follows: \eeq \subsubsection{Limitations} -There is one important limitation: +There is one theoretical limitation: strictly speaking, constraint forces can only be calculated between groups that are not connected by constraints to the rest of the system. If a group contains part of a molecule of which the bond lengths diff --git a/share/html/online/mdp_opt.html b/share/html/online/mdp_opt.html index 97c9e67fc4..0678ad1976 100644 --- a/share/html/online/mdp_opt.html +++ b/share/html/online/mdp_opt.html @@ -1283,7 +1283,7 @@ the minimum is 2. Ewald summation can only be used with nwall=2, where one should use ewald-geometry=3dc. The empty layer in the box serves to decrease the unphysical Coulomb -interaction between periodic images. +interaction between periodic images.
@@ -1334,10 +1334,7 @@ the relative weights are one, between pull-r1 and 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. -
position
-
Pull to the position of the reference group plus -pull-init + time*pull-rate*pull-vec.
+from the COM of the pull group has both a radial and an axial component.
pull-dim: (Y Y Y)
the distance components to be used with geometry distance @@ -1356,62 +1353,72 @@ to the output files
yes
add the COM distance of the starting conformation to pull-init
+
pull-print-reference: (10)
+
+
no
+
do not print the COM of the first group in each pull coordinate
+
yes
+
print the COM of the first group in each pull coordinate
+
pull-nstxout: (10)
frequency for writing out the COMs of all the pull group
pull-nstfout: (1)
frequency for writing out the force of all the pulled group
pull-ngroups: (1)
-
The number of pull groups, not including the reference group. -If there is only one group, there is no difference in treatment -of the reference and pulled group (except with the cylinder geometry). -Below only the pull options for the reference group (ending on 0) -and the first group (ending on 1) are given, -further groups work analogously, but with the number 1 replaced -by the group number.
-
pull-group0:
-
The name of the reference group. When this is empty an absolute reference -of (0,0,0) is used. With an absolute reference the system is no longer -translation invariant and one should think about what to do with -the center of mass motion.
-
pull-weights0:
-
see pull-weights1
-
pull-pbcatom0: (0)
-
see pull-pbcatom1
-
pull-group1:
-
The name of the pull group.
-
pull-weights1:
+
The number of pull groups, not including the absolute reference group, +when used. Pull groups can be reused in multiple pull coordinates. +Below only the pull options for group 1 are given, further groups simply +increase the group index number.
+
pull-ncoords: (1)
+
The number of pull coordinates. Below only the pull options for +coordinate 1 are given, further coordinates simply increase the coordinate +index number.
+ +
pull-group1-name:
+
The name of the pull group, is looked up in the index file +or in the default groups to obtain the atoms involved.
+
pull-group1-weights:
Optional relative weights which are multiplied with the masses of the atoms to give the total weight for the COM. The number should be 0, meaning all 1, or the number of atoms in the pull group.
-
pull-pbcatom1: (0)
+
pull-group1-pbcatom: (0)
The reference atom for the treatment of periodic boundary conditions inside the group (this has no effect on the treatment of the pbc between groups). This option is only important when the diameter of the pull group is larger than half the shortest box vector. For determining the COM, all atoms in the group are put at their periodic image -which is closest to pull-pbcatom1. +which is closest to pull-group1-pbcatom. A value of 0 means that the middle atom (number wise) is used. This parameter is not used with geometry cylinder. A value of -1 turns on cosine weighting, which is useful for a group of molecules in a periodic system, e.g. a water slab (see Engin et al. J. Chem. Phys. B 2010).
-
pull-vec1: (0.0 0.0 0.0)
+ +
pull-coord1-groups:
+
The two groups indices should be given on which this pull coordinate +will operate. The first index can be 0, in which case an absolute reference +of pull-coord1-origin is used. With an absolute reference the system +is no longer translation invariant and one should think about what to do with +the center of mass motion.
+
pull-coord1-origin: (0.0 0.0 0.0)
+
The pull reference position for use with an absolute reference.
+
pull-coord1-vec: (0.0 0.0 0.0)
The pull direction. grompp normalizes the vector.
-
pull-init1: (0.0) / (0.0 0.0 0.0) [nm]
-
The reference distance at t=0. This is a single value, -except for geometry position which uses a vector.
-
pull-rate1: (0) [nm/ps]
+
pull-coord1-init: (0.0) [nm]
+
The reference distance at t=0.
+
pull-coord1-rate: (0) [nm/ps]
The rate of change of the reference position.
-
pull-k1: (0) [kJ mol-1 nm-2] / [kJ mol-1 nm-1]
+
pull-coord1-k: (0) [kJ mol-1 nm-2] / [kJ mol-1 nm-1]
The force constant. For umbrella pulling this is the harmonic force constant in [kJ mol-1 nm-2]. For constant force pulling this is the force constant of the linear potential, and thus minus (!) the constant force in [kJ mol-1 nm-1].
-
pull-kB1: (pull-k1) [kJ mol-1 nm-2] / [kJ mol-1 nm-1]
-
As pull-k1, but for state B. This is only used when +
pull-coord1-kB: (pull-k1) [kJ mol-1 nm-2] / [kJ mol-1 nm-1]
+
As pull-coord1-k, but for state B. This is only used when free-energy is turned on. -The force constant is then (1 - lambda)*pull-k1 + lambda*pull-kB1. +The force constant is then (1 - lambda)*pull-coord1-k + lambda*pull-coord1-kB. +
diff --git a/src/gromacs/fileio/tpxio.c b/src/gromacs/fileio/tpxio.c index 67d7ba117b..55a0cd6922 100644 --- a/src/gromacs/fileio/tpxio.c +++ b/src/gromacs/fileio/tpxio.c @@ -70,7 +70,7 @@ static const char *tpx_tag = TPX_TAG_RELEASE; /* This number should be increased whenever the file format changes! */ -static const int tpx_version = 94; +static const int tpx_version = 95; /* This number should only be increased when you edit the TOPOLOGY section * or the HEADER of the tpx format. @@ -233,10 +233,14 @@ static void _do_section(t_fileio *fio, int key, gmx_bool bRead, const char *src, * Now the higer level routines that do io of the structures and arrays * **************************************************************/ -static void do_pullgrp(t_fileio *fio, t_pullgrp *pgrp, gmx_bool bRead, - int file_version) +static void do_pullgrp_tpx_pre95(t_fileio *fio, + t_pull_group *pgrp, + t_pull_coord *pcrd, + gmx_bool bRead, + int file_version) { - int i; + int i; + rvec tmp; gmx_fio_do_int(fio, pgrp->nat); if (bRead) @@ -251,18 +255,53 @@ static void do_pullgrp(t_fileio *fio, t_pullgrp *pgrp, gmx_bool bRead, } gmx_fio_ndo_real(fio, pgrp->weight, pgrp->nweight); gmx_fio_do_int(fio, pgrp->pbcatom); - gmx_fio_do_rvec(fio, pgrp->vec); - gmx_fio_do_rvec(fio, pgrp->init); - gmx_fio_do_real(fio, pgrp->rate); - gmx_fio_do_real(fio, pgrp->k); + gmx_fio_do_rvec(fio, pcrd->vec); + clear_rvec(pcrd->origin); + gmx_fio_do_rvec(fio, tmp); + pcrd->init = tmp[0]; + gmx_fio_do_real(fio, pcrd->rate); + gmx_fio_do_real(fio, pcrd->k); if (file_version >= 56) { - gmx_fio_do_real(fio, pgrp->kB); + gmx_fio_do_real(fio, pcrd->kB); } else { - pgrp->kB = pgrp->k; + pcrd->kB = pcrd->k; + } +} + +static void do_pull_group(t_fileio *fio, t_pull_group *pgrp, gmx_bool bRead) +{ + int i; + + gmx_fio_do_int(fio, pgrp->nat); + if (bRead) + { + snew(pgrp->ind, pgrp->nat); + } + gmx_fio_ndo_int(fio, pgrp->ind, pgrp->nat); + gmx_fio_do_int(fio, pgrp->nweight); + if (bRead) + { + snew(pgrp->weight, pgrp->nweight); } + gmx_fio_ndo_real(fio, pgrp->weight, pgrp->nweight); + gmx_fio_do_int(fio, pgrp->pbcatom); +} + +static void do_pull_coord(t_fileio *fio, t_pull_coord *pcrd) +{ + int i; + + gmx_fio_do_int(fio, pcrd->group[0]); + gmx_fio_do_int(fio, pcrd->group[1]); + gmx_fio_do_rvec(fio, pcrd->origin); + gmx_fio_do_rvec(fio, pcrd->vec); + gmx_fio_do_real(fio, pcrd->init); + gmx_fio_do_real(fio, pcrd->rate); + gmx_fio_do_real(fio, pcrd->k); + gmx_fio_do_real(fio, pcrd->kB); } static void do_expandedvals(t_fileio *fio, t_expanded *expand, t_lambda *fepvals, gmx_bool bRead, int file_version) @@ -561,21 +600,67 @@ static void do_pull(t_fileio *fio, t_pull *pull, gmx_bool bRead, int file_versio { int g; - gmx_fio_do_int(fio, pull->ngrp); + if (file_version >= 95) + { + gmx_fio_do_int(fio, pull->ngroup); + } + gmx_fio_do_int(fio, pull->ncoord); + if (file_version < 95) + { + pull->ngroup = pull->ncoord + 1; + } gmx_fio_do_int(fio, pull->eGeom); gmx_fio_do_ivec(fio, pull->dim); gmx_fio_do_real(fio, pull->cyl_r1); gmx_fio_do_real(fio, pull->cyl_r0); gmx_fio_do_real(fio, pull->constr_tol); + if (file_version >= 95) + { + gmx_fio_do_int(fio, pull->bPrintRef); + } gmx_fio_do_int(fio, pull->nstxout); gmx_fio_do_int(fio, pull->nstfout); if (bRead) { - snew(pull->grp, pull->ngrp+1); + snew(pull->group, pull->ngroup); + snew(pull->coord, pull->ncoord); } - for (g = 0; g < pull->ngrp+1; g++) + if (file_version < 95) { - do_pullgrp(fio, &pull->grp[g], bRead, file_version); + /* epullgPOS for position pulling, before epullgDIRPBC was removed */ + if (pull->eGeom == epullgDIRPBC) + { + gmx_fatal(FARGS, "pull-geometry=position is no longer supported"); + } + if (pull->eGeom > epullgDIRPBC) + { + pull->eGeom -= 1; + } + + for (g = 0; g < pull->ngroup; g++) + { + /* We read and ignore a pull coordinate for group 0 */ + do_pullgrp_tpx_pre95(fio, &pull->group[g], &pull->coord[max(g-1,0)], + bRead, file_version); + if (g > 0) + { + pull->coord[g-1].group[0] = 0; + pull->coord[g-1].group[1] = g; + } + } + + pull->bPrintRef = (pull->group[0].nat > 0); + } + else + { + for (g = 0; g < pull->ngroup; g++) + { + do_pull_group(fio, &pull->group[g], bRead); + } + for (g = 0; g < pull->ncoord; g++) + { + do_pull_coord(fio, &pull->coord[g]); + } } } diff --git a/src/gromacs/gmxana/gmx_wham.cpp b/src/gromacs/gmxana/gmx_wham.cpp index e124e24f35..b6acc2d17e 100644 --- a/src/gromacs/gmxana/gmx_wham.cpp +++ b/src/gromacs/gmxana/gmx_wham.cpp @@ -111,12 +111,13 @@ typedef struct * \name Using umbrella pull code since gromacs 4.x */ /*!\{*/ - int npullgrps; //!< nr of pull groups in tpr file - int pull_geometry; //!< such as distance, position + int npullcrds; //!< nr of pull coordinates in tpr file + int pull_geometry; //!< such as distance, direction ivec pull_dim; //!< pull dimension with geometry distance int pull_ndim; //!< nr of pull_dim != 0 + gmx_bool bPrintRef; //!< Coordinates of reference groups written to pullx.xvg ? real *k; //!< force constants in tpr file - rvec *init_dist; //!< reference displacements + real *init_dist; //!< reference displacements real *umbInitDist; //!< reference displacement in umbrella direction /*!\}*/ /*! @@ -1982,7 +1983,7 @@ void read_pdo_files(char **fn, int nfiles, t_UmbrellaHeader* header, void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions *opt) { t_inputrec ir; - int i, ngrp, d; + int i, ncrd; t_state state; static int first = 1; @@ -1996,26 +1997,20 @@ void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions } /* nr of pull groups */ - ngrp = ir.pull->ngrp; - if (ngrp < 1) + ncrd = ir.pull->ncoord; + if (ncrd < 1) { - gmx_fatal(FARGS, "This is not a tpr of umbrella simulation. Found only %d pull groups\n", ngrp); + gmx_fatal(FARGS, "This is not a tpr of umbrella simulation. Found only %d pull coordinates\n", ncrd); } - header->npullgrps = ir.pull->ngrp; + header->npullcrds = ir.pull->ncoord; header->pull_geometry = ir.pull->eGeom; + header->bPrintRef = ir.pull->bPrintRef; copy_ivec(ir.pull->dim, header->pull_dim); header->pull_ndim = header->pull_dim[0]+header->pull_dim[1]+header->pull_dim[2]; - if (header->pull_geometry == epullgPOS && header->pull_ndim > 1) - { - gmx_fatal(FARGS, "Found pull geometry 'position' and more than 1 pull dimension (%d).\n" - "Hence, the pull potential does not correspond to a one-dimensional umbrella potential.\n" - "If you have some special umbrella setup you may want to write your own pdo files\n" - "and feed them into g_wham. Check g_wham -h !\n", header->pull_ndim); - } - snew(header->k, ngrp); - snew(header->init_dist, ngrp); - snew(header->umbInitDist, ngrp); + snew(header->k, ncrd); + snew(header->init_dist, ncrd); + snew(header->umbInitDist, ncrd); /* only z-direction with epullgCYL? */ if (header->pull_geometry == epullgCYL) @@ -2029,35 +2024,26 @@ void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions } } - for (i = 0; i < ngrp; i++) + for (i = 0; i < ncrd; i++) { - header->k[i] = ir.pull->grp[i+1].k; + header->k[i] = ir.pull->coord[i].k; if (header->k[i] == 0.0) { - gmx_fatal(FARGS, "Pull group %d has force constant of of 0.0 in %s.\n" + gmx_fatal(FARGS, "Pull coordinate %d has force constant of of 0.0 in %s.\n" "That doesn't seem to be an Umbrella tpr.\n", i, fn); } - copy_rvec(ir.pull->grp[i+1].init, header->init_dist[i]); + header->init_dist[i] = ir.pull->coord[i].init; /* initial distance to reference */ switch (header->pull_geometry) { - case epullgPOS: - for (d = 0; d < DIM; d++) - { - if (header->pull_dim[d]) - { - header->umbInitDist[i] = header->init_dist[i][d]; - } - } - break; case epullgCYL: - /* umbrella distance stored in init_dist[i][0] for geometry cylinder (not in ...[i][ZZ]) */ + /* umbrella distance stored in init_dist[i] for geometry cylinder (not in ...[i][ZZ]) */ case epullgDIST: case epullgDIR: case epullgDIRPBC: - header->umbInitDist[i] = header->init_dist[i][0]; + header->umbInitDist[i] = header->init_dist[i]; break; default: gmx_fatal(FARGS, "Pull geometry %s not supported\n", epullg_names[header->pull_geometry]); @@ -2066,18 +2052,19 @@ void read_tpr_header(const char *fn, t_UmbrellaHeader* header, t_UmbrellaOptions if (opt->verbose || first) { - printf("File %s, %d groups, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n", - fn, header->npullgrps, epullg_names[header->pull_geometry], + printf("File %s, %d coordinates, geometry \"%s\", dimensions [%s %s %s], (%d dimensions)\n" + "\tPull group coordinates%s expected in pullx files.\n", + fn, header->npullcrds, epullg_names[header->pull_geometry], int2YN(header->pull_dim[0]), int2YN(header->pull_dim[1]), int2YN(header->pull_dim[2]), - header->pull_ndim); - for (i = 0; i < ngrp; i++) + header->pull_ndim, (header->bPrintRef ? "" : " not")); + for (i = 0; i < ncrd; i++) { - printf("\tgrp %d) k = %-5g position = %g\n", i, header->k[i], header->umbInitDist[i]); + printf("\tcrd %d) k = %-5g position = %g\n", i, header->k[i], header->umbInitDist[i]); } } if (!opt->verbose && first) { - printf("\tUse option -v to see this output for all input tpr files\n"); + printf("\tUse option -v to see this output for all input tpr files\n\n"); } first = 0; @@ -2103,7 +2090,7 @@ void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header, t_groupselection *groupsel) { double **y = 0, pos = 0., t, force, time0 = 0., dt; - int ny, nt, bins, ibin, i, g, gUsed, dstep = 1, nColPerGrp, nColRefOnce, nColRefEachGrp, nColExpect, ntot; + int ny, nt, bins, ibin, i, g, gUsed, dstep = 1, nColPerCrd, nColRefEachCrd, nColExpect, ntot, column; real min, max, minfound = 1e20, maxfound = -1e20; gmx_bool dt_ok, timeok, bHaveForce; const char *quantity; @@ -2112,44 +2099,33 @@ void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header, static gmx_bool bFirst = TRUE; /* - in force output pullf.xvg: - No reference, one column per pull group - in position output pullx.xvg (not cylinder) - ndim reference, ndim columns per pull group - in position output pullx.xvg (in geometry cylinder): - ndim*2 columns per pull group (ndim for ref, ndim for group) + * Data columns in pull output: + * - in force output pullf.xvg: + * No reference columns, one column per pull coordinate + * + * - in position output pullx.xvg + * bPrintRef == TRUE: for each pull coordinate: ndim reference columns, and ndim dx columns + * bPrintRef == FALSE: for each pull coordinate: no reference columns, but ndim dx columns */ - nColPerGrp = opt->bPullx ? header->pull_ndim : 1; + nColPerCrd = opt->bPullx ? header->pull_ndim : 1; quantity = opt->bPullx ? "position" : "force"; - if (opt->bPullx) + if (opt->bPullx && header->bPrintRef) { - if (header->pull_geometry == epullgCYL) - { - /* Geometry cylinder -> reference group before each pull group */ - nColRefEachGrp = header->pull_ndim; - nColRefOnce = 0; - } - else - { - /* Geometry NOT cylinder -> reference group only once after time column */ - nColRefEachGrp = 0; - nColRefOnce = header->pull_ndim; - } + nColRefEachCrd = header->pull_ndim; } - else /* read forces, no reference groups */ + else { - nColRefEachGrp = 0; - nColRefOnce = 0; + nColRefEachCrd = 0; } - nColExpect = 1 + nColRefOnce + header->npullgrps*(nColRefEachGrp+nColPerGrp); + nColExpect = 1 + header->npullcrds*(nColRefEachCrd+nColPerCrd); bHaveForce = opt->bPullf; /* With geometry "distance" or "distance_periodic", only force reading is supported so far. - That avoids the somewhat tedious extraction of the right columns from the pullx files - to compute the distances projection on the vector. Sorry for the laziness. */ + That avoids the somewhat tedious extraction of the reaction coordinate from the pullx files. + Sorry for the laziness, this is a To-do. */ if ( (header->pull_geometry == epullgDIR || header->pull_geometry == epullgDIRPBC) && opt->bPullx) { @@ -2173,22 +2149,21 @@ void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header, bHaveForce ? "force" : "position", epullg_names[header->pull_geometry], header->pull_ndim); printf("Expecting these columns in pull file:\n" - "\t%d reference columns for all pull groups together\n" - "\t%d reference columns for each individual pull group\n" - "\t%d data columns for each pull group\n", nColRefOnce, nColRefEachGrp, nColPerGrp); - printf("With %d pull groups, expect %d columns (including the time column)\n", header->npullgrps, nColExpect); + "\t%d reference columns for each individual pull coordinate\n" + "\t%d data columns for each pull coordinate\n", nColRefEachCrd, nColPerCrd); + printf("With %d pull groups, expect %d columns (including the time column)\n", header->npullcrds, nColExpect); bFirst = FALSE; } if (ny != nColExpect) { - gmx_fatal(FARGS, "Found %d pull groups in %s,\n but %d data columns in %s (expected %d)\n" + gmx_fatal(FARGS, "Found %d pull coodinates in %s,\n but %d data columns in %s (expected %d)\n" "\nMaybe you confused options -ix and -if ?\n", - header->npullgrps, fntpr, ny-1, fn, nColExpect-1); + header->npullcrds, fntpr, ny-1, fn, nColExpect-1); } if (opt->verbose) { - printf("Found %d times and %d %s sets %s\n", nt, (ny-1)/nColPerGrp, quantity, fn); + printf("Found %d times and %d %s sets %s\n", nt, (ny-1)/nColPerCrd, quantity, fn); } if (!bGetMinMax) @@ -2210,16 +2185,16 @@ void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header, if (groupsel) { /* Use only groups selected with option -is file */ - if (header->npullgrps != groupsel->n) + if (header->npullcrds != groupsel->n) { gmx_fatal(FARGS, "tpr file contains %d pull groups, but expected %d from group selection file\n", - header->npullgrps, groupsel->n); + header->npullcrds, groupsel->n); } window->nPull = groupsel->nUse; } else { - window->nPull = header->npullgrps; + window->nPull = header->npullcrds; } window->nBin = bins; @@ -2259,7 +2234,7 @@ void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header, /* Copying umbrella center and force const is more involved since not all pull groups from header (tpr file) may be used in window variable */ - for (g = 0, gUsed = 0; g < header->npullgrps; ++g) + for (g = 0, gUsed = 0; g < header->npullcrds; ++g) { if (groupsel && (groupsel->bUse[g] == FALSE)) { @@ -2310,7 +2285,6 @@ void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header, /*if (opt->verbose) printf(" time = %f, (tmin,tmax)=(%e,%e), dt_ok=%d timeok=%d\n", t,opt->tmin, opt->tmax, dt_ok,timeok); */ - if (timeok) { /* Note: if groupsel == NULL: @@ -2319,7 +2293,7 @@ void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header, * only groups with groupsel.bUse[g]==TRUE are stored. gUsed is not always equal g */ gUsed = -1; - for (g = 0; g < header->npullgrps; ++g) + for (g = 0; g < header->npullcrds; ++g) { /* was this group selected for application in WHAM? */ if (groupsel && (groupsel->bUse[g] == FALSE)) @@ -2331,41 +2305,30 @@ void read_pull_xf(const char *fn, const char *fntpr, t_UmbrellaHeader * header, if (bHaveForce) { - /* y has 1 time column y[0] and one column per force y[1],...,y[nGrps] */ + /* y has 1 time column y[0] and one column per force y[1],...,y[nCrds] */ force = y[g+1][i]; pos = -force/header->k[g] + header->umbInitDist[g]; } else { + /* pick the right column from: + * time ref1[ndim] coord1[ndim] ref2[ndim] coord2[ndim] ... + */ + column = 1 + nColRefEachCrd + g*(nColRefEachCrd+nColPerCrd); switch (header->pull_geometry) { case epullgDIST: - /* y has 1 time column y[0] and nColPerGrps columns per pull group; - Distance to reference: */ - /* pos=dist_ndim(y+1+nColRef+g*nColPerGrp,header->pull_ndim,i); gmx 4.0 */ - pos = dist_ndim(y + 1 + nColRefOnce + g*nColPerGrp, header->pull_ndim, i); + pos = dist_ndim(y + column, header->pull_ndim, i); break; - case epullgPOS: - /* Columns - Time ref[ndim] group1[ndim] group2[ndim] ... */ case epullgCYL: - /* Columns - Time ref1[ndim] group1[ndim] ref2[ndim] group2[ndim] ... */ - - /* * with geometry==position, we have the reference once (nColRefOnce==ndim), but - no extra reference group columns before each group (nColRefEachGrp==0) - - * with geometry==cylinder, we have no initial ref group column (nColRefOnce==0), - but ndim ref group colums before every group (nColRefEachGrp==ndim) - Distance to reference: */ - pos = y[1 + nColRefOnce + nColRefEachGrp + g*(nColRefEachGrp+nColPerGrp)][i]; + pos = y[column][i]; break; default: gmx_fatal(FARGS, "Bad error, this error should have been catched before. Ups.\n"); } } - /* printf("grp %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */ + /* printf("crd %d dpos %f poseq %f pos %f \n",g,dpos,poseq,pos); */ if (bGetMinMax) { if (pos < minfound) diff --git a/src/gromacs/gmxlib/mvdata.c b/src/gromacs/gmxlib/mvdata.c index dcbcc03491..193dd02a1f 100644 --- a/src/gromacs/gmxlib/mvdata.c +++ b/src/gromacs/gmxlib/mvdata.c @@ -486,7 +486,7 @@ static void bc_cosines(const t_commrec *cr, t_cosines *cs) } } -static void bc_pullgrp(const t_commrec *cr, t_pullgrp *pgrp) +static void bc_pull_group(const t_commrec *cr, t_pull_group *pgrp) { block_bc(cr, *pgrp); if (pgrp->nat > 0) @@ -506,11 +506,13 @@ static void bc_pull(const t_commrec *cr, t_pull *pull) int g; block_bc(cr, *pull); - snew_bc(cr, pull->grp, pull->ngrp+1); - for (g = 0; g < pull->ngrp+1; g++) + snew_bc(cr, pull->group, pull->ngroup); + for (g = 0; g < pull->ngroup; g++) { - bc_pullgrp(cr, &pull->grp[g]); + bc_pull_group(cr, &pull->group[g]); } + snew_bc(cr, pull->coord, pull->ncoord); + nblock_bc(cr, pull->ncoord, pull->coord); } static void bc_rotgrp(const t_commrec *cr, t_rotgrp *rotg) diff --git a/src/gromacs/gmxlib/names.c b/src/gromacs/gmxlib/names.c index 9b7e6b6cea..6b5f527dbf 100644 --- a/src/gromacs/gmxlib/names.c +++ b/src/gromacs/gmxlib/names.c @@ -220,7 +220,7 @@ const char *epull_names[epullNR+1] = { }; const char *epullg_names[epullgNR+1] = { - "distance", "direction", "cylinder", "position", "direction-periodic", NULL + "distance", "direction", "cylinder", "direction-periodic", NULL }; const char *erotg_names[erotgNR+1] = { diff --git a/src/gromacs/gmxlib/txtdump.c b/src/gromacs/gmxlib/txtdump.c index b9d73ed725..b1e3acff6f 100644 --- a/src/gromacs/gmxlib/txtdump.c +++ b/src/gromacs/gmxlib/txtdump.c @@ -545,19 +545,28 @@ static void pr_cosine(FILE *fp, int indent, const char *title, t_cosines *cos, #define PR(t, s) pr_real(fp, indent, t, s) #define PD(t, s) pr_double(fp, indent, t, s) -static void pr_pullgrp(FILE *fp, int indent, int g, t_pullgrp *pg) +static void pr_pull_group(FILE *fp, int indent, int g, t_pull_group *pgrp) { pr_indent(fp, indent); fprintf(fp, "pull-group %d:\n", g); indent += 2; - pr_ivec_block(fp, indent, "atom", pg->ind, pg->nat, TRUE); - pr_rvec(fp, indent, "weight", pg->weight, pg->nweight, TRUE); - PI("pbcatom", pg->pbcatom); - pr_rvec(fp, indent, "vec", pg->vec, DIM, TRUE); - pr_rvec(fp, indent, "init", pg->init, DIM, TRUE); - PR("rate", pg->rate); - PR("k", pg->k); - PR("kB", pg->kB); + pr_ivec_block(fp, indent, "atom", pgrp->ind, pgrp->nat, TRUE); + pr_rvec(fp, indent, "weight", pgrp->weight, pgrp->nweight, TRUE); + PI("pbcatom", pgrp->pbcatom); +} + +static void pr_pull_coord(FILE *fp, int indent, int c, t_pull_coord *pcrd) +{ + pr_indent(fp, indent); + fprintf(fp, "pull-coord %d:\n", c); + PI("group[0]", pcrd->group[0]); + PI("group[1]", pcrd->group[1]); + pr_rvec(fp, indent, "origin", pcrd->origin, DIM, TRUE); + pr_rvec(fp, indent, "vec", pcrd->vec, DIM, TRUE); + PR("init", pcrd->init); + PR("rate", pcrd->rate); + PR("k", pcrd->k); + PR("kB", pcrd->kB); } static void pr_simtempvals(FILE *fp, int indent, t_simtemp *simtemp, int n_lambda) @@ -678,12 +687,18 @@ static void pr_pull(FILE *fp, int indent, t_pull *pull) PR("pull-r1", pull->cyl_r1); PR("pull-r0", pull->cyl_r0); PR("pull-constr-tol", pull->constr_tol); + PS("pull-bPrintRef", EBOOL(pull->bPrintRef)); PI("pull-nstxout", pull->nstxout); PI("pull-nstfout", pull->nstfout); - PI("pull-ngrp", pull->ngrp); - for (g = 0; g < pull->ngrp+1; g++) + PI("pull-ngroup", pull->ngroup); + for (g = 0; g < pull->ngroup; g++) + { + pr_pull_group(fp, indent, g, &pull->group[g]); + } + PI("pull-ncoord", pull->ncoord); + for (g = 0; g < pull->ncoord; g++) { - pr_pullgrp(fp, indent, g, &pull->grp[g]); + pr_pull_coord(fp, indent, g, &pull->coord[g]); } } diff --git a/src/gromacs/gmxlib/typedefs.c b/src/gromacs/gmxlib/typedefs.c index 040e814c41..9fe3cb08ea 100644 --- a/src/gromacs/gmxlib/typedefs.c +++ b/src/gromacs/gmxlib/typedefs.c @@ -387,22 +387,25 @@ void done_top(t_topology *top) done_blocka(&(top->excls)); } -static void done_pullgrp(t_pullgrp *pgrp) +static void done_pull_group(t_pull_group *pgrp) { - sfree(pgrp->ind); - sfree(pgrp->ind_loc); - sfree(pgrp->weight); - sfree(pgrp->weight_loc); + if (pgrp->nat > 0) + { + sfree(pgrp->ind); + sfree(pgrp->ind_loc); + sfree(pgrp->weight); + sfree(pgrp->weight_loc); + } } static void done_pull(t_pull *pull) { int i; - for (i = 0; i < pull->ngrp+1; i++) + for (i = 0; i < pull->ngroup+1; i++) { - done_pullgrp(pull->grp); - done_pullgrp(pull->dyna); + done_pull_group(pull->group); + done_pull_group(pull->dyna); } } diff --git a/src/gromacs/gmxpreprocess/readir.c b/src/gromacs/gmxpreprocess/readir.c index 3903173308..becdc4eca5 100644 --- a/src/gromacs/gmxpreprocess/readir.c +++ b/src/gromacs/gmxpreprocess/readir.c @@ -2307,7 +2307,7 @@ void get_ir(const char *mdparin, const char *mdparout, sfree(dumstr[1]); } -static int search_QMstring(char *s, int ng, const char *gn[]) +static int search_QMstring(const char *s, int ng, const char *gn[]) { /* same as normal search_string, but this one searches QM strings */ int i; @@ -2326,8 +2326,8 @@ static int search_QMstring(char *s, int ng, const char *gn[]) } /* search_QMstring */ - -int search_string(char *s, int ng, char *gn[]) +/* We would like gn to be const as well, but C doesn't allow this */ +int search_string(const char *s, int ng, char *gn[]) { int i; @@ -2497,7 +2497,7 @@ static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames) t_grpopts *opts; gmx_groups_t *groups; t_pull *pull; - int natoms, ai, aj, i, j, d, g, imin, jmin, nc; + int natoms, ai, aj, i, j, d, g, imin, jmin; t_iatom *ia; int *nrdf2, *na_vcm, na_tot; double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0; @@ -2641,43 +2641,34 @@ static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames) * to determine the optimal nrdf assignment. */ pull = ir->pull; - if (pull->eGeom == epullgPOS) + + for (i = 0; i < pull->ncoord; i++) { - nc = 0; - for (i = 0; i < DIM; i++) + imin = 1; + + for (j = 0; j < 2; j++) { - if (pull->dim[i]) + const t_pull_group *pgrp; + + pgrp = &pull->group[pull->coord[i].group[j]]; + + if (pgrp->nat > 0) { - nc++; + /* Subtract 1/2 dof from each group */ + ai = pgrp->ind[0]; + nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin; + nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin; + if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0) + { + gmx_fatal(FARGS, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups->grps[egcTC].nm_ind[ggrpnr(groups, egcTC, ai)]]); + } } - } - } - else - { - nc = 1; - } - for (i = 0; i < pull->ngrp; i++) - { - imin = 2*nc; - if (pull->grp[0].nat > 0) - { - /* Subtract 1/2 dof from the reference group */ - ai = pull->grp[0].ind[0]; - if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1) + else { - nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5; - nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5; - imin--; + /* We need to subtract the whole DOF from group j=1 */ + imin += 1; } } - /* Subtract 1/2 dof from the pulled group */ - ai = pull->grp[1+i].ind[0]; - nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin; - nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin; - if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0) - { - gmx_fatal(FARGS, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups->grps[egcTC].nm_ind[ggrpnr(groups, egcTC, ai)]]); - } } } @@ -3174,6 +3165,8 @@ void do_index(const char* mdparin, const char *ndx, if (ir->ePull != epullNO) { make_pull_groups(ir->pull, pull_grp, grps, gnames); + + make_pull_coords(ir->pull); } if (ir->bRot) @@ -3700,7 +3693,7 @@ void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys, warninp_t wi) { char err_buf[256]; - int i, m, g, nmol, npct; + int i, m, c, nmol, npct; gmx_bool bCharge, bAcc; real gdt_max, *mgrp, mt; rvec acc; @@ -3855,7 +3848,16 @@ void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys, if (ir->ePull != epullNO) { - if (ir->pull->grp[0].nat == 0) + gmx_bool bPullAbsoluteRef; + + bPullAbsoluteRef = FALSE; + for (i = 0; i < ir->pull->ncoord; i++) + { + bPullAbsoluteRef = bPullAbsoluteRef || + ir->pull->coord[i].group[0] == 0 || + ir->pull->coord[i].group[1] == 0; + } + if (bPullAbsoluteRef) { absolute_reference(ir, sys, FALSE, AbsRef); for (m = 0; m < DIM; m++) @@ -3877,9 +3879,9 @@ void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys, if ((ir->epc != epcNO && ir->compress[i][m] != 0) || ir->deform[i][m] != 0) { - for (g = 1; g < ir->pull->ngrp; g++) + for (c = 0; c < ir->pull->ncoord; c++) { - if (ir->pull->grp[g].vec[m] != 0) + if (ir->pull->coord[c].vec[m] != 0) { gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m); } diff --git a/src/gromacs/gmxpreprocess/readir.h b/src/gromacs/gmxpreprocess/readir.h index 613d9b878d..8dcaa6ed5d 100644 --- a/src/gromacs/gmxpreprocess/readir.h +++ b/src/gromacs/gmxpreprocess/readir.h @@ -77,7 +77,7 @@ extern void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts, /* Validate inputrec data. * Fatal errors will be added to nerror. */ -extern int search_string(char *s, int ng, char *gn[]); +extern int search_string(const char *s, int ng, char *gn[]); /* Returns the index of string s in the index groups */ extern void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, @@ -120,9 +120,13 @@ extern char **read_pullparams(int *ninp_p, t_inpfile **inp, warninp_t wi); /* Reads the pull parameters, returns a list of the pull group names */ -extern void make_pull_groups(t_pull *pull, char **pgnames, - t_blocka *grps, char **gnames); -/* Process the pull parameters after reading the index groups */ +extern void make_pull_groups(t_pull *pull, + char **pgnames, + const t_blocka *grps, char **gnames); +/* Process the pull group parameters after reading the index groups */ + +extern void make_pull_coords(t_pull *pull); +/* Process the pull coordinates after reading the pull groups */ extern void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real lambda, const output_env_t oenv, gmx_bool bStart); diff --git a/src/gromacs/gmxpreprocess/readpull.c b/src/gromacs/gmxpreprocess/readpull.c index 0b2c2f985f..55f7989172 100644 --- a/src/gromacs/gmxpreprocess/readpull.c +++ b/src/gromacs/gmxpreprocess/readpull.c @@ -60,20 +60,21 @@ static char pulldim[STRLEN]; -static void string2dvec(char buf[], dvec nums) +static void string2dvec(const char buf[], dvec nums) { - if (sscanf(buf, "%lf%lf%lf", &nums[0], &nums[1], &nums[2]) != 3) + double dum; + + if (sscanf(buf, "%lf%lf%lf%lf", &nums[0], &nums[1], &nums[2], &dum) != 3) { gmx_fatal(FARGS, "Expected three numbers at input line %s", buf); } } -static void init_pullgrp(t_pullgrp *pg, char *wbuf, - gmx_bool bRef, int eGeom, char *s_vec) +static void init_pull_group(t_pull_group *pg, + const char *wbuf) { double d; int n, m; - dvec vec; pg->nweight = 0; while (sscanf(wbuf, "%lf %n", &d, &n) == 1) @@ -85,124 +86,146 @@ static void init_pullgrp(t_pullgrp *pg, char *wbuf, pg->weight[pg->nweight++] = d; wbuf += n; } - if (!bRef) +} + +static void init_pull_coord(t_pull_coord *pcrd, int eGeom, + const char *origin_buf, const char *vec_buf) +{ + int m; + dvec origin, vec; + + string2dvec(origin_buf, origin); + if (pcrd->group[0] != 0 && dnorm(origin) > 0) { - if (eGeom == epullgDIST) - { - clear_dvec(vec); - } - else - { - string2dvec(s_vec, vec); - if (eGeom == epullgDIR || eGeom == epullgCYL || - (eGeom == epullgPOS && dnorm(vec) != 0)) - { - /* Normalize the direction vector */ - dsvmul(1/dnorm(vec), vec, vec); - } - } - for (m = 0; m < DIM; m++) + gmx_fatal(FARGS,"The pull origin can only be set with an absolute reference"); + } + + if (eGeom == epullgDIST) + { + clear_dvec(vec); + } + else + { + string2dvec(vec_buf, vec); + if (eGeom == epullgDIR || eGeom == epullgCYL) { - pg->vec[m] = vec[m]; + /* Normalize the direction vector */ + dsvmul(1/dnorm(vec), vec, vec); } } + for (m = 0; m < DIM; m++) + { + pcrd->origin[m] = origin[m]; + pcrd->vec[m] = vec[m]; + } } char **read_pullparams(int *ninp_p, t_inpfile **inp_p, t_pull *pull, gmx_bool *bStart, warninp_t wi) { - int ninp, nerror = 0, i, nchar, ndim, nscan, m; - t_inpfile *inp; - const char *tmp; - char **grpbuf; - char dummy[STRLEN], buf[STRLEN], init[STRLEN]; - const char *init_def1 = "0.0", *init_def3 = "0.0 0.0 0.0"; - char wbuf[STRLEN], VecTemp[STRLEN]; - dvec vec; + int ninp, nerror = 0, i, nchar, nscan, m, idum; + t_inpfile *inp; + const char *tmp; + char **grpbuf; + char dummy[STRLEN], buf[STRLEN], groups[STRLEN], init[STRLEN]; + const char *init_def1 = "0.0", *init_def3 = "0.0 0.0 0.0"; + char wbuf[STRLEN], origin_buf[STRLEN], vec_buf[STRLEN]; - t_pullgrp *pgrp; + t_pull_group *pgrp; + t_pull_coord *pcrd; ninp = *ninp_p; inp = *inp_p; /* read pull parameters */ - CTYPE("Pull geometry: distance, direction, cylinder or position"); - EETYPE("pull_geometry", pull->eGeom, epullg_names); + CTYPE("Pull geometry: distance, direction, direction-periodic or cylinder"); + EETYPE("pull-geometry", pull->eGeom, epullg_names); CTYPE("Select components for the pull vector. default: Y Y Y"); - STYPE("pull_dim", pulldim, "Y Y Y"); + STYPE("pull-dim", pulldim, "Y Y Y"); CTYPE("Cylinder radius for dynamic reaction force groups (nm)"); - RTYPE("pull_r1", pull->cyl_r1, 1.0); + RTYPE("pull-r1", pull->cyl_r1, 1.0); CTYPE("Switch from r1 to r0 in case of dynamic reaction force"); - RTYPE("pull_r0", pull->cyl_r0, 1.5); + RTYPE("pull-r0", pull->cyl_r0, 1.5); RTYPE("pull_constr_tol", pull->constr_tol, 1E-6); - EETYPE("pull_start", *bStart, yesno_names); - ITYPE("pull_nstxout", pull->nstxout, 10); - ITYPE("pull_nstfout", pull->nstfout, 1); + EETYPE("pull-start", *bStart, yesno_names); + EETYPE("pull-print-reference", pull->bPrintRef, yesno_names); + ITYPE("pull-nstxout", pull->nstxout, 10); + ITYPE("pull-nstfout", pull->nstfout, 1); CTYPE("Number of pull groups"); - ITYPE("pull_ngroups", pull->ngrp, 1); + ITYPE("pull-ngroups", pull->ngroup, 1); + CTYPE("Number of pull coordinates"); + ITYPE("pull-ncoords", pull->ncoord, 1); if (pull->cyl_r1 > pull->cyl_r0) { - warning_error(wi, "pull_r1 > pull_r0"); + warning_error(wi, "pull-r1 > pull_r0"); } - if (pull->ngrp < 1) + if (pull->ngroup < 1) { - gmx_fatal(FARGS, "pull_ngroups should be >= 1"); + gmx_fatal(FARGS, "pull-ngroups should be >= 1"); } + /* We always add an absolute reference group (index 0), even if not used */ + pull->ngroup += 1; - snew(pull->grp, pull->ngrp+1); - - if (pull->eGeom == epullgPOS) - { - ndim = 3; - } - else + if (pull->ncoord < 1) { - ndim = 1; + gmx_fatal(FARGS, "pull-ncoords should be >= 1"); } + snew(pull->group, pull->ngroup); + + snew(pull->coord, pull->ncoord); + /* pull group options */ CTYPE("Group name, weight (default all 1), vector, init, rate (nm/ps), kJ/(mol*nm^2)"); + /* Read the pull groups */ - snew(grpbuf, pull->ngrp+1); - for (i = 0; i < pull->ngrp+1; i++) + snew(grpbuf, pull->ngroup); + /* Group 0 is the absolute reference, we don't read anything for 0 */ + for (i = 1; i < pull->ngroup; i++) { - pgrp = &pull->grp[i]; + pgrp = &pull->group[i]; snew(grpbuf[i], STRLEN); - sprintf(buf, "pull_group%d", i); + sprintf(buf, "pull-group%d-name", i); STYPE(buf, grpbuf[i], ""); - sprintf(buf, "pull_weights%d", i); + sprintf(buf, "pull-group%d-weights", i); STYPE(buf, wbuf, ""); - sprintf(buf, "pull_pbcatom%d", i); + sprintf(buf, "pull-group%d-pbcatom", i); ITYPE(buf, pgrp->pbcatom, 0); - if (i > 0) - { - sprintf(buf, "pull_vec%d", i); - STYPE(buf, VecTemp, "0.0 0.0 0.0"); - sprintf(buf, "pull_init%d", i); - STYPE(buf, init, ndim == 1 ? init_def1 : init_def3); - nscan = sscanf(init, "%lf %lf %lf", &vec[0], &vec[1], &vec[2]); - if (nscan != ndim) - { - fprintf(stderr, "ERROR: %s should have %d components\n", buf, ndim); - nerror++; - } - for (m = 0; m < DIM; m++) - { - pgrp->init[m] = (m < ndim ? vec[m] : 0.0); - } - sprintf(buf, "pull_rate%d", i); - RTYPE(buf, pgrp->rate, 0.0); - sprintf(buf, "pull_k%d", i); - RTYPE(buf, pgrp->k, 0.0); - sprintf(buf, "pull_kB%d", i); - RTYPE(buf, pgrp->kB, pgrp->k); - } /* Initialize the pull group */ - init_pullgrp(pgrp, wbuf, i == 0, pull->eGeom, VecTemp); + init_pull_group(pgrp, wbuf); + } + + /* Read the pull coordinates */ + for (i = 1; i < pull->ncoord + 1; i++) + { + pcrd = &pull->coord[i-1]; + sprintf(buf, "pull-coord%d-groups", i); + STYPE(buf, groups, ""); + nscan = sscanf(groups, "%d %d %d", &pcrd->group[0], &pcrd->group[1], &idum); + if (nscan != 2) + { + fprintf(stderr, "ERROR: %s should have %d components\n", buf, 2); + nerror++; + } + sprintf(buf, "pull-coord%d-origin", i); + STYPE(buf, origin_buf, "0.0 0.0 0.0"); + sprintf(buf, "pull-coord%d-vec", i); + STYPE(buf, vec_buf, "0.0 0.0 0.0"); + sprintf(buf, "pull-coord%d-init", i); + RTYPE(buf, pcrd->init, 0.0); + sprintf(buf, "pull-coord%d-rate", i); + RTYPE(buf, pcrd->rate, 0.0); + sprintf(buf, "pull-coord%d-k", i); + RTYPE(buf, pcrd->k, 0.0); + sprintf(buf, "pull-coord%d-kB", i); + RTYPE(buf, pcrd->kB, pcrd->k); + + /* Initialize the pull coordinate */ + init_pull_coord(pcrd, pull->eGeom, origin_buf, vec_buf); } *ninp_p = ninp; @@ -211,14 +234,81 @@ char **read_pullparams(int *ninp_p, t_inpfile **inp_p, return grpbuf; } -void make_pull_groups(t_pull *pull, char **pgnames, t_blocka *grps, char **gnames) +void make_pull_groups(t_pull *pull, + char **pgnames, + const t_blocka *grps, char **gnames) { - int d, nchar, g, ig = -1, i; - char *ptr, pulldim1[STRLEN]; - t_pullgrp *pgrp; + int g, ig = -1, i; + t_pull_group *pgrp; - ptr = pulldim; - i = 0; + /* Absolute reference group (might not be used) is special */ + pgrp = &pull->group[0]; + pgrp->nat = 0; + pgrp->pbcatom = -1; + + for (g = 1; g < pull->ngroup; g++) + { + pgrp = &pull->group[g]; + + ig = search_string(pgnames[g], grps->nr, gnames); + pgrp->nat = grps->index[ig+1] - grps->index[ig]; + + fprintf(stderr, "Pull group %d '%s' has %d atoms\n", + g, pgnames[g], pgrp->nat); + + if (pgrp->nat == 0) + { + gmx_fatal(FARGS, "Pull group %d '%s' is empty", g, pgnames[g]); + } + + snew(pgrp->ind, pgrp->nat); + for (i = 0; i < pgrp->nat; i++) + { + pgrp->ind[i] = grps->a[grps->index[ig]+i]; + } + + if (pull->eGeom == epullgCYL && g == 1 && pgrp->nweight > 0) + { + gmx_fatal(FARGS, "Weights are not supported for the reference group with cylinder pulling"); + } + if (pgrp->nweight > 0 && pgrp->nweight != pgrp->nat) + { + gmx_fatal(FARGS, "Number of weights (%d) for pull group %d '%s' does not match the number of atoms (%d)", + pgrp->nweight, g, pgnames[g], pgrp->nat); + } + + if (pgrp->nat == 1) + { + /* No pbc is required for this group */ + pgrp->pbcatom = -1; + } + else + { + if (pgrp->pbcatom > 0) + { + pgrp->pbcatom -= 1; + } + else if (pgrp->pbcatom == 0) + { + pgrp->pbcatom = pgrp->ind[(pgrp->nat-1)/2]; + } + else + { + /* Use cosine weighting */ + pgrp->pbcatom = -1; + } + } + } +} + +void make_pull_coords(t_pull *pull) +{ + int ndim, d, nchar, c; + char *ptr, pulldim1[STRLEN]; + t_pull_coord *pcrd; + + ptr = pulldim; + ndim = 0; for (d = 0; d < DIM; d++) { if (sscanf(ptr, "%s%n", pulldim1, &nchar) != 1) @@ -234,7 +324,7 @@ void make_pull_groups(t_pull *pull, char **pgnames, t_blocka *grps, char **gname else if (gmx_strncasecmp(pulldim1, "Y", 1) == 0) { pull->dim[d] = 1; - i++; + ndim++; } else { @@ -243,103 +333,47 @@ void make_pull_groups(t_pull *pull, char **pgnames, t_blocka *grps, char **gname } ptr += nchar; } - if (i == 0) + if (ndim == 0) { gmx_fatal(FARGS, "All entries in pull_dim are N"); } - for (g = 0; g < pull->ngrp+1; g++) + for (c = 0; c < pull->ncoord; c++) { - pgrp = &pull->grp[g]; - if (g == 0 && strcmp(pgnames[g], "") == 0) + pcrd = &pull->coord[c]; + + if (pcrd->group[0] < 0 || pcrd->group[0] >= pull->ngroup || + pcrd->group[1] < 0 || pcrd->group[1] >= pull->ngroup) { - pgrp->nat = 0; + gmx_fatal(FARGS, "Pull group index in pull-coord%d-groups out of range, should be between %d and %d", c+1, 0, pull->ngroup+1); } - else + + if (pcrd->group[0] == pcrd->group[1]) { - ig = search_string(pgnames[g], grps->nr, gnames); - pgrp->nat = grps->index[ig+1] - grps->index[ig]; + gmx_fatal(FARGS, "Identical pull group indices in pull-coord%d-groups", c+1); } - if (pgrp->nat > 0) - { - fprintf(stderr, "Pull group %d '%s' has %d atoms\n", - g, pgnames[g], pgrp->nat); - snew(pgrp->ind, pgrp->nat); - for (i = 0; i < pgrp->nat; i++) - { - pgrp->ind[i] = grps->a[grps->index[ig]+i]; - } - - if (pull->eGeom == epullgCYL && g == 0 && pgrp->nweight > 0) - { - gmx_fatal(FARGS, "Weights are not supported for the reference group with cylinder pulling"); - } - if (pgrp->nweight > 0 && pgrp->nweight != pgrp->nat) - { - gmx_fatal(FARGS, "Number of weights (%d) for pull group %d '%s' does not match the number of atoms (%d)", - pgrp->nweight, g, pgnames[g], pgrp->nat); - } - if (pgrp->nat == 1) - { - /* No pbc is required for this group */ - pgrp->pbcatom = -1; - } - else - { - if (pgrp->pbcatom > 0) - { - pgrp->pbcatom -= 1; - } - else if (pgrp->pbcatom == 0) - { - pgrp->pbcatom = pgrp->ind[(pgrp->nat-1)/2]; - } - else - { - /* Use cosine weighting */ - pgrp->pbcatom = -1; - } - } + if (pull->eGeom == epullgCYL && pcrd->group[0] != 1) + { + gmx_fatal(FARGS, "With pull geometry %s, the first pull group should always be 1", EPULLGEOM(pull->eGeom)); + } - if (g > 0 && pull->eGeom != epullgDIST) + if (pull->eGeom != epullgDIST) + { + for (d = 0; d < DIM; d++) { - for (d = 0; d < DIM; d++) + if (pcrd->vec[d] != 0 && pull->dim[d] == 0) { - if (pgrp->vec[d] != 0 && pull->dim[d] == 0) - { - gmx_fatal(FARGS, "ERROR: pull_vec%d has non-zero %c-component while pull_dim in N\n", g, 'x'+d); - } + gmx_fatal(FARGS, "ERROR: pull-group%d-vec has non-zero %c-component while pull_dim for the %c-dimension is N\n", c+1, 'x'+d, 'x'+d); } } - - if ((pull->eGeom == epullgDIR || pull->eGeom == epullgCYL) && - g > 0 && norm2(pgrp->vec) == 0) - { - gmx_fatal(FARGS, "pull_vec%d can not be zero with geometry %s", - g, EPULLGEOM(pull->eGeom)); - } - if ((pull->eGeom == epullgPOS) && pgrp->rate != 0 && - g > 0 && norm2(pgrp->vec) == 0) - { - gmx_fatal(FARGS, "pull_vec%d can not be zero with geometry %s and non-zero rate", - g, EPULLGEOM(pull->eGeom)); - } } - else + + if ((pull->eGeom == epullgDIR || pull->eGeom == epullgCYL) && + norm2(pcrd->vec) == 0) { - if (g == 0) - { - if (pull->eGeom == epullgCYL) - { - gmx_fatal(FARGS, "Absolute reference groups are not supported with geometry %s", EPULLGEOM(pull->eGeom)); - } - } - else - { - gmx_fatal(FARGS, "Pull group %d '%s' is empty", g, pgnames[g]); - } - pgrp->pbcatom = -1; + gmx_fatal(FARGS, "pull-group%d-vec can not be zero with geometry %s", + c+1, EPULLGEOM(pull->eGeom)); } } } @@ -347,14 +381,16 @@ void make_pull_groups(t_pull *pull, char **pgnames, t_blocka *grps, char **gname void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real lambda, const output_env_t oenv, gmx_bool bStart) { - t_mdatoms *md; - t_pull *pull; - t_pullgrp *pgrp; - t_pbc pbc; - int ndim, g, m; - double t_start, tinvrate; - rvec init; - dvec dr, dev; + t_mdatoms *md; + t_pull *pull; + t_pull_coord *pcrd; + t_pull_group *pgrp0, *pgrp1; + t_pbc pbc; + int c, m; + double t_start, tinvrate; + real init; + dvec dr; + double dev; init_pull(NULL, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0); md = init_mdatoms(NULL, mtop, ir->efep); @@ -364,14 +400,6 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l update_mdatoms(md, lambda); } pull = ir->pull; - if (pull->eGeom == epullgPOS) - { - ndim = 3; - } - else - { - ndim = 1; - } set_pbc(&pbc, ir->ePBC, box); @@ -380,56 +408,39 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l pull_calc_coms(NULL, pull, md, &pbc, t_start, x, NULL); fprintf(stderr, "Pull group natoms pbc atom distance at start reference at t=0\n"); - for (g = 0; g < pull->ngrp+1; g++) + for (c = 0; c < pull->ncoord; c++) { - pgrp = &pull->grp[g]; - fprintf(stderr, "%8d %8d %8d ", g, pgrp->nat, pgrp->pbcatom+1); - copy_rvec(pgrp->init, init); - clear_rvec(pgrp->init); - if (g > 0) + pcrd = &pull->coord[c]; + + pgrp0 = &pull->group[pcrd->group[0]]; + pgrp1 = &pull->group[pcrd->group[1]]; + fprintf(stderr, "%8d %8d %8d\n", + pcrd->group[0], pgrp0->nat, pgrp0->pbcatom+1); + fprintf(stderr, "%8d %8d %8d ", + pcrd->group[1], pgrp1->nat, pgrp1->pbcatom+1); + + init = pcrd->init; + pcrd->init = 0; + + if (pcrd->rate == 0) { - if (pgrp->rate == 0) - { - tinvrate = 0; - } - else - { - tinvrate = t_start/pgrp->rate; - } - get_pullgrp_distance(pull, &pbc, g, 0, dr, dev); - for (m = 0; m < DIM; m++) - { - if (m < ndim) - { - fprintf(stderr, " %6.3f", dev[m]); - } - else - { - fprintf(stderr, " "); - } - } - fprintf(stderr, " "); - for (m = 0; m < DIM; m++) - { - if (bStart) - { - pgrp->init[m] = init[m] + dev[m] - - tinvrate*(pull->eGeom == epullgPOS ? pgrp->vec[m] : 1); - } - else - { - pgrp->init[m] = init[m]; - } - if (m < ndim) - { - fprintf(stderr, " %6.3f", pgrp->init[m]); - } - else - { - fprintf(stderr, " "); - } - } + tinvrate = 0; + } + else + { + tinvrate = t_start/pcrd->rate; + } + get_pull_coord_distance(pull, c, &pbc, 0, dr, &dev); + fprintf(stderr, " %6.3f ", dev); + + if (bStart) + { + pcrd->init = init + dev - tinvrate; + } + else + { + pcrd->init = init; } - fprintf(stderr, "\n"); + fprintf(stderr, " %6.3f\n", pcrd->init); } } diff --git a/src/gromacs/legacyheaders/pull.h b/src/gromacs/legacyheaders/pull.h index f917e7c5b4..39b985b062 100644 --- a/src/gromacs/legacyheaders/pull.h +++ b/src/gromacs/legacyheaders/pull.h @@ -47,9 +47,11 @@ extern "C" { /* This file contains datatypes and function declarations necessary for mdrun to interface with the pull code */ -/* Get the distance to the reference and deviation for pull group g */ -void get_pullgrp_distance(t_pull *pull, t_pbc *pbc, int g, double t, - dvec dr, dvec dev); +/* Get the distance to the reference and deviation for pull coord coord_ind */ +void get_pull_coord_distance(const t_pull *pull, + int coord_ind, + const t_pbc *pbc, double t, + dvec dr, double *dev); /* Set the all the pull forces to zero */ void clear_pull_forces(t_pull *pull); diff --git a/src/gromacs/legacyheaders/types/enums.h b/src/gromacs/legacyheaders/types/enums.h index 3b3dc7e5f8..0e57e3c344 100644 --- a/src/gromacs/legacyheaders/types/enums.h +++ b/src/gromacs/legacyheaders/types/enums.h @@ -300,7 +300,7 @@ enum { }; enum { - epullgDIST, epullgDIR, epullgCYL, epullgPOS, epullgDIRPBC, epullgNR + epullgDIST, epullgDIR, epullgCYL, epullgDIRPBC, epullgNR }; #define PULL_CYL(pull) ((pull)->eGeom == epullgCYL) diff --git a/src/gromacs/legacyheaders/types/inputrec.h b/src/gromacs/legacyheaders/types/inputrec.h index f66c2c16e1..c773cbb547 100644 --- a/src/gromacs/legacyheaders/types/inputrec.h +++ b/src/gromacs/legacyheaders/types/inputrec.h @@ -111,19 +111,28 @@ typedef struct { real *weight_loc; /* Weights for the local indices */ int epgrppbc; /* The type of pbc for this pull group, see enum above */ atom_id pbcatom; /* The reference atom for pbc (global number) */ - rvec vec; /* The pull vector, direction or position */ - rvec init; /* Initial reference displacement */ - real rate; /* Rate of motion (nm/ps) */ - real k; /* force constant */ - real kB; /* force constant for state B */ + + /* Variables not present in mdp, but used at run time */ real wscale; /* scaling factor for the weights: sum w m/sum w w m */ real invtm; /* inverse total mass of the group: 1/wscale sum w m */ dvec x; /* center of mass before update */ dvec xp; /* center of mass after update before constraining */ +} t_pull_group; + +typedef struct { + int group[2]; /* The pull groups, index in group in t_pull */ + rvec origin; /* The origin for the absolute reference */ + rvec vec; /* The pull vector, direction or position */ + real init; /* Initial reference displacement */ + real rate; /* Rate of motion (nm/ps) */ + real k; /* force constant */ + real kB; /* force constant for state B */ + + /* Variables not present in mdp, but used at run time */ dvec dr; /* The distance from the reference group */ double f_scal; /* Scalar force for directional pulling */ dvec f; /* force due to the pulling/constraining */ -} t_pullgrp; +} t_pull_coord; typedef struct { int eSimTempScale; /* simulated temperature scaling; linear or exponential */ @@ -195,12 +204,14 @@ typedef struct { } t_expanded; typedef struct { - int ngrp; /* number of groups */ + int ngroup; /* number of pull groups */ + int ncoord; /* number of pull coordinates */ int eGeom; /* pull geometry */ ivec dim; /* used to select components for constraint */ real cyl_r1; /* radius of cylinder for dynamic COM */ real cyl_r0; /* radius of cylinder including switch length */ real constr_tol; /* absolute tolerance for constraints in (nm) */ + gmx_bool bPrintRef; /* Print coordinates of the first group in each coord */ int nstxout; /* Output frequency for pull x */ int nstfout; /* Output frequency for pull f */ int ePBC; /* the boundary conditions */ @@ -208,8 +219,11 @@ typedef struct { gmx_bool bRefAt; /* do we need reference atoms for a group COM ? */ int cosdim; /* dimension for cosine weighting, -1 if none */ gmx_bool bVirial; /* do we need to add the pull virial? */ - t_pullgrp *grp; /* groups to pull/restrain/etc/ */ - t_pullgrp *dyna; /* dynamic groups for use with local constraints */ + t_pull_group *group; /* groups to pull/restrain/etc/ */ + t_pull_coord *coord; /* the pull coordinates */ + + /* Variables not present in mdp, but used at run time */ + t_pull_group *dyna; /* dynamic groups for use with local constraints */ rvec *rbuf; /* COM calculation buffer */ dvec *dbuf; /* COM calculation buffer */ double *dbuf_cyl; /* cylinder ref. groups COM calculation buffer */ diff --git a/src/gromacs/mdlib/pull.c b/src/gromacs/mdlib/pull.c index d5e97bccdb..633c093925 100644 --- a/src/gromacs/mdlib/pull.c +++ b/src/gromacs/mdlib/pull.c @@ -62,7 +62,7 @@ #include "copyrite.h" #include "macros.h" -static void pull_print_x_grp(FILE *out, gmx_bool bRef, ivec dim, t_pullgrp *pgrp) +static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp) { int m; @@ -70,60 +70,60 @@ static void pull_print_x_grp(FILE *out, gmx_bool bRef, ivec dim, t_pullgrp *pgrp { if (dim[m]) { - fprintf(out, "\t%g", bRef ? pgrp->x[m] : pgrp->dr[m]); + fprintf(out, "\t%g", pgrp->x[m]); } } } -static void pull_print_x(FILE *out, t_pull *pull, double t) +static void pull_print_coord_dr(FILE *out, ivec dim, const t_pull_coord *pcrd) { - int g; - - fprintf(out, "%.4f", t); + int m; - if (PULL_CYL(pull)) + for (m = 0; m < DIM; m++) { - for (g = 1; g < 1+pull->ngrp; g++) + if (dim[m]) { - pull_print_x_grp(out, TRUE, pull->dim, &pull->dyna[g]); - pull_print_x_grp(out, FALSE, pull->dim, &pull->grp[g]); + fprintf(out, "\t%g", pcrd->dr[m]); } } - else +} + +static void pull_print_x(FILE *out, t_pull *pull, double t) +{ + int c; + const t_pull_coord *pcrd; + + fprintf(out, "%.4f", t); + + for (c = 0; c < pull->ncoord; c++) { - for (g = 0; g < 1+pull->ngrp; g++) + pcrd = &pull->coord[c]; + + if (pull->bPrintRef) { - if (pull->grp[g].nat > 0) + if (PULL_CYL(pull)) { - pull_print_x_grp(out, g == 0, pull->dim, &pull->grp[g]); + pull_print_group_x(out, pull->dim, &pull->dyna[c]); + } + else + { + pull_print_group_x(out, pull->dim, &pull->group[pcrd->group[0]]); } } + pull_print_coord_dr(out, pull->dim, pcrd); } fprintf(out, "\n"); } static void pull_print_f(FILE *out, t_pull *pull, double t) { - int g, d; + int c, d; fprintf(out, "%.4f", t); - for (g = 1; g < 1+pull->ngrp; g++) + for (c = 0; c < pull->ncoord; c++) { - if (pull->eGeom == epullgPOS) - { - for (d = 0; d < DIM; d++) - { - if (pull->dim[d]) - { - fprintf(out, "\t%g", pull->grp[g].f[d]); - } - } - } - else - { - fprintf(out, "\t%g", pull->grp[g].f_scal); - } + fprintf(out, "\t%g", pull->coord[c].f_scal); } fprintf(out, "\n"); } @@ -145,7 +145,7 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv gmx_bool bCoord, unsigned long Flags) { FILE *fp; - int nsets, g, m; + int nsets, c, m; char **setname, buf[10]; if (Flags & MD_APPENDFILES) @@ -166,53 +166,48 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv exvggtXNY, oenv); } - snew(setname, (1+pull->ngrp)*DIM); + snew(setname, 2*pull->ncoord*DIM); nsets = 0; - for (g = 0; g < 1+pull->ngrp; g++) + for (c = 0; c < pull->ncoord; c++) { - if (pull->grp[g].nat > 0 && - (g > 0 || (bCoord && !PULL_CYL(pull)))) + if (bCoord) { - if (bCoord || pull->eGeom == epullgPOS) + if (pull->bPrintRef) { - if (PULL_CYL(pull)) - { - for (m = 0; m < DIM; m++) - { - if (pull->dim[m]) - { - sprintf(buf, "%d %s%c", g, "c", 'X'+m); - setname[nsets] = strdup(buf); - nsets++; - } - } - } for (m = 0; m < DIM; m++) { if (pull->dim[m]) { - sprintf(buf, "%d %s%c", - g, (bCoord && g > 0) ? "d" : "", 'X'+m); + sprintf(buf, "%d %s%c", c+1, "c", 'X'+m); setname[nsets] = strdup(buf); nsets++; } } } - else + for (m = 0; m < DIM; m++) { - sprintf(buf, "%d", g); - setname[nsets] = strdup(buf); - nsets++; + if (pull->dim[m]) + { + sprintf(buf, "%d %s%c", c+1, "d", 'X'+m); + setname[nsets] = strdup(buf); + nsets++; + } } } + else + { + sprintf(buf, "%d", c+1); + setname[nsets] = strdup(buf); + nsets++; + } } - if (bCoord || nsets > 1) + if (nsets > 1) { xvgr_legend(fp, nsets, (const char**)setname, oenv); } - for (g = 0; g < nsets; g++) + for (c = 0; c < nsets; c++) { - sfree(setname[g]); + sfree(setname[c]); } sfree(setname); } @@ -221,8 +216,8 @@ static FILE *open_pull_out(const char *fn, t_pull *pull, const output_env_t oenv } /* Apply forces in a mass weighted fashion */ -static void apply_forces_grp(t_pullgrp *pgrp, t_mdatoms * md, - dvec f_pull, int sign, rvec *f) +static void apply_forces_grp(const t_pull_group *pgrp, const t_mdatoms *md, + const dvec f_pull, int sign, rvec *f) { int i, ii, m, start, end; double wmass, inv_wm; @@ -251,24 +246,25 @@ static void apply_forces_grp(t_pullgrp *pgrp, t_mdatoms * md, /* Apply forces in a mass weighted fashion */ static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f) { - int i; - t_pullgrp *pgrp; + int c; + const t_pull_coord *pcrd; - for (i = 1; i < pull->ngrp+1; i++) + for (c = 0; c < pull->ncoord; c++) { - pgrp = &(pull->grp[i]); - apply_forces_grp(pgrp, md, pgrp->f, 1, f); - if (pull->grp[0].nat) + pcrd = &pull->coord[c]; + + if (PULL_CYL(pull)) { - if (PULL_CYL(pull)) - { - apply_forces_grp(&(pull->dyna[i]), md, pgrp->f, -1, f); - } - else + apply_forces_grp(&pull->dyna[c], md, pcrd->f, -1, f); + } + else + { + if (pull->group[pcrd->group[0]].nat > 0) { - apply_forces_grp(&(pull->grp[0]), md, pgrp->f, -1, f); + apply_forces_grp(&pull->group[pcrd->group[0]], md, pcrd->f, -1, f); } } + apply_forces_grp(&pull->group[pcrd->group[1]], md, pcrd->f, 1, f); } } @@ -293,16 +289,28 @@ static double max_pull_distance2(const t_pull *pull, const t_pbc *pbc) return 0.25*max_d2; } -static void get_pullgrps_dr(const t_pull *pull, const t_pbc *pbc, int g, double t, - dvec xg, dvec xref, double max_dist2, - dvec dr) +static void low_get_pull_coord_dr(const t_pull *pull, + const t_pull_coord *pcrd, + const t_pbc *pbc, double t, + dvec xg, dvec xref, double max_dist2, + dvec dr) { - t_pullgrp *pref, *pgrp; - int m; - dvec xrefr, dref = {0, 0, 0}; - double dr2; + const t_pull_group *pgrp0, *pgrp1; + int m; + dvec xrefr, dref = {0, 0, 0}; + double dr2; + + pgrp0 = &pull->group[pcrd->group[0]]; + pgrp1 = &pull->group[pcrd->group[1]]; - pgrp = &pull->grp[g]; + /* Only the first group can be an absolute reference, in that case nat=0 */ + if (pgrp0->nat == 0) + { + for (m = 0; m < DIM; m++) + { + xref[m] = pcrd->origin[m]; + } + } copy_dvec(xref, xrefr); @@ -310,7 +318,7 @@ static void get_pullgrps_dr(const t_pull *pull, const t_pbc *pbc, int g, double { for (m = 0; m < DIM; m++) { - dref[m] = (pgrp->init[0] + pgrp->rate*t)*pull->grp[g].vec[m]; + dref[m] = (pcrd->init + pcrd->rate*t)*pcrd->vec[m]; } /* Add the reference position, so we use the correct periodic image */ dvec_inc(xrefr, dref); @@ -325,7 +333,8 @@ static void get_pullgrps_dr(const t_pull *pull, const t_pbc *pbc, int g, double } if (max_dist2 >= 0 && dr2 > 0.98*0.98*max_dist2) { - gmx_fatal(FARGS, "Distance of pull group %d (%f nm) is larger than 0.49 times the box size (%f)", g, sqrt(dr2), sqrt(max_dist2)); + gmx_fatal(FARGS, "Distance between pull groups %d and %d (%f nm) is larger than 0.49 times the box size (%f)", + pcrd->group[0], pcrd->group[1], sqrt(dr2), sqrt(max_dist2)); } if (pull->eGeom == epullgDIRPBC) @@ -334,10 +343,13 @@ static void get_pullgrps_dr(const t_pull *pull, const t_pbc *pbc, int g, double } } -static void get_pullgrp_dr(const t_pull *pull, const t_pbc *pbc, int g, double t, - dvec dr) +static void get_pull_coord_dr(const t_pull *pull, + int coord_ind, + const t_pbc *pbc, double t, + dvec dr) { - double md2; + double md2; + const t_pull_coord *pcrd; if (pull->eGeom == epullgDIRPBC) { @@ -348,46 +360,39 @@ static void get_pullgrp_dr(const t_pull *pull, const t_pbc *pbc, int g, double t md2 = max_pull_distance2(pull, pbc); } - get_pullgrps_dr(pull, pbc, g, t, - pull->grp[g].x, - PULL_CYL(pull) ? pull->dyna[g].x : pull->grp[0].x, - md2, - dr); + pcrd = &pull->coord[coord_ind]; + + low_get_pull_coord_dr(pull, pcrd, pbc, t, + pull->group[pcrd->group[1]].x, + PULL_CYL(pull) ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x, + md2, + dr); } -void get_pullgrp_distance(t_pull *pull, t_pbc *pbc, int g, double t, - dvec dr, dvec dev) +void get_pull_coord_distance(const t_pull *pull, + int coord_ind, + const t_pbc *pbc, double t, + dvec dr, double *dev) { static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety, but is fairly benign */ - t_pullgrp *pgrp; - int m; - dvec ref; - double drs, inpr; + const t_pull_coord *pcrd; + int m; + double ref, drs, inpr; - pgrp = &pull->grp[g]; + pcrd = &pull->coord[coord_ind]; - get_pullgrp_dr(pull, pbc, g, t, dr); + get_pull_coord_dr(pull, coord_ind, pbc, t, dr); - if (pull->eGeom == epullgPOS) - { - for (m = 0; m < DIM; m++) - { - ref[m] = pgrp->init[m] + pgrp->rate*t*pgrp->vec[m]; - } - } - else - { - ref[0] = pgrp->init[0] + pgrp->rate*t; - } + ref = pcrd->init + pcrd->rate*t; switch (pull->eGeom) { case epullgDIST: /* Pull along the vector between the com's */ - if (ref[0] < 0 && !bWarned) + if (ref < 0 && !bWarned) { - fprintf(stderr, "\nPull reference distance for group %d is negative (%f)\n", g, ref[0]); + fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, ref); bWarned = TRUE; } drs = dnorm(dr); @@ -396,12 +401,12 @@ void get_pullgrp_distance(t_pull *pull, t_pbc *pbc, int g, double t, /* With no vector we can not determine the direction for the force, * so we set the force to zero. */ - dev[0] = 0; + *dev = 0; } else { /* Determine the deviation */ - dev[0] = drs - ref[0]; + *dev = drs - ref; } break; case epullgDIR: @@ -411,16 +416,9 @@ void get_pullgrp_distance(t_pull *pull, t_pbc *pbc, int g, double t, inpr = 0; for (m = 0; m < DIM; m++) { - inpr += pgrp->vec[m]*dr[m]; - } - dev[0] = inpr - ref[0]; - break; - case epullgPOS: - /* Determine the difference of dr and ref along each dimension */ - for (m = 0; m < DIM; m++) - { - dev[m] = (dr[m] - ref[m])*pull->dim[m]; + inpr += pcrd->vec[m]*dr[m]; } + *dev = inpr - ref; break; } } @@ -433,10 +431,10 @@ void clear_pull_forces(t_pull *pull) * It can happen that multiple constraint steps need to be applied * and therefore the constraint forces need to be accumulated. */ - for (i = 0; i < 1+pull->ngrp; i++) + for (i = 0; i < pull->ncoord; i++) { - clear_dvec(pull->grp[i].f); - pull->grp[i].f_scal = 0; + clear_dvec(pull->coord[i].f); + pull->coord[i].f_scal = 0; } } @@ -449,63 +447,48 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, dvec *r_ij; /* x[i] com of i in prev. step. Obeys constr. -> r_ij[i] */ dvec unc_ij; /* xp[i] com of i this step, before constr. -> unc_ij */ - - dvec *rinew; /* current 'new' position of group i */ - dvec *rjnew; /* current 'new' position of group j */ - dvec ref, vec; + dvec *rnew; /* current 'new' positions of the groups */ + double *dr_tot; /* the total update of the coords */ + double ref; + dvec vec; double d0, inpr; double lambda, rm, mass, invdt = 0; gmx_bool bConverged_all, bConverged = FALSE; - int niter = 0, g, ii, j, m, max_iter = 100; - double q, a, b, c; /* for solving the quadratic equation, - see Num. Recipes in C ed 2 p. 184 */ - dvec *dr; /* correction for group i */ - dvec ref_dr; /* correction for group j */ + int niter = 0, g, c, ii, j, m, max_iter = 100; + double a; dvec f; /* the pull force */ dvec tmp, tmp3; - t_pullgrp *pdyna, *pgrp, *pref; + t_pull_group *pdyna, *pgrp0, *pgrp1; + t_pull_coord *pcrd; - snew(r_ij, pull->ngrp+1); - if (PULL_CYL(pull)) - { - snew(rjnew, pull->ngrp+1); - } - else - { - snew(rjnew, 1); - } - snew(dr, pull->ngrp+1); - snew(rinew, pull->ngrp+1); + snew(r_ij, pull->ncoord); + snew(dr_tot, pull->ncoord); + + snew(rnew, pull->ngroup); /* copy the current unconstrained positions for use in iterations. We iterate until rinew[i] and rjnew[j] obey the constraints. Then rinew - pull.x_unc[i] is the correction dr to group i */ - for (g = 1; g < 1+pull->ngrp; g++) + for (g = 0; g < pull->ngroup; g++) { - copy_dvec(pull->grp[g].xp, rinew[g]); + copy_dvec(pull->group[g].xp, rnew[g]); } if (PULL_CYL(pull)) { - for (g = 1; g < 1+pull->ngrp; g++) - { - copy_dvec(pull->dyna[g].xp, rjnew[g]); - } - } - else - { - copy_dvec(pull->grp[0].xp, rjnew[0]); + /* There is only one pull coordinate and reference group */ + copy_dvec(pull->dyna[0].xp, rnew[pull->coord[0].group[0]]); } /* Determine the constraint directions from the old positions */ - for (g = 1; g < 1+pull->ngrp; g++) + for (c = 0; c < pull->ncoord; c++) { - get_pullgrp_dr(pull, pbc, g, t, r_ij[g]); + get_pull_coord_dr(pull, c, pbc, t, r_ij[c]); /* Store the difference vector at time t for printing */ - copy_dvec(r_ij[g], pull->grp[g].dr); + copy_dvec(r_ij[c], pull->coord[c].dr); if (debug) { - fprintf(debug, "Pull group %d dr %f %f %f\n", - g, r_ij[g][XX], r_ij[g][YY], r_ij[g][ZZ]); + fprintf(debug, "Pull coord %d dr %f %f %f\n", + c, r_ij[c][XX], r_ij[c][YY], r_ij[c][ZZ]); } if (pull->eGeom == epullgDIR || pull->eGeom == epullgDIRPBC) @@ -514,11 +497,11 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, a = 0; for (m = 0; m < DIM; m++) { - a += pull->grp[g].vec[m]*r_ij[g][m]; + a += pull->coord[c].vec[m]*r_ij[c][m]; } for (m = 0; m < DIM; m++) { - r_ij[g][m] = a*pull->grp[g].vec[m]; + r_ij[c][m] = a*pull->coord[c].vec[m]; } } } @@ -529,77 +512,67 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, bConverged_all = TRUE; /* loop over all constraints */ - for (g = 1; g < 1+pull->ngrp; g++) + for (c = 0; c < pull->ncoord; c++) { - pgrp = &pull->grp[g]; - if (PULL_CYL(pull)) - { - pref = &pull->dyna[g]; - } - else - { - pref = &pull->grp[0]; - } + dvec dr0, dr1; + + pcrd = &pull->coord[c]; + pgrp0 = &pull->group[pcrd->group[0]]; + pgrp1 = &pull->group[pcrd->group[1]]; /* Get the current difference vector */ - get_pullgrps_dr(pull, pbc, g, t, rinew[g], rjnew[PULL_CYL(pull) ? g : 0], - -1, unc_ij); + low_get_pull_coord_dr(pull, pcrd, pbc, t, + rnew[pcrd->group[1]], + rnew[pcrd->group[0]], + -1, unc_ij); - if (pull->eGeom == epullgPOS) - { - for (m = 0; m < DIM; m++) - { - ref[m] = pgrp->init[m] + pgrp->rate*t*pgrp->vec[m]; - } - } - else - { - ref[0] = pgrp->init[0] + pgrp->rate*t; - /* Keep the compiler happy */ - ref[1] = 0; - ref[2] = 0; - } + ref = pcrd->init + pcrd->rate*t; if (debug) { - fprintf(debug, "Pull group %d, iteration %d\n", g, niter); + fprintf(debug, "Pull coord %d, iteration %d\n", c, niter); } - rm = 1.0/(pull->grp[g].invtm + pref->invtm); + rm = 1.0/(pgrp0->invtm + pgrp1->invtm); switch (pull->eGeom) { case epullgDIST: - if (ref[0] <= 0) + if (ref <= 0) { - gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", g, ref[0]); + gmx_fatal(FARGS, "The pull constraint reference distance for group %d is <= 0 (%f)", c, ref); } - a = diprod(r_ij[g], r_ij[g]); - b = diprod(unc_ij, r_ij[g])*2; - c = diprod(unc_ij, unc_ij) - dsqr(ref[0]); - - if (b < 0) - { - q = -0.5*(b - sqrt(b*b - 4*a*c)); - lambda = -q/a; - } - else { - q = -0.5*(b + sqrt(b*b - 4*a*c)); - lambda = -c/q; - } + double q, c_a, c_b, c_c; - if (debug) - { - fprintf(debug, - "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", - a, b, c, lambda); + c_a = diprod(r_ij[c], r_ij[c]); + c_b = diprod(unc_ij, r_ij[c])*2; + c_c = diprod(unc_ij, unc_ij) - dsqr(ref); + + if (c_b < 0) + { + q = -0.5*(c_b - sqrt(c_b*c_b - 4*c_a*c_c)); + lambda = -q/c_a; + } + else + { + q = -0.5*(c_b + sqrt(c_b*c_b - 4*c_a*c_c)); + lambda = -c_c/q; + } + + if (debug) + { + fprintf(debug, + "Pull ax^2+bx+c=0: a=%e b=%e c=%e lambda=%e\n", + c_a, c_b, c_c, lambda); + } } /* The position corrections dr due to the constraints */ - dsvmul(-lambda*rm*pgrp->invtm, r_ij[g], dr[g]); - dsvmul( lambda*rm*pref->invtm, r_ij[g], ref_dr); + dsvmul(-lambda*rm*pgrp1->invtm, r_ij[c], dr1); + dsvmul( lambda*rm*pgrp0->invtm, r_ij[c], dr0); + dr_tot[c] += -lambda*dnorm(r_ij[c]); break; case epullgDIR: case epullgDIRPBC: @@ -608,112 +581,78 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, a = 0; for (m = 0; m < DIM; m++) { - vec[m] = pgrp->vec[m]; + vec[m] = pcrd->vec[m]; a += unc_ij[m]*vec[m]; } /* Select only the component along the vector */ dsvmul(a, vec, unc_ij); - lambda = a - ref[0]; + lambda = a - ref; if (debug) { fprintf(debug, "Pull inpr %e lambda: %e\n", a, lambda); } /* The position corrections dr due to the constraints */ - dsvmul(-lambda*rm*pull->grp[g].invtm, vec, dr[g]); - dsvmul( lambda*rm* pref->invtm, vec, ref_dr); - break; - case epullgPOS: - for (m = 0; m < DIM; m++) - { - if (pull->dim[m]) - { - lambda = r_ij[g][m] - ref[m]; - /* The position corrections dr due to the constraints */ - dr[g][m] = -lambda*rm*pull->grp[g].invtm; - ref_dr[m] = lambda*rm*pref->invtm; - } - else - { - dr[g][m] = 0; - ref_dr[m] = 0; - } - } + dsvmul(-lambda*rm*pgrp1->invtm, vec, dr1); + dsvmul( lambda*rm*pgrp0->invtm, vec, dr0); + dr_tot[c] += -lambda; break; } /* DEBUG */ if (debug) { - j = (PULL_CYL(pull) ? g : 0); - get_pullgrps_dr(pull, pbc, g, t, rinew[g], rjnew[j], -1, tmp); - get_pullgrps_dr(pull, pbc, g, t, dr[g], ref_dr, -1, tmp3); + int g0, g1; + + g0 = pcrd->group[0]; + g1 = pcrd->group[1]; + low_get_pull_coord_dr(pull, pcrd, pbc, t, rnew[g1], rnew[g0], -1, tmp); + low_get_pull_coord_dr(pull, pcrd, pbc, t, dr1, dr0, -1, tmp3); fprintf(debug, "Pull cur %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", - rinew[g][0], rinew[g][1], rinew[g][2], - rjnew[j][0], rjnew[j][1], rjnew[j][2], dnorm(tmp)); - if (pull->eGeom == epullgPOS) - { - fprintf(debug, - "Pull ref %8.5f %8.5f %8.5f\n", - pgrp->vec[0], pgrp->vec[1], pgrp->vec[2]); - } - else - { - fprintf(debug, - "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f %8.5f %8.5f\n", - "", "", "", "", "", "", ref[0], ref[1], ref[2]); - } + rnew[g0][0], rnew[g0][1], rnew[g0][2], + rnew[g1][0], rnew[g1][1], rnew[g1][2], dnorm(tmp)); + fprintf(debug, + "Pull ref %8s %8s %8s %8s %8s %8s d: %8.5f\n", + "", "", "", "", "", "", ref); fprintf(debug, "Pull cor %8.5f %8.5f %8.5f j:%8.5f %8.5f %8.5f d: %8.5f\n", - dr[g][0], dr[g][1], dr[g][2], - ref_dr[0], ref_dr[1], ref_dr[2], + dr0[0], dr0[1], dr0[2], + dr1[0], dr1[1], dr1[2], dnorm(tmp3)); - fprintf(debug, - "Pull cor %10.7f %10.7f %10.7f\n", - dr[g][0], dr[g][1], dr[g][2]); } /* END DEBUG */ /* Update the COMs with dr */ - dvec_inc(rinew[g], dr[g]); - dvec_inc(rjnew[PULL_CYL(pull) ? g : 0], ref_dr); + dvec_inc(rnew[pcrd->group[1]], dr1); + dvec_inc(rnew[pcrd->group[0]], dr0); } /* Check if all constraints are fullfilled now */ - for (g = 1; g < 1+pull->ngrp; g++) + for (c = 0; c < pull->ncoord; c++) { - pgrp = &pull->grp[g]; + pcrd = &pull->coord[c]; - get_pullgrps_dr(pull, pbc, g, t, rinew[g], rjnew[PULL_CYL(pull) ? g : 0], - -1, unc_ij); + low_get_pull_coord_dr(pull, pcrd, pbc, t, + rnew[pcrd->group[1]], + rnew[pcrd->group[0]], + -1, unc_ij); switch (pull->eGeom) { case epullgDIST: - bConverged = fabs(dnorm(unc_ij) - ref[0]) < pull->constr_tol; + bConverged = fabs(dnorm(unc_ij) - ref) < pull->constr_tol; break; case epullgDIR: case epullgDIRPBC: case epullgCYL: for (m = 0; m < DIM; m++) { - vec[m] = pgrp->vec[m]; + vec[m] = pcrd->vec[m]; } inpr = diprod(unc_ij, vec); dsvmul(inpr, vec, unc_ij); bConverged = - fabs(diprod(unc_ij, vec) - ref[0]) < pull->constr_tol; - break; - case epullgPOS: - bConverged = TRUE; - for (m = 0; m < DIM; m++) - { - if (pull->dim[m] && - fabs(unc_ij[m] - ref[m]) >= pull->constr_tol) - { - bConverged = FALSE; - } - } + fabs(diprod(unc_ij, vec) - ref) < pull->constr_tol; break; } @@ -722,8 +661,8 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, if (debug) { fprintf(debug, "NOT CONVERGED YET: Group %d:" - "d_ref = %f %f %f, current d = %f\n", - g, ref[0], ref[1], ref[2], dnorm(unc_ij)); + "d_ref = %f, current d = %f\n", + g, ref, dnorm(unc_ij)); } bConverged_all = FALSE; @@ -746,59 +685,37 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, invdt = 1/dt; } - /* update the normal groups */ - for (g = 1; g < 1+pull->ngrp; g++) + /* update atoms in the groups */ + for (g = 0; g < pull->ngroup; g++) { - pgrp = &pull->grp[g]; - /* get the final dr and constraint force for group i */ - dvec_sub(rinew[g], pgrp->xp, dr[g]); - /* select components of dr */ - for (m = 0; m < DIM; m++) + const t_pull_group *pgrp; + dvec dr; + + if (PULL_CYL(pull) && g == pull->coord[0].group[0]) { - dr[g][m] *= pull->dim[m]; + pgrp = &pull->dyna[0]; } - dsvmul(1.0/(pgrp->invtm*dt*dt), dr[g], f); - dvec_inc(pgrp->f, f); - switch (pull->eGeom) + else { - case epullgDIST: - for (m = 0; m < DIM; m++) - { - pgrp->f_scal += r_ij[g][m]*f[m]/dnorm(r_ij[g]); - } - break; - case epullgDIR: - case epullgDIRPBC: - case epullgCYL: - for (m = 0; m < DIM; m++) - { - pgrp->f_scal += pgrp->vec[m]*f[m]; - } - break; - case epullgPOS: - break; + pgrp = &pull->group[g]; } - if (vir && bMaster) + /* get the final constraint displacement dr for group g */ + dvec_sub(rnew[g], pgrp->xp, dr); + /* select components of dr */ + for (m = 0; m < DIM; m++) { - /* Add the pull contribution to the virial */ - for (j = 0; j < DIM; j++) - { - for (m = 0; m < DIM; m++) - { - vir[j][m] -= 0.5*f[j]*r_ij[g][m]; - } - } + dr[m] *= pull->dim[m]; } /* update the atom positions */ - copy_dvec(dr[g], tmp); + copy_dvec(dr, tmp); for (j = 0; j < pgrp->nat_loc; j++) { ii = pgrp->ind_loc[j]; if (pgrp->weight_loc) { - dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr[g], tmp); + dsvmul(pgrp->wscale*pgrp->weight_loc[j], dr, tmp); } for (m = 0; m < DIM; m++) { @@ -814,67 +731,24 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, } } - /* update the reference groups */ - if (PULL_CYL(pull)) + /* calculate the constraint forces, used for output and virial only */ + for (c = 0; c < pull->ncoord; c++) { - /* update the dynamic reference groups */ - for (g = 1; g < 1+pull->ngrp; g++) - { - pdyna = &pull->dyna[g]; - dvec_sub(rjnew[g], pdyna->xp, ref_dr); - /* select components of ref_dr */ - for (m = 0; m < DIM; m++) - { - ref_dr[m] *= pull->dim[m]; - } + pcrd = &pull->coord[c]; + pcrd->f_scal = dr_tot[c]/((pull->group[pcrd->group[0]].invtm + pull->group[pcrd->group[1]].invtm)*dt*dt); - for (j = 0; j < pdyna->nat_loc; j++) - { - /* reset the atoms with dr, weighted by w_i */ - dsvmul(pdyna->wscale*pdyna->weight_loc[j], ref_dr, tmp); - ii = pdyna->ind_loc[j]; - for (m = 0; m < DIM; m++) - { - x[ii][m] += tmp[m]; - } - if (v) - { - for (m = 0; m < DIM; m++) - { - v[ii][m] += invdt*tmp[m]; - } - } - } - } - } - else - { - pgrp = &pull->grp[0]; - /* update the reference group */ - dvec_sub(rjnew[0], pgrp->xp, ref_dr); - /* select components of ref_dr */ - for (m = 0; m < DIM; m++) + if (vir && bMaster) { - ref_dr[m] *= pull->dim[m]; - } + double f_invr; - copy_dvec(ref_dr, tmp); - for (j = 0; j < pgrp->nat_loc; j++) - { - ii = pgrp->ind_loc[j]; - if (pgrp->weight_loc) - { - dsvmul(pgrp->wscale*pgrp->weight_loc[j], ref_dr, tmp); - } - for (m = 0; m < DIM; m++) - { - x[ii][m] += tmp[m]; - } - if (v) + /* Add the pull contribution to the virial */ + f_invr = pcrd->f_scal/dnorm(r_ij[c]); + + for (j = 0; j < DIM; j++) { for (m = 0; m < DIM; m++) { - v[ii][m] += invdt*tmp[m]; + vir[j][m] -= 0.5*f_invr*r_ij[c][j]*r_ij[c][m]; } } } @@ -882,9 +756,8 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, /* finished! I hope. Give back some memory */ sfree(r_ij); - sfree(rinew); - sfree(rjnew); - sfree(dr); + sfree(dr_tot); + sfree(rnew); } /* Pulling with a harmonic umbrella potential or constant force */ @@ -892,43 +765,43 @@ static void do_pull_pot(int ePull, t_pull *pull, t_pbc *pbc, double t, real lambda, real *V, tensor vir, real *dVdl) { - int g, j, m; - dvec dev; - double ndr, invdr; - real k, dkdl; - t_pullgrp *pgrp; + int c, j, m; + double dev, ndr, invdr; + real k, dkdl; + t_pull_coord *pcrd; - /* loop over the groups that are being pulled */ + /* loop over the pull coordinates */ *V = 0; *dVdl = 0; - for (g = 1; g < 1+pull->ngrp; g++) + for (c = 0; c < pull->ncoord; c++) { - pgrp = &pull->grp[g]; - get_pullgrp_distance(pull, pbc, g, t, pgrp->dr, dev); + pcrd = &pull->coord[c]; + + get_pull_coord_distance(pull, c, pbc, t, pcrd->dr, &dev); - k = (1.0 - lambda)*pgrp->k + lambda*pgrp->kB; - dkdl = pgrp->kB - pgrp->k; + k = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB; + dkdl = pcrd->kB - pcrd->k; switch (pull->eGeom) { case epullgDIST: - ndr = dnorm(pgrp->dr); + ndr = dnorm(pcrd->dr); invdr = 1/ndr; if (ePull == epullUMBRELLA) { - pgrp->f_scal = -k*dev[0]; - *V += 0.5* k*dsqr(dev[0]); - *dVdl += 0.5*dkdl*dsqr(dev[0]); + pcrd->f_scal = -k*dev; + *V += 0.5* k*dsqr(dev); + *dVdl += 0.5*dkdl*dsqr(dev); } else { - pgrp->f_scal = -k; + pcrd->f_scal = -k; *V += k*ndr; *dVdl += dkdl*ndr; } for (m = 0; m < DIM; m++) { - pgrp->f[m] = pgrp->f_scal*pgrp->dr[m]*invdr; + pcrd->f[m] = pcrd->f_scal*pcrd->dr[m]*invdr; } break; case epullgDIR: @@ -936,41 +809,24 @@ static void do_pull_pot(int ePull, case epullgCYL: if (ePull == epullUMBRELLA) { - pgrp->f_scal = -k*dev[0]; - *V += 0.5* k*dsqr(dev[0]); - *dVdl += 0.5*dkdl*dsqr(dev[0]); + pcrd->f_scal = -k*dev; + *V += 0.5* k*dsqr(dev); + *dVdl += 0.5*dkdl*dsqr(dev); } else { ndr = 0; for (m = 0; m < DIM; m++) { - ndr += pgrp->vec[m]*pgrp->dr[m]; + ndr += pcrd->vec[m]*pcrd->dr[m]; } - pgrp->f_scal = -k; + pcrd->f_scal = -k; *V += k*ndr; *dVdl += dkdl*ndr; } for (m = 0; m < DIM; m++) { - pgrp->f[m] = pgrp->f_scal*pgrp->vec[m]; - } - break; - case epullgPOS: - for (m = 0; m < DIM; m++) - { - if (ePull == epullUMBRELLA) - { - pgrp->f[m] = -k*dev[m]; - *V += 0.5* k*dsqr(dev[m]); - *dVdl += 0.5*dkdl*dsqr(dev[m]); - } - else - { - pgrp->f[m] = -k*pull->dim[m]; - *V += k*pgrp->dr[m]*pull->dim[m]; - *dVdl += dkdl*pgrp->dr[m]*pull->dim[m]; - } + pcrd->f[m] = pcrd->f_scal*pcrd->vec[m]; } break; } @@ -982,7 +838,7 @@ static void do_pull_pot(int ePull, { for (m = 0; m < DIM; m++) { - vir[j][m] -= 0.5*pgrp->f[j]*pgrp->dr[m]; + vir[j][m] -= 0.5*pcrd->f[j]*pcrd->dr[m]; } } } @@ -1021,7 +877,7 @@ void pull_constraint(t_pull *pull, t_mdatoms *md, t_pbc *pbc, } static void make_local_pull_group(gmx_ga2la_t ga2la, - t_pullgrp *pg, int start, int end) + t_pull_group *pg, int start, int end) { int i, ii; @@ -1072,19 +928,16 @@ void dd_make_local_pull_groups(gmx_domdec_t *dd, t_pull *pull, t_mdatoms *md) ga2la = NULL; } - if (pull->grp[0].nat > 0) + for (g = 0; g < pull->ngroup; g++) { - make_local_pull_group(ga2la, &pull->grp[0], md->start, md->start+md->homenr); - } - for (g = 1; g < 1+pull->ngrp; g++) - { - make_local_pull_group(ga2la, &pull->grp[g], md->start, md->start+md->homenr); + make_local_pull_group(ga2la, &pull->group[g], + md->start, md->start+md->homenr); } } static void init_pull_group_index(FILE *fplog, t_commrec *cr, int start, int end, - int g, t_pullgrp *pg, ivec pulldims, + int g, t_pull_group *pg, ivec pulldims, gmx_mtop_t *mtop, t_inputrec *ir, real lambda) { int i, ii, d, nfrozen, ndim; @@ -1263,9 +1116,8 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], gmx_bool bOutFile, unsigned long Flags) { t_pull *pull; - t_pullgrp *pgrp; - int g, start = 0, end = 0, m; - gmx_bool bCite; + t_pull_group *pgrp; + int c, g, start = 0, end = 0, m; pull = ir->pull; @@ -1279,30 +1131,39 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], if (fplog) { - fprintf(fplog, "\nWill apply %s COM pulling in geometry '%s'\n", - EPULLTYPE(ir->ePull), EPULLGEOM(pull->eGeom)); - if (pull->grp[0].nat > 0) + gmx_bool bAbs, bCos; + + bAbs = FALSE; + for (c = 0; c < pull->ncoord; c++) { - fprintf(fplog, "between a reference group and %d group%s\n", - pull->ngrp, pull->ngrp == 1 ? "" : "s"); + if (pull->group[pull->coord[c].group[0]].nat == 0 || + pull->group[pull->coord[c].group[1]].nat == 0) + { + bAbs = TRUE; + } } - else + + fprintf(fplog, "\nWill apply %s COM pulling in geometry '%s'\n", + EPULLTYPE(ir->ePull), EPULLGEOM(pull->eGeom)); + fprintf(fplog, "with %d pull coordinate%s and %d group%s\n", + pull->ncoord, pull->ncoord == 1 ? "" : "s", + pull->ngroup, pull->ngroup == 1 ? "" : "s"); + if (bAbs) { - fprintf(fplog, "with an absolute reference on %d group%s\n", - pull->ngrp, pull->ngrp == 1 ? "" : "s"); + fprintf(fplog, "with an absolute reference\n"); } - bCite = FALSE; - for (g = 0; g < pull->ngrp+1; g++) + bCos = FALSE; + for (g = 0; g < pull->ngroup; g++) { - if (pull->grp[g].nat > 1 && - pull->grp[g].pbcatom < 0) + if (pull->group[g].nat > 1 && + pull->group[g].pbcatom < 0) { /* We are using cosine weighting */ fprintf(fplog, "Cosine weighting is used for group %d\n", g); - bCite = TRUE; + bCos = TRUE; } } - if (bCite) + if (bCos) { please_cite(fplog, "Engin2010"); } @@ -1330,9 +1191,9 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], pull->dbuf_cyl = NULL; pull->bRefAt = FALSE; pull->cosdim = -1; - for (g = 0; g < pull->ngrp+1; g++) + for (g = 0; g < pull->ngroup; g++) { - pgrp = &pull->grp[g]; + pgrp = &pull->group[g]; pgrp->epgrppbc = epgrppbcNONE; if (pgrp->nat > 0) { @@ -1381,11 +1242,23 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], /* if we use dynamic reference groups, do some initialising for them */ if (PULL_CYL(pull)) { - if (pull->grp[0].nat == 0) + if (ir->ePull == epullCONSTRAINT && pull->ncoord > 1) { - gmx_fatal(FARGS, "Dynamic reference groups are not supported when using absolute reference!\n"); + /* We can't easily update the single reference group with multiple + * constraints. This would require recalculating COMs. + */ + gmx_fatal(FARGS, "Constraint COM pulling supports only one coordinate with geometry=cylinder, you can use umbrella pulling with multiple coordinates"); } - snew(pull->dyna, pull->ngrp+1); + + for (c = 0; c < pull->ncoord; c++) + { + if (pull->group[pull->coord[c].group[0]].nat == 0) + { + gmx_fatal(FARGS, "Dynamic reference groups are not supported when using absolute reference!\n"); + } + } + + snew(pull->dyna, pull->ncoord); } /* Only do I/O when we are doing dynamics and if we are the MASTER */ diff --git a/src/gromacs/mdlib/pullutil.c b/src/gromacs/mdlib/pullutil.c index 1b50e5ffaa..cb6d2dd42b 100644 --- a/src/gromacs/mdlib/pullutil.c +++ b/src/gromacs/mdlib/pullutil.c @@ -56,7 +56,7 @@ #include "pull.h" #include "gmx_ga2la.h" -static void pull_set_pbcatom(t_commrec *cr, t_pullgrp *pg, +static void pull_set_pbcatom(t_commrec *cr, t_pull_group *pgrp, t_mdatoms *md, rvec *x, rvec x_pbc) { @@ -66,14 +66,14 @@ static void pull_set_pbcatom(t_commrec *cr, t_pullgrp *pg, { if (DOMAINDECOMP(cr)) { - if (!ga2la_get_home(cr->dd->ga2la, pg->pbcatom, &a)) + if (!ga2la_get_home(cr->dd->ga2la, pgrp->pbcatom, &a)) { a = -1; } } else { - a = pg->pbcatom; + a = pgrp->pbcatom; } if (a >= md->start && a < md->start+md->homenr) @@ -87,7 +87,7 @@ static void pull_set_pbcatom(t_commrec *cr, t_pullgrp *pg, } else { - copy_rvec(x[pg->pbcatom], x_pbc); + copy_rvec(x[pgrp->pbcatom], x_pbc); } } @@ -98,15 +98,15 @@ static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull, int g, n, m; n = 0; - for (g = 0; g < 1+pull->ngrp; g++) + for (g = 0; g < pull->ngroup; g++) { - if ((g == 0 && PULL_CYL(pull)) || pull->grp[g].pbcatom == -1) + if ((g == 0 && PULL_CYL(pull)) || pull->group[g].pbcatom == -1) { clear_rvec(x_pbc[g]); } else { - pull_set_pbcatom(cr, &pull->grp[g], md, x, x_pbc[g]); + pull_set_pbcatom(cr, &pull->group[g], md, x, x_pbc[g]); for (m = 0; m < DIM; m++) { if (pull->dim[m] == 0) @@ -121,7 +121,7 @@ static void pull_set_pbcatoms(t_commrec *cr, t_pull *pull, if (cr && PAR(cr) && n > 0) { /* Sum over the nodes to get x_pbc from the home node of pbcatom */ - gmx_sum((1+pull->ngrp)*DIM, x_pbc[0], cr); + gmx_sum(pull->ngroup*DIM, x_pbc[0], cr); } } @@ -149,15 +149,16 @@ static real get_weight(real x, real r1, real r0) static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, t_pbc *pbc, double t, rvec *x, rvec *xp) { - int g, i, ii, m, start, end; - rvec g_x, dx, dir; - double r0_2, sum_a, sum_ap, dr2, mass, weight, wmass, wwmass, inp; - t_pullgrp *pref, *pgrp, *pdyna; - gmx_ga2la_t ga2la = NULL; + int c, i, ii, m, start, end; + rvec g_x, dx, dir; + double r0_2, sum_a, sum_ap, dr2, mass, weight, wmass, wwmass, inp; + t_pull_coord *pcrd; + t_pull_group *pref, *pgrp, *pdyna; + gmx_ga2la_t ga2la = NULL; if (pull->dbuf_cyl == NULL) { - snew(pull->dbuf_cyl, pull->ngrp*4); + snew(pull->dbuf_cyl, pull->ncoord*4); } if (cr && DOMAINDECOMP(cr)) @@ -171,12 +172,15 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, r0_2 = dsqr(pull->cyl_r0); /* loop over all groups to make a reference group for each*/ - pref = &pull->grp[0]; - for (g = 1; g < 1+pull->ngrp; g++) + for (c = 0; c < pull->ncoord; c++) { - pgrp = &pull->grp[g]; - pdyna = &pull->dyna[g]; - copy_rvec(pgrp->vec, dir); + pcrd = &pull->coord[c]; + + /* pref will be the same group for all pull coordinates */ + pref = &pull->group[pcrd->group[0]]; + pgrp = &pull->group[pcrd->group[1]]; + pdyna = &pull->dyna[c]; + copy_rvec(pcrd->vec, dir); sum_a = 0; sum_ap = 0; wmass = 0; @@ -185,13 +189,13 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, for (m = 0; m < DIM; m++) { - g_x[m] = pgrp->x[m] - pgrp->vec[m]*(pgrp->init[0] + pgrp->rate*t); + g_x[m] = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t); } /* loop over all atoms in the main ref group */ for (i = 0; i < pref->nat; i++) { - ii = pull->grp[0].ind[i]; + ii = pref->ind[i]; if (ga2la) { if (!ga2la_get_home(ga2la, pref->ind[i], &ii)) @@ -235,42 +239,44 @@ static void make_cyl_refgrps(t_commrec *cr, t_pull *pull, t_mdatoms *md, } } } - pull->dbuf_cyl[(g-1)*4+0] = wmass; - pull->dbuf_cyl[(g-1)*4+1] = wwmass; - pull->dbuf_cyl[(g-1)*4+2] = sum_a; - pull->dbuf_cyl[(g-1)*4+3] = sum_ap; + pull->dbuf_cyl[c*4+0] = wmass; + pull->dbuf_cyl[c*4+1] = wwmass; + pull->dbuf_cyl[c*4+2] = sum_a; + pull->dbuf_cyl[c*4+3] = sum_ap; } if (cr && PAR(cr)) { /* Sum the contributions over the nodes */ - gmx_sumd(pull->ngrp*4, pull->dbuf_cyl, cr); + gmx_sumd(pull->ncoord*4, pull->dbuf_cyl, cr); } - for (g = 1; g < 1+pull->ngrp; g++) + for (c = 0; c < pull->ncoord; c++) { - pgrp = &pull->grp[g]; - pdyna = &pull->dyna[g]; + pcrd = &pull->coord[c]; + + pdyna = &pull->dyna[c]; + pgrp = &pull->group[pcrd->group[1]]; - wmass = pull->dbuf_cyl[(g-1)*4+0]; - wwmass = pull->dbuf_cyl[(g-1)*4+1]; + wmass = pull->dbuf_cyl[c*4+0]; + wwmass = pull->dbuf_cyl[c*4+1]; pdyna->wscale = wmass/wwmass; pdyna->invtm = 1.0/(pdyna->wscale*wmass); for (m = 0; m < DIM; m++) { - g_x[m] = pgrp->x[m] - pgrp->vec[m]*(pgrp->init[0] + pgrp->rate*t); - pdyna->x[m] = g_x[m] + pgrp->vec[m]*pull->dbuf_cyl[(g-1)*4+2]/wmass; + g_x[m] = pgrp->x[m] - pcrd->vec[m]*(pcrd->init + pcrd->rate*t); + pdyna->x[m] = g_x[m] + pcrd->vec[m]*pull->dbuf_cyl[c*4+2]/wmass; if (xp) { - pdyna->xp[m] = g_x[m] + pgrp->vec[m]*pull->dbuf_cyl[(g-1)*4+3]/wmass; + pdyna->xp[m] = g_x[m] + pcrd->vec[m]*pull->dbuf_cyl[c*4+3]/wmass; } } if (debug) { fprintf(debug, "Pull cylinder group %d:%8.3f%8.3f%8.3f m:%8.3f\n", - g, pdyna->x[0], pdyna->x[1], + c, pdyna->x[0], pdyna->x[1], pdyna->x[2], 1.0/pdyna->invtm); } } @@ -293,21 +299,21 @@ void pull_calc_coms(t_commrec *cr, t_pull *pull, t_mdatoms *md, t_pbc *pbc, double t, rvec x[], rvec *xp) { - int g, i, ii, m; - real mass, w, wm, twopi_box = 0; - double wmass, wwmass, invwmass; - dvec com, comp; - double cm, sm, cmp, smp, ccm, csm, ssm, csw, snw; - rvec *xx[2], x_pbc = {0, 0, 0}, dx; - t_pullgrp *pgrp; + int g, i, ii, m; + real mass, w, wm, twopi_box = 0; + double wmass, wwmass, invwmass; + dvec com, comp; + double cm, sm, cmp, smp, ccm, csm, ssm, csw, snw; + rvec *xx[2], x_pbc = {0, 0, 0}, dx; + t_pull_group *pgrp; if (pull->rbuf == NULL) { - snew(pull->rbuf, 1+pull->ngrp); + snew(pull->rbuf, pull->ngroup); } if (pull->dbuf == NULL) { - snew(pull->dbuf, 3*(1+pull->ngrp)); + snew(pull->dbuf, 3*pull->ngroup); } if (pull->bRefAt) @@ -327,9 +333,9 @@ void pull_calc_coms(t_commrec *cr, twopi_box = 2.0*M_PI/pbc->box[pull->cosdim][pull->cosdim]; } - for (g = 0; g < 1+pull->ngrp; g++) + for (g = 0; g < pull->ngroup; g++) { - pgrp = &pull->grp[g]; + pgrp = &pull->group[g]; clear_dvec(com); clear_dvec(comp); wmass = 0; @@ -448,12 +454,12 @@ void pull_calc_coms(t_commrec *cr, if (cr && PAR(cr)) { /* Sum the contributions over the nodes */ - gmx_sumd((1+pull->ngrp)*3*DIM, pull->dbuf[0], cr); + gmx_sumd(pull->ngroup*3*DIM, pull->dbuf[0], cr); } - for (g = 0; g < 1+pull->ngrp; g++) + for (g = 0; g < pull->ngroup; g++) { - pgrp = &pull->grp[g]; + pgrp = &pull->group[g]; if (pgrp->nat > 0 && !(g == 0 && PULL_CYL(pull))) { if (pgrp->epgrppbc != epgrppbcCOS) diff --git a/src/programs/gmx/tpbcmp.c b/src/programs/gmx/tpbcmp.c index d4a7130f29..0c92274be4 100644 --- a/src/programs/gmx/tpbcmp.c +++ b/src/programs/gmx/tpbcmp.c @@ -900,10 +900,10 @@ static void comp_pull_AB(FILE *fp, t_pull *pull, real ftol, real abstol) { int i; - for (i = 0; i < pull->ngrp+1; i++) + for (i = 0; i < pull->ncoord; i++) { - fprintf(fp, "comparing pull group %d\n", i); - cmp_real(fp, "pullgrp->k", -1, pull->grp[i].k, pull->grp[i].kB, ftol, abstol); + fprintf(fp, "comparing pull coord %d\n", i); + cmp_real(fp, "pull-coord->k", -1, pull->coord[i].k, pull->coord[i].kB, ftol, abstol); } } -- 2.22.0