2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
5 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
38 * Implements functions in swapcoords.h.
40 * \author Carsten Kutzner <ckutzne@gwdg.de>
41 * \ingroup module_swap
45 #include "gromacs/utility/enumerationhelpers.h"
46 #include "swapcoords.h"
56 #include "gromacs/domdec/domdec_struct.h"
57 #include "gromacs/domdec/localatomset.h"
58 #include "gromacs/domdec/localatomsetmanager.h"
59 #include "gromacs/fileio/confio.h"
60 #include "gromacs/fileio/gmxfio.h"
61 #include "gromacs/fileio/xvgr.h"
62 #include "gromacs/gmxlib/network.h"
63 #include "gromacs/math/vec.h"
64 #include "gromacs/mdlib/groupcoord.h"
65 #include "gromacs/mdrunutility/handlerestart.h"
66 #include "gromacs/mdtypes/commrec.h"
67 #include "gromacs/mdtypes/imdmodule.h"
68 #include "gromacs/mdtypes/inputrec.h"
69 #include "gromacs/mdtypes/md_enums.h"
70 #include "gromacs/mdtypes/mdrunoptions.h"
71 #include "gromacs/mdtypes/observableshistory.h"
72 #include "gromacs/mdtypes/state.h"
73 #include "gromacs/mdtypes/swaphistory.h"
74 #include "gromacs/pbcutil/pbc.h"
75 #include "gromacs/timing/wallcycle.h"
76 #include "gromacs/topology/mtop_lookup.h"
77 #include "gromacs/topology/topology.h"
78 #include "gromacs/utility/cstringutil.h"
79 #include "gromacs/utility/fatalerror.h"
80 #include "gromacs/utility/pleasecite.h"
81 #include "gromacs/utility/smalloc.h"
82 #include "gromacs/utility/snprintf.h"
84 static const std::string SwS = { "SWAP:" }; /**< For output that comes from the swap module */
85 static const std::string SwSEmpty = { " " }; /**< Placeholder for multi-line output */
86 static constexpr gmx::EnumerationArray<Compartment, const char*> CompStr = { "A", "B" }; /**< Compartment name */
87 static constexpr gmx::EnumerationArray<SwapType, const char*> SwapStr = { "", "X-", "Y-", "Z-" }; /**< Name for the swap types. */
88 static const char* const DimStr[DIM + 1] = { "X", "Y", "Z", nullptr }; /**< Name for the swap dimension. */
90 /** Keep track of through which channel the ions have passed */
91 enum class ChannelHistory : int
98 static constexpr gmx::EnumerationArray<ChannelHistory, const char*> ChannelString = { "none",
100 "channel1" }; /**< Name for the channels */
102 /*! \brief Domain identifier.
104 * Keeps track of from which compartment the ions came before passing the
107 enum class Domain : int
114 static constexpr gmx::EnumerationArray<Domain, const char*> DomainString = { "not_assigned",
116 "Domain_B" }; /**< Name for the domains */
121 extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
124 * \brief Implement Computational Electrophysiology swapping.
126 class SwapCoordinates final : public IMDModule
129 IMdpOptionProvider* mdpOptionProvider() override { return nullptr; }
130 IMDOutputProvider* outputProvider() override { return nullptr; }
131 void initForceProviders(ForceProviders* /* forceProviders */) override {}
132 void subscribeToSimulationSetupNotifications(MDModulesNotifiers* /* notifiers */) override {}
133 void subscribeToPreProcessingNotifications(MDModulesNotifiers* /* notifiers */) override {}
136 std::unique_ptr<IMDModule> createSwapCoordinatesModule()
138 return std::make_unique<SwapCoordinates>();
145 * Structure containing compartment-specific data.
147 typedef struct swap_compartment
149 int nMol; /**< Number of ion or water molecules detected
150 in this compartment. */
151 int nMolBefore; /**< Number of molecules before swapping. */
152 int nMolReq; /**< Requested number of molecules in compartment. */
153 real nMolAv; /**< Time-averaged number of molecules matching
154 the compartment conditions. */
155 int* nMolPast; /**< Past molecule counts for time-averaging. */
156 int* ind; /**< Indices to collective array of atoms. */
157 real* dist; /**< Distance of atom to bulk layer, which is
158 normally the center layer of the compartment */
159 int nalloc; /**< Allocation size for ind array. */
160 int inflow_net; /**< Net inflow of ions into this compartment. */
165 * This structure contains data needed for the groups involved in swapping:
166 * split group 0, split group 1, solvent group, ion groups.
168 typedef struct swap_group
170 /*!\brief Construct a swap group given the managed swap atoms.
172 * \param[in] atomset Managed indices of atoms that are part of the swap group.
174 swap_group(const gmx::LocalAtomSet& atomset);
175 char* molname = nullptr; /**< Name of the group or ion type */
176 int apm = 0; /**< Number of atoms in each molecule */
177 gmx::LocalAtomSet atomset; /**< The atom indices in the swap group */
178 rvec* xc = nullptr; /**< Collective array of group atom positions (size nat) */
179 ivec* xc_shifts = nullptr; /**< Current (collective) shifts (size nat) */
180 ivec* xc_eshifts = nullptr; /**< Extra shifts since last DD step (size nat) */
181 rvec* xc_old = nullptr; /**< Old (collective) positions (size nat) */
182 real q = 0.; /**< Total charge of one molecule of this group */
183 real* m = nullptr; /**< Masses (can be omitted, size apm) */
184 Domain* comp_from = nullptr; /**< (Collective) Stores from which compartment this
185 molecule has come. This way we keep track of
186 through which channel an ion permeates
187 (size nMol = nat/apm) */
188 Domain* comp_now = nullptr; /**< In which compartment this ion is now (size nMol) */
189 ChannelHistory* channel_label = nullptr; /**< Which channel was passed at last by this ion?
191 rvec center; /**< Center of the group; COM if masses are used */
192 gmx::EnumerationArray<Compartment, t_compartment> comp; /**< Distribution of particles of this
193 group across the two compartments */
194 gmx::EnumerationArray<Compartment, real> vacancy; /**< How many molecules need to be swapped in? */
195 gmx::EnumerationArray<Channel, int> fluxfromAtoB; /**< Net flux of ions per channel */
196 gmx::EnumerationArray<Channel, int> nCyl; /**< Number of ions residing in a channel */
197 int nCylBoth = 0; /**< Ions assigned to cyl0 and cyl1. Not good. */
200 t_swapgrp::swap_group(const gmx::LocalAtomSet& atomset) : atomset{ atomset }
205 for (auto compartment : keysOf(comp))
207 comp[compartment] = {};
208 vacancy[compartment] = 0;
210 for (auto channel : keysOf(fluxfromAtoB))
212 fluxfromAtoB[channel] = 0;
218 * Main (private) data structure for the position swapping protocol.
222 int swapdim; /**< One of XX, YY, ZZ */
223 t_pbc* pbc; /**< Needed to make molecules whole. */
224 FILE* fpout; /**< Output file. */
225 int ngrp; /**< Number of t_swapgrp groups */
226 std::vector<t_swapgrp> group; /**< Separate groups for channels, solvent, ions */
227 int fluxleak; /**< Flux not going through any of the channels. */
228 real deltaQ; /**< The charge imbalance between the compartments. */
232 /*! \brief Check whether point is in channel.
234 * A channel is a cylinder defined by a disc
235 * with radius r around its center c. The thickness of the cylinder is
242 * <---------c--------->
248 * \param[in] point The position (xyz) under consideration.
249 * \param[in] center The center of the cylinder.
250 * \param[in] d_up The upper extension of the cylinder.
251 * \param[in] d_down The lower extension.
252 * \param[in] r_cyl2 Cylinder radius squared.
253 * \param[in] pbc Structure with info about periodic boundary conditions.
254 * \param[in] normal The membrane normal direction is typically 3, i.e. z, but can be x or y also.
256 * \returns Whether the point is inside the defined cylindric channel.
258 static gmx_bool is_in_channel(rvec point, rvec center, real d_up, real d_down, real r_cyl2, t_pbc* pbc, int normal)
261 int plane1, plane2; /* Directions tangential to membrane */
264 plane1 = (normal + 1) % 3; /* typically 0, i.e. XX */
265 plane2 = (normal + 2) % 3; /* typically 1, i.e. YY */
267 /* Get the distance vector dr between the point and the center of the cylinder */
268 pbc_dx(pbc, point, center, dr); /* This puts center in the origin */
270 /* Check vertical direction */
271 if ((dr[normal] > d_up) || (dr[normal] < -d_down))
276 /* Check radial direction */
277 if ((dr[plane1] * dr[plane1] + dr[plane2] * dr[plane2]) > r_cyl2)
282 /* All check passed, this point is in the cylinder */
287 /*! \brief Prints output to CompEL output file.
289 * Prints to swap output file how many ions are in each compartment,
290 * where the centers of the split groups are, and how many ions of each type
291 * passed the channels.
293 static void print_ionlist(t_swap* s, double time, const char comment[])
296 fprintf(s->fpout, "%12.5e", time);
298 // Output number of molecules and difference to reference counts for each
299 // compartment and ion type
300 for (auto iComp : gmx::EnumerationWrapper<Compartment>{})
302 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
304 t_compartment* comp = &s->group[ig].comp[iComp];
306 fprintf(s->fpout, "%10d%10.1f%10d", comp->nMol, comp->nMolAv - comp->nMolReq, comp->inflow_net);
310 // Output center of split groups
313 s->group[static_cast<int>(SwapGroupSplittingType::Split0)].center[s->swapdim],
314 s->group[static_cast<int>(SwapGroupSplittingType::Split1)].center[s->swapdim]);
316 // Output ion flux for each channel and ion type
317 for (auto iChan : gmx::EnumerationWrapper<Channel>{})
319 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
321 t_swapgrp* g = &s->group[ig];
322 fprintf(s->fpout, "%10d", g->fluxfromAtoB[iChan]);
326 /* Output the number of molecules that leaked from A to B */
327 fprintf(s->fpout, "%10d", s->fluxleak);
329 fprintf(s->fpout, "%s\n", comment);
333 /*! \brief Get the center of a group of nat atoms.
335 * Since with PBC an atom group might not be whole, use the first atom as the
336 * reference atom and determine the center with respect to this reference.
338 static void get_molecule_center(rvec x[], int nat, const real* weights, rvec center, t_pbc* pbc)
341 rvec weightedPBCimage;
343 rvec reference, correctPBCimage, dx;
346 /* Use the first atom as the reference and put other atoms near that one */
347 /* This does not work for large molecules that span > half of the box! */
348 copy_rvec(x[0], reference);
350 /* Calculate either the weighted center or simply the center of geometry */
353 for (i = 0; i < nat; i++)
355 /* PBC distance between position and reference */
356 pbc_dx(pbc, x[i], reference, dx);
358 /* Add PBC distance to reference */
359 rvec_add(reference, dx, correctPBCimage);
361 /* Take weight into account */
362 if (nullptr == weights)
371 svmul(wi, correctPBCimage, weightedPBCimage);
374 rvec_inc(center, weightedPBCimage);
378 svmul(1.0 / wsum, center, center);
382 /*! \brief Return TRUE if position x of ion (or water) is found in the compartment,
383 * i.e. between w1 and w2.
385 * One can define and additional offset "b" if one wants to exchange ions/water
386 * to or from a plane not directly in the middle of w1 and w2. The offset can be
387 * in ]-1.0, ..., +1.0 [.
388 * A bulkOffset of 0.0 means 'no offset', so the swap-layer is directly in the
389 * middle between w1 and w2. Offsets -1.0 < b < 0.0 will yield swaps nearer to w1,
390 * whereas offsets 0.0 < 0 < +1.0 will yield swaps nearer to w2.
394 * ||--------------+-------------|-------------+------------------------||
395 * w1 ? ? ? ? ? ? ? ? ? ? ? w2
396 * ||--------------+-------------|----b--------+------------------------||
401 * \param[in] w1 Position of 'wall' atom 1.
402 * \param[in] w2 Position of 'wall' atom 2.
403 * \param[in] x Position of the ion or the water molecule under consideration.
404 * \param[in] l Length of the box, from || to || in the sketch.
405 * \param[in] bulkOffset Where is the bulk layer "b" to be found between w1 and w2?
406 * \param[out] distance_from_b Distance of x to the bulk layer "b".
408 * \returns TRUE if x is between w1 and w2.
410 * Also computes the distance of x to the compartment center (the layer that is
411 * normally situated in the middle of w1 and w2 that would be considered as having
412 * the bulk concentration of ions).
414 static gmx_bool compartment_contains_atom(real w1, real w2, real x, real l, real bulkOffset, real* distance_from_b)
420 /* First set the origin in the middle of w1 and w2 */
427 /* Now choose the PBC image of x that is closest to the origin: */
438 *distance_from_b = static_cast<real>(fabs(x - bulkOffset * 0.5 * width));
440 /* Return TRUE if we now are in area "????" */
441 return (x >= w1) && (x < w2);
445 /*! \brief Updates the time-averaged number of ions in a compartment. */
446 static void update_time_window(t_compartment* comp, int values, int replace)
452 /* Put in the new value */
455 comp->nMolPast[replace] = comp->nMol;
458 /* Compute the new time-average */
460 for (i = 0; i < values; i++)
462 average += comp->nMolPast[i];
465 comp->nMolAv = average;
469 /*! \brief Add the atom with collective index ci to the atom list in compartment 'comp'.
471 * \param[in] ci Index of this ion in the collective xc array.
472 * \param[inout] comp Compartment to add this atom to.
473 * \param[in] distance Shortest distance of this atom to the bulk layer,
474 * from which ion/water pairs are selected for swapping.
476 static void add_to_list(int ci, t_compartment* comp, real distance)
480 if (nr >= comp->nalloc)
482 comp->nalloc = over_alloc_dd(nr + 1);
483 srenew(comp->ind, comp->nalloc);
484 srenew(comp->dist, comp->nalloc);
487 comp->dist[nr] = distance;
492 /*! \brief Determine the compartment boundaries from the channel centers. */
493 static void get_compartment_boundaries(Compartment c, t_swap* s, const matrix box, real* left, real* right)
496 real leftpos, rightpos, leftpos_orig;
499 if (c >= Compartment::Count)
501 gmx_fatal(FARGS, "Compartment out of range");
504 pos0 = s->group[static_cast<int>(SwapGroupSplittingType::Split0)].center[s->swapdim];
505 pos1 = s->group[static_cast<int>(SwapGroupSplittingType::Split1)].center[s->swapdim];
518 /* This gets us the other compartment: */
519 if (c == Compartment::B)
521 leftpos_orig = leftpos;
523 rightpos = leftpos_orig + box[s->swapdim][s->swapdim];
531 /*! \brief Determine the per-channel ion flux.
533 * To determine the flux through the individual channels, we
534 * remember the compartment and channel history of each ion. An ion can be
535 * either in channel0 or channel1, or in the remaining volume of compartment
539 * +-----------------+
542 * ||||||||||0|||||||| bilayer with channel 0
547 * |||||1||||||||||||| bilayer with channel 1
550 * +-----------------+
554 static void detect_flux_per_channel(t_swapgrp* g,
560 ChannelHistory* channel_label,
561 const t_swapcoords* sc,
570 gmx_bool in_cyl0, in_cyl1;
576 /* Check whether ion is inside any of the channels */
577 in_cyl0 = is_in_channel(atomPosition,
578 s->group[static_cast<int>(SwapGroupSplittingType::Split0)].center,
584 in_cyl1 = is_in_channel(atomPosition,
585 s->group[static_cast<int>(SwapGroupSplittingType::Split1)].center,
592 if (in_cyl0 && in_cyl1)
594 /* Ion appears to be in both channels. Something is severely wrong! */
596 *comp_now = Domain::Notset;
597 *comp_from = Domain::Notset;
598 *channel_label = ChannelHistory::None;
602 /* Ion is in channel 0 now */
603 *channel_label = ChannelHistory::Ch0;
604 *comp_now = Domain::Notset;
605 g->nCyl[Channel::Zero]++;
609 /* Ion is in channel 1 now */
610 *channel_label = ChannelHistory::Ch1;
611 *comp_now = Domain::Notset;
612 g->nCyl[Channel::One]++;
616 /* Ion is not in any of the channels, so it must be in domain A or B */
617 if (Compartment::A == comp)
619 *comp_now = Domain::A;
623 *comp_now = Domain::B;
627 /* Only take action, if ion is now in domain A or B, and was before
628 * in the other domain!
630 if (Domain::Notset == *comp_from)
632 /* Maybe we can set the domain now */
633 *comp_from = *comp_now; /* Could still be Domain::Notset, though */
635 else if ((*comp_now != Domain::Notset) /* if in channel */
636 && (*comp_from != *comp_now))
638 /* Obviously the ion changed its domain.
639 * Count this for the channel through which it has passed. */
640 switch (*channel_label)
642 case ChannelHistory::None:
646 " %s Warning! Step %s, ion %d moved from %s to %s\n",
648 gmx_step_str(step, buf),
650 DomainString[*comp_from],
651 DomainString[*comp_now]);
654 fprintf(stderr, ", possibly due to a swap in the original simulation.\n");
659 "but did not pass cyl0 or cyl1 as defined in the .mdp file.\n"
660 "Do you have an ion somewhere within the membrane?\n");
661 /* Write this info to the CompEL output file: */
663 " # Warning: step %s, ion %d moved from %s to %s (probably through the "
665 gmx_step_str(step, buf),
667 DomainString[*comp_from],
668 DomainString[*comp_now]);
671 case ChannelHistory::Ch0:
672 case ChannelHistory::Ch1:
673 if (*channel_label == ChannelHistory::Ch0)
682 if (Domain::A == *comp_from)
684 g->fluxfromAtoB[chan_nr]++;
688 g->fluxfromAtoB[chan_nr]--;
690 fprintf(fpout, "# Atom nr. %d finished passing %s.\n", iAtom, ChannelString[*channel_label]);
693 gmx_fatal(FARGS, "%s Unknown channel history entry for ion type '%s'\n", SwS.c_str(), g->molname);
696 /* This ion has moved to the _other_ compartment ... */
697 *comp_from = *comp_now;
698 /* ... and it did not pass any channel yet */
699 *channel_label = ChannelHistory::None;
704 /*! \brief Determines which ions or solvent molecules are in compartment A and B */
705 static void sortMoleculesIntoCompartments(t_swapgrp* g,
707 const t_swapcoords* sc,
715 gmx::EnumerationArray<Compartment, int> nMolNotInComp; /* consistency check */
716 real cyl0_r2 = sc->cyl0r * sc->cyl0r;
717 real cyl1_r2 = sc->cyl1r * sc->cyl1r;
719 /* Get us a counter that cycles in the range of [0 ... sc->nAverage[ */
720 int replace = (step / sc->nstswap) % sc->nAverage;
722 for (auto comp : gmx::EnumerationWrapper<Compartment>{})
726 /* Get lists of atoms that match criteria for this compartment */
727 get_compartment_boundaries(comp, s, box, &left, &right);
729 /* First clear the ion molecule lists */
730 g->comp[comp].nMol = 0;
731 nMolNotInComp[comp] = 0; /* consistency check */
733 /* Loop over the molecules and atoms of this group */
734 for (int iMol = 0, iAtom = 0; iAtom < static_cast<int>(g->atomset.numAtomsGlobal());
735 iAtom += g->apm, iMol++)
740 /* Is this first atom of the molecule in the compartment that we look at? */
741 if (compartment_contains_atom(
742 left, right, g->xc[iAtom][sd], box[sd][sd], sc->bulkOffset[comp], &dist))
744 /* Add the first atom of this molecule to the list of molecules in this compartment */
745 add_to_list(iAtom, &g->comp[comp], dist);
747 /* Master also checks for ion groups through which channel each ion has passed */
748 if (MASTER(cr) && (g->comp_now != nullptr) && !bIsSolvent)
750 int globalAtomNr = g->atomset.globalIndex()[iAtom] + 1; /* PDB index starts at 1 ... */
751 detect_flux_per_channel(g,
757 &g->channel_label[iMol],
769 nMolNotInComp[comp]++;
772 /* Correct the time-averaged number of ions in the compartment */
775 update_time_window(&g->comp[comp], sc->nAverage, replace);
779 /* Flux detection warnings */
780 if (MASTER(cr) && !bIsSolvent)
786 "%s Warning: %d atoms were detected as being in both channels! Probably your "
788 "%s cylinder is way too large, or one compartment has collapsed (step "
795 fprintf(s->fpout, "Warning: %d atoms were assigned to both channels!\n", g->nCylBoth);
801 if (bIsSolvent && nullptr != fpout)
804 "# Solv. molecules in comp.%s: %d comp.%s: %d\n",
805 CompStr[Compartment::A],
806 g->comp[Compartment::A].nMol,
807 CompStr[Compartment::B],
808 g->comp[Compartment::B].nMol);
811 /* Consistency checks */
812 const auto numMolecules = static_cast<int>(g->atomset.numAtomsGlobal() / g->apm);
813 if (nMolNotInComp[Compartment::A] + nMolNotInComp[Compartment::B] != numMolecules)
816 "%s Warning: Inconsistency while assigning '%s' molecules to compartments. !inA: "
817 "%d, !inB: %d, total molecules %d\n",
820 nMolNotInComp[Compartment::A],
821 nMolNotInComp[Compartment::B],
825 int sum = g->comp[Compartment::A].nMol + g->comp[Compartment::B].nMol;
826 if (sum != numMolecules)
829 "%s Warning: %d molecules are in group '%s', but altogether %d have been assigned "
830 "to the compartments.\n",
839 /*! \brief Find out how many group atoms are in the compartments initially */
840 static void get_initial_ioncounts(const t_inputrec* ir,
842 const rvec x[], /* the initial positions */
852 /* Loop over the user-defined (ion) groups */
853 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
857 /* Copy the initial positions of the atoms in the group
858 * to the collective array so that we can compartmentalize */
859 for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
861 int ind = g->atomset.globalIndex()[i];
862 copy_rvec(x[ind], g->xc[i]);
865 /* Set up the compartments and get lists of atoms in each compartment */
866 sortMoleculesIntoCompartments(g, cr, sc, s, box, 0, s->fpout, bRerun, FALSE);
868 /* Set initial molecule counts if requested (as signaled by "-1" value) */
869 for (auto ic : gmx::EnumerationWrapper<Compartment>{})
871 int requested = sc->grp[ig].nmolReq[ic];
874 g->comp[ic].nMolReq = g->comp[ic].nMol;
878 g->comp[ic].nMolReq = requested;
882 /* Check whether the number of requested molecules adds up to the total number */
883 int req = g->comp[Compartment::A].nMolReq + g->comp[Compartment::B].nMolReq;
884 int tot = g->comp[Compartment::A].nMol + g->comp[Compartment::B].nMol;
889 "Mismatch of the number of %s ions summed over both compartments.\n"
890 "You requested a total of %d ions (%d in A and %d in B),\n"
891 "but there are a total of %d ions of this type in the system.\n",
894 g->comp[Compartment::A].nMolReq,
895 g->comp[Compartment::B].nMolReq,
899 /* Initialize time-averaging:
900 * Write initial concentrations to all time bins to start with */
901 for (auto ic : gmx::EnumerationWrapper<Compartment>{})
903 g->comp[ic].nMolAv = g->comp[ic].nMol;
904 for (int i = 0; i < sc->nAverage; i++)
906 g->comp[ic].nMolPast[i] = g->comp[ic].nMol;
913 /*! \brief Copy history of ion counts from checkpoint file.
915 * When called, the checkpoint file has already been read in. Here we copy
916 * over the values from .cpt file to the swap data structure.
918 static void get_initial_ioncounts_from_cpt(const t_inputrec* ir,
920 swaphistory_t* swapstate,
932 /* Copy the past values from the checkpoint values that have been read in already */
935 fprintf(stderr, "%s Copying values from checkpoint\n", SwS.c_str());
938 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
941 gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
943 for (auto ic : gmx::EnumerationWrapper<Compartment>{})
945 g->comp[ic].nMolReq = gs->nMolReq[ic];
946 g->comp[ic].inflow_net = gs->inflow_net[ic];
951 "%s ... Influx netto: %d Requested: %d Past values: ",
953 g->comp[ic].inflow_net,
954 g->comp[ic].nMolReq);
957 for (int j = 0; j < sc->nAverage; j++)
959 g->comp[ic].nMolPast[j] = gs->nMolPast[ic][j];
962 fprintf(stderr, "%d ", g->comp[ic].nMolPast[j]);
967 fprintf(stderr, "\n");
975 /*! \brief The master lets all others know about the initial ion counts. */
976 static void bc_initial_concentrations(t_commrec* cr, t_swapcoords* swap, t_swap* s)
978 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
980 t_swapgrp* g = &s->group[ig];
982 for (auto ic : gmx::EnumerationWrapper<Compartment>{})
984 gmx_bcast(sizeof(g->comp[ic].nMolReq), &(g->comp[ic].nMolReq), cr->mpi_comm_mygroup);
985 gmx_bcast(sizeof(g->comp[ic].nMol), &(g->comp[ic].nMol), cr->mpi_comm_mygroup);
986 gmx_bcast(swap->nAverage * sizeof(g->comp[ic].nMolPast[0]), g->comp[ic].nMolPast, cr->mpi_comm_mygroup);
992 /*! \brief Ensure that each atom belongs to at most one of the swap groups. */
993 static void check_swap_groups(t_swap* s, int nat, gmx_bool bVerbose)
995 int* nGroup = nullptr; /* This array counts for each atom in the MD system to
996 how many swap groups it belongs (should be 0 or 1!) */
998 int nMultiple = 0; /* Number of atoms belonging to multiple groups */
1003 fprintf(stderr, "%s Making sure each atom belongs to at most one of the swap groups.\n", SwS.c_str());
1006 /* Add one to the group count of atoms belonging to a swap group: */
1008 for (int i = 0; i < s->ngrp; i++)
1010 t_swapgrp* g = &s->group[i];
1011 for (size_t j = 0; j < g->atomset.numAtomsGlobal(); j++)
1013 /* Get the global index of this atom of this group: */
1014 ind = g->atomset.globalIndex()[j];
1018 /* Make sure each atom belongs to at most one of the groups: */
1019 for (int i = 0; i < nat; i++)
1031 "%s Cannot perform swapping since %d atom%s allocated to more than one swap "
1033 "%s Each atom must be allocated to at most one of the split groups, the swap "
1034 "groups, or the solvent.\n"
1035 "%s Check the .mdp file settings regarding the swap index groups or the index "
1036 "groups themselves.\n",
1039 (1 == nMultiple) ? " is" : "s are",
1046 /*! \brief Get the number of atoms per molecule for this group.
1048 * Also ensure that all the molecules in this group have this number of atoms.
1050 static int get_group_apm_check(int igroup, t_swap* s, gmx_bool bVerbose, const gmx_mtop_t& mtop)
1052 t_swapgrp* g = &s->group[igroup];
1053 const int* ind = s->group[igroup].atomset.globalIndex().data();
1054 int nat = s->group[igroup].atomset.numAtomsGlobal();
1056 /* Determine the number of solvent atoms per solvent molecule from the
1057 * first solvent atom: */
1059 mtopGetMolblockIndex(mtop, ind[0], &molb, nullptr, nullptr);
1060 int apm = mtop.moleculeBlockIndices[molb].numAtomsPerMolecule;
1065 "%s Checking whether all %s molecules consist of %d atom%s\n",
1069 apm > 1 ? "s" : "");
1072 /* Check whether this is also true for all other solvent atoms */
1073 for (int i = 1; i < nat; i++)
1075 mtopGetMolblockIndex(mtop, ind[i], &molb, nullptr, nullptr);
1076 if (apm != mtop.moleculeBlockIndices[molb].numAtomsPerMolecule)
1078 gmx_fatal(FARGS, "Not all molecules of swap group %d consist of %d atoms.", igroup, apm);
1082 // TODO: check whether charges and masses of each molecule are identical!
1087 /*! \brief Print the legend to the swap output file.
1089 * Also print the initial values of ion counts and position of split groups.
1091 static void print_ionlist_legend(const t_inputrec* ir, t_swap* s, const gmx_output_env_t* oenv)
1093 const char** legend;
1097 int nIonTypes = ir->swap->ngrp - static_cast<int>(SwapGroupSplittingType::Count);
1099 static_cast<int>(Compartment::Count) * nIonTypes * 3 + 2
1100 + static_cast<int>(Channel::Count) * nIonTypes + 1);
1102 // Number of molecules and difference to reference counts for each
1103 // compartment and ion type
1104 for (auto ic : gmx::EnumerationWrapper<Compartment>{})
1106 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1108 t_swapGroup* g = &ir->swap->grp[ig];
1109 real q = s->group[ig].q;
1111 snprintf(buf, STRLEN, "%s %s ions (charge %s%g)", CompStr[ic], g->molname, q > 0 ? "+" : "", q);
1112 legend[count++] = gmx_strdup(buf);
1116 "%s av. mismatch to %d %s ions",
1118 s->group[ig].comp[ic].nMolReq,
1120 legend[count++] = gmx_strdup(buf);
1122 snprintf(buf, STRLEN, "%s net %s ion influx", CompStr[ic], g->molname);
1123 legend[count++] = gmx_strdup(buf);
1127 // Center of split groups
1130 "%scenter of %s of split group 0",
1131 SwapStr[ir->eSwapCoords],
1132 (nullptr != s->group[static_cast<int>(SwapGroupSplittingType::Split0)].m) ? "mass" : "geometry");
1133 legend[count++] = gmx_strdup(buf);
1136 "%scenter of %s of split group 1",
1137 SwapStr[ir->eSwapCoords],
1138 (nullptr != s->group[static_cast<int>(SwapGroupSplittingType::Split1)].m) ? "mass" : "geometry");
1139 legend[count++] = gmx_strdup(buf);
1141 // Ion flux for each channel and ion type
1142 for (auto ic : gmx::EnumerationWrapper<Channel>{})
1144 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1146 t_swapGroup* g = &ir->swap->grp[ig];
1147 snprintf(buf, STRLEN, "A->ch%d->B %s permeations", static_cast<int>(ic), g->molname);
1148 legend[count++] = gmx_strdup(buf);
1152 // Number of molecules that leaked from A to B
1153 snprintf(buf, STRLEN, "leakage");
1154 legend[count++] = gmx_strdup(buf);
1156 xvgr_legend(s->fpout, count, legend, oenv);
1159 "# Instantaneous ion counts and time-averaged differences to requested numbers\n");
1161 // We add a simple text legend helping to identify the columns with xvgr legend strings
1162 fprintf(s->fpout, "# time (ps)");
1163 for (int i = 0; i < count; i++)
1165 snprintf(buf, STRLEN, "s%d", i);
1166 fprintf(s->fpout, "%10s", buf);
1168 fprintf(s->fpout, "\n");
1173 /*! \brief Initialize channel ion flux detection routine.
1175 * Initialize arrays that keep track of where the ions come from and where
1178 static void detect_flux_per_channel_init(t_swap* s, swaphistory_t* swapstate, const bool isRestart)
1181 swapstateIons_t* gs;
1183 /* All these flux detection routines run on the master only */
1184 if (swapstate == nullptr)
1189 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1192 gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1194 /******************************************************/
1195 /* Channel and domain history for the individual ions */
1196 /******************************************************/
1197 if (isRestart) /* set the pointers right */
1199 g->comp_from = gs->comp_from;
1200 g->channel_label = gs->channel_label;
1202 else /* allocate memory for molecule counts */
1204 snew(g->comp_from, g->atomset.numAtomsGlobal() / g->apm);
1205 gs->comp_from = g->comp_from;
1206 snew(g->channel_label, g->atomset.numAtomsGlobal() / g->apm);
1207 gs->channel_label = g->channel_label;
1209 snew(g->comp_now, g->atomset.numAtomsGlobal() / g->apm);
1211 /* Initialize the channel and domain history counters */
1212 for (size_t i = 0; i < g->atomset.numAtomsGlobal() / g->apm; i++)
1214 g->comp_now[i] = Domain::Notset;
1217 g->comp_from[i] = Domain::Notset;
1218 g->channel_label[i] = ChannelHistory::None;
1222 /************************************/
1223 /* Channel fluxes for both channels */
1224 /************************************/
1225 g->nCyl[Channel::Zero] = 0;
1226 g->nCyl[Channel::One] = 0;
1232 fprintf(stderr, "%s Copying channel fluxes from checkpoint file data\n", SwS.c_str());
1236 // Loop over ion types (and both channels)
1237 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1240 gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1242 for (auto ic : gmx::EnumerationWrapper<Channel>{})
1245 "%s Channel %d flux history for ion type %s (charge %g): ",
1247 static_cast<int>(ic),
1252 g->fluxfromAtoB[ic] = gs->fluxfromAtoB[ic];
1256 g->fluxfromAtoB[ic] = 0;
1259 fprintf(stderr, "%d molecule%s", g->fluxfromAtoB[ic], g->fluxfromAtoB[ic] == 1 ? "" : "s");
1260 fprintf(stderr, "\n");
1264 /* Set pointers for checkpoint writing */
1265 swapstate->fluxleak_p = &s->fluxleak;
1266 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1269 gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1271 for (auto ic : gmx::EnumerationWrapper<Channel>{})
1273 gs->fluxfromAtoB_p[ic] = &g->fluxfromAtoB[ic];
1279 /*! \brief Outputs the initial structure to PDB file for debugging reasons.
1281 * Output the starting structure so that in case of multimeric channels
1282 * the user can check whether we have the correct PBC image for all atoms.
1283 * If this is not correct, the ion counts per channel will be very likely
1286 static void outputStartStructureIfWanted(const gmx_mtop_t& mtop, rvec* x, PbcType pbcType, const matrix box)
1288 char* env = getenv("GMX_COMPELDUMP");
1293 "\n%s Found env.var. GMX_COMPELDUMP, will output CompEL starting structure made "
1295 "%s In case of multimeric channels, please check whether they have the correct PBC "
1296 "representation.\n",
1300 write_sto_conf_mtop(
1301 "CompELAssumedWholeConfiguration.pdb", *mtop.name, mtop, x, nullptr, pbcType, box);
1306 /*! \brief Initialize the swapstate structure, used for checkpoint writing.
1308 * The swapstate struct stores the information we need to make the channels
1309 * whole again after restarts from a checkpoint file. Here we do the following:
1310 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
1311 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
1312 * c) in any case, for subsequent checkpoint writing, we set the pointers in
1313 * swapstate to the x_old arrays, which contain the correct PBC representation of
1314 * multimeric channels at the last time step.
1316 static void init_swapstate(swaphistory_t* swapstate,
1319 const gmx_mtop_t& mtop,
1320 const rvec* x, /* the initial positions */
1322 const t_inputrec* ir)
1324 rvec* x_pbc = nullptr; /* positions of the whole MD system with molecules made whole */
1328 /* We always need the last whole positions such that
1329 * in the next time step we can make the channels whole again in PBC */
1330 if (swapstate->bFromCpt)
1332 /* Copy the last whole positions of each channel from .cpt */
1333 g = &(s->group[static_cast<int>(SwapGroupSplittingType::Split0)]);
1334 for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1336 copy_rvec(swapstate->xc_old_whole[Channel::Zero][i], g->xc_old[i]);
1338 g = &(s->group[static_cast<int>(SwapGroupSplittingType::Split1)]);
1339 for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1341 copy_rvec(swapstate->xc_old_whole[Channel::One][i], g->xc_old[i]);
1346 swapstate->eSwapCoords = ir->eSwapCoords;
1348 /* Set the number of ion types and allocate memory for checkpointing */
1349 swapstate->nIonTypes = s->ngrp - static_cast<int>(SwapGroupSplittingType::Count);
1350 snew(swapstate->ionType, swapstate->nIonTypes);
1352 /* Store the total number of ions of each type in the swapstateIons
1353 * structure that is accessible during checkpoint writing */
1354 for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1356 swapstateIons_t* gs = &swapstate->ionType[ii];
1357 gs->nMol = sc->grp[ii + static_cast<int>(SwapGroupSplittingType::Count)].nat;
1360 /* Extract the initial split group positions. */
1362 /* Remove pbc, make molecule whole. */
1363 snew(x_pbc, mtop.natoms);
1364 copy_rvecn(x, x_pbc, 0, mtop.natoms);
1366 /* This can only make individual molecules whole, not multimers */
1367 do_pbc_mtop(ir->pbcType, box, &mtop, x_pbc);
1369 /* Output the starting structure? */
1370 outputStartStructureIfWanted(mtop, x_pbc, ir->pbcType, box);
1372 /* If this is the first run (i.e. no checkpoint present) we assume
1373 * that the starting positions give us the correct PBC representation */
1374 for (int ig = static_cast<int>(SwapGroupSplittingType::Split0);
1375 ig <= static_cast<int>(SwapGroupSplittingType::Split1);
1378 g = &(s->group[ig]);
1379 for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1381 copy_rvec(x_pbc[g->atomset.globalIndex()[i]], g->xc_old[i]);
1386 /* Prepare swapstate arrays for later checkpoint writing */
1387 swapstate->nat[Channel::Zero] =
1388 s->group[static_cast<int>(SwapGroupSplittingType::Split0)].atomset.numAtomsGlobal();
1389 swapstate->nat[Channel::One] =
1390 s->group[static_cast<int>(SwapGroupSplittingType::Split1)].atomset.numAtomsGlobal();
1393 /* For subsequent checkpoint writing, set the swapstate pointers to the xc_old
1394 * arrays that get updated at every swapping step */
1395 swapstate->xc_old_whole_p[Channel::Zero] =
1396 &s->group[static_cast<int>(SwapGroupSplittingType::Split0)].xc_old;
1397 swapstate->xc_old_whole_p[Channel::One] =
1398 &s->group[static_cast<int>(SwapGroupSplittingType::Split1)].xc_old;
1401 /*! \brief Determine the total charge imbalance resulting from the swap groups */
1402 static real getRequestedChargeImbalance(t_swap* s)
1407 real particle_charge;
1408 gmx::EnumerationArray<Compartment, real> particle_number;
1410 for (ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1414 particle_charge = g->q;
1415 particle_number[Compartment::A] = g->comp[Compartment::A].nMolReq;
1416 particle_number[Compartment::B] = g->comp[Compartment::B].nMolReq;
1418 DeltaQ += particle_charge * (particle_number[Compartment::A] - particle_number[Compartment::B]);
1425 /*! \brief Sorts anions and cations into two separate groups
1427 * This routine should be called for the 'anions' and 'cations' group,
1428 * of which the indices were lumped together in the older version of the code.
1430 static void copyIndicesToGroup(const int* indIons, int nIons, t_swapGroup* g, t_commrec* cr)
1434 /* If explicit ion counts were requested in the .mdp file
1435 * (by setting positive values for the number of ions),
1436 * we can make an additional consistency check here */
1437 if ((g->nmolReq[Compartment::A] < 0) && (g->nmolReq[Compartment::B] < 0))
1439 if (g->nat != (g->nmolReq[Compartment::A] + g->nmolReq[Compartment::B]))
1441 gmx_fatal_collective(FARGS,
1444 "%s Inconsistency while importing swap-related data from an old "
1445 "input file version.\n"
1446 "%s The requested ion counts in compartments A (%d) and B (%d)\n"
1447 "%s do not add up to the number of ions (%d) of this type for the "
1451 g->nmolReq[Compartment::A],
1452 g->nmolReq[Compartment::B],
1459 srenew(g->ind, g->nat);
1460 for (int i = 0; i < g->nat; i++)
1462 g->ind[i] = indIons[i];
1467 /*! \brief Converts old .tpr file CompEL contents to new data layout.
1469 * If we have read an old .tpr file (tpxv <= tpxv_CompElPolyatomicIonsAndMultipleIonTypes),
1470 * anions and cations are stored together in group #3. In the new
1471 * format we store each ion type in a separate group.
1472 * The 'classic' groups are:
1473 * #0 split group 0 - OK
1474 * #1 split group 1 - OK
1476 * #3 anions - contains also cations, needs to be converted
1477 * #4 cations - empty before conversion
1480 static void convertOldToNewGroupFormat(t_swapcoords* sc, const gmx_mtop_t& mtop, gmx_bool bVerbose, t_commrec* cr)
1482 t_swapGroup* g = &sc->grp[3];
1484 /* Loop through the atom indices of group #3 (anions) and put all indices
1485 * that belong to cations into the cation group.
1489 int* indAnions = nullptr;
1490 int* indCations = nullptr;
1491 snew(indAnions, g->nat);
1492 snew(indCations, g->nat);
1495 for (int i = 0; i < g->nat; i++)
1497 const t_atom& atom = mtopGetAtomParameters(mtop, g->ind[i], &molb);
1500 // This is an anion, add it to the list of anions
1501 indAnions[nAnions++] = g->ind[i];
1505 // This is a cation, add it to the list of cations
1506 indCations[nCations++] = g->ind[i];
1513 "%s Sorted %d ions into separate groups of %d anions and %d cations.\n",
1521 /* Now we have the correct lists of anions and cations.
1522 * Copy it to the right groups.
1524 copyIndicesToGroup(indAnions, nAnions, g, cr);
1526 copyIndicesToGroup(indCations, nCations, g, cr);
1532 /*! \brief Returns TRUE if we started from an old .tpr
1534 * Then we need to re-sort anions and cations into separate groups */
1535 static gmx_bool bConvertFromOldTpr(t_swapcoords* sc)
1537 // If the last group has no atoms it means we need to convert!
1538 return (sc->ngrp >= 5) && (0 == sc->grp[4].nat);
1542 t_swap* init_swapcoords(FILE* fplog,
1543 const t_inputrec* ir,
1545 const gmx_mtop_t& mtop,
1546 const t_state* globalState,
1547 ObservablesHistory* oh,
1549 gmx::LocalAtomSetManager* atomSets,
1550 const gmx_output_env_t* oenv,
1551 const gmx::MdrunOptions& mdrunOptions,
1552 const gmx::StartingBehavior startingBehavior)
1555 swapstateIons_t* gs;
1556 swaphistory_t* swapstate = nullptr;
1558 if ((PAR(cr)) && !DOMAINDECOMP(cr))
1560 gmx_fatal(FARGS, "Position swapping is only implemented for domain decomposition!");
1563 auto* sc = ir->swap;
1564 auto* s = new t_swap();
1566 if (mdrunOptions.rerun)
1571 "%s This module does not support reruns in parallel\nPlease request a serial "
1572 "run with -nt 1 / -np 1\n",
1576 fprintf(stderr, "%s Rerun - using every available frame\n", SwS.c_str());
1578 sc->nAverage = 1; /* averaging makes no sense for reruns */
1581 if (MASTER(cr) && startingBehavior == gmx::StartingBehavior::NewSimulation)
1583 fprintf(fplog, "\nInitializing ion/water position exchanges\n");
1584 please_cite(fplog, "Kutzner2011b");
1587 switch (ir->eSwapCoords)
1589 case SwapType::X: s->swapdim = XX; break;
1590 case SwapType::Y: s->swapdim = YY; break;
1591 case SwapType::Z: s->swapdim = ZZ; break;
1592 default: s->swapdim = -1; break;
1595 const gmx_bool bVerbose = mdrunOptions.verbose;
1597 // For compatibility with old .tpr files
1598 if (bConvertFromOldTpr(sc))
1600 convertOldToNewGroupFormat(sc, mtop, bVerbose && MASTER(cr), cr);
1603 /* Copy some data and pointers to the group structures for convenience */
1604 /* Number of atoms in the group */
1606 for (int i = 0; i < s->ngrp; i++)
1608 s->group.emplace_back(atomSets->add(
1609 gmx::ArrayRef<const int>(sc->grp[i].ind, sc->grp[i].ind + sc->grp[i].nat)));
1610 s->group[i].molname = sc->grp[i].molname;
1613 /* Check for overlapping atoms */
1614 check_swap_groups(s, mtop.natoms, bVerbose && MASTER(cr));
1616 /* Allocate space for the collective arrays for all groups */
1617 /* For the collective position array */
1618 for (int i = 0; i < s->ngrp; i++)
1621 snew(g->xc, g->atomset.numAtomsGlobal());
1623 /* For the split groups (the channels) we need some extra memory to
1624 * be able to make the molecules whole even if they span more than
1625 * half of the box size. */
1626 if ((i == static_cast<int>(SwapGroupSplittingType::Split0))
1627 || (i == static_cast<int>(SwapGroupSplittingType::Split1)))
1629 snew(g->xc_shifts, g->atomset.numAtomsGlobal());
1630 snew(g->xc_eshifts, g->atomset.numAtomsGlobal());
1631 snew(g->xc_old, g->atomset.numAtomsGlobal());
1637 if (oh->swapHistory == nullptr)
1639 oh->swapHistory = std::make_unique<swaphistory_t>(swaphistory_t{});
1641 swapstate = oh->swapHistory.get();
1643 init_swapstate(swapstate, sc, s, mtop, globalState->x.rvec_array(), globalState->box, ir);
1646 /* After init_swapstate we have a set of (old) whole positions for our
1647 * channels. Now transfer that to all nodes */
1650 for (int ig = static_cast<int>(SwapGroupSplittingType::Split0);
1651 ig <= static_cast<int>(SwapGroupSplittingType::Split1);
1654 g = &(s->group[ig]);
1655 gmx_bcast((g->atomset.numAtomsGlobal()) * sizeof((g->xc_old)[0]), g->xc_old, cr->mpi_comm_mygroup);
1659 /* Make sure that all molecules in the solvent and ion groups contain the
1660 * same number of atoms each */
1661 for (int ig = static_cast<int>(SwapGroupSplittingType::Solvent); ig < s->ngrp; ig++)
1665 g = &(s->group[ig]);
1666 g->apm = get_group_apm_check(ig, s, MASTER(cr) && bVerbose, mtop);
1668 /* Since all molecules of a group are equal, we only need enough space
1669 * to determine properties of a single molecule at at time */
1670 snew(g->m, g->apm); /* For the center of mass */
1671 charge = 0; /* To determine the charge imbalance */
1673 for (int j = 0; j < g->apm; j++)
1675 const t_atom& atom = mtopGetAtomParameters(mtop, g->atomset.globalIndex()[j], &molb);
1679 /* Total charge of one molecule of this group: */
1684 /* Need mass-weighted center of split group? */
1685 for (int j = static_cast<int>(SwapGroupSplittingType::Split0);
1686 j <= static_cast<int>(SwapGroupSplittingType::Split1);
1690 if (sc->massw_split[j])
1692 /* Save the split group masses if mass-weighting is requested */
1693 snew(g->m, g->atomset.numAtomsGlobal());
1695 for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1697 g->m[i] = mtopGetAtomMass(mtop, g->atomset.globalIndex()[i], &molb);
1702 /* Make a t_pbc struct on all nodes so that the molecules
1703 * chosen for an exchange can be made whole. */
1706 bool restartWithAppending = (startingBehavior == gmx::StartingBehavior::RestartWithAppending);
1712 "%s Opening output file %s%s\n",
1715 restartWithAppending ? " for appending" : "");
1718 s->fpout = gmx_fio_fopen(fn, restartWithAppending ? "a" : "w");
1720 if (!restartWithAppending)
1722 xvgr_header(s->fpout, "Molecule counts", "Time (ps)", "counts", exvggtXNY, oenv);
1724 for (int ig = 0; ig < s->ngrp; ig++)
1726 auto enumValue = static_cast<SwapGroupSplittingType>(ig);
1727 g = &(s->group[ig]);
1729 "# %s group '%s' contains %d atom%s",
1730 ig < static_cast<int>(SwapGroupSplittingType::Count) ? enumValueToString(enumValue)
1733 static_cast<int>(g->atomset.numAtomsGlobal()),
1734 (g->atomset.numAtomsGlobal() > 1) ? "s" : "");
1735 if (!(SwapGroupSplittingType::Split0 == enumValue
1736 || SwapGroupSplittingType::Split1 == enumValue))
1739 " with %d atom%s in each molecule of charge %g",
1741 (g->apm > 1) ? "s" : "",
1744 fprintf(s->fpout, ".\n");
1747 fprintf(s->fpout, "#\n# Initial positions of split groups:\n");
1750 for (int j = static_cast<int>(SwapGroupSplittingType::Split0);
1751 j <= static_cast<int>(SwapGroupSplittingType::Split1);
1754 auto enumValue = static_cast<SwapGroupSplittingType>(j);
1756 for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1758 copy_rvec(globalState->x[sc->grp[j].ind[i]], g->xc[i]);
1760 /* xc has the correct PBC representation for the two channels, so we do
1761 * not need to correct for that */
1762 get_center(g->xc, g->m, g->atomset.numAtomsGlobal(), g->center);
1763 if (!restartWithAppending)
1766 "# %s group %s-center %5f nm\n",
1767 enumValueToString(enumValue),
1769 g->center[s->swapdim]);
1773 if (!restartWithAppending)
1775 if ((0 != sc->bulkOffset[Compartment::A]) || (0 != sc->bulkOffset[Compartment::B]))
1777 fprintf(s->fpout, "#\n");
1779 "# You provided an offset for the position of the bulk layer(s).\n");
1781 "# That means the layers to/from which ions and water molecules are "
1784 "# are not midway (= at 0.0) between the compartment-defining layers (at "
1786 fprintf(s->fpout, "# bulk-offsetA = %g\n", sc->bulkOffset[Compartment::A]);
1787 fprintf(s->fpout, "# bulk-offsetB = %g\n", sc->bulkOffset[Compartment::B]);
1790 fprintf(s->fpout, "#\n");
1792 "# Split0 cylinder radius %f nm, up %f nm, down %f nm\n",
1797 "# Split1 cylinder radius %f nm, up %f nm, down %f nm\n",
1802 fprintf(s->fpout, "#\n");
1803 if (!mdrunOptions.rerun)
1806 "# Coupling constant (number of swap attempt steps to average over): %d "
1807 "(translates to %f ps).\n",
1809 sc->nAverage * sc->nstswap * ir->delta_t);
1810 fprintf(s->fpout, "# Threshold is %f\n", sc->threshold);
1811 fprintf(s->fpout, "#\n");
1813 "# Remarks about which atoms passed which channel use global atoms numbers "
1814 "starting at one.\n");
1823 /* Allocate memory to remember the past particle counts for time averaging */
1824 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1826 g = &(s->group[ig]);
1827 for (auto ic : keysOf(g->comp))
1829 snew(g->comp[ic].nMolPast, sc->nAverage);
1833 /* Get the initial particle concentrations and let the other nodes know */
1836 if (startingBehavior != gmx::StartingBehavior::NewSimulation)
1838 get_initial_ioncounts_from_cpt(ir, s, swapstate, cr, bVerbose);
1842 fprintf(stderr, "%s Determining initial numbers of ions per compartment.\n", SwS.c_str());
1843 get_initial_ioncounts(
1844 ir, s, globalState->x.rvec_array(), globalState->box, cr, mdrunOptions.rerun);
1847 /* Prepare (further) checkpoint writes ... */
1848 if (startingBehavior != gmx::StartingBehavior::NewSimulation)
1850 /* Consistency check */
1851 if (swapstate->nAverage != sc->nAverage)
1854 "%s Ion count averaging steps mismatch! checkpoint: %d, tpr: %d",
1856 swapstate->nAverage,
1862 swapstate->nAverage = sc->nAverage;
1864 fprintf(stderr, "%s Setting pointers for checkpoint writing\n", SwS.c_str());
1865 for (auto ic : gmx::EnumerationWrapper<Compartment>{})
1867 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1870 gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1872 gs->nMolReq_p[ic] = &(g->comp[ic].nMolReq);
1873 gs->nMolPast_p[ic] = &(g->comp[ic].nMolPast[0]);
1874 gs->inflow_net_p[ic] = &(g->comp[ic].inflow_net);
1878 /* Determine the total charge imbalance */
1879 s->deltaQ = getRequestedChargeImbalance(s);
1883 fprintf(stderr, "%s Requested charge imbalance is Q(A) - Q(B) = %g e.\n", SwS.c_str(), s->deltaQ);
1885 if (!restartWithAppending)
1887 fprintf(s->fpout, "# Requested charge imbalance is Q(A)-Q(B) = %g e.\n", s->deltaQ);
1893 bc_initial_concentrations(cr, ir->swap, s);
1896 /* Update the time-averaged number of molecules for all groups and compartments */
1897 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < sc->ngrp; ig++)
1900 for (auto ic : keysOf(g->comp))
1902 update_time_window(&g->comp[ic], sc->nAverage, -1);
1906 /* Initialize arrays that keep track of through which channel the ions go */
1907 detect_flux_per_channel_init(s, swapstate, startingBehavior != gmx::StartingBehavior::NewSimulation);
1909 /* We need to print the legend if we open this file for the first time. */
1910 if (MASTER(cr) && !restartWithAppending)
1912 print_ionlist_legend(ir, s, oenv);
1918 void finish_swapcoords(t_swap* s)
1926 // Close the swap output file
1927 gmx_fio_fclose(s->fpout);
1931 /*! \brief Do we need to swap a molecule in any of the ion groups with a water molecule at this step?
1933 * From the requested and average molecule counts we determine whether a swap is needed
1934 * at this time step.
1936 static gmx_bool need_swap(const t_swapcoords* sc, t_swap* s)
1938 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < sc->ngrp; ig++)
1940 t_swapgrp* g = &s->group[ig];
1942 for (auto ic : keysOf(g->comp))
1944 if (g->comp[ic].nMolReq - g->comp[ic].nMolAv >= sc->threshold)
1954 /*! \brief Return the index of an atom or molecule suitable for swapping.
1956 * Returns the index of an atom that is far off the compartment boundaries,
1957 * that is near to the bulk layer to/from which the swaps take place.
1958 * Other atoms of the molecule (if any) will directly follow the returned index.
1960 * \param[in] comp Structure containing compartment-specific data.
1961 * \param[in] molname Name of the molecule.
1963 * \returns Index of the first atom of the molecule chosen for a position exchange.
1965 static int get_index_of_distant_atom(t_compartment* comp, const char molname[])
1968 real d = GMX_REAL_MAX;
1971 /* comp->nat contains the original number of atoms in this compartment
1972 * prior to doing any swaps. Some of these atoms may already have been
1973 * swapped out, but then they are marked with a distance of GMX_REAL_MAX
1975 for (int iMol = 0; iMol < comp->nMolBefore; iMol++)
1977 if (comp->dist[iMol] < d)
1980 d = comp->dist[ibest];
1987 "Could not get index of %s atom. Compartment contains %d %s molecules before "
1994 /* Set the distance of this index to infinity such that it won't get selected again in
1997 comp->dist[ibest] = GMX_REAL_MAX;
1999 return comp->ind[ibest];
2003 /*! \brief Swaps centers of mass and makes molecule whole if broken */
2004 static void translate_positions(rvec* x, int apm, rvec old_com, rvec new_com, t_pbc* pbc)
2007 rvec reference, dx, correctPBCimage;
2010 /* Use the first atom as the reference for PBC */
2011 copy_rvec(x[0], reference);
2013 for (i = 0; i < apm; i++)
2015 /* PBC distance between position and reference */
2016 pbc_dx(pbc, x[i], reference, dx);
2018 /* Add PBC distance to reference */
2019 rvec_add(reference, dx, correctPBCimage);
2021 /* Subtract old_com from correct image and add new_com */
2022 rvec_dec(correctPBCimage, old_com);
2023 rvec_inc(correctPBCimage, new_com);
2025 copy_rvec(correctPBCimage, x[i]);
2030 /*! \brief Write back the modified local positions from the collective array to the official positions. */
2031 static void apply_modified_positions(swap_group* g, rvec x[])
2033 auto collectiveIndex = g->atomset.collectiveIndex().begin();
2034 for (const auto localIndex : g->atomset.localIndex())
2036 /* Copy the possibly modified position */
2037 copy_rvec(g->xc[*collectiveIndex], x[localIndex]);
2043 gmx_bool do_swapcoords(t_commrec* cr,
2046 const t_inputrec* ir,
2048 gmx_wallcycle* wcycle,
2054 const t_swapcoords* sc = ir->swap;
2055 gmx_bool bSwap = FALSE;
2058 rvec com_solvent, com_particle; /* solvent and swap molecule's center of mass */
2061 wallcycle_start(wcycle, WallCycleCounter::Swap);
2063 set_pbc(s->pbc, ir->pbcType, box);
2065 /* Assemble the positions of the split groups, i.e. the channels.
2066 * Here we also pass a shifts array to communicate_group_positions(), so that it can make
2067 * the molecules whole even in cases where they span more than half of the box in
2069 for (int ig = static_cast<int>(SwapGroupSplittingType::Split0);
2070 ig <= static_cast<int>(SwapGroupSplittingType::Split1);
2073 t_swapgrp* g = &(s->group[ig]);
2074 communicate_group_positions(cr,
2080 g->atomset.numAtomsGlobal(),
2081 g->atomset.numAtomsLocal(),
2082 g->atomset.localIndex().data(),
2083 g->atomset.collectiveIndex().data(),
2087 get_center(g->xc, g->m, g->atomset.numAtomsGlobal(), g->center); /* center of split groups == channels */
2090 /* Assemble the positions of the ions (ig = 3, 4, ...). These molecules should
2091 * be small and we can always make them whole with a simple distance check.
2092 * Therefore we pass NULL as third argument. */
2093 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
2095 t_swapgrp* g = &(s->group[ig]);
2096 communicate_group_positions(cr,
2102 g->atomset.numAtomsGlobal(),
2103 g->atomset.numAtomsLocal(),
2104 g->atomset.localIndex().data(),
2105 g->atomset.collectiveIndex().data(),
2109 /* Determine how many ions of this type each compartment contains */
2110 sortMoleculesIntoCompartments(g, cr, sc, s, box, step, s->fpout, bRerun, FALSE);
2113 /* Output how many ions are in the compartments */
2116 print_ionlist(s, t, "");
2119 /* If we are doing a rerun, we are finished here, since we cannot perform
2126 /* Do we have to perform a swap? */
2127 bSwap = need_swap(sc, s);
2130 /* Since we here know that we have to perform ion/water position exchanges,
2131 * we now assemble the solvent positions */
2132 t_swapgrp* g = &(s->group[static_cast<int>(SwapGroupSplittingType::Solvent)]);
2133 communicate_group_positions(cr,
2139 g->atomset.numAtomsGlobal(),
2140 g->atomset.numAtomsLocal(),
2141 g->atomset.localIndex().data(),
2142 g->atomset.collectiveIndex().data(),
2146 /* Determine how many molecules of solvent each compartment contains */
2147 sortMoleculesIntoCompartments(g, cr, sc, s, box, step, s->fpout, bRerun, TRUE);
2149 /* Save number of solvent molecules per compartment prior to any swaps */
2150 g->comp[Compartment::A].nMolBefore = g->comp[Compartment::A].nMol;
2151 g->comp[Compartment::B].nMolBefore = g->comp[Compartment::B].nMol;
2153 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
2155 t_swapgrp* g = &(s->group[ig]);
2157 for (auto ic : keysOf(g->comp))
2159 /* Determine in which compartment ions are missing and where they are too many */
2160 g->vacancy[ic] = g->comp[ic].nMolReq - g->comp[ic].nMolAv;
2162 /* Save number of ions per compartment prior to swaps */
2163 g->comp[ic].nMolBefore = g->comp[ic].nMol;
2167 /* Now actually perform the particle exchanges, one swap group after another */
2168 gsol = &s->group[static_cast<int>(SwapGroupSplittingType::Solvent)];
2169 for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
2172 t_swapgrp* g = &s->group[ig];
2173 for (auto thisC : gmx::EnumerationWrapper<Compartment>{})
2175 /* Index to the other compartment */
2176 auto otherC = thisC == Compartment::A ? Compartment::B : Compartment::A;
2178 while (g->vacancy[thisC] >= sc->threshold)
2180 /* Swap in an ion */
2182 /* Get the xc-index of the first atom of a solvent molecule of this compartment */
2183 isol = get_index_of_distant_atom(&gsol->comp[thisC], gsol->molname);
2185 /* Get the xc-index of a particle from the other compartment */
2186 iion = get_index_of_distant_atom(&g->comp[otherC], g->molname);
2188 get_molecule_center(&gsol->xc[isol], gsol->apm, gsol->m, com_solvent, s->pbc);
2189 get_molecule_center(&g->xc[iion], g->apm, g->m, com_particle, s->pbc);
2191 /* Subtract solvent molecule's center of mass and add swap particle's center of mass */
2192 translate_positions(&gsol->xc[isol], gsol->apm, com_solvent, com_particle, s->pbc);
2193 /* Similarly for the swap particle, subtract com_particle and add com_solvent */
2194 translate_positions(&g->xc[iion], g->apm, com_particle, com_solvent, s->pbc);
2196 /* Keep track of the changes */
2197 g->vacancy[thisC]--;
2198 g->vacancy[otherC]++;
2199 g->comp[thisC].nMol++;
2200 g->comp[otherC].nMol--;
2201 g->comp[thisC].inflow_net++;
2202 g->comp[otherC].inflow_net--;
2203 /* Correct the past time window to still get the right averages from now on */
2204 g->comp[thisC].nMolAv++;
2205 g->comp[otherC].nMolAv--;
2206 for (int j = 0; j < sc->nAverage; j++)
2208 g->comp[thisC].nMolPast[j]++;
2209 g->comp[otherC].nMolPast[j]--;
2211 /* Clear ion history */
2214 int iMol = iion / g->apm;
2215 g->channel_label[iMol] = ChannelHistory::None;
2216 g->comp_from[iMol] = Domain::Notset;
2218 /* That was the swap */
2223 if (nswaps && bVerbose)
2226 "%s Performed %d swap%s in step %" PRId64 " for iontype %s.\n",
2229 nswaps > 1 ? "s" : "",
2235 if (s->fpout != nullptr)
2237 print_ionlist(s, t, " # after swap");
2240 /* For the solvent and user-defined swap groups, each rank writes back its
2241 * (possibly modified) local positions to the official position array. */
2242 for (int ig = static_cast<int>(SwapGroupSplittingType::Solvent); ig < s->ngrp; ig++)
2244 t_swapgrp* g = &s->group[ig];
2245 apply_modified_positions(g, x);
2248 } /* end of if(bSwap) */
2250 wallcycle_stop(wcycle, WallCycleCounter::Swap);