2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * Implements functions in swapcoords.h.
39 * \author Carsten Kutzner <ckutzne@gwdg.de>
40 * \ingroup module_swap
44 #include "swapcoords.h"
51 #include "gromacs/fileio/confio.h"
52 #include "gromacs/fileio/gmxfio.h"
53 #include "gromacs/fileio/xvgr.h"
54 #include "gromacs/legacyheaders/copyrite.h"
55 #include "gromacs/legacyheaders/macros.h"
56 #include "gromacs/legacyheaders/mdrun.h"
57 #include "gromacs/legacyheaders/names.h"
58 #include "gromacs/legacyheaders/network.h"
59 #include "gromacs/legacyheaders/typedefs.h"
60 #include "gromacs/math/vec.h"
61 #include "gromacs/mdlib/groupcoord.h"
62 #include "gromacs/pbcutil/pbc.h"
63 #include "gromacs/timing/wallcycle.h"
64 #include "gromacs/topology/mtop_util.h"
65 #include "gromacs/utility/cstringutil.h"
66 #include "gromacs/utility/smalloc.h"
68 static const char *SwS = {"SWAP:"}; /**< For output that comes from the swap module */
69 static const char *SwSEmpty = {" "}; /**< Placeholder for multi-line output */
70 static const char* IonString[eIonNR] = {"anion", "cation" }; /**< Type of ion, used for verbose output */
71 static const char* IonStr[eIonNR] = {"-", "+" }; /**< Type of ion, used for short output */
72 static const char* CompStr[eCompNR] = {"A", "B" }; /**< Compartment name */
73 static const char *SwapStr[eSwapTypesNR+1] = { "", "X-", "Y-", "Z-", NULL}; /**< Name for the swap types. */
74 static const char *DimStr[DIM+1] = { "X", "Y", "Z", NULL}; /**< Name for the swap dimension. */
76 /* eGrpSplit0 and eGrpSplit1 _must_ be neighbors in this list because
77 * we sometimes loop from eGrpSplit0 to eGrpSplit1 */
79 eGrpIons, eGrpSplit0, eGrpSplit1, eGrpSolvent, eGrpNr
80 }; /**< Group identifier */
81 static const char* GrpString[eGrpNr] = { "ion", "split0", "split1", "solvent" }; /**< Group name */
83 /** Keep track of through which channel the ions have passed */
84 enum eChannelHistory {
85 eChHistPassedNone, eChHistPassedCh0, eChHistPassedCh1, eChHistNr
87 static const char* ChannelString[eChHistNr] = { "none", "channel0", "channel1" }; /**< Name for the channels */
89 /*! \brief Domain identifier.
91 * Keeps track of from which compartment the ions came before passing the
95 eDomainNotset, eDomainA, eDomainB, eDomainNr
97 static const char* DomainString[eDomainNr] = { "not_assigned", "Domain_A", "Domain_B" }; /**< Name for the domains */
102 * Structure containing compartment-specific data.
104 typedef struct swap_compartment
106 int nat; /**< Number of atoms matching the
107 compartment conditions. */
108 int nat_old; /**< Number of atoms before swapping. */
109 int nat_req; /**< Requested number of atoms. */
110 real nat_av; /**< Time-averaged number of atoms matching
111 the compartment conditions. */
112 int *nat_past; /**< Past ion counts for time-averaging. */
113 int *ind; /**< Indices to coll array of atoms. */
114 real *dist; /**< Distance of atom to compartment center. */
115 int nalloc; /**< Allocation size for ind array. */
116 int inflow_netto; /**< Net inflow of ions into this compartment. */
121 * This structure contains data needed for each of the groups involved in swapping: ions, water,
124 typedef struct swap_group
126 int nat; /**< Number of atoms in the group */
127 int apm; /**< Number of atoms in each molecule */
128 atom_id *ind; /**< Global atom indices of the group */
129 atom_id *ind_loc; /**< Local atom indices of the group */
130 int nat_loc; /**< Number of local group atoms */
131 int nalloc_loc; /**< Allocation size for ind_loc */
132 rvec *xc; /**< Collective array of group atom positions */
133 ivec *xc_shifts; /**< Current (collective) shifts */
134 ivec *xc_eshifts; /**< Extra shifts since last DD step */
135 rvec *xc_old; /**< Old (collective) positions */
136 real *qc; /**< Collective array of charges */
137 int *c_ind_loc; /**< Position of local atoms in the
138 collective array, [0..nat_loc] */
139 real *m; /**< Masses (can be omitted) */
140 unsigned char *comp_from; /**< (Collective) Stores from which compartment this
141 atom has come. This way we keep track of through
142 which channel an ion permeates (only used for
144 unsigned char *comp_now; /**< In which compartment this ion is now */
145 unsigned char *channel_label; /**< Which channel was passed at last by this ion? */
146 rvec center; /**< Center of the group; COM if masses are used */
151 * Main (private) data structure for the position swapping protocol.
153 typedef struct t_swap
155 int swapdim; /**< One of XX, YY, ZZ */
156 t_pbc *pbc; /**< Needed to make molecules whole. */
157 FILE *fpout; /**< Output file. */
158 t_group group[eGrpNr]; /**< Ions, solvent or channels? */
159 t_compartment comp[eCompNR][eIonNR]; /**< Data for a specific compartment and ion type. */
160 t_compartment compsol[eCompNR]; /**< Solvent compartments. */
161 int fluxfromAtoB[eChanNR][eIonNR]; /**< Net flux per channels and ion type. */
162 int ncyl0ions; /**< Number of ions residing in channel 0. */
163 int ncyl1ions; /**< Same for channel 1. */
164 int cyl0and1; /**< Ions assigned to cyl0 and cyl1. Not good. */
165 int *fluxleak; /**< Pointer to a single int value holding the
166 flux not going through any of the channels. */
167 real deltaQ; /**< The charge imbalance between the compartments. */
172 /*! \brief Check whether point is in channel.
174 * A channel is a cylinder defined by a disc
175 * with radius r around its center c. The thickness of the cylinder is
182 * <---------c--------->
188 static gmx_bool is_in_channel(
189 rvec point, /* Point under consideration */
190 rvec center, /* 'Center' of cylinder */
191 real d_up, /* Upper extension */
192 real d_down, /* Lower extensions */
193 real r_cyl2, /* Cylinder radius squared */
195 int normal) /* The membrane normal direction is typically 3, i.e. ZZ, but can be X or Y also */
198 int plane1, plane2; /* Directions tangential to membrane */
201 plane1 = (normal + 1) % 3; /* typically 0, i.e. XX */
202 plane2 = (normal + 2) % 3; /* typically 1, i.e. YY */
204 /* Get the distance vector dr between the point and the center of the cylinder */
205 pbc_dx(pbc, point, center, dr); /* This puts center in the origin */
207 /* Check vertical direction */
208 if ( (dr[normal] > d_up) || (dr[normal] < -d_down) )
213 /* Check radial direction */
214 if ( (dr[plane1]*dr[plane1] + dr[plane2]*dr[plane2]) > r_cyl2)
219 /* All check passed, this point is in the cylinder */
224 /*! \brief Prints to swap output file which ions are in which compartment. */
225 static void print_ionlist(
228 const char comment[])
230 int itype, icomp, i, j;
234 fprintf(s->fpout, "%12.5e", time);
235 for (icomp = 0; icomp < eCompNR; icomp++)
237 for (itype = 0; itype < eIonNR; itype++)
239 comp = &(s->comp[icomp][itype]);
240 fprintf(s->fpout, "%7d%7.1f%7d", comp->nat, comp->nat_av-comp->nat_req, comp->inflow_netto);
243 fprintf(s->fpout, "%12.3e%12.3e",
244 s->group[eGrpSplit0].center[s->swapdim],
245 s->group[eGrpSplit1].center[s->swapdim]);
247 for (i = 0; i < eChanNR; i++)
249 for (j = 0; j < eIonNR; j++)
251 fprintf(s->fpout, "%12d", s->fluxfromAtoB[i][j]);
255 /* Also print the number of ions that leaked from A to B: */
256 fprintf(s->fpout, "%12d", *s->fluxleak);
258 fprintf(s->fpout, "%s\n", comment);
262 /*! \brief Get the center of a group of nat atoms.
264 * Since with PBC an atom group might not be whole, use the first atom as the
265 * reference atom and determine the center with respect to this reference.
267 static void get_molecule_center(
275 rvec weightedPBCimage;
277 rvec reference, correctPBCimage, dx;
280 /* Use the first atom as the reference and put other atoms near that one */
281 /* This does not work for large molecules that span > half of the box! */
282 copy_rvec(x[0], reference);
284 /* Calculate either the weighted center or simply the center of geometry */
287 for (i = 0; i < nat; i++)
289 /* PBC distance between position and reference */
290 pbc_dx(pbc, x[i], reference, dx);
292 /* Add PBC distance to reference */
293 rvec_add(reference, dx, correctPBCimage);
295 /* Take weight into account */
305 svmul(wi, correctPBCimage, weightedPBCimage);
308 rvec_inc(center, weightedPBCimage);
312 svmul(1.0/wsum, center, center);
317 /*! \brief Return TRUE if ion is found in the compartment.
319 * Returns TRUE if x is between (w1+gap) and (w2-gap)
323 * ||-----------|--+--|----------o----------|--+--|---------------------||
324 * w1 ????????????????????? w2
328 static gmx_bool compartment_contains_atom(
329 real w1, /* position of wall atom 1 */
330 real w2, /* position of wall atom 2 */
333 real l, /* length of the box, from || to || in the sketch */
334 real *distance_from_center)
339 /* First set the origin in the middle of w1 and w2 */
345 /* Now choose the PBC image of x that is closest to the origin: */
356 *distance_from_center = (real)fabs(x);
358 /* Return TRUE if we now are in area "????" */
359 if ( (x >= (w1+gap)) && (x < (w2-gap)) )
370 /*! \brief Updates the time-averaged number of ions in a compartment. */
371 static void update_time_window(t_compartment *comp, int values, int replace)
377 /* Put in the new value */
380 comp->nat_past[replace] = comp->nat;
383 /* Compute the new time-average */
385 for (i = 0; i < values; i++)
387 average += comp->nat_past[i];
390 comp->nat_av = average;
394 /*! \brief Add atom with collective index ci to the list 'comp'. */
395 static void add_to_list(
396 int ci, /* index of this ion in the collective array xc, qc */
397 t_compartment *comp, /* Compartment to add this atom to */
398 real distance) /* Shortest distance of this atom to the compartment center */
405 if (nr >= comp->nalloc)
407 comp->nalloc = over_alloc_dd(nr+1);
408 srenew(comp->ind, comp->nalloc);
409 srenew(comp->dist, comp->nalloc);
412 comp->dist[nr] = distance;
417 /*! \brief Determine the compartment boundaries from the channel centers. */
418 static void get_compartment_boundaries(
422 real *left, real *right)
425 real leftpos, rightpos, leftpos_orig;
430 gmx_fatal(FARGS, "No compartment %d.", c);
433 pos0 = s->group[eGrpSplit0].center[s->swapdim];
434 pos1 = s->group[eGrpSplit1].center[s->swapdim];
447 /* This gets us the other compartment: */
450 leftpos_orig = leftpos;
452 rightpos = leftpos_orig + box[s->swapdim][s->swapdim];
460 /*! \brief Determine the per-channel ion flux.
462 * To determine the flux through the individual channels, we
463 * remember the compartment and channel history of each ion. An ion can be
464 * either in channel0 or channel1, or in the remaining volume of compartment
468 * +-----------------+
471 * ||||||||||0|||||||| bilayer with channel 0
476 * |||||1||||||||||||| bilayer with channel 1
479 * +-----------------+
483 static void detect_flux_per_channel(
488 unsigned char *comp_now,
489 unsigned char *comp_from,
490 unsigned char *channel_label,
500 gmx_bool in_cyl0, in_cyl1;
507 /* Check whether ion is inside any of the channels */
508 in_cyl0 = is_in_channel(ion_pos, s->group[eGrpSplit0].center, sc->cyl0u, sc->cyl0l, cyl0_r2, s->pbc, sd);
509 in_cyl1 = is_in_channel(ion_pos, s->group[eGrpSplit1].center, sc->cyl1u, sc->cyl1l, cyl1_r2, s->pbc, sd);
511 if (in_cyl0 && in_cyl1)
513 /* Ion appears to be in both channels. Something is severely wrong! */
515 *comp_now = eDomainNotset;
516 *comp_from = eDomainNotset;
517 *channel_label = eChHistPassedNone;
521 /* Ion is in channel 0 now */
522 *channel_label = eChHistPassedCh0;
523 *comp_now = eDomainNotset;
528 /* Ion is in channel 1 now */
529 *channel_label = eChHistPassedCh1;
530 *comp_now = eDomainNotset;
535 /* Ion is not in any of the channels, so it must be in domain A or B */
538 *comp_now = eDomainA;
542 *comp_now = eDomainB;
546 /* Only take action, if ion is now in domain A or B, and was before
547 * in the other domain!
549 if (eDomainNotset == *comp_from)
551 /* Maybe we can set the domain now */
552 *comp_from = *comp_now; /* Could still be eDomainNotset, though */
554 else if ( (*comp_now != eDomainNotset ) /* if in channel */
555 && (*comp_from != *comp_now) )
557 /* Obviously the ion changed its domain.
558 * Count this for the channel through which it has passed. */
559 switch (*channel_label)
561 case eChHistPassedNone:
562 *s->fluxleak = *s->fluxleak + 1;
564 fprintf(stderr, " %s Warning! Step %s, ion %d (%s) moved from %s to %s\n",
565 SwS, gmx_step_str(step, buf), iion, IonStr[iontype], DomainString[*comp_from], DomainString[*comp_now]);
568 fprintf(stderr, ", possibly due to a swap in the original simulation.\n");
572 fprintf(stderr, "but did not pass cyl0 or cyl1 as defined in the .mdp file.\n"
573 "Do you have an ion somewhere within the membrane?\n");
574 /* Write this info to the CompEL output file: */
575 fprintf(s->fpout, " # Warning: step %s, ion %d (%s) moved from %s to %s (probably through the membrane)\n",
576 gmx_step_str(step, buf), iion, IonStr[iontype],
577 DomainString[*comp_from], DomainString[*comp_now]);
581 case eChHistPassedCh0:
582 case eChHistPassedCh1:
583 if (*channel_label == eChHistPassedCh0)
592 if (eDomainA == *comp_from)
594 s->fluxfromAtoB[chan_nr][iontype]++;
598 s->fluxfromAtoB[chan_nr][iontype]--;
600 fprintf(fpout, "# Atom nr. %d finished passing %s.\n", iion, ChannelString[*channel_label]);
603 gmx_fatal(FARGS, "%s Unknown channel history entry!\n", SwS);
607 /* This ion has moved to the _other_ compartment ... */
608 *comp_from = *comp_now;
609 /* ... and it did not pass any channel yet */
610 *channel_label = eChHistPassedNone;
615 /*! \brief Get the lists of ions for the two compartments */
616 static void compartmentalize_ions(
629 real cyl0_r2, cyl1_r2;
631 int sum, not_in_comp[eCompNR]; /* consistency check */
636 iong = &s->group[eGrpIons];
639 cyl0_r2 = sc->cyl0r * sc->cyl0r;
640 cyl1_r2 = sc->cyl1r * sc->cyl1r;
643 /* Get us a counter that cycles in the range of [0 ... sc->nAverage[ */
644 replace = (step/sc->nstswap) % sc->nAverage;
646 for (comp = eCompA; comp <= eCompB; comp++)
648 /* Get lists of atoms that match criteria for this compartment */
649 get_compartment_boundaries(comp, sc->si_priv, box, &left, &right);
651 /* First clear the ion lists */
652 s->comp[comp][eIonNEG].nat = 0;
653 s->comp[comp][eIonPOS].nat = 0;
654 not_in_comp[comp] = 0; /* consistency check */
656 /* Loop over the IONS */
657 for (i = 0; i < iong->nat; i++)
659 /* Anion or cation? */
660 type = iong->qc[i] < 0 ? eIonNEG : eIonPOS;
662 /* Is this ion in the compartment that we look at? */
663 if (compartment_contains_atom(left, right, 0, iong->xc[i][sd], box[sd][sd], &dist) )
665 /* Now put it into the list containing only ions of its type */
666 add_to_list(i, &s->comp[comp][type], dist);
668 /* Master also checks through which channel each ion has passed */
669 if (MASTER(cr) && (iong->comp_now != NULL))
671 ion_nr_global = iong->ind[i] + 1; /* PDB index starts at 1 ... */
672 detect_flux_per_channel(ion_nr_global, comp, type, iong->xc[i],
673 &iong->comp_now[i], &iong->comp_from[i], &iong->channel_label[i],
674 sc, cyl0_r2, cyl1_r2, step, bRerun, fpout);
679 not_in_comp[comp] += 1;
682 /* Correct the time-averaged number of ions in both compartments */
683 update_time_window(&s->comp[comp][eIonNEG], sc->nAverage, replace);
684 update_time_window(&s->comp[comp][eIonPOS], sc->nAverage, replace);
687 /* Flux detection warnings */
693 "%s Warning: %d atoms were detected as being in both channels! Probably your split\n"
694 "%s cylinder is way too large, or one compartment has collapsed (step %" GMX_PRId64 ")\n",
695 SwS, s->cyl0and1, SwS, step);
697 fprintf(s->fpout, "Warning: %d atoms were assigned to both channels!\n", s->cyl0and1);
704 /* Consistency checks */
705 if (not_in_comp[eCompA] + not_in_comp[eCompB] != iong->nat)
709 fprintf(fpout, "# Warning: Inconsistency during ion compartmentalization. !inA: %d, !inB: %d, total ions %d\n",
710 not_in_comp[eCompA], not_in_comp[eCompB], iong->nat);
715 fprintf(stderr, "%s rank %d: Inconsistency during ion compartmentalization. !inA: %d, !inB: %d, total ions %d\n",
716 SwS, cr->nodeid, not_in_comp[eCompA], not_in_comp[eCompB], iong->nat);
720 sum = s->comp[eCompA][eIonNEG].nat + s->comp[eCompA][eIonPOS].nat
721 + s->comp[eCompB][eIonNEG].nat + s->comp[eCompB][eIonPOS].nat;
722 if (sum != iong->nat)
726 fprintf(fpout, "# Warning: %d atoms are in the ion group, but altogether %d have been assigned to the compartments.\n",
732 fprintf(stderr, "%s rank %d: %d atoms are in the ion group, but altogether %d have been assigned to the compartments.\n",
733 SwS, cr->nodeid, iong->nat, sum);
741 /*! \brief Set up the compartments and get lists of solvent atoms in each compartment */
742 static void compartmentalize_solvent(
754 int sum, not_in_comp[eCompNR]; /* consistency check */
758 solg = &s->group[eGrpSolvent];
762 for (comp = eCompA; comp <= eCompB; comp++)
764 /* Get lists of atoms that match criteria for this compartment */
765 get_compartment_boundaries(comp, sc->si_priv, box, &left, &right);
767 /* First clear the solvent molecule lists */
768 s->compsol[comp].nat = 0;
769 not_in_comp[comp] = 0; /* consistency check */
771 /* Loop over the solvent MOLECULES */
772 for (i = 0; i < sc->nat_sol; i += apm)
774 if (compartment_contains_atom(left, right, 0, solg->xc[i][sd], box[sd][sd], &dist))
776 /* Add the whole molecule to the list */
777 for (j = 0; j < apm; j++)
779 add_to_list(i+j, &s->compsol[comp], dist);
784 not_in_comp[comp] += apm;
791 fprintf(fpout, "# Solv. molecules in comp.%s: %d comp.%s: %d\n",
792 CompStr[eCompA], s->compsol[eCompA].nat/apm,
793 CompStr[eCompB], s->compsol[eCompB].nat/apm);
796 /* Consistency checks */
797 if (not_in_comp[eCompA] + not_in_comp[eCompB] != solg->nat)
801 fprintf(fpout, "# Warning: Inconsistency during solvent compartmentalization. !inA: %d, !inB: %d, solvent atoms %d\n",
802 not_in_comp[eCompA], not_in_comp[eCompB], solg->nat);
807 fprintf(stderr, "%s rank %d: Inconsistency during solvent compartmentalization. !inA: %d, !inB: %d, solvent atoms %d\n",
808 SwS, cr->nodeid, not_in_comp[eCompA], not_in_comp[eCompB], solg->nat);
811 sum = s->compsol[eCompA].nat + s->compsol[eCompB].nat;
812 if (sum != solg->nat)
816 fprintf(fpout, "# Warning: %d atoms in solvent group, but %d have been assigned to the compartments.\n",
822 fprintf(stderr, "%s rank %d: %d atoms in solvent group, but %d have been assigned to the compartments.\n",
823 SwS, cr->nodeid, solg->nat, sum);
829 /*! \brief Find out how many group atoms are in the compartments initially */
830 static void get_initial_ioncounts(
832 rvec x[], /* the initial positions */
846 /* Copy the initial swap group positions to the collective array so
847 * that we can compartmentalize */
848 for (i = 0; i < sc->nat; i++)
851 copy_rvec(x[ind], s->group[eGrpIons].xc[i]);
854 /* Set up the compartments and get lists of atoms in each compartment */
855 compartmentalize_ions(cr, sc, box, 0, s->fpout, bRerun);
857 /* Set initial concentrations if requested */
858 for (ic = 0; ic < eCompNR; ic++)
860 s->comp[ic][eIonPOS].nat_req = sc->ncations[ic];
861 s->comp[ic][eIonNEG].nat_req = sc->nanions[ic];
863 for (ic = 0; ic < eCompNR; ic++)
865 for (ii = 0; ii < eIonNR; ii++)
867 if (s->comp[ic][ii].nat_req < 0)
869 s->comp[ic][ii].nat_req = s->comp[ic][ii].nat;
874 /* Check whether the number of requested ions adds up to the total number of ions */
875 for (ii = 0; ii < eIonNR; ii++)
877 req[ii] = s->comp[eCompA][ii].nat_req + s->comp[eCompB][ii].nat_req;
878 tot[ii] = s->comp[eCompA][ii].nat + s->comp[eCompB][ii].nat;
880 if ( (req[eCompA] != tot[eCompA]) || (req[eCompB] != tot[eCompB ]) )
882 gmx_fatal(FARGS, "Mismatch of the number of ions summed over both compartments.\n"
883 "You requested a total of %d anions and %d cations,\n"
884 "but there are a total of %d anions and %d cations in the system.\n",
885 req[eIonNEG], req[eIonPOS],
886 tot[eIonNEG], tot[eIonPOS]);
889 /* Initialize time-averaging:
890 * Write initial concentrations to all time bins to start with */
891 for (ic = 0; ic < eCompNR; ic++)
893 for (ii = 0; ii < eIonNR; ii++)
895 s->comp[ic][ii].nat_av = s->comp[ic][ii].nat;
896 for (i = 0; i < sc->nAverage; i++)
898 s->comp[ic][ii].nat_past[i] = s->comp[ic][ii].nat;
905 /*! \brief Copy history of ion counts from checkpoint file.
907 * When called, the checkpoint file has already been read in. Here we copy
908 * over the values from .cpt file to the swap data structure.
910 static void get_initial_ioncounts_from_cpt(
911 t_inputrec *ir, swapstate_t *swapstate,
912 t_commrec *cr, gmx_bool bVerbose)
923 /* Copy the past values from the checkpoint values that have been read in already */
926 fprintf(stderr, "%s Copying values from checkpoint\n", SwS);
929 for (ic = 0; ic < eCompNR; ic++)
931 for (ii = 0; ii < eIonNR; ii++)
933 s->comp[ic][ii].nat_req = swapstate->nat_req[ic][ii];
934 s->comp[ic][ii].inflow_netto = swapstate->inflow_netto[ic][ii];
938 fprintf(stderr, "%s ... Influx netto: %d Requested: %d Past values: ", SwS,
939 s->comp[ic][ii].inflow_netto, s->comp[ic][ii].nat_req);
942 for (j = 0; j < sc->nAverage; j++)
944 s->comp[ic][ii].nat_past[j] = swapstate->nat_past[ic][ii][j];
947 fprintf(stderr, "%d ", s->comp[ic][ii].nat_past[j]);
952 fprintf(stderr, "\n");
960 /*! \brief The master lets all others know about the initial ion counts. */
961 static void bc_initial_concentrations(
970 for (ic = 0; ic < eCompNR; ic++)
972 for (ii = 0; ii < eIonNR; ii++)
974 gmx_bcast(sizeof(s->comp[ic][ii].nat_req), &(s->comp[ic][ii].nat_req), cr);
975 gmx_bcast(sizeof(s->comp[ic][ii].nat ), &(s->comp[ic][ii].nat ), cr);
976 gmx_bcast( swap->nAverage * sizeof(s->comp[ic][ii].nat_past[0]), s->comp[ic][ii].nat_past, cr);
982 /*! \brief Ensure that each atom belongs to at most one of the swap groups. */
983 static void check_swap_groups(t_swap *s, int nat, gmx_bool bVerbose)
987 atom_id *nGroup = NULL; /* This array counts for each atom in the MD system to
988 how many swap groups it belongs (should be 0 or 1!) */
990 int nMultiple = 0; /* Number of atoms belonging to multiple groups */
995 fprintf(stderr, "%s Making sure each atom belongs to at most one of the swap groups.\n", SwS);
998 /* Add one to the group count of atoms belonging to a swap group: */
1000 for (i = 0; i < eGrpNr; i++)
1003 for (j = 0; j < g->nat; j++)
1005 /* Get the global index of this atom of this group: */
1010 /* Make sure each atom belongs to at most one swap group: */
1011 for (j = 0; j < g->nat; j++)
1022 gmx_fatal(FARGS, "%s Cannot perform swapping since %d atom%s allocated to more than one swap index group.\n"
1023 "%s Each atom must be allocated to at most one of the split groups, the swap group, or the solvent.\n"
1024 "%s Check the .mdp file settings regarding the swap index groups or the index groups themselves.\n",
1025 SwS, nMultiple, (1 == nMultiple) ? " is" : "s are", SwSEmpty, SwSEmpty);
1030 /*! \brief Get the number of atoms per molecule for this group.
1032 * Also ensure that all the molecules in this group have this number of atoms.
1034 static int get_group_apm_check(
1038 const gmx_mtop_atomlookup_t alook,
1043 int molb, molnr, atnr_mol;
1046 ind = s->group[group].ind;
1047 nat = s->group[group].nat;
1049 /* Determine the number of solvent atoms per solvent molecule from the
1050 * first solvent atom: */
1052 gmx_mtop_atomnr_to_molblock_ind(alook, ind[i], &molb, &molnr, &atnr_mol);
1053 apm = mtop->molblock[molb].natoms_mol;
1057 fprintf(stderr, "%s Checking whether all %s molecules consist of %d atom%s\n",
1058 SwS, GrpString[group], apm, apm > 1 ? "s" : "");
1061 /* Check whether this is also true for all other solvent atoms */
1062 for (i = 1; i < nat; i++)
1064 gmx_mtop_atomnr_to_molblock_ind(alook, ind[i], &molb, &molnr, &atnr_mol);
1065 if (apm != mtop->molblock[molb].natoms_mol)
1067 gmx_fatal(FARGS, "Not all %s group molecules consist of %d atoms.",
1068 GrpString[group], apm);
1076 /*! \brief Print the legend to the swap output file.
1078 * Also print the initial ion counts
1080 static void print_ionlist_legend(t_inputrec *ir, const output_env_t oenv)
1082 const char **legend;
1088 s = ir->swap->si_priv;
1090 snew(legend, eCompNR*eIonNR*3 + 2 + eChanNR*eIonNR + 1);
1091 for (ic = count = 0; ic < eCompNR; ic++)
1093 for (ii = 0; ii < eIonNR; ii++)
1095 sprintf(buf, "%s %ss", CompStr[ic], IonString[ii]);
1096 legend[count++] = gmx_strdup(buf);
1097 sprintf(buf, "%s av. mismatch to %d%s",
1098 CompStr[ic], s->comp[ic][ii].nat_req, IonStr[ii]);
1099 legend[count++] = gmx_strdup(buf);
1100 sprintf(buf, "%s netto %s influx", CompStr[ic], IonString[ii]);
1101 legend[count++] = gmx_strdup(buf);
1104 sprintf(buf, "%scenter of %s of split group 0", SwapStr[ir->eSwapCoords], (NULL != s->group[eGrpSplit0].m) ? "mass" : "geometry");
1105 legend[count++] = gmx_strdup(buf);
1106 sprintf(buf, "%scenter of %s of split group 1", SwapStr[ir->eSwapCoords], (NULL != s->group[eGrpSplit1].m) ? "mass" : "geometry");
1107 legend[count++] = gmx_strdup(buf);
1109 for (ic = 0; ic < eChanNR; ic++)
1111 for (ii = 0; ii < eIonNR; ii++)
1113 sprintf(buf, "A->ch%d->B %s permeations", ic, IonString[ii]);
1114 legend[count++] = gmx_strdup(buf);
1118 sprintf(buf, "leakage");
1119 legend[count++] = gmx_strdup(buf);
1121 xvgr_legend(s->fpout, count, legend, oenv);
1123 fprintf(s->fpout, "# Instantaneous ion counts and time-averaged differences to requested numbers\n");
1124 fprintf(s->fpout, "# time[ps] A_an diff t_in A_cat diff t_in B_an diff t_in B_cat diff t_in ");
1125 fprintf(s->fpout, " %s-Split0 %s-Split1", DimStr[s->swapdim], DimStr[s->swapdim]);
1126 fprintf(s->fpout, " A-ch0-B_an A-ch0-B_cat A-ch1-B_an A-ch1-B_cat ion_leakage\n");
1131 /*! \brief Initialize channel ion flux detection routine.
1133 * Initialize arrays that keep track of where the ions come from and where
1136 static void detect_flux_per_channel_init(
1139 swapstate_t *swapstate,
1140 gmx_bool bStartFromCpt)
1146 g = &(s->group[eGrpIons]);
1148 /* All these flux detection routines run on the master only */
1152 g->comp_from = NULL;
1153 g->channel_label = NULL;
1158 /******************************************************/
1159 /* Channel and domain history for the individual ions */
1160 /******************************************************/
1161 if (bStartFromCpt) /* set the pointers right */
1163 g->comp_from = swapstate->comp_from;
1164 g->channel_label = swapstate->channel_label;
1166 else /* allocate memory */
1168 snew(g->comp_from, g->nat);
1169 swapstate->comp_from = g->comp_from;
1170 snew(g->channel_label, g->nat);
1171 swapstate->channel_label = g->channel_label;
1173 snew(g->comp_now, g->nat);
1175 /* Initialize the channel and domain history counters */
1176 for (i = 0; i < g->nat; i++)
1178 g->comp_now[i] = eDomainNotset;
1181 g->comp_from[i] = eDomainNotset;
1182 g->channel_label[i] = eChHistPassedNone;
1186 /************************************/
1187 /* Channel fluxes for both channels */
1188 /************************************/
1195 fprintf(stderr, "%s Copying channel fluxes from checkpoint file data\n", SwS);
1198 for (ic = 0; ic < eChanNR; ic++)
1200 fprintf(stderr, "%s Channel %d flux history: ", SwS, ic);
1201 for (ii = 0; ii < eIonNR; ii++)
1205 s->fluxfromAtoB[ic][ii] = swapstate->fluxfromAtoB[ic][ii];
1209 s->fluxfromAtoB[ic][ii] = 0;
1212 fprintf(stderr, "%d %s%s ", s->fluxfromAtoB[ic][ii], IonString[ii], s->fluxfromAtoB[ic][ii] == 1 ? "" : "s");
1214 fprintf(stderr, "\n");
1218 s->fluxleak = swapstate->fluxleak;
1222 snew(s->fluxleak, 1);
1224 /* Set pointer for checkpoint writing */
1225 swapstate->fluxleak = s->fluxleak;
1228 /* Set pointers for checkpoint writing */
1229 for (ic = 0; ic < eChanNR; ic++)
1231 for (ii = 0; ii < eIonNR; ii++)
1233 swapstate->fluxfromAtoB_p[ic][ii] = &(s->fluxfromAtoB[ic][ii]);
1239 /*! \brief Outputs the initial structure to PDB file for debugging reasons.
1241 * Output the starting structure so that in case of multimeric channels
1242 * the user can check whether we have the correct PBC image for all atoms.
1243 * If this is not correct, the ion counts per channel will be very likely
1246 static void outputStartStructureIfWanted(gmx_mtop_t *mtop, rvec *x, int ePBC, matrix box)
1248 char *env = getenv("GMX_COMPELDUMP");
1252 fprintf(stderr, "\n%s Found env.var. GMX_COMPELDUMP, will output CompEL starting structure made whole.\n"
1253 "%s In case of multimeric channels, please check whether they have the correct PBC representation.\n",
1256 write_sto_conf_mtop("CompELAssumedWholeConfiguration.pdb", *mtop->name, mtop, x, NULL, ePBC, box);
1261 /*! \brief Initialize the swapstate structure, used for checkpoint writing.
1263 * The swapstate struct stores the information we need to make the channels
1264 * whole again after restarts from a checkpoint file. Here we do the following:\n
1265 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,\n
1266 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,\n
1267 * c) in any case, for subsequent checkpoint writing, we set the pointers in\n
1268 * swapstate to the x_old arrays, which contain the correct PBC representation of
1269 * multimeric channels at the last time step.
1271 static void init_swapstate(
1272 swapstate_t *swapstate,
1275 rvec x[], /* the initial positions */
1280 rvec *x_pbc = NULL; /* positions of the whole MD system with molecules made whole */
1287 /* We always need the last whole positions such that
1288 * in the next time step we can make the channels whole again in PBC */
1289 if (swapstate->bFromCpt)
1291 /* Copy the last whole positions of each channel from .cpt */
1292 g = &(s->group[eGrpSplit0]);
1293 for (i = 0; i < g->nat; i++)
1295 copy_rvec(swapstate->xc_old_whole[eChan0][i], g->xc_old[i]);
1297 g = &(s->group[eGrpSplit1]);
1298 for (i = 0; i < g->nat; i++)
1300 copy_rvec(swapstate->xc_old_whole[eChan1][i], g->xc_old[i]);
1305 /* Extract the initial split group positions. */
1307 /* Remove pbc, make molecule whole. */
1308 snew(x_pbc, mtop->natoms);
1309 m_rveccopy(mtop->natoms, x, x_pbc);
1311 /* This can only make individual molecules whole, not multimers */
1312 do_pbc_mtop(NULL, ePBC, box, mtop, x_pbc);
1314 /* Output the starting structure? */
1315 outputStartStructureIfWanted(mtop, x_pbc, ePBC, box);
1317 /* If this is the first run (i.e. no checkpoint present) we assume
1318 * that the starting positions give us the correct PBC representation */
1319 for (ig = eGrpSplit0; ig <= eGrpSplit1; ig++)
1321 g = &(s->group[ig]);
1322 for (i = 0; i < g->nat; i++)
1324 copy_rvec(x_pbc[g->ind[i]], g->xc_old[i]);
1329 /* Prepare swapstate arrays for later checkpoint writing */
1330 swapstate->nat[eChan0] = s->group[eGrpSplit0].nat;
1331 swapstate->nat[eChan1] = s->group[eGrpSplit1].nat;
1334 /* For subsequent checkpoint writing, set the swapstate pointers to the xc_old
1335 * arrays that get updated at every swapping step */
1336 swapstate->xc_old_whole_p[eChan0] = &s->group[eGrpSplit0].xc_old;
1337 swapstate->xc_old_whole_p[eChan1] = &s->group[eGrpSplit1].xc_old;
1341 extern void init_swapcoords(
1349 swapstate_t *swapstate,
1351 const output_env_t oenv,
1352 unsigned long Flags)
1354 int i, ic, ig, ii, j;
1359 gmx_bool bAppend, bStartFromCpt, bRerun;
1360 gmx_mtop_atomlookup_t alook = NULL;
1363 alook = gmx_mtop_atomlookup_init(mtop);
1365 if ( (PAR(cr)) && !DOMAINDECOMP(cr) )
1367 gmx_fatal(FARGS, "Position swapping is only implemented for domain decomposition!");
1370 bAppend = Flags & MD_APPENDFILES;
1371 bStartFromCpt = Flags & MD_STARTFROMCPT;
1372 bRerun = Flags & MD_RERUN;
1375 snew(sc->si_priv, 1);
1382 gmx_fatal(FARGS, "%s This module does not support reruns in parallel\nPlease request a serial run with -nt 1 / -np 1\n", SwS);
1385 fprintf(stderr, "%s Rerun - using every available frame\n", SwS);
1387 sc->nAverage = 1; /* averaging makes no sense for reruns */
1390 if (MASTER(cr) && !bAppend)
1392 fprintf(fplog, "\nInitializing ion/water position exchanges\n");
1393 please_cite(fplog, "Kutzner2011b");
1396 switch (ir->eSwapCoords)
1412 /* Copy some data to the group structures for convenience */
1413 /* Number of atoms in the group */
1414 s->group[eGrpIons ].nat = sc->nat;
1415 s->group[eGrpSplit0 ].nat = sc->nat_split[0];
1416 s->group[eGrpSplit1 ].nat = sc->nat_split[1];
1417 s->group[eGrpSolvent].nat = sc->nat_sol;
1418 /* Pointer to the indices */
1419 s->group[eGrpIons ].ind = sc->ind;
1420 s->group[eGrpSplit0 ].ind = sc->ind_split[0];
1421 s->group[eGrpSplit1 ].ind = sc->ind_split[1];
1422 s->group[eGrpSolvent].ind = sc->ind_sol;
1424 check_swap_groups(s, mtop->natoms, bVerbose && MASTER(cr));
1426 /* Allocate space for the collective arrays for all groups */
1427 for (ig = 0; ig < eGrpNr; ig++)
1429 g = &(s->group[ig]);
1430 snew(g->xc, g->nat);
1431 snew(g->c_ind_loc, g->nat);
1432 /* For the split groups (the channels) we need some extra memory to
1433 * be able to make the molecules whole even if they span more than
1434 * half of the box size. */
1435 if (eGrpSplit0 == ig || eGrpSplit1 == ig)
1437 snew(g->xc_shifts, g->nat);
1438 snew(g->xc_eshifts, g->nat);
1439 snew(g->xc_old, g->nat);
1445 init_swapstate(swapstate, sc, mtop, x, box, ir->ePBC);
1448 /* After init_swapstate we have a set of (old) whole positions for our
1449 * channels. Now transfer that to all nodes */
1452 for (ig = eGrpSplit0; ig <= eGrpSplit1; ig++)
1454 g = &(s->group[ig]);
1455 gmx_bcast((g->nat)*sizeof((g->xc_old)[0]), g->xc_old, (cr));
1459 /* Make sure that all molecules in the ion and solvent groups contain the
1460 * same number of atoms each */
1461 s->group[eGrpIons ].apm = get_group_apm_check(eGrpIons, s, MASTER(cr) && bVerbose, alook, mtop);
1462 s->group[eGrpSolvent].apm = get_group_apm_check(eGrpSolvent, s, MASTER(cr) && bVerbose, alook, mtop);
1464 /* Save masses where needed */
1465 s->group[eGrpIons ].m = NULL;
1466 /* We only need enough space to determine a single solvent molecule's
1467 * center at at time */
1468 g = &(s->group[eGrpSolvent]);
1471 /* Need mass-weighted center of split group? */
1472 for (j = 0, ig = eGrpSplit0; j < eChanNR; ig++, j++)
1474 g = &(s->group[ig]);
1475 if (TRUE == sc->massw_split[j])
1477 /* Save the split group charges if mass-weighting is requested */
1479 for (i = 0; i < g->nat; i++)
1481 gmx_mtop_atomnr_to_atom(alook, g->ind[i], &atom);
1491 /* Save the ionic charges */
1492 g = &(s->group[eGrpIons]);
1493 snew(g->qc, g->nat);
1494 for (i = 0; i < g->nat; i++)
1496 gmx_mtop_atomnr_to_atom(alook, g->ind[i], &atom);
1501 set_pbc(s->pbc, -1, box);
1508 fprintf(stderr, "%s Opening output file %s%s\n", SwS, fn, bAppend ? " for appending" : "");
1511 s->fpout = gmx_fio_fopen(fn, bAppend ? "a" : "w" );
1515 xvgr_header(s->fpout, "Ion counts", "Time (ps)", "counts", exvggtXNY, oenv);
1517 for (ig = 0; ig < eGrpNr; ig++)
1519 g = &(s->group[ig]);
1520 fprintf(s->fpout, "# %s group contains %d atom%s", GrpString[ig], g->nat, (g->nat > 1) ? "s" : "");
1521 if (eGrpSolvent == ig || eGrpIons == ig)
1523 fprintf(s->fpout, " with %d atom%s in each molecule", g->apm, (g->apm > 1) ? "s" : "");
1525 fprintf(s->fpout, ".\n");
1528 fprintf(s->fpout, "#\n# Initial positions of split groups:\n");
1531 for (j = 0, ig = eGrpSplit0; j < eChanNR; j++, ig++)
1533 g = &(s->group[ig]);
1534 for (i = 0; i < g->nat; i++)
1536 copy_rvec(x[sc->ind_split[j][i]], g->xc[i]);
1538 if (eGrpSplit0 == ig || eGrpSplit1 == ig)
1540 /* xc has the correct PBC representation for the two channels, so we do
1541 * not need to correct for that */
1542 get_center(g->xc, g->m, g->nat, g->center);
1546 /* For the water molecules, we need to make the molecules whole */
1547 get_molecule_center(g->xc, g->nat, g->m, g->center, s->pbc);
1551 fprintf(s->fpout, "# %s group %s-center %5f nm\n", GrpString[ig],
1552 DimStr[s->swapdim], g->center[s->swapdim]);
1558 fprintf(s->fpout, "#\n");
1559 fprintf(s->fpout, "# split0 cylinder radius %f nm, up %f nm, down %f nm\n",
1560 sc->cyl0r, sc->cyl0u, sc->cyl0l);
1561 fprintf(s->fpout, "# split1 cylinder radius %f nm, up %f nm, down %f nm\n",
1562 sc->cyl1r, sc->cyl1u, sc->cyl1l);
1567 fprintf(s->fpout, "#\n");
1570 fprintf(s->fpout, "# Coupling constant (number of swap attempt steps to average over): %d (translates to %f ps).\n",
1571 sc->nAverage, sc->nAverage*sc->nstswap*ir->delta_t);
1572 fprintf(s->fpout, "# Threshold is %f\n", sc->threshold);
1573 fprintf(s->fpout, "#\n");
1574 fprintf(s->fpout, "# Remarks about which atoms passed which channel use global atoms numbers starting at one.\n");
1583 /* Prepare for parallel or serial run */
1586 for (ig = 0; ig < eGrpNr; ig++)
1588 g = &(s->group[ig]);
1596 for (ig = 0; ig < eGrpNr; ig++)
1598 g = &(s->group[ig]);
1599 g->nat_loc = g->nat;
1600 g->ind_loc = g->ind;
1601 /* c_ind_loc needs to be set to identity in the serial case */
1602 for (i = 0; i < g->nat; i++)
1604 g->c_ind_loc[i] = i;
1609 /* Allocate memory for the ion counts time window */
1610 for (ic = 0; ic < eCompNR; ic++)
1612 for (ii = 0; ii < eIonNR; ii++)
1614 snew(s->comp[ic][ii].nat_past, sc->nAverage);
1618 /* Get the initial ion concentrations and let the other nodes know */
1621 swapstate->nions = s->group[eGrpIons].nat;
1625 get_initial_ioncounts_from_cpt(ir, swapstate, cr, bVerbose);
1629 fprintf(stderr, "%s Determining initial ion counts.\n", SwS);
1630 get_initial_ioncounts(ir, x, box, cr, bRerun);
1633 /* Prepare (further) checkpoint writes ... */
1636 /* Consistency check */
1637 if (swapstate->nAverage != sc->nAverage)
1639 gmx_fatal(FARGS, "%s Ion count averaging steps mismatch! checkpoint: %d, tpr: %d",
1640 SwS, swapstate->nAverage, sc->nAverage);
1645 swapstate->nAverage = sc->nAverage;
1647 fprintf(stderr, "%s Setting pointers for checkpoint writing\n", SwS);
1648 for (ic = 0; ic < eCompNR; ic++)
1650 for (ii = 0; ii < eIonNR; ii++)
1652 swapstate->nat_req_p[ic][ii] = &(s->comp[ic][ii].nat_req);
1653 swapstate->nat_past_p[ic][ii] = &(s->comp[ic][ii].nat_past[0]);
1654 swapstate->inflow_netto_p[ic][ii] = &(s->comp[ic][ii].inflow_netto);
1658 /* Determine the total charge imbalance */
1659 s->deltaQ = ( (-1) * s->comp[eCompA][eIonNEG].nat_req + s->comp[eCompA][eIonPOS].nat_req )
1660 - ( (-1) * s->comp[eCompB][eIonNEG].nat_req + s->comp[eCompB][eIonPOS].nat_req );
1664 fprintf(stderr, "%s Requested charge imbalance is Q(A) - Q(B) = %gz.\n", SwS, s->deltaQ);
1668 fprintf(s->fpout, "# Requested charge imbalance is Q(A)-Q(B) = %gz.\n", s->deltaQ);
1674 bc_initial_concentrations(cr, ir->swap);
1677 /* Put the time-averaged number of ions for all compartments */
1678 for (ic = 0; ic < eCompNR; ic++)
1680 for (ii = 0; ii < eIonNR; ii++)
1682 update_time_window(&(s->comp[ic][ii]), sc->nAverage, -1);
1686 /* Initialize arrays that keep track of through which channel the ions go */
1687 detect_flux_per_channel_init(cr, s, swapstate, bStartFromCpt);
1689 /* We need to print the legend if we open this file for the first time. */
1690 if (MASTER(cr) && !bAppend)
1692 print_ionlist_legend(ir, oenv);
1697 extern void dd_make_local_swap_groups(gmx_domdec_t *dd, t_swapcoords *sc)
1703 /* Make ion group, split groups and solvent group */
1704 for (ig = 0; ig < eGrpNr; ig++)
1706 g = &(sc->si_priv->group[ig]);
1707 dd_make_local_group_indices(dd->ga2la, g->nat, g->ind,
1708 &(g->nat_loc), &(g->ind_loc), &(g->nalloc_loc), g->c_ind_loc);
1713 /*! \brief Do we need to swap ions with water molecules at this step?
1715 * From the requested and average ion counts we determine whether a swap is needed
1716 * at this time step.
1718 static gmx_bool need_swap(t_swapcoords *sc)
1725 for (ic = 0; ic < eCompNR; ic++)
1727 for (ii = 0; ii < eIonNR; ii++)
1729 if (s->comp[ic][ii].nat_req - s->comp[ic][ii].nat_av >= sc->threshold)
1739 /*! \brief Return index of atom that we can use for swapping.
1741 * Returns the index of an atom that is far off the compartment boundaries.
1742 * Other atoms of the molecule (if any) will directly follow the returned index
1744 static int get_index_of_distant_atom(
1745 t_compartment *comp,
1746 int apm) /* Atoms per molecule - just return the first atom index of a molecule */
1749 real d = GMX_REAL_MAX;
1752 /* comp->nat contains the original number of atoms in this compartment
1753 * prior to doing any swaps. Some of these atoms may already have been
1754 * swapped out, but then they are marked with a distance of GMX_REAL_MAX
1756 for (i = 0; i < comp->nat_old; i += apm)
1758 if (comp->dist[i] < d)
1761 d = comp->dist[ibest];
1767 gmx_fatal(FARGS, "Could not get index of swap atom. Compartment atoms %d before swaps, atoms per molecule %d.",
1768 comp->nat_old, apm);
1771 /* Set the distance of this index to infinity such that it won't get selected again in
1774 comp->dist[ibest] = GMX_REAL_MAX;
1776 return comp->ind[ibest];
1780 /*! \brief Swaps centers of mass and makes molecule whole if broken */
1781 static void translate_positions(
1789 rvec reference, dx, correctPBCimage;
1792 /* Use the first atom as the reference for PBC */
1793 copy_rvec(x[0], reference);
1795 for (i = 0; i < apm; i++)
1797 /* PBC distance between position and reference */
1798 pbc_dx(pbc, x[i], reference, dx);
1800 /* Add PBC distance to reference */
1801 rvec_add(reference, dx, correctPBCimage);
1803 /* Subtract old_com from correct image and add new_com */
1804 rvec_dec(correctPBCimage, old_com);
1805 rvec_inc(correctPBCimage, new_com);
1807 copy_rvec(correctPBCimage, x[i]);
1812 /*! \brief Write back the the modified local positions from the collective array to the official positions. */
1813 static void apply_modified_positions(
1820 for (l = 0; l < g->nat_loc; l++)
1822 /* Get the right local index to write to */
1824 /* Where is the local atom in the collective array? */
1825 cind = g->c_ind_loc[l];
1827 /* Copy the possibly modified position */
1828 copy_rvec(g->xc[cind], x[ii]);
1833 extern gmx_bool do_swapcoords(
1838 gmx_wallcycle_t wcycle,
1847 int j, ii, ic, ig, im, gmax, nswaps;
1848 gmx_bool bSwap = FALSE;
1850 real vacancy[eCompNR][eIonNR];
1852 rvec solvent_center, ion_center;
1854 gmx_mtop_atomlookup_t alook = NULL;
1857 wallcycle_start(wcycle, ewcSWAP);
1862 /* Assemble all the positions of the swap group (ig = 0), the split groups
1863 * (ig = 1,2), and possibly the solvent group (ig = 3) */
1866 for (ig = 0; ig < gmax; ig++)
1868 g = &(s->group[ig]);
1869 if (eGrpSplit0 == ig || eGrpSplit1 == ig)
1871 /* The split groups, i.e. the channels. Here we need the full
1872 * communicate_group_positions(), so that we can make the molecules
1873 * whole even in cases where they span more than half of the box in
1875 communicate_group_positions(cr, g->xc, g->xc_shifts, g->xc_eshifts, TRUE,
1876 x, g->nat, g->nat_loc, g->ind_loc, g->c_ind_loc, g->xc_old, box);
1878 get_center(g->xc, g->m, g->nat, g->center); /* center of split groups == channels */
1882 /* Swap group (ions), and solvent group. These molecules are small
1883 * and we can always make them whole with a simple distance check.
1884 * Therefore we pass NULL as third argument. */
1885 communicate_group_positions(cr, g->xc, NULL, NULL, FALSE,
1886 x, g->nat, g->nat_loc, g->ind_loc, g->c_ind_loc, NULL, NULL);
1890 /* Set up the compartments and get lists of atoms in each compartment,
1891 * determine how many ions each compartment contains */
1892 compartmentalize_ions(cr, sc, box, step, s->fpout, bRerun);
1894 /* Output how many ions are in the compartments */
1897 print_ionlist(s, t, "");
1900 /* If we are doing a rerun, we are finished here, since we cannot perform
1907 /* Do we have to perform a swap? */
1908 bSwap = need_swap(sc);
1911 g = &(s->group[eGrpSolvent]);
1912 communicate_group_positions(cr, g->xc, NULL, NULL, FALSE,
1913 x, g->nat, g->nat_loc, g->ind_loc, g->c_ind_loc, NULL, NULL);
1915 compartmentalize_solvent(cr, sc, box, s->fpout);
1917 /* Determine where ions are missing and where ions are too many */
1918 for (ic = 0; ic < eCompNR; ic++)
1920 for (ii = 0; ii < eIonNR; ii++)
1922 vacancy[ic][ii] = s->comp[ic][ii].nat_req - s->comp[ic][ii].nat_av;
1926 /* Remember the original number of ions per compartment */
1927 for (ic = 0; ic < eCompNR; ic++)
1929 s->compsol[ic].nat_old = s->compsol[ic].nat;
1930 for (ii = 0; ii < eIonNR; ii++)
1932 s->comp[ic][ii].nat_old = s->comp[ic][ii].nat;
1936 /* Now actually correct the number of ions */
1938 alook = gmx_mtop_atomlookup_init(mtop);
1939 for (ic = 0; ic < eCompNR; ic++)
1941 for (ii = 0; ii < eIonNR; ii++)
1943 while (vacancy[ic][ii] >= sc->threshold)
1945 /* Swap in an ion */
1947 /* Get the xc-index of the first atom of a solvent molecule of this compartment */
1948 isol = get_index_of_distant_atom(&(s->compsol[ic]), s->group[eGrpSolvent].apm );
1950 /* Get the xc-index of an ion from the other compartment */
1951 iion = get_index_of_distant_atom(&(s->comp[(ic+1)%eCompNR][ii]), s->group[eGrpIons].apm );
1953 /* Get the solvent molecule's center of mass */
1954 for (im = 0; im < s->group[eGrpSolvent].apm; im++)
1956 gmx_mtop_atomnr_to_atom(alook, s->group[eGrpSolvent].ind[isol+im], &atom);
1957 s->group[eGrpSolvent].m[im] = atom->m;
1959 get_molecule_center(&(s->group[eGrpSolvent].xc[isol]), s->group[eGrpSolvent].apm, s->group[eGrpSolvent].m, solvent_center, s->pbc);
1960 get_molecule_center(&(s->group[eGrpIons ].xc[iion]), s->group[eGrpIons ].apm, NULL, ion_center, s->pbc);
1962 /* subtract com_solvent and add com_ion */
1963 translate_positions(&(s->group[eGrpSolvent].xc[isol]), s->group[eGrpSolvent].apm, solvent_center, ion_center, s->pbc);
1964 /* For the ion, subtract com_ion and add com_solvent */
1965 translate_positions(&(s->group[eGrpIons ].xc[iion]), s->group[eGrpIons ].apm, ion_center, solvent_center, s->pbc);
1968 vacancy[(ic+1) % eCompNR][ii]++;
1970 /* Keep track of the changes */
1971 s->comp[ic ][ii].nat++;
1972 s->comp[(ic+1) % eCompNR][ii].nat--;
1973 s->comp[ic ][ii].inflow_netto++;
1974 s->comp[(ic+1) % eCompNR][ii].inflow_netto--;
1975 /* Correct the past time window to still get the right averages from now on */
1976 s->comp[ic ][ii].nat_av++;
1977 s->comp[(ic+1) % eCompNR][ii].nat_av--;
1978 for (j = 0; j < sc->nAverage; j++)
1980 s->comp[ic ][ii].nat_past[j]++;
1981 s->comp[(ic+1) % eCompNR][ii].nat_past[j]--;
1983 /* Clear ion history */
1986 s->group[eGrpIons].channel_label[iion] = eChHistPassedNone;
1987 s->group[eGrpIons].comp_from[iion] = eDomainNotset;
1989 /* That was the swap */
1994 gmx_mtop_atomlookup_destroy(alook);
1998 fprintf(stderr, "%s Performed %d swap%s in step %" GMX_PRId64 ".\n", SwS, nswaps, nswaps > 1 ? "s" : "", step);
2000 if (s->fpout != NULL)
2002 print_ionlist(s, t, " # after swap");
2005 /* Write back the the modified local positions from the collective array to the official coordinates */
2006 apply_modified_positions(&(s->group[eGrpIons ]), x);
2007 apply_modified_positions(&(s->group[eGrpSolvent]), x);
2008 } /* end of if(bSwap) */
2010 wallcycle_stop(wcycle, ewcSWAP);