2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015,2016,2017, 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"
52 #include "gromacs/domdec/domdec_struct.h"
53 #include "gromacs/fileio/confio.h"
54 #include "gromacs/fileio/gmxfio.h"
55 #include "gromacs/fileio/xvgr.h"
56 #include "gromacs/gmxlib/network.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/mdlib/groupcoord.h"
59 #include "gromacs/mdlib/mdrun.h"
60 #include "gromacs/mdlib/sim_util.h"
61 #include "gromacs/mdtypes/inputrec.h"
62 #include "gromacs/mdtypes/md_enums.h"
63 #include "gromacs/pbcutil/pbc.h"
64 #include "gromacs/timing/wallcycle.h"
65 #include "gromacs/topology/mtop_util.h"
66 #include "gromacs/topology/topology.h"
67 #include "gromacs/utility/cstringutil.h"
68 #include "gromacs/utility/fatalerror.h"
69 #include "gromacs/utility/pleasecite.h"
70 #include "gromacs/utility/smalloc.h"
71 #include "gromacs/utility/snprintf.h"
73 static const char *SwS = {"SWAP:"}; /**< For output that comes from the swap module */
74 static const char *SwSEmpty = {" "}; /**< Placeholder for multi-line output */
75 static const char* CompStr[eCompNR] = {"A", "B" }; /**< Compartment name */
76 static const char *SwapStr[eSwapTypesNR+1] = { "", "X-", "Y-", "Z-", NULL}; /**< Name for the swap types. */
77 static const char *DimStr[DIM+1] = { "X", "Y", "Z", NULL}; /**< Name for the swap dimension. */
79 /** Keep track of through which channel the ions have passed */
80 enum eChannelHistory {
81 eChHistPassedNone, eChHistPassedCh0, eChHistPassedCh1, eChHistNr
83 static const char* ChannelString[eChHistNr] = { "none", "channel0", "channel1" }; /**< Name for the channels */
85 /*! \brief Domain identifier.
87 * Keeps track of from which compartment the ions came before passing the
91 eDomainNotset, eDomainA, eDomainB, eDomainNr
93 static const char* DomainString[eDomainNr] = { "not_assigned", "Domain_A", "Domain_B" }; /**< Name for the domains */
98 * Structure containing compartment-specific data.
100 typedef struct swap_compartment
102 int nMol; /**< Number of ion or water molecules detected
103 in this compartment. */
104 int nMolBefore; /**< Number of molecules before swapping. */
105 int nMolReq; /**< Requested number of molecules in compartment. */
106 real nMolAv; /**< Time-averaged number of molecules matching
107 the compartment conditions. */
108 int *nMolPast; /**< Past molecule counts for time-averaging. */
109 int *ind; /**< Indices to collective array of atoms. */
110 real *dist; /**< Distance of atom to bulk layer, which is
111 normally the center layer of the compartment */
112 int nalloc; /**< Allocation size for ind array. */
113 int inflow_net; /**< Net inflow of ions into this compartment. */
118 * This structure contains data needed for the groups involved in swapping:
119 * split group 0, split group 1, solvent group, ion groups.
121 typedef struct swap_group
123 char *molname; /**< Name of the group or ion type */
124 int nat; /**< Number of atoms in the group */
125 int apm; /**< Number of atoms in each molecule */
126 int *ind; /**< Global atom indices of the group (size nat) */
127 int *ind_loc; /**< Local atom indices of the group */
128 int nat_loc; /**< Number of local group atoms */
129 int nalloc_loc; /**< Allocation size for ind_loc */
130 rvec *xc; /**< Collective array of group atom positions (size nat) */
131 ivec *xc_shifts; /**< Current (collective) shifts (size nat) */
132 ivec *xc_eshifts; /**< Extra shifts since last DD step (size nat) */
133 rvec *xc_old; /**< Old (collective) positions (size nat) */
134 real q; /**< Total charge of one molecule of this group */
135 int *c_ind_loc; /**< Position of local atoms in the
136 collective array, [0..nat_loc] */
137 real *m; /**< Masses (can be omitted, size apm) */
138 unsigned char *comp_from; /**< (Collective) Stores from which compartment this
139 molecule has come. This way we keep track of
140 through which channel an ion permeates
141 (size nMol = nat/apm) */
142 unsigned char *comp_now; /**< In which compartment this ion is now (size nMol) */
143 unsigned char *channel_label; /**< Which channel was passed at last by this ion?
145 rvec center; /**< Center of the group; COM if masses are used */
146 t_compartment comp[eCompNR]; /**< Distribution of particles of this group across
147 the two compartments */
148 real vacancy[eCompNR]; /**< How many molecules need to be swapped in? */
149 int fluxfromAtoB[eChanNR]; /**< Net flux of ions per channel */
150 int nCyl[eChanNR]; /**< Number of ions residing in a channel */
151 int nCylBoth; /**< Ions assigned to cyl0 and cyl1. Not good. */
156 * Main (private) data structure for the position swapping protocol.
158 typedef struct t_swap
160 int swapdim; /**< One of XX, YY, ZZ */
161 t_pbc *pbc; /**< Needed to make molecules whole. */
162 FILE *fpout; /**< Output file. */
163 int ngrp; /**< Number of t_swapgrp groups */
164 t_swapgrp *group; /**< Separate groups for channels, solvent, ions */
165 int fluxleak; /**< Flux not going through any of the channels. */
166 real deltaQ; /**< The charge imbalance between the compartments. */
171 /*! \brief Check whether point is in channel.
173 * A channel is a cylinder defined by a disc
174 * with radius r around its center c. The thickness of the cylinder is
181 * <---------c--------->
187 * \param[in] point The position (xyz) under consideration.
188 * \param[in] center The center of the cylinder.
189 * \param[in] d_up The upper extension of the cylinder.
190 * \param[in] d_down The lower extension.
191 * \param[in] r_cyl2 Cylinder radius squared.
192 * \param[in] pbc Structure with info about periodic boundary conditions.
193 * \param[in] normal The membrane normal direction is typically 3, i.e. z, but can be x or y also.
195 * \returns Whether the point is inside the defined cylindric channel.
197 static gmx_bool is_in_channel(
207 int plane1, plane2; /* Directions tangential to membrane */
210 plane1 = (normal + 1) % 3; /* typically 0, i.e. XX */
211 plane2 = (normal + 2) % 3; /* typically 1, i.e. YY */
213 /* Get the distance vector dr between the point and the center of the cylinder */
214 pbc_dx(pbc, point, center, dr); /* This puts center in the origin */
216 /* Check vertical direction */
217 if ( (dr[normal] > d_up) || (dr[normal] < -d_down) )
222 /* Check radial direction */
223 if ( (dr[plane1]*dr[plane1] + dr[plane2]*dr[plane2]) > r_cyl2)
228 /* All check passed, this point is in the cylinder */
233 /*! \brief Prints output to CompEL output file.
235 * Prints to swap output file how many ions are in each compartment,
236 * where the centers of the split groups are, and how many ions of each type
237 * passed the channels.
239 static void print_ionlist(
242 const char comment[])
245 fprintf(s->fpout, "%12.5e", time);
247 // Output number of molecules and difference to reference counts for each
248 // compartment and ion type
249 for (int iComp = 0; iComp < eCompNR; iComp++)
251 for (int ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
253 t_compartment *comp = &s->group[ig].comp[iComp];
255 fprintf(s->fpout, "%10d%10.1f%10d", comp->nMol, comp->nMolAv - comp->nMolReq, comp->inflow_net);
259 // Output center of split groups
260 fprintf(s->fpout, "%10g%10g",
261 s->group[eGrpSplit0].center[s->swapdim],
262 s->group[eGrpSplit1].center[s->swapdim]);
264 // Output ion flux for each channel and ion type
265 for (int iChan = 0; iChan < eChanNR; iChan++)
267 for (int ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
269 t_swapgrp *g = &s->group[ig];
270 fprintf(s->fpout, "%10d", g->fluxfromAtoB[iChan]);
274 /* Output the number of molecules that leaked from A to B */
275 fprintf(s->fpout, "%10d", s->fluxleak);
277 fprintf(s->fpout, "%s\n", comment);
281 /*! \brief Get the center of a group of nat atoms.
283 * Since with PBC an atom group might not be whole, use the first atom as the
284 * reference atom and determine the center with respect to this reference.
286 static void get_molecule_center(
294 rvec weightedPBCimage;
296 rvec reference, correctPBCimage, dx;
299 /* Use the first atom as the reference and put other atoms near that one */
300 /* This does not work for large molecules that span > half of the box! */
301 copy_rvec(x[0], reference);
303 /* Calculate either the weighted center or simply the center of geometry */
306 for (i = 0; i < nat; i++)
308 /* PBC distance between position and reference */
309 pbc_dx(pbc, x[i], reference, dx);
311 /* Add PBC distance to reference */
312 rvec_add(reference, dx, correctPBCimage);
314 /* Take weight into account */
324 svmul(wi, correctPBCimage, weightedPBCimage);
327 rvec_inc(center, weightedPBCimage);
331 svmul(1.0/wsum, center, center);
336 /*! \brief Return TRUE if position x of ion (or water) is found in the compartment,
337 * i.e. between w1 and w2.
339 * One can define and additional offset "b" if one wants to exchange ions/water
340 * to or from a plane not directly in the middle of w1 and w2. The offset can be
341 * in ]-1.0, ..., +1.0 [.
342 * A bulkOffset of 0.0 means 'no offset', so the swap-layer is directly in the
343 * middle between w1 and w2. Offsets -1.0 < b < 0.0 will yield swaps nearer to w1,
344 * whereas offsets 0.0 < 0 < +1.0 will yield swaps nearer to w2.
348 * ||--------------+-------------|-------------+------------------------||
349 * w1 ? ? ? ? ? ? ? ? ? ? ? w2
350 * ||--------------+-------------|----b--------+------------------------||
355 * \param[in] w1 Position of 'wall' atom 1.
356 * \param[in] w2 Position of 'wall' atom 2.
357 * \param[in] x Position of the ion or the water molecule under consideration.
358 * \param[in] l Length of the box, from || to || in the sketch.
359 * \param[in] bulkOffset Where is the bulk layer "b" to be found between w1 and w2?
360 * \param[out] distance_from_b Distance of x to the bulk layer "b".
362 * \returns TRUE if x is between w1 and w2.
364 * Also computes the distance of x to the compartment center (the layer that is
365 * normally situated in the middle of w1 and w2 that would be considered as having
366 * the bulk concentration of ions).
368 static gmx_bool compartment_contains_atom(
374 real *distance_from_b)
380 /* First set the origin in the middle of w1 and w2 */
387 /* Now choose the PBC image of x that is closest to the origin: */
398 *distance_from_b = (real)fabs(x - bulkOffset*0.5*width);
400 /* Return TRUE if we now are in area "????" */
401 if ( (x >= w1) && (x < w2) )
412 /*! \brief Updates the time-averaged number of ions in a compartment. */
413 static void update_time_window(t_compartment *comp, int values, int replace)
419 /* Put in the new value */
422 comp->nMolPast[replace] = comp->nMol;
425 /* Compute the new time-average */
427 for (i = 0; i < values; i++)
429 average += comp->nMolPast[i];
432 comp->nMolAv = average;
436 /*! \brief Add the atom with collective index ci to the atom list in compartment 'comp'.
438 * \param[in] ci Index of this ion in the collective xc array.
439 * \param[inout] comp Compartment to add this atom to.
440 * \param[in] distance Shortest distance of this atom to the bulk layer,
441 * from which ion/water pairs are selected for swapping.
443 static void add_to_list(
450 if (nr >= comp->nalloc)
452 comp->nalloc = over_alloc_dd(nr+1);
453 srenew(comp->ind, comp->nalloc);
454 srenew(comp->dist, comp->nalloc);
457 comp->dist[nr] = distance;
462 /*! \brief Determine the compartment boundaries from the channel centers. */
463 static void get_compartment_boundaries(
467 real *left, real *right)
470 real leftpos, rightpos, leftpos_orig;
475 gmx_fatal(FARGS, "No compartment %c.", c+'A');
478 pos0 = s->group[eGrpSplit0].center[s->swapdim];
479 pos1 = s->group[eGrpSplit1].center[s->swapdim];
492 /* This gets us the other compartment: */
495 leftpos_orig = leftpos;
497 rightpos = leftpos_orig + box[s->swapdim][s->swapdim];
505 /*! \brief Determine the per-channel ion flux.
507 * To determine the flux through the individual channels, we
508 * remember the compartment and channel history of each ion. An ion can be
509 * either in channel0 or channel1, or in the remaining volume of compartment
513 * +-----------------+
516 * ||||||||||0|||||||| bilayer with channel 0
521 * |||||1||||||||||||| bilayer with channel 1
524 * +-----------------+
528 static void detect_flux_per_channel(
533 unsigned char *comp_now,
534 unsigned char *comp_from,
535 unsigned char *channel_label,
545 gmx_bool in_cyl0, in_cyl1;
552 /* Check whether ion is inside any of the channels */
553 in_cyl0 = is_in_channel(atomPosition, s->group[eGrpSplit0].center, sc->cyl0u, sc->cyl0l, cyl0_r2, s->pbc, sd);
554 in_cyl1 = is_in_channel(atomPosition, s->group[eGrpSplit1].center, sc->cyl1u, sc->cyl1l, cyl1_r2, s->pbc, sd);
556 if (in_cyl0 && in_cyl1)
558 /* Ion appears to be in both channels. Something is severely wrong! */
560 *comp_now = eDomainNotset;
561 *comp_from = eDomainNotset;
562 *channel_label = eChHistPassedNone;
566 /* Ion is in channel 0 now */
567 *channel_label = eChHistPassedCh0;
568 *comp_now = eDomainNotset;
573 /* Ion is in channel 1 now */
574 *channel_label = eChHistPassedCh1;
575 *comp_now = eDomainNotset;
580 /* Ion is not in any of the channels, so it must be in domain A or B */
583 *comp_now = eDomainA;
587 *comp_now = eDomainB;
591 /* Only take action, if ion is now in domain A or B, and was before
592 * in the other domain!
594 if (eDomainNotset == *comp_from)
596 /* Maybe we can set the domain now */
597 *comp_from = *comp_now; /* Could still be eDomainNotset, though */
599 else if ( (*comp_now != eDomainNotset ) /* if in channel */
600 && (*comp_from != *comp_now) )
602 /* Obviously the ion changed its domain.
603 * Count this for the channel through which it has passed. */
604 switch (*channel_label)
606 case eChHistPassedNone:
609 fprintf(stderr, " %s Warning! Step %s, ion %d moved from %s to %s\n",
610 SwS, gmx_step_str(step, buf), iAtom, DomainString[*comp_from], DomainString[*comp_now]);
613 fprintf(stderr, ", possibly due to a swap in the original simulation.\n");
617 fprintf(stderr, "but did not pass cyl0 or cyl1 as defined in the .mdp file.\n"
618 "Do you have an ion somewhere within the membrane?\n");
619 /* Write this info to the CompEL output file: */
620 fprintf(s->fpout, " # Warning: step %s, ion %d moved from %s to %s (probably through the membrane)\n",
621 gmx_step_str(step, buf), iAtom,
622 DomainString[*comp_from], DomainString[*comp_now]);
626 case eChHistPassedCh0:
627 case eChHistPassedCh1:
628 if (*channel_label == eChHistPassedCh0)
637 if (eDomainA == *comp_from)
639 g->fluxfromAtoB[chan_nr]++;
643 g->fluxfromAtoB[chan_nr]--;
645 fprintf(fpout, "# Atom nr. %d finished passing %s.\n", iAtom, ChannelString[*channel_label]);
648 gmx_fatal(FARGS, "%s Unknown channel history entry for ion type '%s'\n",
653 /* This ion has moved to the _other_ compartment ... */
654 *comp_from = *comp_now;
655 /* ... and it did not pass any channel yet */
656 *channel_label = eChHistPassedNone;
661 /*! \brief Determines which ions or solvent molecules are in compartment A and B */
662 static void sortMoleculesIntoCompartments(
672 gmx_swapcoords_t s = sc->si_priv;
673 int nMolNotInComp[eCompNR]; /* consistency check */
674 real cyl0_r2 = sc->cyl0r * sc->cyl0r;
675 real cyl1_r2 = sc->cyl1r * sc->cyl1r;
677 /* Get us a counter that cycles in the range of [0 ... sc->nAverage[ */
678 int replace = (step/sc->nstswap) % sc->nAverage;
680 for (int comp = eCompA; comp <= eCompB; comp++)
684 /* Get lists of atoms that match criteria for this compartment */
685 get_compartment_boundaries(comp, sc->si_priv, box, &left, &right);
687 /* First clear the ion molecule lists */
688 g->comp[comp].nMol = 0;
689 nMolNotInComp[comp] = 0; /* consistency check */
691 /* Loop over the molecules and atoms of this group */
692 for (int iMol = 0, iAtom = 0; iAtom < g->nat; iAtom += g->apm, iMol++)
697 /* Is this first atom of the molecule in the compartment that we look at? */
698 if (compartment_contains_atom(left, right, g->xc[iAtom][sd], box[sd][sd], sc->bulkOffset[comp], &dist) )
700 /* Add the first atom of this molecule to the list of molecules in this compartment */
701 add_to_list(iAtom, &g->comp[comp], dist);
703 /* Master also checks for ion groups through which channel each ion has passed */
704 if (MASTER(cr) && (g->comp_now != NULL) && !bIsSolvent)
706 int globalAtomNr = g->ind[iAtom] + 1; /* PDB index starts at 1 ... */
707 detect_flux_per_channel(g, globalAtomNr, comp, g->xc[iAtom],
708 &g->comp_now[iMol], &g->comp_from[iMol], &g->channel_label[iMol],
709 sc, cyl0_r2, cyl1_r2, step, bRerun, fpout);
714 nMolNotInComp[comp]++;
717 /* Correct the time-averaged number of ions in the compartment */
720 update_time_window(&g->comp[comp], sc->nAverage, replace);
724 /* Flux detection warnings */
725 if (MASTER(cr) && !bIsSolvent)
730 "%s Warning: %d atoms were detected as being in both channels! Probably your split\n"
731 "%s cylinder is way too large, or one compartment has collapsed (step %" GMX_PRId64 ")\n",
732 SwS, g->nCylBoth, SwS, step);
734 fprintf(s->fpout, "Warning: %d atoms were assigned to both channels!\n", g->nCylBoth);
740 if (bIsSolvent && NULL != fpout)
742 fprintf(fpout, "# Solv. molecules in comp.%s: %d comp.%s: %d\n",
743 CompStr[eCompA], g->comp[eCompA].nMol,
744 CompStr[eCompB], g->comp[eCompB].nMol);
747 /* Consistency checks */
748 if (nMolNotInComp[eCompA] + nMolNotInComp[eCompB] != g->nat / g->apm)
750 fprintf(stderr, "%s Warning: Inconsistency while assigning '%s' molecules to compartments. !inA: %d, !inB: %d, total molecules %d\n",
751 SwS, g->molname, nMolNotInComp[eCompA], nMolNotInComp[eCompB], g->nat/g->apm);
754 int sum = g->comp[eCompA].nMol + g->comp[eCompB].nMol;
755 if (sum != g->nat/g->apm)
757 fprintf(stderr, "%s Warning: %d molecules are in group '%s', but altogether %d have been assigned to the compartments.\n",
758 SwS, g->nat/g->apm, g->molname, sum);
763 /*! \brief Find out how many group atoms are in the compartments initially */
764 static void get_initial_ioncounts(
766 rvec x[], /* the initial positions */
782 /* Loop over the user-defined (ion) groups */
783 for (ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
787 /* Copy the initial positions of the atoms in the group
788 * to the collective array so that we can compartmentalize */
789 for (i = 0; i < g->nat; i++)
792 copy_rvec(x[ind], g->xc[i]);
795 /* Set up the compartments and get lists of atoms in each compartment */
796 sortMoleculesIntoCompartments(g, cr, sc, box, 0, s->fpout, bRerun, FALSE);
798 /* Set initial molecule counts if requested (as signaled by "-1" value) */
799 for (ic = 0; ic < eCompNR; ic++)
801 int requested = sc->grp[ig].nmolReq[ic];
804 g->comp[ic].nMolReq = g->comp[ic].nMol;
808 g->comp[ic].nMolReq = requested;
812 /* Check whether the number of requested molecules adds up to the total number */
813 req = g->comp[eCompA].nMolReq + g->comp[eCompB].nMolReq;
814 tot = g->comp[eCompA].nMol + g->comp[eCompB].nMol;
818 gmx_fatal(FARGS, "Mismatch of the number of %s ions summed over both compartments.\n"
819 "You requested a total of %d ions (%d in A and %d in B),\n"
820 "but there are a total of %d ions of this type in the system.\n",
821 g->molname, req, g->comp[eCompA].nMolReq,
822 g->comp[eCompB].nMolReq, tot);
825 /* Initialize time-averaging:
826 * Write initial concentrations to all time bins to start with */
827 for (ic = 0; ic < eCompNR; ic++)
829 g->comp[ic].nMolAv = g->comp[ic].nMol;
830 for (i = 0; i < sc->nAverage; i++)
832 g->comp[ic].nMolPast[i] = g->comp[ic].nMol;
839 /*! \brief Copy history of ion counts from checkpoint file.
841 * When called, the checkpoint file has already been read in. Here we copy
842 * over the values from .cpt file to the swap data structure.
844 static void get_initial_ioncounts_from_cpt(
845 t_inputrec *ir, swapstate_t *swapstate,
846 t_commrec *cr, gmx_bool bVerbose)
858 /* Copy the past values from the checkpoint values that have been read in already */
861 fprintf(stderr, "%s Copying values from checkpoint\n", SwS);
864 for (int ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
867 gs = &swapstate->ionType[ig - eSwapFixedGrpNR];
869 for (int ic = 0; ic < eCompNR; ic++)
871 g->comp[ic].nMolReq = gs->nMolReq[ic];
872 g->comp[ic].inflow_net = gs->inflow_net[ic];
876 fprintf(stderr, "%s ... Influx netto: %d Requested: %d Past values: ", SwS,
877 g->comp[ic].inflow_net, g->comp[ic].nMolReq);
880 for (int j = 0; j < sc->nAverage; j++)
882 g->comp[ic].nMolPast[j] = gs->nMolPast[ic][j];
885 fprintf(stderr, "%d ", g->comp[ic].nMolPast[j]);
890 fprintf(stderr, "\n");
898 /*! \brief The master lets all others know about the initial ion counts. */
899 static void bc_initial_concentrations(
910 for (ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
914 for (ic = 0; ic < eCompNR; ic++)
916 gmx_bcast(sizeof(g->comp[ic].nMolReq), &(g->comp[ic].nMolReq), cr);
917 gmx_bcast(sizeof(g->comp[ic].nMol ), &(g->comp[ic].nMol ), cr);
918 gmx_bcast( swap->nAverage * sizeof(g->comp[ic].nMolPast[0]), g->comp[ic].nMolPast, cr);
924 /*! \brief Ensure that each atom belongs to at most one of the swap groups. */
925 static void check_swap_groups(t_swap *s, int nat, gmx_bool bVerbose)
927 int *nGroup = NULL; /* This array counts for each atom in the MD system to
928 how many swap groups it belongs (should be 0 or 1!) */
930 int nMultiple = 0; /* Number of atoms belonging to multiple groups */
935 fprintf(stderr, "%s Making sure each atom belongs to at most one of the swap groups.\n", SwS);
938 /* Add one to the group count of atoms belonging to a swap group: */
940 for (int i = 0; i < s->ngrp; i++)
942 t_swapgrp *g = &s->group[i];
943 for (int j = 0; j < g->nat; j++)
945 /* Get the global index of this atom of this group: */
950 /* Make sure each atom belongs to at most one of the groups: */
951 for (int i = 0; i < nat; i++)
962 gmx_fatal(FARGS, "%s Cannot perform swapping since %d atom%s allocated to more than one swap index group.\n"
963 "%s Each atom must be allocated to at most one of the split groups, the swap groups, or the solvent.\n"
964 "%s Check the .mdp file settings regarding the swap index groups or the index groups themselves.\n",
965 SwS, nMultiple, (1 == nMultiple) ? " is" : "s are", SwSEmpty, SwSEmpty);
970 /*! \brief Get the number of atoms per molecule for this group.
972 * Also ensure that all the molecules in this group have this number of atoms.
974 static int get_group_apm_check(
978 const gmx_mtop_atomlookup_t alook,
981 int molb, molnr, atnr_mol;
982 t_swapgrp *g = &s->group[igroup];
983 int *ind = s->group[igroup].ind;
984 int nat = s->group[igroup].nat;
986 /* Determine the number of solvent atoms per solvent molecule from the
987 * first solvent atom: */
988 gmx_mtop_atomnr_to_molblock_ind(alook, ind[0], &molb, &molnr, &atnr_mol);
989 int apm = mtop->molblock[molb].natoms_mol;
993 fprintf(stderr, "%s Checking whether all %s molecules consist of %d atom%s\n", SwS,
994 g->molname, apm, apm > 1 ? "s" : "");
997 /* Check whether this is also true for all other solvent atoms */
998 for (int i = 1; i < nat; i++)
1000 gmx_mtop_atomnr_to_molblock_ind(alook, ind[i], &molb, &molnr, &atnr_mol);
1001 if (apm != mtop->molblock[molb].natoms_mol)
1003 gmx_fatal(FARGS, "Not all molecules of swap group %d consist of %d atoms.",
1008 //TODO: check whether charges and masses of each molecule are identical!
1013 /*! \brief Print the legend to the swap output file.
1015 * Also print the initial values of ion counts and position of split groups.
1017 static void print_ionlist_legend(t_inputrec *ir,
1018 const gmx_output_env_t *oenv)
1020 const char **legend;
1024 t_swap *s = ir->swap->si_priv;
1025 int nIonTypes = ir->swap->ngrp - eSwapFixedGrpNR;
1026 snew(legend, eCompNR*nIonTypes*3 + 2 + eChanNR*nIonTypes + 1);
1028 // Number of molecules and difference to reference counts for each
1029 // compartment and ion type
1030 for (int ic = count = 0; ic < eCompNR; ic++)
1032 for (int ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
1034 t_swapGroup *g = &ir->swap->grp[ig];
1035 real q = s->group[ig].q;
1037 snprintf(buf, STRLEN, "%s %s ions (charge %s%g)", CompStr[ic], g->molname, q > 0 ? "+" : "", q);
1038 legend[count++] = gmx_strdup(buf);
1040 snprintf(buf, STRLEN, "%s av. mismatch to %d %s ions",
1041 CompStr[ic], s->group[ig].comp[ic].nMolReq, g->molname);
1042 legend[count++] = gmx_strdup(buf);
1044 snprintf(buf, STRLEN, "%s net %s ion influx", CompStr[ic], g->molname);
1045 legend[count++] = gmx_strdup(buf);
1049 // Center of split groups
1050 snprintf(buf, STRLEN, "%scenter of %s of split group 0", SwapStr[ir->eSwapCoords], (NULL != s->group[eGrpSplit0].m) ? "mass" : "geometry");
1051 legend[count++] = gmx_strdup(buf);
1052 snprintf(buf, STRLEN, "%scenter of %s of split group 1", SwapStr[ir->eSwapCoords], (NULL != s->group[eGrpSplit1].m) ? "mass" : "geometry");
1053 legend[count++] = gmx_strdup(buf);
1055 // Ion flux for each channel and ion type
1056 for (int ic = 0; ic < eChanNR; ic++)
1058 for (int ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
1060 t_swapGroup *g = &ir->swap->grp[ig];
1061 snprintf(buf, STRLEN, "A->ch%d->B %s permeations", ic, g->molname);
1062 legend[count++] = gmx_strdup(buf);
1066 // Number of molecules that leaked from A to B
1067 snprintf(buf, STRLEN, "leakage");
1068 legend[count++] = gmx_strdup(buf);
1070 xvgr_legend(s->fpout, count, legend, oenv);
1072 fprintf(s->fpout, "# Instantaneous ion counts and time-averaged differences to requested numbers\n");
1074 // We add a simple text legend helping to identify the columns with xvgr legend strings
1075 fprintf(s->fpout, "# time (ps)");
1076 for (int i = 0; i < count; i++)
1078 snprintf(buf, STRLEN, "s%d", i);
1079 fprintf(s->fpout, "%10s", buf);
1081 fprintf(s->fpout, "\n");
1086 /*! \brief Initialize channel ion flux detection routine.
1088 * Initialize arrays that keep track of where the ions come from and where
1091 static void detect_flux_per_channel_init(
1094 swapstate_t *swapstate,
1095 gmx_bool bStartFromCpt)
1098 swapstateIons_t *gs;
1101 for (int ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
1104 gs = &swapstate->ionType[ig - eSwapFixedGrpNR];
1106 /* All these flux detection routines run on the master only */
1110 g->comp_from = NULL;
1111 g->channel_label = NULL;
1116 /******************************************************/
1117 /* Channel and domain history for the individual ions */
1118 /******************************************************/
1119 if (bStartFromCpt) /* set the pointers right */
1121 g->comp_from = gs->comp_from;
1122 g->channel_label = gs->channel_label;
1124 else /* allocate memory for molecule counts */
1126 snew(g->comp_from, g->nat/g->apm);
1127 gs->comp_from = g->comp_from;
1128 snew(g->channel_label, g->nat/g->apm);
1129 gs->channel_label = g->channel_label;
1131 snew(g->comp_now, g->nat/g->apm);
1133 /* Initialize the channel and domain history counters */
1134 for (int i = 0; i < g->nat/g->apm; i++)
1136 g->comp_now[i] = eDomainNotset;
1139 g->comp_from[i] = eDomainNotset;
1140 g->channel_label[i] = eChHistPassedNone;
1144 /************************************/
1145 /* Channel fluxes for both channels */
1146 /************************************/
1147 g->nCyl[eChan0] = 0;
1148 g->nCyl[eChan1] = 0;
1154 fprintf(stderr, "%s Copying channel fluxes from checkpoint file data\n", SwS);
1158 // Loop over ion types (and both channels)
1159 for (int ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
1162 gs = &swapstate->ionType[ig - eSwapFixedGrpNR];
1164 for (int ic = 0; ic < eChanNR; ic++)
1166 fprintf(stderr, "%s Channel %d flux history for ion type %s (charge %g): ", SwS, ic, g->molname, g->q);
1169 g->fluxfromAtoB[ic] = gs->fluxfromAtoB[ic];
1173 g->fluxfromAtoB[ic] = 0;
1176 fprintf(stderr, "%d molecule%s",
1177 g->fluxfromAtoB[ic], g->fluxfromAtoB[ic] == 1 ? "" : "s");
1178 fprintf(stderr, "\n");
1182 /* Set pointers for checkpoint writing */
1183 swapstate->fluxleak_p = &s->fluxleak;
1184 for (int ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
1187 gs = &swapstate->ionType[ig - eSwapFixedGrpNR];
1189 for (int ic = 0; ic < eChanNR; ic++)
1191 gs->fluxfromAtoB_p[ic] = &g->fluxfromAtoB[ic];
1197 /*! \brief Outputs the initial structure to PDB file for debugging reasons.
1199 * Output the starting structure so that in case of multimeric channels
1200 * the user can check whether we have the correct PBC image for all atoms.
1201 * If this is not correct, the ion counts per channel will be very likely
1204 static void outputStartStructureIfWanted(gmx_mtop_t *mtop, rvec *x, int ePBC, matrix box)
1206 char *env = getenv("GMX_COMPELDUMP");
1210 fprintf(stderr, "\n%s Found env.var. GMX_COMPELDUMP, will output CompEL starting structure made whole.\n"
1211 "%s In case of multimeric channels, please check whether they have the correct PBC representation.\n",
1214 write_sto_conf_mtop("CompELAssumedWholeConfiguration.pdb", *mtop->name, mtop, x, NULL, ePBC, box);
1219 /*! \brief Initialize the swapstate structure, used for checkpoint writing.
1221 * The swapstate struct stores the information we need to make the channels
1222 * whole again after restarts from a checkpoint file. Here we do the following:
1223 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
1224 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
1225 * c) in any case, for subsequent checkpoint writing, we set the pointers in
1226 * swapstate to the x_old arrays, which contain the correct PBC representation of
1227 * multimeric channels at the last time step.
1229 static void init_swapstate(
1230 swapstate_t *swapstate,
1233 rvec x[], /* the initial positions */
1237 rvec *x_pbc = NULL; /* positions of the whole MD system with molecules made whole */
1244 /* We always need the last whole positions such that
1245 * in the next time step we can make the channels whole again in PBC */
1246 if (swapstate->bFromCpt)
1248 /* Copy the last whole positions of each channel from .cpt */
1249 g = &(s->group[eGrpSplit0]);
1250 for (int i = 0; i < g->nat; i++)
1252 copy_rvec(swapstate->xc_old_whole[eChan0][i], g->xc_old[i]);
1254 g = &(s->group[eGrpSplit1]);
1255 for (int i = 0; i < g->nat; i++)
1257 copy_rvec(swapstate->xc_old_whole[eChan1][i], g->xc_old[i]);
1262 /* Set the number of ion types and allocate memory for checkpointing */
1263 swapstate->nIonTypes = s->ngrp - eSwapFixedGrpNR;
1264 snew(swapstate->ionType, swapstate->nIonTypes);
1266 /* Store the total number of ions of each type in the swapstateIons
1267 * structure that is accessible during checkpoint writing */
1268 for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1270 swapstateIons_t *gs = &swapstate->ionType[ii];
1271 gs->nMol = sc->grp[ii + eSwapFixedGrpNR].nat;
1274 /* Extract the initial split group positions. */
1276 /* Remove pbc, make molecule whole. */
1277 snew(x_pbc, mtop->natoms);
1278 copy_rvecn(x, x_pbc, 0, mtop->natoms);
1280 /* This can only make individual molecules whole, not multimers */
1281 do_pbc_mtop(NULL, ePBC, box, mtop, x_pbc);
1283 /* Output the starting structure? */
1284 outputStartStructureIfWanted(mtop, x_pbc, ePBC, box);
1286 /* If this is the first run (i.e. no checkpoint present) we assume
1287 * that the starting positions give us the correct PBC representation */
1288 for (int ig = eGrpSplit0; ig <= eGrpSplit1; ig++)
1290 g = &(s->group[ig]);
1291 for (int i = 0; i < g->nat; i++)
1293 copy_rvec(x_pbc[g->ind[i]], g->xc_old[i]);
1298 /* Prepare swapstate arrays for later checkpoint writing */
1299 swapstate->nat[eChan0] = s->group[eGrpSplit0].nat;
1300 swapstate->nat[eChan1] = s->group[eGrpSplit1].nat;
1303 /* For subsequent checkpoint writing, set the swapstate pointers to the xc_old
1304 * arrays that get updated at every swapping step */
1305 swapstate->xc_old_whole_p[eChan0] = &s->group[eGrpSplit0].xc_old;
1306 swapstate->xc_old_whole_p[eChan1] = &s->group[eGrpSplit1].xc_old;
1309 /*! \brief Determine the total charge imbalance resulting from the swap groups */
1310 static real getRequestedChargeImbalance(t_swap *s)
1315 real particle_charge;
1316 real particle_number[eCompNR];
1318 // s->deltaQ = ( (-1) * s->comp[eCompA][eIonNEG].nat_req + s->comp[eCompA][eIonPOS].nat_req )
1319 // - ( (-1) * s->comp[eCompB][eIonNEG].nat_req + s->comp[eCompB][eIonPOS].nat_req );
1321 for (ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
1325 particle_charge = g->q;
1326 particle_number[eCompA] = g->comp[eCompA].nMolReq;
1327 particle_number[eCompB] = g->comp[eCompB].nMolReq;
1329 DeltaQ += particle_charge * (particle_number[eCompA] - particle_number[eCompB]);
1336 /*! \brief Sorts anions and cations into two separate groups
1338 * This routine should be called for the 'anions' and 'cations' group,
1339 * of which the indices were lumped together in the older version of the code.
1341 void copyIndicesToGroup(
1349 /* If explicit ion counts were requested in the .mdp file
1350 * (by setting positive values for the number of ions),
1351 * we can make an additional consistency check here */
1352 if ( (g->nmolReq[eCompA] < 0) && (g->nmolReq[eCompB] < 0) )
1354 if (g->nat != (g->nmolReq[eCompA] + g->nmolReq[eCompB]) )
1356 gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
1357 "%s Inconsistency while importing swap-related data from an old input file version.\n"
1358 "%s The requested ion counts in compartments A (%d) and B (%d)\n"
1359 "%s do not add up to the number of ions (%d) of this type for the group '%s'.\n",
1360 SwS, SwSEmpty, g->nmolReq[eCompA], g->nmolReq[eCompB], SwSEmpty, g->nat, g->molname);
1364 srenew(g->ind, g->nat);
1365 for (int i = 0; i < g->nat; i++)
1367 g->ind[i] = indIons[i];
1372 /*! \brief Converts old .tpr file CompEL contents to new data layout.
1374 * If we have read an old .tpr file (tpxv <= tpxv_CompElPolyatomicIonsAndMultipleIonTypes),
1375 * anions and cations are stored together in group #3. In the new
1376 * format we store each ion type in a separate group.
1377 * The 'classic' groups are:
1378 * #0 split group 0 - OK
1379 * #1 split group 1 - OK
1381 * #3 anions - contains also cations, needs to be converted
1382 * #4 cations - empty before conversion
1385 void convertOldToNewGroupFormat(
1392 gmx_mtop_atomlookup_t alook = gmx_mtop_atomlookup_init(mtop);
1393 t_swapGroup *g = &sc->grp[3];
1395 /* Loop through the atom indices of group #3 (anions) and put all indices
1396 * that belong to cations into the cation group.
1400 int *indAnions = NULL;
1401 int *indCations = NULL;
1402 snew(indAnions, g->nat);
1403 snew(indCations, g->nat);
1405 for (int i = 0; i < g->nat; i++)
1407 gmx_mtop_atomnr_to_atom(alook, g->ind[i], &atom);
1410 // This is an anion, add it to the list of anions
1411 indAnions[nAnions++] = g->ind[i];
1415 // This is a cation, add it to the list of cations
1416 indCations[nCations++] = g->ind[i];
1422 fprintf(stdout, "%s Sorted %d ions into separate groups of %d anions and %d cations.\n",
1423 SwS, g->nat, nAnions, nCations);
1427 /* Now we have the correct lists of anions and cations.
1428 * Copy it to the right groups.
1430 copyIndicesToGroup(indAnions, nAnions, g, cr);
1432 copyIndicesToGroup(indCations, nCations, g, cr);
1440 /*! \brief Returns TRUE if we started from an old .tpr
1442 * Then we need to re-sort anions and cations into separate groups */
1443 gmx_bool bConvertFromOldTpr(t_swapcoords *sc)
1445 // If the last group has no atoms it means we need to convert!
1446 if ( (sc->ngrp >= 5) && (0 == sc->grp[4].nat) )
1454 void init_swapcoords(
1462 swapstate_t *swapstate,
1464 const gmx_output_env_t *oenv,
1465 unsigned long Flags)
1471 swapstateIons_t *gs;
1472 gmx_bool bAppend, bStartFromCpt, bRerun;
1473 gmx_mtop_atomlookup_t alook = NULL;
1475 alook = gmx_mtop_atomlookup_init(mtop);
1477 if ( (PAR(cr)) && !DOMAINDECOMP(cr) )
1479 gmx_fatal(FARGS, "Position swapping is only implemented for domain decomposition!");
1482 bAppend = Flags & MD_APPENDFILES;
1483 bStartFromCpt = Flags & MD_STARTFROMCPT;
1484 bRerun = Flags & MD_RERUN;
1487 snew(sc->si_priv, 1);
1494 gmx_fatal(FARGS, "%s This module does not support reruns in parallel\nPlease request a serial run with -nt 1 / -np 1\n", SwS);
1497 fprintf(stderr, "%s Rerun - using every available frame\n", SwS);
1499 sc->nAverage = 1; /* averaging makes no sense for reruns */
1502 if (MASTER(cr) && !bAppend)
1504 fprintf(fplog, "\nInitializing ion/water position exchanges\n");
1505 please_cite(fplog, "Kutzner2011b");
1508 switch (ir->eSwapCoords)
1524 // For compatibility with old .tpr files
1525 if (bConvertFromOldTpr(sc) )
1527 convertOldToNewGroupFormat(sc, mtop, bVerbose && MASTER(cr), cr);
1530 /* Copy some data and pointers to the group structures for convenience */
1531 /* Number of atoms in the group */
1533 snew(s->group, s->ngrp);
1534 for (int i = 0; i < s->ngrp; i++)
1536 s->group[i].nat = sc->grp[i].nat;
1537 s->group[i].ind = sc->grp[i].ind;
1538 s->group[i].molname = sc->grp[i].molname;
1541 /* Check for overlapping atoms */
1542 check_swap_groups(s, mtop->natoms, bVerbose && MASTER(cr));
1544 /* Allocate space for the collective arrays for all groups */
1545 /* For the collective position array */
1546 for (int i = 0; i < s->ngrp; i++)
1549 snew(g->xc, g->nat);
1550 snew(g->c_ind_loc, g->nat);
1552 /* For the split groups (the channels) we need some extra memory to
1553 * be able to make the molecules whole even if they span more than
1554 * half of the box size. */
1555 if ( (i == eGrpSplit0) || (i == eGrpSplit1) )
1557 snew(g->xc_shifts, g->nat);
1558 snew(g->xc_eshifts, g->nat);
1559 snew(g->xc_old, g->nat);
1565 init_swapstate(swapstate, sc, mtop, x, box, ir->ePBC);
1568 /* After init_swapstate we have a set of (old) whole positions for our
1569 * channels. Now transfer that to all nodes */
1572 for (int ig = eGrpSplit0; ig <= eGrpSplit1; ig++)
1574 g = &(s->group[ig]);
1575 gmx_bcast((g->nat)*sizeof((g->xc_old)[0]), g->xc_old, (cr));
1579 /* Make sure that all molecules in the solvent and ion groups contain the
1580 * same number of atoms each */
1581 for (int ig = eGrpSolvent; ig < s->ngrp; ig++)
1585 g = &(s->group[ig]);
1586 g->apm = get_group_apm_check(ig, s, MASTER(cr) && bVerbose, alook, mtop);
1588 /* Since all molecules of a group are equal, we only need enough space
1589 * to determine properties of a single molecule at at time */
1590 snew(g->m, g->apm); /* For the center of mass */
1591 charge = 0; /* To determine the charge imbalance */
1592 for (int j = 0; j < g->apm; j++)
1594 gmx_mtop_atomnr_to_atom(alook, g->ind[j], &atom);
1598 /* Total charge of one molecule of this group: */
1603 /* Need mass-weighted center of split group? */
1604 for (int j = eGrpSplit0; j <= eGrpSplit1; j++)
1607 if (TRUE == sc->massw_split[j])
1609 /* Save the split group masses if mass-weighting is requested */
1611 for (int i = 0; i < g->nat; i++)
1613 gmx_mtop_atomnr_to_atom(alook, g->ind[i], &atom);
1619 /* Make a t_pbc struct on all nodes so that the molecules
1620 * chosen for an exchange can be made whole. */
1627 fprintf(stderr, "%s Opening output file %s%s\n", SwS, fn, bAppend ? " for appending" : "");
1630 s->fpout = gmx_fio_fopen(fn, bAppend ? "a" : "w" );
1634 xvgr_header(s->fpout, "Molecule counts", "Time (ps)", "counts", exvggtXNY, oenv);
1636 for (int ig = 0; ig < s->ngrp; ig++)
1638 g = &(s->group[ig]);
1639 fprintf(s->fpout, "# %s group '%s' contains %d atom%s",
1640 ig < eSwapFixedGrpNR ? eSwapFixedGrp_names[ig] : "Ion",
1641 g->molname, g->nat, (g->nat > 1) ? "s" : "");
1642 if (!(eGrpSplit0 == ig || eGrpSplit1 == ig) )
1644 fprintf(s->fpout, " with %d atom%s in each molecule of charge %g",
1645 g->apm, (g->apm > 1) ? "s" : "", g->q);
1647 fprintf(s->fpout, ".\n");
1650 fprintf(s->fpout, "#\n# Initial positions of split groups:\n");
1653 for (int j = eGrpSplit0; j <= eGrpSplit1; j++)
1656 for (int i = 0; i < g->nat; i++)
1658 copy_rvec(x[sc->grp[j].ind[i]], g->xc[i]);
1660 /* xc has the correct PBC representation for the two channels, so we do
1661 * not need to correct for that */
1662 get_center(g->xc, g->m, g->nat, g->center);
1665 fprintf(s->fpout, "# %s group %s-center %5f nm\n", eSwapFixedGrp_names[j],
1666 DimStr[s->swapdim], g->center[s->swapdim]);
1672 if ( (0 != sc->bulkOffset[eCompA]) || (0 != sc->bulkOffset[eCompB]) )
1674 fprintf(s->fpout, "#\n");
1675 fprintf(s->fpout, "# You provided an offset for the position of the bulk layer(s).\n");
1676 fprintf(s->fpout, "# That means the layers to/from which ions and water molecules are swapped\n");
1677 fprintf(s->fpout, "# are not midway (= at 0.0) between the compartment-defining layers (at +/- 1.0).\n");
1678 fprintf(s->fpout, "# bulk-offsetA = %g\n", sc->bulkOffset[eCompA]);
1679 fprintf(s->fpout, "# bulk-offsetB = %g\n", sc->bulkOffset[eCompB]);
1682 fprintf(s->fpout, "#\n");
1683 fprintf(s->fpout, "# Split0 cylinder radius %f nm, up %f nm, down %f nm\n",
1684 sc->cyl0r, sc->cyl0u, sc->cyl0l);
1685 fprintf(s->fpout, "# Split1 cylinder radius %f nm, up %f nm, down %f nm\n",
1686 sc->cyl1r, sc->cyl1u, sc->cyl1l);
1688 fprintf(s->fpout, "#\n");
1691 fprintf(s->fpout, "# Coupling constant (number of swap attempt steps to average over): %d (translates to %f ps).\n",
1692 sc->nAverage, sc->nAverage*sc->nstswap*ir->delta_t);
1693 fprintf(s->fpout, "# Threshold is %f\n", sc->threshold);
1694 fprintf(s->fpout, "#\n");
1695 fprintf(s->fpout, "# Remarks about which atoms passed which channel use global atoms numbers starting at one.\n");
1704 /* Prepare for parallel or serial run */
1707 for (int ig = 0; ig < s->ngrp; ig++)
1709 g = &(s->group[ig]);
1717 for (int ig = 0; ig < s->ngrp; ig++)
1719 g = &(s->group[ig]);
1720 g->nat_loc = g->nat;
1721 g->ind_loc = g->ind;
1722 /* c_ind_loc needs to be set to identity in the serial case */
1723 for (int i = 0; i < g->nat; i++)
1725 g->c_ind_loc[i] = i;
1730 /* Allocate memory to remember the past particle counts for time averaging */
1731 for (int ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
1733 g = &(s->group[ig]);
1734 for (int ic = 0; ic < eCompNR; ic++)
1736 snew(g->comp[ic].nMolPast, sc->nAverage);
1740 /* Get the initial particle concentrations and let the other nodes know */
1745 get_initial_ioncounts_from_cpt(ir, swapstate, cr, bVerbose);
1749 fprintf(stderr, "%s Determining initial numbers of ions per compartment.\n", SwS);
1750 get_initial_ioncounts(ir, x, box, cr, bRerun);
1753 /* Prepare (further) checkpoint writes ... */
1756 /* Consistency check */
1757 if (swapstate->nAverage != sc->nAverage)
1759 gmx_fatal(FARGS, "%s Ion count averaging steps mismatch! checkpoint: %d, tpr: %d",
1760 SwS, swapstate->nAverage, sc->nAverage);
1765 swapstate->nAverage = sc->nAverage;
1767 fprintf(stderr, "%s Setting pointers for checkpoint writing\n", SwS);
1768 for (int ic = 0; ic < eCompNR; ic++)
1770 for (int ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
1773 gs = &swapstate->ionType[ig - eSwapFixedGrpNR];
1775 gs->nMolReq_p[ic] = &(g->comp[ic].nMolReq);
1776 gs->nMolPast_p[ic] = &(g->comp[ic].nMolPast[0]);
1777 gs->inflow_net_p[ic] = &(g->comp[ic].inflow_net);
1781 /* Determine the total charge imbalance */
1782 s->deltaQ = getRequestedChargeImbalance(s);
1786 fprintf(stderr, "%s Requested charge imbalance is Q(A) - Q(B) = %g e.\n", SwS, s->deltaQ);
1790 fprintf(s->fpout, "# Requested charge imbalance is Q(A)-Q(B) = %g e.\n", s->deltaQ);
1796 bc_initial_concentrations(cr, ir->swap);
1799 /* Update the time-averaged number of molecules for all groups and compartments */
1800 for (int ig = eSwapFixedGrpNR; ig < sc->ngrp; ig++)
1803 for (int ic = 0; ic < eCompNR; ic++)
1805 update_time_window(&g->comp[ic], sc->nAverage, -1);
1809 /* Initialize arrays that keep track of through which channel the ions go */
1810 detect_flux_per_channel_init(cr, s, swapstate, bStartFromCpt);
1812 /* We need to print the legend if we open this file for the first time. */
1813 if (MASTER(cr) && !bAppend)
1815 print_ionlist_legend(ir, oenv);
1820 void finish_swapcoords(t_swapcoords *sc)
1822 if (sc->si_priv->fpout)
1824 // Close the swap output file
1825 gmx_fio_fclose(sc->si_priv->fpout);
1830 void dd_make_local_swap_groups(gmx_domdec_t *dd, t_swapcoords *sc)
1836 /* Make split groups, solvent group, and user-defined groups of particles
1837 * under control of the swap protocol */
1838 for (ig = 0; ig < sc->ngrp; ig++)
1840 g = &(sc->si_priv->group[ig]);
1841 dd_make_local_group_indices(dd->ga2la, g->nat, g->ind,
1842 &(g->nat_loc), &(g->ind_loc), &(g->nalloc_loc), g->c_ind_loc);
1847 /*! \brief Do we need to swap a molecule in any of the ion groups with a water molecule at this step?
1849 * From the requested and average molecule counts we determine whether a swap is needed
1850 * at this time step.
1852 static gmx_bool need_swap(t_swapcoords *sc)
1860 for (ig = eSwapFixedGrpNR; ig < sc->ngrp; ig++)
1864 for (ic = 0; ic < eCompNR; ic++)
1866 if (g->comp[ic].nMolReq - g->comp[ic].nMolAv >= sc->threshold)
1876 /*! \brief Return the index of an atom or molecule suitable for swapping.
1878 * Returns the index of an atom that is far off the compartment boundaries,
1879 * that is near to the bulk layer to/from which the swaps take place.
1880 * Other atoms of the molecule (if any) will directly follow the returned index.
1882 * \param[in] comp Structure containing compartment-specific data.
1883 * \param[in] molname Name of the molecule.
1885 * \returns Index of the first atom of the molecule chosen for a position exchange.
1887 static int get_index_of_distant_atom(
1888 t_compartment *comp,
1889 const char molname[])
1892 real d = GMX_REAL_MAX;
1895 /* comp->nat contains the original number of atoms in this compartment
1896 * prior to doing any swaps. Some of these atoms may already have been
1897 * swapped out, but then they are marked with a distance of GMX_REAL_MAX
1899 for (int iMol = 0; iMol < comp->nMolBefore; iMol++)
1901 if (comp->dist[iMol] < d)
1904 d = comp->dist[ibest];
1910 gmx_fatal(FARGS, "Could not get index of %s atom. Compartment contains %d %s molecules before swaps.",
1911 molname, comp->nMolBefore, molname);
1914 /* Set the distance of this index to infinity such that it won't get selected again in
1917 comp->dist[ibest] = GMX_REAL_MAX;
1919 return comp->ind[ibest];
1923 /*! \brief Swaps centers of mass and makes molecule whole if broken */
1924 static void translate_positions(
1932 rvec reference, dx, correctPBCimage;
1935 /* Use the first atom as the reference for PBC */
1936 copy_rvec(x[0], reference);
1938 for (i = 0; i < apm; i++)
1940 /* PBC distance between position and reference */
1941 pbc_dx(pbc, x[i], reference, dx);
1943 /* Add PBC distance to reference */
1944 rvec_add(reference, dx, correctPBCimage);
1946 /* Subtract old_com from correct image and add new_com */
1947 rvec_dec(correctPBCimage, old_com);
1948 rvec_inc(correctPBCimage, new_com);
1950 copy_rvec(correctPBCimage, x[i]);
1955 /*! \brief Write back the the modified local positions from the collective array to the official positions. */
1956 static void apply_modified_positions(
1963 for (l = 0; l < g->nat_loc; l++)
1965 /* Get the right local index to write to */
1967 /* Where is the local atom in the collective array? */
1968 cind = g->c_ind_loc[l];
1970 /* Copy the possibly modified position */
1971 copy_rvec(g->xc[cind], x[ii]);
1976 gmx_bool do_swapcoords(
1981 gmx_wallcycle_t wcycle,
1990 int j, ic, ig, nswaps;
1991 int thisC, otherC; /* Index into this compartment and the other one */
1992 gmx_bool bSwap = FALSE;
1993 t_swapgrp *g, *gsol;
1995 rvec com_solvent, com_particle; /* solvent and swap molecule's center of mass */
1996 gmx_mtop_atomlookup_t alook = NULL;
1999 wallcycle_start(wcycle, ewcSWAP);
2004 set_pbc(s->pbc, ir->ePBC, box);
2006 /* Assemble the positions of the split groups, i.e. the channels.
2007 * Here we also pass a shifts array to communicate_group_positions(), so that it can make
2008 * the molecules whole even in cases where they span more than half of the box in
2010 for (ig = eGrpSplit0; ig <= eGrpSplit1; ig++)
2012 g = &(s->group[ig]);
2013 communicate_group_positions(cr, g->xc, g->xc_shifts, g->xc_eshifts, TRUE,
2014 x, g->nat, g->nat_loc, g->ind_loc, g->c_ind_loc, g->xc_old, box);
2016 get_center(g->xc, g->m, g->nat, g->center); /* center of split groups == channels */
2019 /* Assemble the positions of the ions (ig = 3, 4, ...). These molecules should
2020 * be small and we can always make them whole with a simple distance check.
2021 * Therefore we pass NULL as third argument. */
2022 for (ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
2024 g = &(s->group[ig]);
2025 communicate_group_positions(cr, g->xc, NULL, NULL, FALSE,
2026 x, g->nat, g->nat_loc, g->ind_loc, g->c_ind_loc, NULL, NULL);
2028 /* Determine how many ions of this type each compartment contains */
2029 sortMoleculesIntoCompartments(g, cr, sc, box, step, s->fpout, bRerun, FALSE);
2032 /* Output how many ions are in the compartments */
2035 print_ionlist(s, t, "");
2038 /* If we are doing a rerun, we are finished here, since we cannot perform
2045 /* Do we have to perform a swap? */
2046 bSwap = need_swap(sc);
2049 /* Since we here know that we have to perform ion/water position exchanges,
2050 * we now assemble the solvent positions */
2051 g = &(s->group[eGrpSolvent]);
2052 communicate_group_positions(cr, g->xc, NULL, NULL, FALSE,
2053 x, g->nat, g->nat_loc, g->ind_loc, g->c_ind_loc, NULL, NULL);
2055 /* Determine how many molecules of solvent each compartment contains */
2056 sortMoleculesIntoCompartments(g, cr, sc, box, step, s->fpout, bRerun, TRUE);
2058 /* Save number of solvent molecules per compartment prior to any swaps */
2059 g->comp[eCompA].nMolBefore = g->comp[eCompA].nMol;
2060 g->comp[eCompB].nMolBefore = g->comp[eCompB].nMol;
2062 for (ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
2064 g = &(s->group[ig]);
2066 for (ic = 0; ic < eCompNR; ic++)
2068 /* Determine in which compartment ions are missing and where they are too many */
2069 g->vacancy[ic] = g->comp[ic].nMolReq - g->comp[ic].nMolAv;
2071 /* Save number of ions per compartment prior to swaps */
2072 g->comp[ic].nMolBefore = g->comp[ic].nMol;
2076 /* Now actually perform the particle exchanges, one swap group after another */
2077 alook = gmx_mtop_atomlookup_init(mtop);
2078 gsol = &s->group[eGrpSolvent];
2079 for (ig = eSwapFixedGrpNR; ig < s->ngrp; ig++)
2083 for (thisC = 0; thisC < eCompNR; thisC++)
2085 /* Index to the other compartment */
2086 otherC = (thisC+1) % eCompNR;
2088 while (g->vacancy[thisC] >= sc->threshold)
2090 /* Swap in an ion */
2092 /* Get the xc-index of the first atom of a solvent molecule of this compartment */
2093 isol = get_index_of_distant_atom(&gsol->comp[thisC], gsol->molname);
2095 /* Get the xc-index of a particle from the other compartment */
2096 iion = get_index_of_distant_atom(&g->comp[otherC], g->molname);
2098 get_molecule_center(&gsol->xc[isol], gsol->apm, gsol->m, com_solvent, s->pbc);
2099 get_molecule_center(&g->xc[iion], g->apm, g->m, com_particle, s->pbc);
2101 /* Subtract solvent molecule's center of mass and add swap particle's center of mass */
2102 translate_positions(&gsol->xc[isol], gsol->apm, com_solvent, com_particle, s->pbc);
2103 /* Similarly for the swap particle, subtract com_particle and add com_solvent */
2104 translate_positions(&g->xc[iion], g->apm, com_particle, com_solvent, s->pbc);
2106 /* Keep track of the changes */
2107 g->vacancy[thisC ]--;
2108 g->vacancy[otherC]++;
2109 g->comp [thisC ].nMol++;
2110 g->comp [otherC].nMol--;
2111 g->comp [thisC ].inflow_net++;
2112 g->comp [otherC].inflow_net--;
2113 /* Correct the past time window to still get the right averages from now on */
2114 g->comp [thisC ].nMolAv++;
2115 g->comp [otherC].nMolAv--;
2116 for (j = 0; j < sc->nAverage; j++)
2118 g->comp[thisC ].nMolPast[j]++;
2119 g->comp[otherC].nMolPast[j]--;
2121 /* Clear ion history */
2124 int iMol = iion / g->apm;
2125 g->channel_label[iMol] = eChHistPassedNone;
2126 g->comp_from[iMol] = eDomainNotset;
2128 /* That was the swap */
2133 if (nswaps && bVerbose)
2135 fprintf(stderr, "%s Performed %d swap%s in step %" GMX_PRId64 " for iontype %s.\n",
2136 SwS, nswaps, nswaps > 1 ? "s" : "", step, g->molname);
2139 gmx_mtop_atomlookup_destroy(alook);
2141 if (s->fpout != NULL)
2143 print_ionlist(s, t, " # after swap");
2146 /* For the solvent and user-defined swap groups, each rank writes back its
2147 * (possibly modified) local positions to the official position array. */
2148 for (ig = eGrpSolvent; ig < s->ngrp; ig++)
2151 apply_modified_positions(g, x);
2154 } /* end of if(bSwap) */
2156 wallcycle_stop(wcycle, ewcSWAP);