\begin{figure}
\centerline{\includegraphics[width=13.5cm]{plots/compelsetup.pdf}}
-\caption{Typical double-membrane setup for CompEL simulations (A, B). Plot (C) shows
-the potential difference $\Delta U$ resulting
+\caption{Typical double-membrane setup for CompEL simulations (A, B).
+Ion\,/\,water molecule exchanges will be performed as needed
+between the two light blue volumes around the dotted black lines (A).
+Plot (C) shows the potential difference $\Delta U$ resulting
from the selected charge imbalance $\Delta q_{ref}$ between the compartments.}
\label{fig:compelsetup}
\end{figure}
{\tt swap-frequency} determines how often a swap attempt will be made.
This step requires that the positions of the ions, solvent, and swap groups are
communicated between the parallel processes, so if chosen too small it can decrease the simulation
-performance.
+performance. The {\tt Position swapping} entry in the cycle and time accounting
+table at the end of the {\tt md.log} file summarizes the amount of runtime spent
+in the swap module.
{\small
\begin{verbatim}
-split-group0 = channel0 ; Defines compartment boundary
-split-group1 = channel1 ; Defines other compartment boundary
+split-group0 = channel0 ; Defines compartment boundary
+split-group1 = channel1 ; Defines other compartment boundary
massw-split0 = no ; use mass-weighted center?
massw-split1 = no
\end{verbatim}}
channel is selected as split group, the center of the channel will define the dividing
plane between the compartments (dashed horizontal line in the figure). All index groups
must be defined in the index file.
+
+If, to restore the requested ion counts, an ion from one compartment has to be exchanged
+with a water molecule from the other compartment, then those molecules are swapped
+which have the largest distance to the compartment-defining boundaries
+(dashed horizontal lines in \figref{compelsetup} A). Depending on the ion concentration,
+this effectively results in exchanges of molecules between the light blue volumes.
+If a channel is very asymmetric in $z$-direction and would extend into one of the
+swap volumes, one can offset the swap exchange plane with the {\tt bulk-offset}
+parameter. A value of 0.0 means no offset $b$, values $-1.0 < b < 0$ move the
+swap exchange plane closer to the lower, values $0 < b < 1.0$ closer to the upper
+membrane. The figure \figref{compelsetup} A (left) depicts that for the {\bf A} compartment.
+
{\small
\begin{verbatim}
swap-group = NA+_CL- ; Ions to be included in exchange
struct t_gmx_IMD *setup; /* Stores non-inputrec IMD data */
} t_IMD;
-/* Abstract types for position swapping only defined in swapcoords.c */
+/* Abstract types for position swapping only defined in swapcoords.cpp */
typedef struct t_swap *gmx_swapcoords_t;
typedef struct t_swapcoords {
- int nstswap; /* Every how many steps a swap is attempted? */
- int nat; /* Number of atoms in the ion group */
- int nat_split[2]; /* Number of atoms in the split group */
- int nat_sol; /* Number of atoms in the solvent group */
- atom_id *ind; /* The global ion group atoms numbers */
- atom_id *ind_split[2]; /* Split groups for compartment partitioning */
- atom_id *ind_sol; /* The global solvent group atom numbers */
- gmx_bool massw_split[2]; /* Use mass-weighted positions in split group? */
- real cyl0r, cyl1r; /* Split cylinders defined by radius, upper and */
- real cyl0u, cyl1u; /* ... lower extension. The split cylinders de- */
- real cyl0l, cyl1l; /* ... fine the channels and are each anchored */
- /* ... in the center of the split group */
- int nanions[eCompNR]; /* Requested number of anions and */
- int nAverage; /* Coupling constant (nr of swap attempt steps) */
- real threshold; /* Ion counts may deviate from the requested
- values by +-threshold before a swap is done */
- int ncations[eCompNR]; /* ... cations for both compartments */
- gmx_swapcoords_t si_priv; /* swap private data accessible in
- * swapcoords.c */
+ int nstswap; /* Every how many steps a swap is attempted? */
+ int nat; /* Number of atoms in the ion group */
+ int nat_split[2]; /* Number of atoms in the split group */
+ int nat_sol; /* Number of atoms in the solvent group */
+ atom_id *ind; /* The global ion group atoms numbers */
+ atom_id *ind_split[2]; /* Split groups for compartment partitioning */
+ atom_id *ind_sol; /* The global solvent group atom numbers */
+ gmx_bool massw_split[2]; /* Use mass-weighted positions in split group? */
+ real cyl0r, cyl1r; /* Split cylinders defined by radius, upper and */
+ real cyl0u, cyl1u; /* ... lower extension. The split cylinders de- */
+ real cyl0l, cyl1l; /* ... fine the channels and are each anchored */
+ /* ... in the center of the split group */
+ int nanions[eCompNR]; /* Requested number of anions and */
+ int nAverage; /* Coupling constant (nr of swap attempt steps) */
+ real threshold; /* Ion counts may deviate from the requested
+ values by +-threshold before a swap is done */
+ int ncations[eCompNR]; /* ... cations for both compartments */
+ real bulkOffset[eCompNR]; /* Offset of the swap layer (='bulk') w.r.t.
+ the compartment-defining layers */
+ gmx_swapcoords_t si_priv; /* swap private data accessible in
+ * swapcoords.cpp */
} t_swapcoords;
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/snprintf.h"
static const char *SwS = {"SWAP:"}; /**< For output that comes from the swap module */
static const char *SwSEmpty = {" "}; /**< Placeholder for multi-line output */
the compartment conditions. */
int *nat_past; /**< Past ion counts for time-averaging. */
int *ind; /**< Indices to coll array of atoms. */
- real *dist; /**< Distance of atom to compartment center. */
+ real *dist; /**< Distance of atom to bulk layer, which is
+ normally the center layer of the compartment */
int nalloc; /**< Allocation size for ind array. */
int inflow_netto; /**< Net inflow of ions into this compartment. */
} t_compartment;
* v d_down
*
* \endcode
+ *
+ * \param[in] point The position (xyz) under consideration.
+ * \param[in] center The center of the cylinder.
+ * \param[in] d_up The upper extension of the cylinder.
+ * \param[in] d_down The lower extension.
+ * \param[in] r_cyl2 Cylinder radius squared.
+ * \param[in] pbc Structure with info about periodic boundary conditions.
+ * \param[in] normal The membrane normal direction is typically 3, i.e. z, but can be x or y also.
+ *
+ * \returns Whether the point is inside the defined cylindric channel.
*/
static gmx_bool is_in_channel(
- rvec point, /* Point under consideration */
- rvec center, /* 'Center' of cylinder */
- real d_up, /* Upper extension */
- real d_down, /* Lower extensions */
- real r_cyl2, /* Cylinder radius squared */
+ rvec point,
+ rvec center,
+ real d_up,
+ real d_down,
+ real r_cyl2,
t_pbc *pbc,
- int normal) /* The membrane normal direction is typically 3, i.e. ZZ, but can be X or Y also */
+ int normal)
{
rvec dr;
int plane1, plane2; /* Directions tangential to membrane */
-/*! \brief Return TRUE if ion is found in the compartment.
+/*! \brief Return TRUE if position x of ion (or water) is found in the compartment,
+ * i.e. between w1 and w2.
*
- * Returns TRUE if x is between (w1+gap) and (w2-gap)
+ * One can define and additional offset "b" if one wants to exchange ions/water
+ * to or from a plane not directly in the middle of w1 and w2. The offset can be
+ * in ]-1.0, ..., +1.0 [.
+ * A bulkOffset of 0.0 means 'no offset', so the swap-layer is directly in the
+ * middle between w1 and w2. Offsets -1.0 < b < 0.0 will yield swaps nearer to w1,
+ * whereas offsets 0.0 < 0 < +1.0 will yield swaps nearer to w2.
*
* \code
*
- * ||-----------|--+--|----------o----------|--+--|---------------------||
- * w1 ????????????????????? w2
+ * ||--------------+-------------|-------------+------------------------||
+ * w1 ? ? ? ? ? ? ? ? ? ? ? w2
+ * ||--------------+-------------|----b--------+------------------------||
+ * -1 0.0 +1
*
* \endcode
+ *
+ * \param[in] w1 Position of 'wall' atom 1.
+ * \param[in] w2 Position of 'wall' atom 2.
+ * \param[in] x Position of the ion or the water molecule under consideration.
+ * \param[in] l Length of the box, from || to || in the sketch.
+ * \param[in] bulkOffset Where is the bulk layer "b" to be found between w1 and w2?
+ * \param[out] distance_from_b Distance of x to the bulk layer "b".
+ *
+ * \returns TRUE if x is between w1 and w2.
+ *
+ * Also computes the distance of x to the compartment center (the layer that is
+ * normally situated in the middle of w1 and w2 that would be considered as having
+ * the bulk concentration of ions).
*/
static gmx_bool compartment_contains_atom(
- real w1, /* position of wall atom 1 */
- real w2, /* position of wall atom 2 */
- real gap,
+ real w1,
+ real w2,
real x,
- real l, /* length of the box, from || to || in the sketch */
- real *distance_from_center)
+ real l,
+ real bulkOffset,
+ real *distance_from_b)
{
real m, l_2;
+ real width;
/* First set the origin in the middle of w1 and w2 */
- m = 0.5 * (w1 + w2);
- w1 -= m;
- w2 -= m;
- x -= m;
+ m = 0.5 * (w1 + w2);
+ w1 -= m;
+ w2 -= m;
+ x -= m;
+ width = w2 - w1;
/* Now choose the PBC image of x that is closest to the origin: */
l_2 = 0.5*l;
x += l;
}
- *distance_from_center = (real)fabs(x);
+ *distance_from_b = (real)fabs(x - bulkOffset*0.5*width);
/* Return TRUE if we now are in area "????" */
- if ( (x >= (w1+gap)) && (x < (w2-gap)) )
+ if ( (x >= w1) && (x < w2) )
{
return TRUE;
}
}
-/*! \brief Add atom with collective index ci to the list 'comp'. */
+/*! \brief Add the atom with collective index ci to the atom list in compartment 'comp'.
+ *
+ * \param[in] ci Index of this ion in the collective xc array.
+ * \param[inout] comp Compartment to add this atom to.
+ * \param[in] distance Shortest distance of this atom to the bulk layer,
+ * from which ion/water pairs are selected for swapping.
+ */
static void add_to_list(
- int ci, /* index of this ion in the collective array xc, qc */
- t_compartment *comp, /* Compartment to add this atom to */
- real distance) /* Shortest distance of this atom to the compartment center */
+ int ci,
+ t_compartment *comp,
+ real distance)
{
int nr;
if (c >= eCompNR)
{
- gmx_fatal(FARGS, "No compartment %d.", c);
+ gmx_fatal(FARGS, "No compartment %c.", c+'A');
}
pos0 = s->group[eGrpSplit0].center[s->swapdim];
type = iong->qc[i] < 0 ? eIonNEG : eIonPOS;
/* Is this ion in the compartment that we look at? */
- if (compartment_contains_atom(left, right, 0, iong->xc[i][sd], box[sd][sd], &dist) )
+ if (compartment_contains_atom(left, right, iong->xc[i][sd], box[sd][sd], sc->bulkOffset[comp], &dist) )
{
/* Now put it into the list containing only ions of its type */
add_to_list(i, &s->comp[comp][type], dist);
/* Loop over the solvent MOLECULES */
for (i = 0; i < sc->nat_sol; i += apm)
{
- if (compartment_contains_atom(left, right, 0, solg->xc[i][sd], box[sd][sd], &dist))
+ if (compartment_contains_atom(left, right, solg->xc[i][sd], box[sd][sd], sc->bulkOffset[comp], &dist))
{
/* Add the whole molecule to the list */
for (j = 0; j < apm; j++)
{
const char **legend;
int ic, count, ii;
- char buf[256];
+ char buf[STRLEN];
t_swap *s;
{
for (ii = 0; ii < eIonNR; ii++)
{
- sprintf(buf, "%s %ss", CompStr[ic], IonString[ii]);
+ snprintf(buf, STRLEN, "%s %ss", CompStr[ic], IonString[ii]);
legend[count++] = gmx_strdup(buf);
- sprintf(buf, "%s av. mismatch to %d%s",
- CompStr[ic], s->comp[ic][ii].nat_req, IonStr[ii]);
+ snprintf(buf, STRLEN, "%s av. mismatch to %d%s",
+ CompStr[ic], s->comp[ic][ii].nat_req, IonStr[ii]);
legend[count++] = gmx_strdup(buf);
- sprintf(buf, "%s netto %s influx", CompStr[ic], IonString[ii]);
+ snprintf(buf, STRLEN, "%s netto %s influx", CompStr[ic], IonString[ii]);
legend[count++] = gmx_strdup(buf);
}
}
- sprintf(buf, "%scenter of %s of split group 0", SwapStr[ir->eSwapCoords], (NULL != s->group[eGrpSplit0].m) ? "mass" : "geometry");
+ snprintf(buf, STRLEN, "%scenter of %s of split group 0", SwapStr[ir->eSwapCoords], (NULL != s->group[eGrpSplit0].m) ? "mass" : "geometry");
legend[count++] = gmx_strdup(buf);
- sprintf(buf, "%scenter of %s of split group 1", SwapStr[ir->eSwapCoords], (NULL != s->group[eGrpSplit1].m) ? "mass" : "geometry");
+ snprintf(buf, STRLEN, "%scenter of %s of split group 1", SwapStr[ir->eSwapCoords], (NULL != s->group[eGrpSplit1].m) ? "mass" : "geometry");
legend[count++] = gmx_strdup(buf);
for (ic = 0; ic < eChanNR; ic++)
{
for (ii = 0; ii < eIonNR; ii++)
{
- sprintf(buf, "A->ch%d->B %s permeations", ic, IonString[ii]);
+ snprintf(buf, STRLEN, "A->ch%d->B %s permeations", ic, IonString[ii]);
legend[count++] = gmx_strdup(buf);
}
}
- sprintf(buf, "leakage");
+ snprintf(buf, STRLEN, "leakage");
legend[count++] = gmx_strdup(buf);
xvgr_legend(s->fpout, count, legend, oenv);
if (!bAppend)
{
+ if ( (0 != sc->bulkOffset[eCompA]) || (0 != sc->bulkOffset[eCompB]) )
+ {
+ fprintf(s->fpout, "#\n");
+ fprintf(s->fpout, "# You provided an offset for the position of the bulk layer(s).\n");
+ fprintf(s->fpout, "# That means the layers to/from which ions and water molecules are swapped\n");
+ fprintf(s->fpout, "# are not midway (= at 0.0) between the compartment-defining layers (at +/- 1.0).\n");
+ fprintf(s->fpout, "# bulk-offsetA = %g\n", sc->bulkOffset[eCompA]);
+ fprintf(s->fpout, "# bulk-offsetB = %g\n", sc->bulkOffset[eCompB]);
+ }
+
fprintf(s->fpout, "#\n");
fprintf(s->fpout, "# split0 cylinder radius %f nm, up %f nm, down %f nm\n",
sc->cyl0r, sc->cyl0u, sc->cyl0l);
fprintf(s->fpout, "# split1 cylinder radius %f nm, up %f nm, down %f nm\n",
sc->cyl1r, sc->cyl1u, sc->cyl1l);
- }
- if (!bAppend)
- {
fprintf(s->fpout, "#\n");
if (!bRerun)
{
if (bVerbose)
{
- fprintf(stderr, "%s Requested charge imbalance is Q(A) - Q(B) = %gz.\n", SwS, s->deltaQ);
+ fprintf(stderr, "%s Requested charge imbalance is Q(A) - Q(B) = %g e.\n", SwS, s->deltaQ);
}
if (!bAppend)
{
- fprintf(s->fpout, "# Requested charge imbalance is Q(A)-Q(B) = %gz.\n", s->deltaQ);
+ fprintf(s->fpout, "# Requested charge imbalance is Q(A)-Q(B) = %g e.\n", s->deltaQ);
}
}
}
-/*! \brief Return index of atom that we can use for swapping.
+/*! \brief Return the index of an atom or molecule suitable for swapping.
+ *
+ * Returns the index of an atom that is far off the compartment boundaries,
+ * that is near to the bulk layer to/from which the swaps take place.
+ * Other atoms of the molecule (if any) will directly follow the returned index.
+ *
+ * \param[in] comp Structure containing compartment-specific data.
+ * \param[in] apm Atoms per molecule - just return the first atom index of a molecule
*
- * Returns the index of an atom that is far off the compartment boundaries.
- * Other atoms of the molecule (if any) will directly follow the returned index
+ * \returns Index of the first atom of the molecule chosen for a position exchange.
*/
static int get_index_of_distant_atom(
t_compartment *comp,
- int apm) /* Atoms per molecule - just return the first atom index of a molecule */
+ int apm)
{
int i, ibest = -1;
real d = GMX_REAL_MAX;
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff