From 3bd1342556994115a8e4bcd3aad8debca607a13c Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Wed, 18 Feb 2015 14:37:18 +0100 Subject: [PATCH] Added pull geometry direction-relative The new pull geometry direction-relative enables to pull along a vector defines by two (additional) pull groups. Change-Id: I210f6bfbd29dcb69664298bab6472919ca88091d --- docs/manual/plots/pulldirrel.eps | 257 +++++++++++++++++++++ docs/manual/plots/pulldirrel.fig | 47 ++++ docs/manual/special.tex | 23 +- docs/user-guide/mdp-options.rst | 17 +- src/gromacs/fileio/tpxio.c | 10 +- src/gromacs/gmxlib/names.c | 2 +- src/gromacs/gmxlib/txtdump.c | 5 + src/gromacs/gmxlib/warninp.c | 8 +- src/gromacs/gmxpreprocess/readpull.c | 34 +-- src/gromacs/legacyheaders/types/enums.h | 2 +- src/gromacs/legacyheaders/types/inputrec.h | 3 +- src/gromacs/pulling/pull.c | 166 ++++++++++--- src/gromacs/pulling/pull.h | 49 ++-- 13 files changed, 532 insertions(+), 91 deletions(-) create mode 100644 docs/manual/plots/pulldirrel.eps create mode 100644 docs/manual/plots/pulldirrel.fig diff --git a/docs/manual/plots/pulldirrel.eps b/docs/manual/plots/pulldirrel.eps new file mode 100644 index 0000000000..a3ecd3f2f1 --- /dev/null +++ b/docs/manual/plots/pulldirrel.eps @@ -0,0 +1,257 @@ +%!PS-Adobe-3.0 EPSF-3.0 +%%Title: pulldirrel.fig +%%Creator: fig2dev Version 3.2 Patchlevel 5e +%%CreationDate: Tue Feb 24 12:19:24 2015 +%%BoundingBox: 0 0 213 142 +%Magnification: 1.0000 +%%EndComments +%%BeginProlog +/$F2psDict 200 dict def +$F2psDict begin +$F2psDict /mtrx matrix put +/col-1 {0 setgray} bind def +/col0 {0.000 0.000 0.000 srgb} bind def +/col1 {0.000 0.000 1.000 srgb} bind def +/col2 {0.000 1.000 0.000 srgb} bind def +/col3 {0.000 1.000 1.000 srgb} bind def +/col4 {1.000 0.000 0.000 srgb} bind def +/col5 {1.000 0.000 1.000 srgb} bind def +/col6 {1.000 1.000 0.000 srgb} bind def +/col7 {1.000 1.000 1.000 srgb} bind def +/col8 {0.000 0.000 0.560 srgb} bind def +/col9 {0.000 0.000 0.690 srgb} bind def +/col10 {0.000 0.000 0.820 srgb} bind def +/col11 {0.530 0.810 1.000 srgb} bind def +/col12 {0.000 0.560 0.000 srgb} bind def +/col13 {0.000 0.690 0.000 srgb} bind def +/col14 {0.000 0.820 0.000 srgb} bind def +/col15 {0.000 0.560 0.560 srgb} bind def +/col16 {0.000 0.690 0.690 srgb} bind def +/col17 {0.000 0.820 0.820 srgb} bind def +/col18 {0.560 0.000 0.000 srgb} bind def +/col19 {0.690 0.000 0.000 srgb} bind def +/col20 {0.820 0.000 0.000 srgb} bind def +/col21 {0.560 0.000 0.560 srgb} bind def +/col22 {0.690 0.000 0.690 srgb} bind def +/col23 {0.820 0.000 0.820 srgb} bind def +/col24 {0.500 0.190 0.000 srgb} bind def +/col25 {0.630 0.250 0.000 srgb} bind def +/col26 {0.750 0.380 0.000 srgb} bind def +/col27 {1.000 0.500 0.500 srgb} bind def +/col28 {1.000 0.630 0.630 srgb} bind def +/col29 {1.000 0.750 0.750 srgb} bind def +/col30 {1.000 0.880 0.880 srgb} bind def +/col31 {1.000 0.840 0.000 srgb} bind def + +end + +/cp {closepath} bind def +/ef {eofill} bind def +/gr {grestore} bind def +/gs {gsave} bind def +/sa {save} bind def +/rs {restore} bind def +/l {lineto} bind def +/m {moveto} bind def +/rm {rmoveto} bind def +/n {newpath} bind def +/s {stroke} bind def +/sh {show} bind def +/slc {setlinecap} bind def +/slj {setlinejoin} bind def +/slw {setlinewidth} bind def +/srgb {setrgbcolor} bind def +/rot {rotate} bind def +/sc {scale} bind def +/sd {setdash} bind def +/ff {findfont} bind def +/sf {setfont} bind def +/scf {scalefont} bind def +/sw {stringwidth} bind def +/tr {translate} bind def +/tnt {dup dup currentrgbcolor + 4 -2 roll dup 1 exch sub 3 -1 roll mul add + 4 -2 roll dup 1 exch sub 3 -1 roll mul add + 4 -2 roll dup 1 exch sub 3 -1 roll mul add srgb} + bind def +/shd {dup dup currentrgbcolor 4 -2 roll mul 4 -2 roll mul + 4 -2 roll mul srgb} bind def + /DrawEllipse { + /endangle exch def + /startangle exch def + /yrad exch def + /xrad exch def + /y exch def + /x exch def + /savematrix mtrx currentmatrix def + x y tr xrad yrad sc 0 0 1 startangle endangle arc + closepath + savematrix setmatrix + } def + +/$F2psBegin {$F2psDict begin /$F2psEnteredState save def} def +/$F2psEnd {$F2psEnteredState restore end} def + +/pageheader { +save +newpath 0 142 moveto 0 0 lineto 213 0 lineto 213 142 lineto closepath clip newpath +-35.7 148.5 translate +1 -1 scale +$F2psBegin +10 setmiterlimit +0 slj 0 slc + 0.06299 0.06299 sc +} bind def +/pagefooter { +$F2psEnd +restore +} bind def +%%EndProlog +pageheader +% +% Fig objects follow +% +% +% here starts figure with depth 51 +% Ellipse +15.000 slw +gs +2025 1800 tr +-90.000 rot +n 0 0 318 318 0 360 DrawEllipse 90.000 rot +gs col0 s gr +gr + +% Ellipse +gs +3600 2025 tr +-90.000 rot +n 0 0 318 318 0 360 DrawEllipse 90.000 rot +gs col0 s gr +gr + +% Ellipse +gs +900 675 tr +-90.000 rot +n 0 0 318 318 0 360 DrawEllipse 90.000 rot +gs col0 s gr +gr + +% Ellipse +gs +2025 450 tr +-90.000 rot +n 0 0 318 318 0 360 DrawEllipse 90.000 rot +gs col0 s gr +gr + +% Polyline +0 slj +0 slc +n 810 855 m 990 945 l + 1080 765 l gs col0 s gr +% Polyline +n 3690 1845 m 3510 1755 l + 3420 1935 l gs col0 s gr +% Polyline + [90] 0 sd +gs clippath +2540 836 m 2598 720 l 2544 693 l 2486 810 l 2486 810 l 2554 743 l 2540 836 l cp +2049 1683 m 1991 1800 l 2045 1826 l 2103 1709 l 2103 1709 l 2036 1777 l 2049 1683 l cp +eoclip +n 2025 1800 m + 2565 720 l gs col0 s gr gr + [] 0 sd +% arrowhead +n 2049 1683 m 2036 1777 l 2103 1709 l col0 s +% arrowhead +n 2540 836 m 2554 743 l 2486 810 l col0 s +% Polyline +n 2340 1170 m 2160 1080 l + 2070 1260 l gs col0 s gr +/Helvetica ff 285.75 scf sf +900 630 m +gs 1 -1 sc (3) dup sw pop 2 div neg 0 rm col0 sh gr +/Helvetica ff 285.75 scf sf +3600 2295 m +gs 1 -1 sc (4) dup sw pop 2 div neg 0 rm col0 sh gr +/Helvetica ff 285.75 scf sf +2520 1170 m +gs 1 -1 sc (d) dup sw pop 2 div neg 0 rm col0 sh gr +/Helvetica ff 222.25 scf sf +2655 1260 m +gs 1 -1 sc (p) dup sw pop 2 div neg 0 rm col0 sh gr +/Helvetica ff 285.75 scf sf +2025 405 m +gs 1 -1 sc (2) dup sw pop 2 div neg 0 rm col0 sh gr +/Helvetica ff 285.75 scf sf +2025 2070 m +gs 1 -1 sc (1) dup sw pop 2 div neg 0 rm col0 sh gr +% Polyline + [90] 0 sd +gs clippath +3401 1993 m 3586 2085 l 3640 1978 l 3454 1885 l 3454 1885 l 3579 2015 l 3401 1993 l cp +eoclip +n 900 675 m + 3600 2025 l gs col0 s gr gr + [] 0 sd +% arrowhead +n 3401 1993 m 3579 2015 l 3454 1885 l 3458 1954 l 3401 1993 l + cp gs 0.00 setgray ef gr col0 s +% Polyline + [90] 0 sd +gs clippath +2085 641 m 2085 435 l 1965 435 l 1965 641 l 1965 641 l 2025 473 l 2085 641 l cp +eoclip +n 2025 450 m + 2025 1800 l gs col0 s gr gr + [] 0 sd +% arrowhead +n 2085 641 m 2025 473 l 1965 641 l 2025 608 l 2085 641 l + cp gs 0.00 setgray ef gr col0 s +% Polyline +gs clippath +694 985 m 628 1118 l 708 1158 l 774 1025 l 774 1025 l 684 1106 l 694 985 l cp +eoclip +n 900 675 m + 675 1125 l gs col1 s gr gr + +% arrowhead +n 694 985 m 684 1106 l 774 1025 l 724 1025 l 694 985 l + cp gs col1 1.00 shd ef gr col1 s +% Polyline +gs clippath +1264 1369 m 1131 1303 l 1091 1383 l 1224 1449 l 1224 1449 l 1144 1359 l 1264 1369 l cp +eoclip +n 2025 1800 m + 1125 1350 l gs col4 s gr gr + +% arrowhead +n 1264 1369 m 1144 1359 l 1224 1449 l 1224 1399 l 1264 1369 l + cp gs col4 1.00 shd ef gr col4 s +% Polyline +gs clippath +2785 880 m 2918 946 l 2958 866 l 2825 800 l 2825 800 l 2906 891 l 2785 880 l cp +eoclip +n 2025 450 m + 2925 900 l gs col4 s gr gr + +% arrowhead +n 2785 880 m 2906 891 l 2825 800 l 2825 850 l 2785 880 l + cp gs col4 1.00 shd ef gr col4 s +% Polyline +gs clippath +3805 1714 m 3871 1581 l 3791 1541 l 3725 1674 l 3725 1674 l 3816 1594 l 3805 1714 l cp +eoclip +n 3600 2025 m + 3825 1575 l gs col1 s gr gr + +% arrowhead +n 3805 1714 m 3816 1594 l 3725 1674 l 3775 1674 l 3805 1714 l + cp gs col1 1.00 shd ef gr col1 s +% here ends figure; +pagefooter +showpage +%%Trailer +%EOF diff --git a/docs/manual/plots/pulldirrel.fig b/docs/manual/plots/pulldirrel.fig new file mode 100644 index 0000000000..39adeba728 --- /dev/null +++ b/docs/manual/plots/pulldirrel.fig @@ -0,0 +1,47 @@ +#FIG 3.2 Produced by xfig version 3.2.5c +Landscape +Center +Metric +A4 +100.00 +Single +-2 +1200 2 +1 3 0 2 0 7 51 -1 -1 0.000 0 1.5708 2025 1800 318 318 2025 1800 2250 2025 +1 3 0 2 0 7 51 -1 -1 0.000 0 1.5708 3600 2025 318 318 3600 2025 3825 2250 +1 3 0 2 0 7 51 -1 -1 0.000 0 1.5708 900 675 318 318 900 675 1125 900 +1 3 0 2 0 7 51 -1 -1 0.000 0 1.5708 2025 450 318 318 2025 450 2250 675 +2 1 0 2 0 7 51 -1 -1 0.000 0 0 -1 0 0 3 + 810 855 990 945 1080 765 +2 1 0 2 0 7 51 -1 -1 0.000 0 0 -1 0 0 3 + 3690 1845 3510 1755 3420 1935 +2 1 1 2 0 7 50 -1 -1 6.000 0 0 -1 1 0 2 + 2 1 2.00 120.00 135.00 + 900 675 3600 2025 +2 1 1 2 0 7 50 -1 -1 6.000 0 0 -1 0 1 2 + 2 1 2.00 120.00 135.00 + 2025 450 2025 1800 +2 1 0 2 1 7 50 -1 -1 0.000 0 0 -1 1 0 2 + 2 1 2.00 90.00 90.00 + 900 675 675 1125 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 -1 1 0 2 + 2 1 2.00 90.00 90.00 + 2025 1800 1125 1350 +2 1 0 2 4 7 50 -1 -1 0.000 0 0 -1 1 0 2 + 2 1 2.00 90.00 90.00 + 2025 450 2925 900 +2 1 0 2 1 7 50 -1 -1 0.000 0 0 -1 1 0 2 + 2 1 2.00 90.00 90.00 + 3600 2025 3825 1575 +2 1 1 2 0 7 51 -1 -1 6.000 0 0 -1 1 1 2 + 0 0 2.00 60.00 90.00 + 0 0 2.00 60.00 90.00 + 2025 1800 2565 720 +2 1 0 2 0 7 51 -1 -1 0.000 0 0 -1 0 0 3 + 2340 1170 2160 1080 2070 1260 +4 1 0 51 -1 16 18 0.0000 4 210 165 900 630 3\001 +4 1 0 51 -1 16 18 0.0000 4 210 165 3600 2295 4\001 +4 1 0 51 -1 16 18 0.0000 4 225 165 2520 1170 d\001 +4 1 0 51 -1 16 14 0.0000 4 165 135 2655 1260 p\001 +4 1 0 51 -1 16 18 0.0000 4 210 165 2025 405 2\001 +4 1 0 51 -1 16 18 0.0000 4 210 165 2025 2070 1\001 diff --git a/docs/manual/special.tex b/docs/manual/special.tex index d79a680e58..045f5a4dbc 100644 --- a/docs/manual/special.tex +++ b/docs/manual/special.tex @@ -180,7 +180,7 @@ a canonical ensemble of the initial state A and $\beta=1/k_B T$. The pull code applies forces or constraints between the centers of mass of one or more pairs of groups of atoms. 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 +on usually two, but sometimes more, 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 @@ -194,7 +194,7 @@ weighting factor can also be used. with umbrella pulling. $V_{rup}$ is the velocity at which the spring is retracted, $Z_{link}$ is the atom to which the spring is attached and $Z_{spring}$ is the location of the spring.} -\label{fi:pull} +\label{fig:pull} \end{figure} Three different types of calculation are supported, @@ -264,7 +264,7 @@ of 0.79 times the radius. reference group applied to interface systems. C is the reference group. The circles represent the center of mass of two groups plus the reference group, $d_c$ is the reference distance.} -\label{fi:pullref} +\label{fig:pullref} \end{figure} For a group of molecules in a periodic system, a plain reference group @@ -306,6 +306,23 @@ as follows: \ve{F}_{\!i} = \frac{w'_i \, m_i}{M} \, \ve{F}_{\!com} \eeq +\subsubsection{Definition of the pull direction} + +The most common setup is to pull along the direction of the vector containing +the two pull groups, this is selected with +{\tt pull-coord?-geometry = distance}. You might want to pull along a certain +vector instead, which is selected with {\tt pull-coord?-geometry = direction}. +But this can cause unwanted torque forces in the system, unless you pull against a reference group with (nearly) fixed orientation, e.g. a membrane protein embedded in a membrane along x/y while pulling along z. If your reference group does not have a fixed orientation, you should probably use +{\tt pull-coord?-geometry = direction-relative}, see \figref{pulldirrel}. +Since the potential now depends on the coordinates of two additional groups defining the orientation, the torque forces will work on these two groups. + +\begin{figure} +\centerline{\includegraphics[width=5cm]{plots/pulldirrel}} +\caption{The pull setup for geometry {\tt direction-relative}. The ``normal'' pull groups are 1 and 2. Groups 3 and 4 define the pull direction and thus the direction of the normal pull forces (red). This leads to reaction forces (blue) on groups 3 and 4, which are perpendicular to the pull direction. Their magnitude is given by the ``normal'' pull force times the ratio of $d_p$ and the distance between groups 3 and 4.} +\label{fig:pulldirrel} +\end{figure} + + \subsubsection{Limitations} There is one theoretical limitation: strictly speaking, constraint forces can only be calculated between diff --git a/docs/user-guide/mdp-options.rst b/docs/user-guide/mdp-options.rst index b1585e6bf7..758de39b48 100644 --- a/docs/user-guide/mdp-options.rst +++ b/docs/user-guide/mdp-options.rst @@ -1690,7 +1690,7 @@ applicable pulling coordinate. their periodic image which is closest to :mdp:`pull-group1-pbcatom`. A value of 0 means that the middle atom (number wise) is used. This parameter is not used with - :mdp:`pull-geometry` cylinder. A value of -1 turns on cosine + :mdp:`pull-group1-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). @@ -1702,7 +1702,8 @@ applicable pulling coordinate. absolute reference of :mdp:`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. + motion. Note that (only) for :mdp:`pull-coord1-geometry` = + :mdp-value:`direction-relative` four groups are required. .. mdp:: pull-coord1-type: @@ -1748,6 +1749,18 @@ applicable pulling coordinate. dynamic (*e.g.* no pressure scaling) in the pull dimensions and the pull force is not added to virial. + .. mdp-value:: direction-relative + + As :mdp-value:`direction`, but the pull vector is the vector + that points from the COM of a third to the COM of a fourth pull + group. This means that 4 groups need to be supplied in + :mdp:`pull-coord1-groups`. Note that the pull force will give + rise to a torque on the pull vector, which is turn leads to + forces perpendicular to the pull vector on the two groups + defining the vector. If you want a pull group to move between + the two groups defining the vector, simply use the union of + these two groups as the reference group. + .. mdp-value:: cylinder Designed for pulling with respect to a layer where the reference diff --git a/src/gromacs/fileio/tpxio.c b/src/gromacs/fileio/tpxio.c index 84cb9d478e..e3defd42f6 100644 --- a/src/gromacs/fileio/tpxio.c +++ b/src/gromacs/fileio/tpxio.c @@ -91,7 +91,8 @@ enum tpxv { tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, /**< potentials for supporting coarse-grained force fields */ tpxv_InteractiveMolecularDynamics, /**< interactive molecular dynamics (IMD) */ tpxv_RemoveObsoleteParameters1, /**< remove optimize_fft, dihre_fc, nstcheckpoint */ - tpxv_PullCoordTypeGeom /**< add pull type and geometry per group and flat-bottom */ + tpxv_PullCoordTypeGeom, /**< add pull type and geometry per group and flat-bottom */ + tpxv_PullGeomDirRel /**< add pull geometry direction-relative */ }; /*! \brief Version number of the file format written to run input @@ -105,7 +106,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_PullCoordTypeGeom; +static const int tpx_version = tpxv_PullGeomDirRel; /* This number should only be increased when you edit the TOPOLOGY section @@ -299,6 +300,11 @@ static void do_pull_coord(t_fileio *fio, t_pull_coord *pcrd, int file_version, { gmx_fio_do_int(fio, pcrd->eType); gmx_fio_do_int(fio, pcrd->eGeom); + if (pcrd->eGeom == epullgDIRRELATIVE) + { + gmx_fio_do_int(fio, pcrd->group[2]); + gmx_fio_do_int(fio, pcrd->group[3]); + } gmx_fio_do_ivec(fio, pcrd->dim); } else diff --git a/src/gromacs/gmxlib/names.c b/src/gromacs/gmxlib/names.c index 3bae62a817..12e94beb9e 100644 --- a/src/gromacs/gmxlib/names.c +++ b/src/gromacs/gmxlib/names.c @@ -223,7 +223,7 @@ const char *epull_names[epullNR+1] = { }; const char *epullg_names[epullgNR+1] = { - "distance", "direction", "cylinder", "direction-periodic", NULL + "distance", "direction", "cylinder", "direction-periodic", "direction-relative", NULL }; const char *erotg_names[erotgNR+1] = { diff --git a/src/gromacs/gmxlib/txtdump.c b/src/gromacs/gmxlib/txtdump.c index ff8bbf07fb..20e3f3d8f3 100644 --- a/src/gromacs/gmxlib/txtdump.c +++ b/src/gromacs/gmxlib/txtdump.c @@ -634,6 +634,11 @@ 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]); + if (pcrd->eGeom == epullgDIRRELATIVE) + { + PI("group[2]", pcrd->group[2]); + PI("group[3]", pcrd->group[3]); + } PS("type", EPULLTYPE(pcrd->eType)); PS("geometry", EPULLGEOM(pcrd->eGeom)); pr_ivec(fp, indent, "dim", pcrd->dim, DIM, TRUE); diff --git a/src/gromacs/gmxlib/warninp.c b/src/gromacs/gmxlib/warninp.c index 1721a47012..62ce4aa2e4 100644 --- a/src/gromacs/gmxlib/warninp.c +++ b/src/gromacs/gmxlib/warninp.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. @@ -74,12 +74,10 @@ warninp_t init_warning(gmx_bool bAllowWarnings, int maxwarning) void set_warning_line(warninp_t wi, const char *s, int line) { - if (s == NULL) + if (s != NULL) { - gmx_incons("Calling set_warning_line with NULL pointer"); + strcpy(wi->filenm, s); } - - strcpy(wi->filenm, s); wi->lineno = line; } diff --git a/src/gromacs/gmxpreprocess/readpull.c b/src/gromacs/gmxpreprocess/readpull.c index dde971b5c5..c92eada292 100644 --- a/src/gromacs/gmxpreprocess/readpull.c +++ b/src/gromacs/gmxpreprocess/readpull.c @@ -129,7 +129,8 @@ static void init_pull_coord(t_pull_coord *pcrd, dvec origin, vec; char buf[STRLEN]; - if (pcrd->eType == epullCONSTRAINT && pcrd->eGeom == epullgCYL) + if (pcrd->eType == epullCONSTRAINT && (pcrd->eGeom == epullgCYL || + pcrd->eGeom == epullgDIRRELATIVE)) { gmx_fatal(FARGS, "Pulling of type %s can not be combined with geometry %s. Consider using pull type %s.", epull_names[pcrd->eType], @@ -145,6 +146,8 @@ static void init_pull_coord(t_pull_coord *pcrd, gmx_fatal(FARGS, "The pull origin can only be set with an absolute reference"); } + /* Check and set the pull vector */ + clear_dvec(vec); if (pcrd->eGeom == epullgDIST) { if (pcrd->init < 0) @@ -160,10 +163,8 @@ static void init_pull_coord(t_pull_coord *pcrd, * be using direction. We will do this later, since an already planned * generalization of the pull code makes pull dim available here. */ - - clear_dvec(vec); } - else + else if (pcrd->eGeom != epullgDIRRELATIVE) { string2dvec(vec_buf, vec); if (dnorm2(vec) == 0) @@ -188,7 +189,7 @@ char **read_pullparams(int *ninp_p, t_inpfile **inp_p, t_pull *pull, warninp_t wi) { - int ninp, nerror = 0, i, nchar, nscan, m, idum; + int ninp, i, nchar, nscan, m, idum; t_inpfile *inp; const char *tmp; char **grpbuf; @@ -258,20 +259,26 @@ char **read_pullparams(int *ninp_p, t_inpfile **inp_p, /* Read the pull coordinates */ for (i = 1; i < pull->ncoord + 1; i++) { + int ngroup; + 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 contain %d pull group indices\n", - 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); + + nscan = sscanf(groups, "%d %d %d %d %d", &pcrd->group[0], &pcrd->group[1], &pcrd->group[2], &pcrd->group[3], &idum); + ngroup = (pcrd->eGeom == epullgDIRRELATIVE) ? 4 : 2; + if (nscan != ngroup) + { + sprintf(wbuf, "%s should contain %d pull group indices with geometry %s", + buf, ngroup, epullg_names[pcrd->eGeom]); + set_warning_line(wi, NULL, -1); + warning_error(wi, wbuf); + } + sprintf(buf, "pull-coord%d-dim", i); STYPE(buf, dim_buf, "Y Y Y"); sprintf(buf, "pull-coord%d-origin", i); @@ -426,7 +433,6 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l int c, m; double t_start; real init = 0; - dvec dr; double dev, value; init_pull(NULL, ir, 0, NULL, mtop, NULL, oenv, lambda, FALSE, 0); @@ -462,7 +468,7 @@ void set_pull_init(t_inputrec *ir, gmx_mtop_t *mtop, rvec *x, matrix box, real l pcrd->init = 0; } - get_pull_coord_distance(pull, c, &pbc, t_start, dr, &dev); + get_pull_coord_distance(pull, c, &pbc, t_start, &dev); /* Calculate the value of the coordinate at t_start */ value = pcrd->init + t_start*pcrd->rate + dev; fprintf(stderr, " %10.3f nm", value); diff --git a/src/gromacs/legacyheaders/types/enums.h b/src/gromacs/legacyheaders/types/enums.h index 92bdcb8237..c86ce031ff 100644 --- a/src/gromacs/legacyheaders/types/enums.h +++ b/src/gromacs/legacyheaders/types/enums.h @@ -295,7 +295,7 @@ enum { }; enum { - epullgDIST, epullgDIR, epullgCYL, epullgDIRPBC, epullgNR + epullgDIST, epullgDIR, epullgCYL, epullgDIRPBC, epullgDIRRELATIVE, epullgNR }; /* Enforced rotation groups */ diff --git a/src/gromacs/legacyheaders/types/inputrec.h b/src/gromacs/legacyheaders/types/inputrec.h index 0bd0c85cd7..01b5451b7b 100644 --- a/src/gromacs/legacyheaders/types/inputrec.h +++ b/src/gromacs/legacyheaders/types/inputrec.h @@ -126,7 +126,7 @@ typedef struct { } t_pull_group; typedef struct { - int group[2]; /* The pull groups, index in group in t_pull */ + int group[4]; /* 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 */ @@ -140,6 +140,7 @@ typedef struct { /* Variables not present in mdp, but used at run time */ dvec dr; /* The distance from the reference group */ + double vec_len; /* Length of vec for direction-relative */ 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 */ diff --git a/src/gromacs/pulling/pull.c b/src/gromacs/pulling/pull.c index 539d015e19..d7cf619295 100644 --- a/src/gromacs/pulling/pull.c +++ b/src/gromacs/pulling/pull.c @@ -366,6 +366,42 @@ static void apply_forces_cyl_grp(const t_pull_group *pgrp, double dv_corr, } } +/* Apply torque forces in a mass weighted fashion to the groups that define + * the pull vector direction for pull coordinate pcrd. + */ +static void apply_forces_vec_torque(const t_pull *pull, + const t_pull_coord *pcrd, + const t_mdatoms *md, + rvec *f) +{ + double inpr; + int m; + dvec f_perp; + + /* The component inpr along the pull vector is accounted for in the usual + * way. Here we account for the component perpendicular to vec. + */ + inpr = 0; + for (m = 0; m < DIM; m++) + { + inpr += pcrd->dr[m]*pcrd->vec[m]; + } + /* The torque force works along the component of the distance vector + * of between the two "usual" pull groups that is perpendicular to + * the pull vector. The magnitude of this force is the "usual" scale force + * multiplied with the ratio of the distance between the two "usual" pull + * groups and the distance between the two groups that define the vector. + */ + for (m = 0; m < DIM; m++) + { + f_perp[m] = (pcrd->dr[m] - inpr*pcrd->vec[m])/pcrd->vec_len*pcrd->f_scal; + } + + /* Apply the force to the groups defining the vector using opposite signs */ + apply_forces_grp(&pull->group[pcrd->group[2]], md, f_perp, -1, f); + apply_forces_grp(&pull->group[pcrd->group[3]], md, f_perp, 1, f); +} + /* Apply forces in a mass weighted fashion */ static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f) { @@ -393,6 +429,14 @@ static void apply_forces(t_pull * pull, t_mdatoms * md, rvec *f) } else { + if (pcrd->eGeom == epullgDIRRELATIVE) + { + /* We need to apply the torque forces to the pull groups + * that define the pull vector. + */ + apply_forces_vec_torque(pull, pcrd, md, f); + } + if (pull->group[pcrd->group[0]].nat > 0) { apply_forces_grp(&pull->group[pcrd->group[0]], md, pcrd->f, -1, f); @@ -420,6 +464,9 @@ static double max_pull_distance2(const t_pull_coord *pcrd, const t_pbc *pbc) return 0.25*max_d2; } +/* This function returns the distance based on coordinates xg and xref. + * Note that the pull coordinate struct pcrd is not modified. + */ static void low_get_pull_coord_dr(const t_pull *pull, const t_pull_coord *pcrd, const t_pbc *pbc, double t, @@ -474,13 +521,15 @@ static void low_get_pull_coord_dr(const t_pull *pull, } } +/* This function returns the distance based on the contents of the pull struct. + * pull->coord[coord_ind].dr, and potentially vector, are updated. + */ static void get_pull_coord_dr(const t_pull *pull, int coord_ind, - const t_pbc *pbc, double t, - dvec dr) + const t_pbc *pbc, double t) { - double md2; - const t_pull_coord *pcrd; + double md2; + t_pull_coord *pcrd; pcrd = &pull->coord[coord_ind]; @@ -493,17 +542,48 @@ static void get_pull_coord_dr(const t_pull *pull, md2 = max_pull_distance2(pcrd, pbc); } + if (pcrd->eGeom == epullgDIRRELATIVE) + { + /* We need to determine the pull vector */ + const t_pull_group *pgrp2, *pgrp3; + dvec vec; + int m; + double invlen; + + pgrp2 = &pull->group[pcrd->group[2]]; + pgrp3 = &pull->group[pcrd->group[3]]; + + pbc_dx_d(pbc, pgrp3->x, pgrp2->x, vec); + + for (m = 0; m < DIM; m++) + { + vec[m] *= pcrd->dim[m]; + } + pcrd->vec_len = dnorm(vec); + for (m = 0; m < DIM; m++) + { + pcrd->vec[m] = vec[m]/pcrd->vec_len; + } + if (debug) + { + fprintf(debug, "pull coord %d vector: %6.3f %6.3f %6.3f normalized: %6.3f %6.3f %6.3f\n", + coord_ind, + vec[XX], vec[YY], vec[ZZ], + pcrd->vec[XX], pcrd->vec[YY], pcrd->vec[ZZ]); + } + } + low_get_pull_coord_dr(pull, pcrd, pbc, t, pull->group[pcrd->group[1]].x, pcrd->eGeom == epullgCYL ? pull->dyna[coord_ind].x : pull->group[pcrd->group[0]].x, md2, - dr); + pcrd->dr); } void get_pull_coord_distance(const t_pull *pull, int coord_ind, const t_pbc *pbc, double t, - dvec dr, double *dev) + double *dev) { static gmx_bool bWarned = FALSE; /* TODO: this should be fixed for thread-safety, but is fairly benign */ @@ -513,7 +593,7 @@ void get_pull_coord_distance(const t_pull *pull, pcrd = &pull->coord[coord_ind]; - get_pull_coord_dr(pull, coord_ind, pbc, t, dr); + get_pull_coord_dr(pull, coord_ind, pbc, t); ref = pcrd->init + pcrd->rate*t; @@ -526,7 +606,7 @@ void get_pull_coord_distance(const t_pull *pull, fprintf(stderr, "\nPull reference distance for coordinate %d is negative (%f)\n", coord_ind+1, ref); bWarned = TRUE; } - drs = dnorm(dr); + drs = dnorm(pcrd->dr); if (drs == 0) { /* With no vector we can not determine the direction for the force, @@ -542,12 +622,13 @@ void get_pull_coord_distance(const t_pull *pull, break; case epullgDIR: case epullgDIRPBC: + case epullgDIRRELATIVE: case epullgCYL: /* Pull along vec */ inpr = 0; for (m = 0; m < DIM; m++) { - inpr += pcrd->vec[m]*dr[m]; + inpr += pcrd->vec[m]*pcrd->dr[m]; } *dev = inpr - ref; break; @@ -615,13 +696,12 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, continue; } - get_pull_coord_dr(pull, c, pbc, t, r_ij[c]); + get_pull_coord_dr(pull, c, pbc, t); /* Store the difference vector at time t for printing */ - 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]); + c, pcrd->dr[XX], pcrd->dr[YY], pcrd->dr[ZZ]); } if (pcrd->eGeom == epullgDIR || pcrd->eGeom == epullgDIRPBC) @@ -630,13 +710,17 @@ static void do_constraint(t_pull *pull, t_pbc *pbc, a = 0; for (m = 0; m < DIM; m++) { - a += pcrd->vec[m]*r_ij[c][m]; + a += pcrd->vec[m]*pcrd->dr[m]; } for (m = 0; m < DIM; m++) { r_ij[c][m] = a*pcrd->vec[m]; } } + else + { + copy_dvec(pcrd->dr, r_ij[c]); + } if (dnorm2(r_ij[c]) == 0) { @@ -933,7 +1017,7 @@ static void do_pull_pot(t_pull *pull, t_pbc *pbc, double t, real lambda, continue; } - get_pull_coord_distance(pull, c, pbc, t, pcrd->dr, &dev); + get_pull_coord_distance(pull, c, pbc, t, &dev); k = (1.0 - lambda)*pcrd->k + lambda*pcrd->kB; dkdl = pcrd->kB - pcrd->k; @@ -1303,11 +1387,22 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], for (c = 0; c < pull->ncoord; c++) { t_pull_coord *pcrd; + int calc_com_start, calc_com_end, g; pcrd = &pull->coord[c]; if (pcrd->eType == epullCONSTRAINT) { + /* Check restrictions of the constraint pull code */ + if (pcrd->eGeom == epullgCYL || + pcrd->eGeom == epullgDIRRELATIVE) + { + 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]); + } + pull->bConstraint = TRUE; } else @@ -1318,28 +1413,16 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], 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) + /* We only need to calculate the plain COM of a group + * when it is not only used as a cylinder group. + */ + calc_com_start = (pcrd->eGeom == epullgCYL ? 1 : 0); + calc_com_end = (pcrd->eGeom == epullgDIRRELATIVE ? 4 : 2); + + for (g = calc_com_start; g <= calc_com_end; g++) { - pull->group[pcrd->group[1]].bCalcCOM = TRUE; + pull->group[pcrd->group[g]].bCalcCOM = TRUE; } } @@ -1409,9 +1492,6 @@ 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. @@ -1420,13 +1500,23 @@ void init_pull(FILE *fplog, t_inputrec *ir, int nfile, const t_filenm fnm[], * But since this is a very exotic case, we don't check for this. * A warning is printed in init_pull_group_index. */ + + gmx_bool bConstraint; + ivec pulldim, pulldim_con; + + /* Loop over all pull coordinates to see along which dimensions + * this group is pulled and if it is involved in constraints. + */ 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) + pull->coord[c].group[1] == g || + (pull->coord[c].eGeom == epullgDIRRELATIVE && + (pull->coord[c].group[2] == g || + pull->coord[c].group[3] == g))) { for (m = 0; m < DIM; m++) { diff --git a/src/gromacs/pulling/pull.h b/src/gromacs/pulling/pull.h index 7b23de006a..72a0f8bfe6 100644 --- a/src/gromacs/pulling/pull.h +++ b/src/gromacs/pulling/pull.h @@ -59,19 +59,20 @@ extern "C" { struct t_pbc; -/*! \brief Get the distance to the reference and deviation for pull coord coord_ind. - * - * \param[in] pull The pull group. - * \param[in] coord_ind Number of the pull coordinate. - * \param[in] pbc Information structure about periodicity. - * \param[in] t Time. - * \param[out] dr The pull coordinate difference vector. - * \param[out] dev The deviation from the reference distance. +/*! \brief Get the deviation for pull coord coord_ind. + * + * This call updates some data in the pull coordinates in \p pull. + * + * \param[in,out] pull The pull struct. + * \param[in] coord_ind Number of the pull coordinate. + * \param[in] pbc Information structure about periodicity. + * \param[in] t Time. + * \param[out] dev The deviation from the reference distance. */ void get_pull_coord_distance(const t_pull *pull, int coord_ind, const struct t_pbc *pbc, double t, - dvec dr, double *dev); + double *dev); /*! \brief Set the all the pull forces to zero. @@ -83,14 +84,14 @@ void clear_pull_forces(t_pull *pull); /*! \brief Determine the COM pull forces and add them to f, return the potential * - * \param[in] pull The pull group. - * \param[in] md All atoms. - * \param[in] pbc Information struct about periodicity. - * \param[in] cr Struct for communication info. - * \param[in] t Time. - * \param[in] lambda The value of lambda in FEP calculations. - * \param[in] x Positions. - * \param[in] f Forces. + * \param[in,out] pull The pull struct. + * \param[in] md All atoms. + * \param[in] pbc Information struct about periodicity. + * \param[in] cr Struct for communication info. + * \param[in] t Time. + * \param[in] lambda The value of lambda in FEP calculations. + * \param[in] x Positions. + * \param[in] f Forces. * \param[in,out] vir The virial, which, if != NULL, gets a pull correction. * \param[out] dvdlambda Pull contribution to dV/d(lambda). * @@ -104,13 +105,13 @@ real pull_potential(t_pull *pull, t_mdatoms *md, struct t_pbc *pbc, /*! \brief Constrain the coordinates xp in the directions in x * and also constrain v when v != NULL. * - * \param[in] pull The pull group. - * \param[in] md All atoms. - * \param[in] pbc Information struct about periodicity. - * \param[in] cr Struct for communication info. - * \param[in] dt The time step length. - * \param[in] t The time. - * \param[in] x Positions. + * \param[in,out] pull The pull data. + * \param[in] md All atoms. + * \param[in] pbc Information struct about periodicity. + * \param[in] cr Struct for communication info. + * \param[in] dt The time step length. + * \param[in] t The time. + * \param[in] x Positions. * \param[in,out] xp Updated x, can be NULL. * \param[in,out] v Velocities, which may get a pull correction. * \param[in,out] vir The virial, which, if != NULL, gets a pull correction. -- 2.22.0