--- /dev/null
+%!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
--- /dev/null
+#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
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
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,
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
\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
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).
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:
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
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
*
* 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
{
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
};
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] = {
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);
*
* 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.
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;
}
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],
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)
* 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)
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;
/* 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);
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);
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);
};
enum {
- epullgDIST, epullgDIR, epullgCYL, epullgDIRPBC, epullgNR
+ epullgDIST, epullgDIR, epullgCYL, epullgDIRPBC, epullgDIRRELATIVE, epullgNR
};
/* Enforced rotation groups */
} 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 */
/* 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 */
}
}
+/* 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)
{
}
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);
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,
}
}
+/* 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];
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 */
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;
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,
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;
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)
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)
{
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;
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
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;
}
}
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.
* 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++)
{
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.
/*! \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).
*
/*! \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.
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff