src/gromacs/swap/swapcoords.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   9  *
  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.
  14  *
  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.
  19  *
  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.
  24  *
  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.
  32  *
  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.
  35  */
  36 /*! \internal \file
  37  * \brief
  38  * Implements functions in swapcoords.h.
  39  *
  40  * \author Carsten Kutzner <ckutzne@gwdg.de>
  41  * \ingroup module_swap
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "gromacs/utility/enumerationhelpers.h"
  46 #include "swapcoords.h"
  47
  48 #include <cstdio>
  49 #include <cstdlib>
  50 #include <ctime>
  51
  52 #include <memory>
  53 #include <string>
  54 #include <vector>
  55
  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"
  83
  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. */
  89
  90 /** Keep track of through which channel the ions have passed */
  91 enum class ChannelHistory : int
  92 {
  93     None,
  94     Ch0,
  95     Ch1,
  96     Count
  97 };
  98 static constexpr gmx::EnumerationArray<ChannelHistory, const char*> ChannelString = { "none",
  99                                                                                       "channel0",
 100                                                                                       "channel1" }; /**< Name for the channels */
 101
 102 /*! \brief Domain identifier.
 103  *
 104  * Keeps track of from which compartment the ions came before passing the
 105  * channel.
 106  */
 107 enum class Domain : int
 108 {
 109     Notset,
 110     A,
 111     B,
 112     Count
 113 };
 114 static constexpr gmx::EnumerationArray<Domain, const char*> DomainString = { "not_assigned",
 115                                                                              "Domain_A",
 116                                                                              "Domain_B" }; /**< Name for the domains */
 117
 118 namespace gmx
 119 {
 120
 121 extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
 122
 123 /*! \internal
 124  * \brief Implement Computational Electrophysiology swapping.
 125  */
 126 class SwapCoordinates final : public IMDModule
 127 {
 128     // From 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 {}
 134 };
 135
 136 std::unique_ptr<IMDModule> createSwapCoordinatesModule()
 137 {
 138     return std::make_unique<SwapCoordinates>();
 139 }
 140
 141
 142 } // namespace gmx
 143
 144 /*! \internal \brief
 145  * Structure containing compartment-specific data.
 146  */
 147 typedef struct swap_compartment
 148 {
 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.     */
 161 } t_compartment;
 162
 163
 164 /*! \internal \brief
 165  * This structure contains data needed for the groups involved in swapping:
 166  * split group 0, split group 1, solvent group, ion groups.
 167  */
 168 typedef struct swap_group
 169 {
 170     /*!\brief Construct a swap group given the managed swap atoms.
 171      *
 172      * \param[in] atomset Managed indices of atoms that are part of the swap group.
 173      */
 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?
 190                                                (size nMol) */
 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.             */
 198 } t_swapgrp;
 199
 200 t_swapgrp::swap_group(const gmx::LocalAtomSet& atomset) : atomset{ atomset }
 201 {
 202     center[0] = 0;
 203     center[1] = 0;
 204     center[2] = 0;
 205     for (auto compartment : keysOf(comp))
 206     {
 207         comp[compartment]    = {};
 208         vacancy[compartment] = 0;
 209     }
 210     for (auto channel : keysOf(fluxfromAtoB))
 211     {
 212         fluxfromAtoB[channel] = 0;
 213         nCyl[channel]         = 0;
 214     }
 215 };
 216
 217 /*! \internal \brief
 218  * Main (private) data structure for the position swapping protocol.
 219  */
 220 struct t_swap
 221 {
 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.  */
 229 };
 230
 231
 232 /*! \brief Check whether point is in channel.
 233  *
 234  * A channel is a cylinder defined by a disc
 235  * with radius r around its center c. The thickness of the cylinder is
 236  * d_up - d_down.
 237  *
 238  * \code
 239  *               ^  d_up
 240  *               |
 241  *     r         |
 242  *     <---------c--------->
 243  *               |
 244  *               v  d_down
 245  *
 246  * \endcode
 247  *
 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.
 255  *
 256  * \returns   Whether the point is inside the defined cylindric channel.
 257  */
 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)
 259 {
 260     rvec dr;
 261     int  plane1, plane2; /* Directions tangential to membrane */
 262
 263
 264     plane1 = (normal + 1) % 3; /* typically 0, i.e. XX */
 265     plane2 = (normal + 2) % 3; /* typically 1, i.e. YY */
 266
 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 */
 269
 270     /* Check vertical direction */
 271     if ((dr[normal] > d_up) || (dr[normal] < -d_down))
 272     {
 273         return FALSE;
 274     }
 275
 276     /* Check radial direction */
 277     if ((dr[plane1] * dr[plane1] + dr[plane2] * dr[plane2]) > r_cyl2)
 278     {
 279         return FALSE;
 280     }
 281
 282     /* All check passed, this point is in the cylinder */
 283     return TRUE;
 284 }
 285
 286
 287 /*! \brief Prints output to CompEL output file.
 288  *
 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.
 292  */
 293 static void print_ionlist(t_swap* s, double time, const char comment[])
 294 {
 295     // Output time
 296     fprintf(s->fpout, "%12.5e", time);
 297
 298     // Output number of molecules and difference to reference counts for each
 299     // compartment and ion type
 300     for (auto iComp : gmx::EnumerationWrapper<Compartment>{})
 301     {
 302         for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
 303         {
 304             t_compartment* comp = &s->group[ig].comp[iComp];
 305
 306             fprintf(s->fpout, "%10d%10.1f%10d", comp->nMol, comp->nMolAv - comp->nMolReq, comp->inflow_net);
 307         }
 308     }
 309
 310     // Output center of split groups
 311     fprintf(s->fpout,
 312             "%10g%10g",
 313             s->group[static_cast<int>(SwapGroupSplittingType::Split0)].center[s->swapdim],
 314             s->group[static_cast<int>(SwapGroupSplittingType::Split1)].center[s->swapdim]);
 315
 316     // Output ion flux for each channel and ion type
 317     for (auto iChan : gmx::EnumerationWrapper<Channel>{})
 318     {
 319         for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
 320         {
 321             t_swapgrp* g = &s->group[ig];
 322             fprintf(s->fpout, "%10d", g->fluxfromAtoB[iChan]);
 323         }
 324     }
 325
 326     /* Output the number of molecules that leaked from A to B */
 327     fprintf(s->fpout, "%10d", s->fluxleak);
 328
 329     fprintf(s->fpout, "%s\n", comment);
 330 }
 331
 332
 333 /*! \brief Get the center of a group of nat atoms.
 334  *
 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.
 337  */
 338 static void get_molecule_center(rvec x[], int nat, const real* weights, rvec center, t_pbc* pbc)
 339 {
 340     int  i;
 341     rvec weightedPBCimage;
 342     real wi, wsum;
 343     rvec reference, correctPBCimage, dx;
 344
 345
 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);
 349
 350     /* Calculate either the weighted center or simply the center of geometry */
 351     wsum = 0.0;
 352     clear_rvec(center);
 353     for (i = 0; i < nat; i++)
 354     {
 355         /* PBC distance between position and reference */
 356         pbc_dx(pbc, x[i], reference, dx);
 357
 358         /* Add PBC distance to reference */
 359         rvec_add(reference, dx, correctPBCimage);
 360
 361         /* Take weight into account */
 362         if (nullptr == weights)
 363         {
 364             wi = 1.0;
 365         }
 366         else
 367         {
 368             wi = weights[i];
 369         }
 370         wsum += wi;
 371         svmul(wi, correctPBCimage, weightedPBCimage);
 372
 373         /* Add to center */
 374         rvec_inc(center, weightedPBCimage);
 375     }
 376
 377     /* Normalize */
 378     svmul(1.0 / wsum, center, center);
 379 }
 380
 381
 382 /*! \brief Return TRUE if position x of ion (or water) is found in the compartment,
 383  * i.e. between w1 and w2.
 384  *
 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.
 391  *
 392  * \code
 393  *
 394  * ||--------------+-------------|-------------+------------------------||
 395  *                w1  ? ? ? ? ? ? ? ? ? ? ?   w2
 396  * ||--------------+-------------|----b--------+------------------------||
 397  *                -1            0.0           +1
 398  *
 399  * \endcode
 400  *
 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".
 407  *
 408  * \returns TRUE if x is between w1 and w2.
 409  *
 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).
 413  */
 414 static gmx_bool compartment_contains_atom(real w1, real w2, real x, real l, real bulkOffset, real* distance_from_b)
 415 {
 416     real m, l_2;
 417     real width;
 418
 419
 420     /* First set the origin in the middle of w1 and w2 */
 421     m = 0.5 * (w1 + w2);
 422     w1 -= m;
 423     w2 -= m;
 424     x -= m;
 425     width = w2 - w1;
 426
 427     /* Now choose the PBC image of x that is closest to the origin: */
 428     l_2 = 0.5 * l;
 429     while (x > l_2)
 430     {
 431         x -= l;
 432     }
 433     while (x <= -l_2)
 434     {
 435         x += l;
 436     }
 437
 438     *distance_from_b = static_cast<real>(fabs(x - bulkOffset * 0.5 * width));
 439
 440     /* Return TRUE if we now are in area "????" */
 441     return (x >= w1) && (x < w2);
 442 }
 443
 444
 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)
 447 {
 448     real average;
 449     int  i;
 450
 451
 452     /* Put in the new value */
 453     if (replace >= 0)
 454     {
 455         comp->nMolPast[replace] = comp->nMol;
 456     }
 457
 458     /* Compute the new time-average */
 459     average = 0.0;
 460     for (i = 0; i < values; i++)
 461     {
 462         average += comp->nMolPast[i];
 463     }
 464     average /= values;
 465     comp->nMolAv = average;
 466 }
 467
 468
 469 /*! \brief Add the atom with collective index ci to the atom list in compartment 'comp'.
 470  *
 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.
 475  */
 476 static void add_to_list(int ci, t_compartment* comp, real distance)
 477 {
 478     int nr = comp->nMol;
 479
 480     if (nr >= comp->nalloc)
 481     {
 482         comp->nalloc = over_alloc_dd(nr + 1);
 483         srenew(comp->ind, comp->nalloc);
 484         srenew(comp->dist, comp->nalloc);
 485     }
 486     comp->ind[nr]  = ci;
 487     comp->dist[nr] = distance;
 488     comp->nMol++;
 489 }
 490
 491
 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)
 494 {
 495     real pos0, pos1;
 496     real leftpos, rightpos, leftpos_orig;
 497
 498
 499     if (c >= Compartment::Count)
 500     {
 501         gmx_fatal(FARGS, "Compartment out of range");
 502     }
 503
 504     pos0 = s->group[static_cast<int>(SwapGroupSplittingType::Split0)].center[s->swapdim];
 505     pos1 = s->group[static_cast<int>(SwapGroupSplittingType::Split1)].center[s->swapdim];
 506
 507     if (pos0 < pos1)
 508     {
 509         leftpos  = pos0;
 510         rightpos = pos1;
 511     }
 512     else
 513     {
 514         leftpos  = pos1;
 515         rightpos = pos0;
 516     }
 517
 518     /* This gets us the other compartment: */
 519     if (c == Compartment::B)
 520     {
 521         leftpos_orig = leftpos;
 522         leftpos      = rightpos;
 523         rightpos     = leftpos_orig + box[s->swapdim][s->swapdim];
 524     }
 525
 526     *left  = leftpos;
 527     *right = rightpos;
 528 }
 529
 530
 531 /*! \brief Determine the per-channel ion flux.
 532  *
 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
 536  * A or B.
 537  *
 538  * \code
 539  *    +-----------------+
 540  *    |                 | B
 541  *    |                 | B compartment
 542  *    ||||||||||0|||||||| bilayer with channel 0
 543  *    |                 | A
 544  *    |                 | A
 545  *    |                 | A compartment
 546  *    |                 | A
 547  *    |||||1||||||||||||| bilayer with channel 1
 548  *    |                 | B
 549  *    |                 | B compartment
 550  *    +-----------------+
 551  *
 552  * \endcode
 553  */
 554 static void detect_flux_per_channel(t_swapgrp*          g,
 555                                     int                 iAtom,
 556                                     Compartment         comp,
 557                                     rvec                atomPosition,
 558                                     Domain*             comp_now,
 559                                     Domain*             comp_from,
 560                                     ChannelHistory*     channel_label,
 561                                     const t_swapcoords* sc,
 562                                     t_swap*             s,
 563                                     real                cyl0_r2,
 564                                     real                cyl1_r2,
 565                                     int64_t             step,
 566                                     gmx_bool            bRerun,
 567                                     FILE*               fpout)
 568 {
 569     int      sd, chan_nr;
 570     gmx_bool in_cyl0, in_cyl1;
 571     char     buf[STRLEN];
 572
 573
 574     sd = s->swapdim;
 575
 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,
 579                             sc->cyl0u,
 580                             sc->cyl0l,
 581                             cyl0_r2,
 582                             s->pbc,
 583                             sd);
 584     in_cyl1 = is_in_channel(atomPosition,
 585                             s->group[static_cast<int>(SwapGroupSplittingType::Split1)].center,
 586                             sc->cyl1u,
 587                             sc->cyl1l,
 588                             cyl1_r2,
 589                             s->pbc,
 590                             sd);
 591
 592     if (in_cyl0 && in_cyl1)
 593     {
 594         /* Ion appears to be in both channels. Something is severely wrong! */
 595         g->nCylBoth++;
 596         *comp_now      = Domain::Notset;
 597         *comp_from     = Domain::Notset;
 598         *channel_label = ChannelHistory::None;
 599     }
 600     else if (in_cyl0)
 601     {
 602         /* Ion is in channel 0 now */
 603         *channel_label = ChannelHistory::Ch0;
 604         *comp_now      = Domain::Notset;
 605         g->nCyl[Channel::Zero]++;
 606     }
 607     else if (in_cyl1)
 608     {
 609         /* Ion is in channel 1 now */
 610         *channel_label = ChannelHistory::Ch1;
 611         *comp_now      = Domain::Notset;
 612         g->nCyl[Channel::One]++;
 613     }
 614     else
 615     {
 616         /* Ion is not in any of the channels, so it must be in domain A or B */
 617         if (Compartment::A == comp)
 618         {
 619             *comp_now = Domain::A;
 620         }
 621         else
 622         {
 623             *comp_now = Domain::B;
 624         }
 625     }
 626
 627     /* Only take action, if ion is now in domain A or B, and was before
 628      * in the other domain!
 629      */
 630     if (Domain::Notset == *comp_from)
 631     {
 632         /* Maybe we can set the domain now */
 633         *comp_from = *comp_now; /* Could still be Domain::Notset, though */
 634     }
 635     else if ((*comp_now != Domain::Notset) /* if in channel */
 636              && (*comp_from != *comp_now))
 637     {
 638         /* Obviously the ion changed its domain.
 639          * Count this for the channel through which it has passed. */
 640         switch (*channel_label)
 641         {
 642             case ChannelHistory::None:
 643                 ++s->fluxleak;
 644
 645                 fprintf(stderr,
 646                         " %s Warning! Step %s, ion %d moved from %s to %s\n",
 647                         SwS.c_str(),
 648                         gmx_step_str(step, buf),
 649                         iAtom,
 650                         DomainString[*comp_from],
 651                         DomainString[*comp_now]);
 652                 if (bRerun)
 653                 {
 654                     fprintf(stderr, ", possibly due to a swap in the original simulation.\n");
 655                 }
 656                 else
 657                 {
 658                     fprintf(stderr,
 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: */
 662                     fprintf(s->fpout,
 663                             " # Warning: step %s, ion %d moved from %s to %s (probably through the "
 664                             "membrane)\n",
 665                             gmx_step_str(step, buf),
 666                             iAtom,
 667                             DomainString[*comp_from],
 668                             DomainString[*comp_now]);
 669                 }
 670                 break;
 671             case ChannelHistory::Ch0:
 672             case ChannelHistory::Ch1:
 673                 if (*channel_label == ChannelHistory::Ch0)
 674                 {
 675                     chan_nr = 0;
 676                 }
 677                 else
 678                 {
 679                     chan_nr = 1;
 680                 }
 681
 682                 if (Domain::A == *comp_from)
 683                 {
 684                     g->fluxfromAtoB[chan_nr]++;
 685                 }
 686                 else
 687                 {
 688                     g->fluxfromAtoB[chan_nr]--;
 689                 }
 690                 fprintf(fpout, "# Atom nr. %d finished passing %s.\n", iAtom, ChannelString[*channel_label]);
 691                 break;
 692             default:
 693                 gmx_fatal(FARGS, "%s Unknown channel history entry for ion type '%s'\n", SwS.c_str(), g->molname);
 694         }
 695
 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;
 700     }
 701 }
 702
 703
 704 /*! \brief Determines which ions or solvent molecules are in compartment A and B */
 705 static void sortMoleculesIntoCompartments(t_swapgrp*          g,
 706                                           t_commrec*          cr,
 707                                           const t_swapcoords* sc,
 708                                           t_swap*             s,
 709                                           const matrix        box,
 710                                           int64_t             step,
 711                                           FILE*               fpout,
 712                                           gmx_bool            bRerun,
 713                                           gmx_bool            bIsSolvent)
 714 {
 715     gmx::EnumerationArray<Compartment, int> nMolNotInComp; /* consistency check */
 716     real                                    cyl0_r2 = sc->cyl0r * sc->cyl0r;
 717     real                                    cyl1_r2 = sc->cyl1r * sc->cyl1r;
 718
 719     /* Get us a counter that cycles in the range of [0 ... sc->nAverage[ */
 720     int replace = (step / sc->nstswap) % sc->nAverage;
 721
 722     for (auto comp : gmx::EnumerationWrapper<Compartment>{})
 723     {
 724         real left, right;
 725
 726         /* Get lists of atoms that match criteria for this compartment */
 727         get_compartment_boundaries(comp, s, box, &left, &right);
 728
 729         /* First clear the ion molecule lists */
 730         g->comp[comp].nMol  = 0;
 731         nMolNotInComp[comp] = 0; /* consistency check */
 732
 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++)
 736         {
 737             real dist;
 738             int  sd = s->swapdim;
 739
 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))
 743             {
 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);
 746
 747                 /* Master also checks for ion groups through which channel each ion has passed */
 748                 if (MASTER(cr) && (g->comp_now != nullptr) && !bIsSolvent)
 749                 {
 750                     int globalAtomNr = g->atomset.globalIndex()[iAtom] + 1; /* PDB index starts at 1 ... */
 751                     detect_flux_per_channel(g,
 752                                             globalAtomNr,
 753                                             comp,
 754                                             g->xc[iAtom],
 755                                             &g->comp_now[iMol],
 756                                             &g->comp_from[iMol],
 757                                             &g->channel_label[iMol],
 758                                             sc,
 759                                             s,
 760                                             cyl0_r2,
 761                                             cyl1_r2,
 762                                             step,
 763                                             bRerun,
 764                                             fpout);
 765                 }
 766             }
 767             else
 768             {
 769                 nMolNotInComp[comp]++;
 770             }
 771         }
 772         /* Correct the time-averaged number of ions in the compartment */
 773         if (!bIsSolvent)
 774         {
 775             update_time_window(&g->comp[comp], sc->nAverage, replace);
 776         }
 777     }
 778
 779     /* Flux detection warnings */
 780     if (MASTER(cr) && !bIsSolvent)
 781     {
 782         if (g->nCylBoth > 0)
 783         {
 784             fprintf(stderr,
 785                     "\n"
 786                     "%s Warning: %d atoms were detected as being in both channels! Probably your "
 787                     "split\n"
 788                     "%s          cylinder is way too large, or one compartment has collapsed (step "
 789                     "%" PRId64 ")\n",
 790                     SwS.c_str(),
 791                     g->nCylBoth,
 792                     SwS.c_str(),
 793                     step);
 794
 795             fprintf(s->fpout, "Warning: %d atoms were assigned to both channels!\n", g->nCylBoth);
 796
 797             g->nCylBoth = 0;
 798         }
 799     }
 800
 801     if (bIsSolvent && nullptr != fpout)
 802     {
 803         fprintf(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);
 809     }
 810
 811     /* Consistency checks */
 812     const auto numMolecules = static_cast<int>(g->atomset.numAtomsGlobal() / g->apm);
 813     if (nMolNotInComp[Compartment::A] + nMolNotInComp[Compartment::B] != numMolecules)
 814     {
 815         fprintf(stderr,
 816                 "%s Warning: Inconsistency while assigning '%s' molecules to compartments. !inA: "
 817                 "%d, !inB: %d, total molecules %d\n",
 818                 SwS.c_str(),
 819                 g->molname,
 820                 nMolNotInComp[Compartment::A],
 821                 nMolNotInComp[Compartment::B],
 822                 numMolecules);
 823     }
 824
 825     int sum = g->comp[Compartment::A].nMol + g->comp[Compartment::B].nMol;
 826     if (sum != numMolecules)
 827     {
 828         fprintf(stderr,
 829                 "%s Warning: %d molecules are in group '%s', but altogether %d have been assigned "
 830                 "to the compartments.\n",
 831                 SwS.c_str(),
 832                 numMolecules,
 833                 g->molname,
 834                 sum);
 835     }
 836 }
 837
 838
 839 /*! \brief Find out how many group atoms are in the compartments initially */
 840 static void get_initial_ioncounts(const t_inputrec* ir,
 841                                   t_swap*           s,
 842                                   const rvec        x[], /* the initial positions */
 843                                   const matrix      box,
 844                                   t_commrec*        cr,
 845                                   gmx_bool          bRerun)
 846 {
 847     t_swapcoords* sc;
 848     t_swapgrp*    g;
 849
 850     sc = ir->swap;
 851
 852     /* Loop over the user-defined (ion) groups */
 853     for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
 854     {
 855         g = &s->group[ig];
 856
 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++)
 860         {
 861             int ind = g->atomset.globalIndex()[i];
 862             copy_rvec(x[ind], g->xc[i]);
 863         }
 864
 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);
 867
 868         /* Set initial molecule counts if requested (as signaled by "-1" value) */
 869         for (auto ic : gmx::EnumerationWrapper<Compartment>{})
 870         {
 871             int requested = sc->grp[ig].nmolReq[ic];
 872             if (requested < 0)
 873             {
 874                 g->comp[ic].nMolReq = g->comp[ic].nMol;
 875             }
 876             else
 877             {
 878                 g->comp[ic].nMolReq = requested;
 879             }
 880         }
 881
 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;
 885
 886         if ((req != tot))
 887         {
 888             gmx_fatal(FARGS,
 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",
 892                       g->molname,
 893                       req,
 894                       g->comp[Compartment::A].nMolReq,
 895                       g->comp[Compartment::B].nMolReq,
 896                       tot);
 897         }
 898
 899         /* Initialize time-averaging:
 900          * Write initial concentrations to all time bins to start with */
 901         for (auto ic : gmx::EnumerationWrapper<Compartment>{})
 902         {
 903             g->comp[ic].nMolAv = g->comp[ic].nMol;
 904             for (int i = 0; i < sc->nAverage; i++)
 905             {
 906                 g->comp[ic].nMolPast[i] = g->comp[ic].nMol;
 907             }
 908         }
 909     }
 910 }
 911
 912
 913 /*! \brief Copy history of ion counts from checkpoint file.
 914  *
 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.
 917  */
 918 static void get_initial_ioncounts_from_cpt(const t_inputrec* ir,
 919                                            t_swap*           s,
 920                                            swaphistory_t*    swapstate,
 921                                            t_commrec*        cr,
 922                                            gmx_bool          bVerbose)
 923 {
 924     t_swapcoords*    sc;
 925     t_swapgrp*       g;
 926     swapstateIons_t* gs;
 927
 928     sc = ir->swap;
 929
 930     if (MASTER(cr))
 931     {
 932         /* Copy the past values from the checkpoint values that have been read in already */
 933         if (bVerbose)
 934         {
 935             fprintf(stderr, "%s Copying values from checkpoint\n", SwS.c_str());
 936         }
 937
 938         for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
 939         {
 940             g  = &s->group[ig];
 941             gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
 942
 943             for (auto ic : gmx::EnumerationWrapper<Compartment>{})
 944             {
 945                 g->comp[ic].nMolReq    = gs->nMolReq[ic];
 946                 g->comp[ic].inflow_net = gs->inflow_net[ic];
 947
 948                 if (bVerbose)
 949                 {
 950                     fprintf(stderr,
 951                             "%s ... Influx netto: %d   Requested: %d   Past values: ",
 952                             SwS.c_str(),
 953                             g->comp[ic].inflow_net,
 954                             g->comp[ic].nMolReq);
 955                 }
 956
 957                 for (int j = 0; j < sc->nAverage; j++)
 958                 {
 959                     g->comp[ic].nMolPast[j] = gs->nMolPast[ic][j];
 960                     if (bVerbose)
 961                     {
 962                         fprintf(stderr, "%d ", g->comp[ic].nMolPast[j]);
 963                     }
 964                 }
 965                 if (bVerbose)
 966                 {
 967                     fprintf(stderr, "\n");
 968                 }
 969             }
 970         }
 971     }
 972 }
 973
 974
 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)
 977 {
 978     for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
 979     {
 980         t_swapgrp* g = &s->group[ig];
 981
 982         for (auto ic : gmx::EnumerationWrapper<Compartment>{})
 983         {
 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);
 987         }
 988     }
 989 }
 990
 991
 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)
 994 {
 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!) */
 997     int ind       = -1;
 998     int nMultiple = 0; /* Number of atoms belonging to multiple groups */
 999
1000
1001     if (bVerbose)
1002     {
1003         fprintf(stderr, "%s Making sure each atom belongs to at most one of the swap groups.\n", SwS.c_str());
1004     }
1005
1006     /* Add one to the group count of atoms belonging to a swap group: */
1007     snew(nGroup, nat);
1008     for (int i = 0; i < s->ngrp; i++)
1009     {
1010         t_swapgrp* g = &s->group[i];
1011         for (size_t j = 0; j < g->atomset.numAtomsGlobal(); j++)
1012         {
1013             /* Get the global index of this atom of this group: */
1014             ind = g->atomset.globalIndex()[j];
1015             nGroup[ind]++;
1016         }
1017     }
1018     /* Make sure each atom belongs to at most one of the groups: */
1019     for (int i = 0; i < nat; i++)
1020     {
1021         if (nGroup[i] > 1)
1022         {
1023             nMultiple++;
1024         }
1025     }
1026     sfree(nGroup);
1027
1028     if (nMultiple)
1029     {
1030         gmx_fatal(FARGS,
1031                   "%s Cannot perform swapping since %d atom%s allocated to more than one swap "
1032                   "index group.\n"
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",
1037                   SwS.c_str(),
1038                   nMultiple,
1039                   (1 == nMultiple) ? " is" : "s are",
1040                   SwSEmpty.c_str(),
1041                   SwSEmpty.c_str());
1042     }
1043 }
1044
1045
1046 /*! \brief Get the number of atoms per molecule for this group.
1047  *
1048  * Also ensure that all the molecules in this group have this number of atoms.
1049  */
1050 static int get_group_apm_check(int igroup, t_swap* s, gmx_bool bVerbose, const gmx_mtop_t& mtop)
1051 {
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();
1055
1056     /* Determine the number of solvent atoms per solvent molecule from the
1057      * first solvent atom: */
1058     int molb = 0;
1059     mtopGetMolblockIndex(mtop, ind[0], &molb, nullptr, nullptr);
1060     int apm = mtop.moleculeBlockIndices[molb].numAtomsPerMolecule;
1061
1062     if (bVerbose)
1063     {
1064         fprintf(stderr,
1065                 "%s Checking whether all %s molecules consist of %d atom%s\n",
1066                 SwS.c_str(),
1067                 g->molname,
1068                 apm,
1069                 apm > 1 ? "s" : "");
1070     }
1071
1072     /* Check whether this is also true for all other solvent atoms */
1073     for (int i = 1; i < nat; i++)
1074     {
1075         mtopGetMolblockIndex(mtop, ind[i], &molb, nullptr, nullptr);
1076         if (apm != mtop.moleculeBlockIndices[molb].numAtomsPerMolecule)
1077         {
1078             gmx_fatal(FARGS, "Not all molecules of swap group %d consist of %d atoms.", igroup, apm);
1079         }
1080     }
1081
1082     // TODO: check whether charges and masses of each molecule are identical!
1083     return apm;
1084 }
1085
1086
1087 /*! \brief Print the legend to the swap output file.
1088  *
1089  * Also print the initial values of ion counts and position of split groups.
1090  */
1091 static void print_ionlist_legend(const t_inputrec* ir, t_swap* s, const gmx_output_env_t* oenv)
1092 {
1093     const char** legend;
1094     int          count = 0;
1095     char         buf[STRLEN];
1096
1097     int nIonTypes = ir->swap->ngrp - static_cast<int>(SwapGroupSplittingType::Count);
1098     snew(legend,
1099          static_cast<int>(Compartment::Count) * nIonTypes * 3 + 2
1100                  + static_cast<int>(Channel::Count) * nIonTypes + 1);
1101
1102     // Number of molecules and difference to reference counts for each
1103     // compartment and ion type
1104     for (auto ic : gmx::EnumerationWrapper<Compartment>{})
1105     {
1106         for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1107         {
1108             t_swapGroup* g = &ir->swap->grp[ig];
1109             real         q = s->group[ig].q;
1110
1111             snprintf(buf, STRLEN, "%s %s ions (charge %s%g)", CompStr[ic], g->molname, q > 0 ? "+" : "", q);
1112             legend[count++] = gmx_strdup(buf);
1113
1114             snprintf(buf,
1115                      STRLEN,
1116                      "%s av. mismatch to %d %s ions",
1117                      CompStr[ic],
1118                      s->group[ig].comp[ic].nMolReq,
1119                      g->molname);
1120             legend[count++] = gmx_strdup(buf);
1121
1122             snprintf(buf, STRLEN, "%s net %s ion influx", CompStr[ic], g->molname);
1123             legend[count++] = gmx_strdup(buf);
1124         }
1125     }
1126
1127     // Center of split groups
1128     snprintf(buf,
1129              STRLEN,
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);
1134     snprintf(buf,
1135              STRLEN,
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);
1140
1141     // Ion flux for each channel and ion type
1142     for (auto ic : gmx::EnumerationWrapper<Channel>{})
1143     {
1144         for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1145         {
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);
1149         }
1150     }
1151
1152     // Number of molecules that leaked from A to B
1153     snprintf(buf, STRLEN, "leakage");
1154     legend[count++] = gmx_strdup(buf);
1155
1156     xvgr_legend(s->fpout, count, legend, oenv);
1157
1158     fprintf(s->fpout,
1159             "# Instantaneous ion counts and time-averaged differences to requested numbers\n");
1160
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++)
1164     {
1165         snprintf(buf, STRLEN, "s%d", i);
1166         fprintf(s->fpout, "%10s", buf);
1167     }
1168     fprintf(s->fpout, "\n");
1169     fflush(s->fpout);
1170 }
1171
1172
1173 /*! \brief Initialize channel ion flux detection routine.
1174  *
1175  * Initialize arrays that keep track of where the ions come from and where
1176  * they go.
1177  */
1178 static void detect_flux_per_channel_init(t_swap* s, swaphistory_t* swapstate, const bool isRestart)
1179 {
1180     t_swapgrp*       g;
1181     swapstateIons_t* gs;
1182
1183     /* All these flux detection routines run on the master only */
1184     if (swapstate == nullptr)
1185     {
1186         return;
1187     }
1188
1189     for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1190     {
1191         g  = &s->group[ig];
1192         gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1193
1194         /******************************************************/
1195         /* Channel and domain history for the individual ions */
1196         /******************************************************/
1197         if (isRestart) /* set the pointers right */
1198         {
1199             g->comp_from     = gs->comp_from;
1200             g->channel_label = gs->channel_label;
1201         }
1202         else /* allocate memory for molecule counts */
1203         {
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;
1208         }
1209         snew(g->comp_now, g->atomset.numAtomsGlobal() / g->apm);
1210
1211         /* Initialize the channel and domain history counters */
1212         for (size_t i = 0; i < g->atomset.numAtomsGlobal() / g->apm; i++)
1213         {
1214             g->comp_now[i] = Domain::Notset;
1215             if (!isRestart)
1216             {
1217                 g->comp_from[i]     = Domain::Notset;
1218                 g->channel_label[i] = ChannelHistory::None;
1219             }
1220         }
1221
1222         /************************************/
1223         /* Channel fluxes for both channels */
1224         /************************************/
1225         g->nCyl[Channel::Zero] = 0;
1226         g->nCyl[Channel::One]  = 0;
1227         g->nCylBoth            = 0;
1228     }
1229
1230     if (isRestart)
1231     {
1232         fprintf(stderr, "%s Copying channel fluxes from checkpoint file data\n", SwS.c_str());
1233     }
1234
1235
1236     // Loop over ion types (and both channels)
1237     for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1238     {
1239         g  = &s->group[ig];
1240         gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1241
1242         for (auto ic : gmx::EnumerationWrapper<Channel>{})
1243         {
1244             fprintf(stderr,
1245                     "%s Channel %d flux history for ion type %s (charge %g): ",
1246                     SwS.c_str(),
1247                     static_cast<int>(ic),
1248                     g->molname,
1249                     g->q);
1250             if (isRestart)
1251             {
1252                 g->fluxfromAtoB[ic] = gs->fluxfromAtoB[ic];
1253             }
1254             else
1255             {
1256                 g->fluxfromAtoB[ic] = 0;
1257             }
1258
1259             fprintf(stderr, "%d molecule%s", g->fluxfromAtoB[ic], g->fluxfromAtoB[ic] == 1 ? "" : "s");
1260             fprintf(stderr, "\n");
1261         }
1262     }
1263
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++)
1267     {
1268         g  = &s->group[ig];
1269         gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1270
1271         for (auto ic : gmx::EnumerationWrapper<Channel>{})
1272         {
1273             gs->fluxfromAtoB_p[ic] = &g->fluxfromAtoB[ic];
1274         }
1275     }
1276 }
1277
1278
1279 /*! \brief Outputs the initial structure to PDB file for debugging reasons.
1280  *
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
1284  * wrong.
1285  */
1286 static void outputStartStructureIfWanted(const gmx_mtop_t& mtop, rvec* x, PbcType pbcType, const matrix box)
1287 {
1288     char* env = getenv("GMX_COMPELDUMP");
1289
1290     if (env != nullptr)
1291     {
1292         fprintf(stderr,
1293                 "\n%s Found env.var. GMX_COMPELDUMP, will output CompEL starting structure made "
1294                 "whole.\n"
1295                 "%s In case of multimeric channels, please check whether they have the correct PBC "
1296                 "representation.\n",
1297                 SwS.c_str(),
1298                 SwSEmpty.c_str());
1299
1300         write_sto_conf_mtop(
1301                 "CompELAssumedWholeConfiguration.pdb", *mtop.name, mtop, x, nullptr, pbcType, box);
1302     }
1303 }
1304
1305
1306 /*! \brief Initialize the swapstate structure, used for checkpoint writing.
1307  *
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.
1315  */
1316 static void init_swapstate(swaphistory_t*    swapstate,
1317                            t_swapcoords*     sc,
1318                            t_swap*           s,
1319                            const gmx_mtop_t& mtop,
1320                            const rvec*       x, /* the initial positions */
1321                            const matrix      box,
1322                            const t_inputrec* ir)
1323 {
1324     rvec*      x_pbc = nullptr; /* positions of the whole MD system with molecules made whole */
1325     t_swapgrp* g;
1326
1327
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)
1331     {
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++)
1335         {
1336             copy_rvec(swapstate->xc_old_whole[Channel::Zero][i], g->xc_old[i]);
1337         }
1338         g = &(s->group[static_cast<int>(SwapGroupSplittingType::Split1)]);
1339         for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1340         {
1341             copy_rvec(swapstate->xc_old_whole[Channel::One][i], g->xc_old[i]);
1342         }
1343     }
1344     else
1345     {
1346         swapstate->eSwapCoords = ir->eSwapCoords;
1347
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);
1351
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++)
1355         {
1356             swapstateIons_t* gs = &swapstate->ionType[ii];
1357             gs->nMol            = sc->grp[ii + static_cast<int>(SwapGroupSplittingType::Count)].nat;
1358         }
1359
1360         /* Extract the initial split group positions. */
1361
1362         /* Remove pbc, make molecule whole. */
1363         snew(x_pbc, mtop.natoms);
1364         copy_rvecn(x, x_pbc, 0, mtop.natoms);
1365
1366         /* This can only make individual molecules whole, not multimers */
1367         do_pbc_mtop(ir->pbcType, box, &mtop, x_pbc);
1368
1369         /* Output the starting structure? */
1370         outputStartStructureIfWanted(mtop, x_pbc, ir->pbcType, box);
1371
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);
1376              ig++)
1377         {
1378             g = &(s->group[ig]);
1379             for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1380             {
1381                 copy_rvec(x_pbc[g->atomset.globalIndex()[i]], g->xc_old[i]);
1382             }
1383         }
1384         sfree(x_pbc);
1385
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();
1391     }
1392
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;
1399 }
1400
1401 /*! \brief Determine the total charge imbalance resulting from the swap groups */
1402 static real getRequestedChargeImbalance(t_swap* s)
1403 {
1404     int                                      ig;
1405     real                                     DeltaQ = 0.0;
1406     t_swapgrp*                               g;
1407     real                                     particle_charge;
1408     gmx::EnumerationArray<Compartment, real> particle_number;
1409
1410     for (ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1411     {
1412         g = &s->group[ig];
1413
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;
1417
1418         DeltaQ += particle_charge * (particle_number[Compartment::A] - particle_number[Compartment::B]);
1419     }
1420
1421     return DeltaQ;
1422 }
1423
1424
1425 /*! \brief Sorts anions and cations into two separate groups
1426  *
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.
1429  */
1430 static void copyIndicesToGroup(const int* indIons, int nIons, t_swapGroup* g, t_commrec* cr)
1431 {
1432     g->nat = nIons;
1433
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))
1438     {
1439         if (g->nat != (g->nmolReq[Compartment::A] + g->nmolReq[Compartment::B]))
1440         {
1441             gmx_fatal_collective(FARGS,
1442                                  cr->mpi_comm_mysim,
1443                                  MASTER(cr),
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 "
1448                                  "group '%s'.\n",
1449                                  SwS.c_str(),
1450                                  SwSEmpty.c_str(),
1451                                  g->nmolReq[Compartment::A],
1452                                  g->nmolReq[Compartment::B],
1453                                  SwSEmpty.c_str(),
1454                                  g->nat,
1455                                  g->molname);
1456         }
1457     }
1458
1459     srenew(g->ind, g->nat);
1460     for (int i = 0; i < g->nat; i++)
1461     {
1462         g->ind[i] = indIons[i];
1463     }
1464 }
1465
1466
1467 /*! \brief Converts old .tpr file CompEL contents to new data layout.
1468  *
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
1475  * #2 solvent        - OK
1476  * #3 anions         - contains also cations, needs to be converted
1477  * #4 cations        - empty before conversion
1478  *
1479  */
1480 static void convertOldToNewGroupFormat(t_swapcoords* sc, const gmx_mtop_t& mtop, gmx_bool bVerbose, t_commrec* cr)
1481 {
1482     t_swapGroup* g = &sc->grp[3];
1483
1484     /* Loop through the atom indices of group #3 (anions) and put all indices
1485      * that belong to cations into the cation group.
1486      */
1487     int  nAnions    = 0;
1488     int  nCations   = 0;
1489     int* indAnions  = nullptr;
1490     int* indCations = nullptr;
1491     snew(indAnions, g->nat);
1492     snew(indCations, g->nat);
1493
1494     int molb = 0;
1495     for (int i = 0; i < g->nat; i++)
1496     {
1497         const t_atom& atom = mtopGetAtomParameters(mtop, g->ind[i], &molb);
1498         if (atom.q < 0)
1499         {
1500             // This is an anion, add it to the list of anions
1501             indAnions[nAnions++] = g->ind[i];
1502         }
1503         else
1504         {
1505             // This is a cation, add it to the list of cations
1506             indCations[nCations++] = g->ind[i];
1507         }
1508     }
1509
1510     if (bVerbose)
1511     {
1512         fprintf(stdout,
1513                 "%s Sorted %d ions into separate groups of %d anions and %d cations.\n",
1514                 SwS.c_str(),
1515                 g->nat,
1516                 nAnions,
1517                 nCations);
1518     }
1519
1520
1521     /* Now we have the correct lists of anions and cations.
1522      * Copy it to the right groups.
1523      */
1524     copyIndicesToGroup(indAnions, nAnions, g, cr);
1525     g = &sc->grp[4];
1526     copyIndicesToGroup(indCations, nCations, g, cr);
1527     sfree(indAnions);
1528     sfree(indCations);
1529 }
1530
1531
1532 /*! \brief Returns TRUE if we started from an old .tpr
1533  *
1534  * Then we need to re-sort anions and cations into separate groups */
1535 static gmx_bool bConvertFromOldTpr(t_swapcoords* sc)
1536 {
1537     // If the last group has no atoms it means we need to convert!
1538     return (sc->ngrp >= 5) && (0 == sc->grp[4].nat);
1539 }
1540
1541
1542 t_swap* init_swapcoords(FILE*                       fplog,
1543                         const t_inputrec*           ir,
1544                         const char*                 fn,
1545                         const gmx_mtop_t&           mtop,
1546                         const t_state*              globalState,
1547                         ObservablesHistory*         oh,
1548                         t_commrec*                  cr,
1549                         gmx::LocalAtomSetManager*   atomSets,
1550                         const gmx_output_env_t*     oenv,
1551                         const gmx::MdrunOptions&    mdrunOptions,
1552                         const gmx::StartingBehavior startingBehavior)
1553 {
1554     t_swapgrp*       g;
1555     swapstateIons_t* gs;
1556     swaphistory_t*   swapstate = nullptr;
1557
1558     if ((PAR(cr)) && !DOMAINDECOMP(cr))
1559     {
1560         gmx_fatal(FARGS, "Position swapping is only implemented for domain decomposition!");
1561     }
1562
1563     auto* sc = ir->swap;
1564     auto* s  = new t_swap();
1565
1566     if (mdrunOptions.rerun)
1567     {
1568         if (PAR(cr))
1569         {
1570             gmx_fatal(FARGS,
1571                       "%s This module does not support reruns in parallel\nPlease request a serial "
1572                       "run with -nt 1 / -np 1\n",
1573                       SwS.c_str());
1574         }
1575
1576         fprintf(stderr, "%s Rerun - using every available frame\n", SwS.c_str());
1577         sc->nstswap  = 1;
1578         sc->nAverage = 1; /* averaging makes no sense for reruns */
1579     }
1580
1581     if (MASTER(cr) && startingBehavior == gmx::StartingBehavior::NewSimulation)
1582     {
1583         fprintf(fplog, "\nInitializing ion/water position exchanges\n");
1584         please_cite(fplog, "Kutzner2011b");
1585     }
1586
1587     switch (ir->eSwapCoords)
1588     {
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;
1593     }
1594
1595     const gmx_bool bVerbose = mdrunOptions.verbose;
1596
1597     // For compatibility with old .tpr files
1598     if (bConvertFromOldTpr(sc))
1599     {
1600         convertOldToNewGroupFormat(sc, mtop, bVerbose && MASTER(cr), cr);
1601     }
1602
1603     /* Copy some data and pointers to the group structures for convenience */
1604     /* Number of atoms in the group */
1605     s->ngrp = sc->ngrp;
1606     for (int i = 0; i < s->ngrp; i++)
1607     {
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;
1611     }
1612
1613     /* Check for overlapping atoms */
1614     check_swap_groups(s, mtop.natoms, bVerbose && MASTER(cr));
1615
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++)
1619     {
1620         g = &s->group[i];
1621         snew(g->xc, g->atomset.numAtomsGlobal());
1622
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)))
1628         {
1629             snew(g->xc_shifts, g->atomset.numAtomsGlobal());
1630             snew(g->xc_eshifts, g->atomset.numAtomsGlobal());
1631             snew(g->xc_old, g->atomset.numAtomsGlobal());
1632         }
1633     }
1634
1635     if (MASTER(cr))
1636     {
1637         if (oh->swapHistory == nullptr)
1638         {
1639             oh->swapHistory = std::make_unique<swaphistory_t>(swaphistory_t{});
1640         }
1641         swapstate = oh->swapHistory.get();
1642
1643         init_swapstate(swapstate, sc, s, mtop, globalState->x.rvec_array(), globalState->box, ir);
1644     }
1645
1646     /* After init_swapstate we have a set of (old) whole positions for our
1647      * channels. Now transfer that to all nodes */
1648     if (PAR(cr))
1649     {
1650         for (int ig = static_cast<int>(SwapGroupSplittingType::Split0);
1651              ig <= static_cast<int>(SwapGroupSplittingType::Split1);
1652              ig++)
1653         {
1654             g = &(s->group[ig]);
1655             gmx_bcast((g->atomset.numAtomsGlobal()) * sizeof((g->xc_old)[0]), g->xc_old, cr->mpi_comm_mygroup);
1656         }
1657     }
1658
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++)
1662     {
1663         real charge;
1664
1665         g      = &(s->group[ig]);
1666         g->apm = get_group_apm_check(ig, s, MASTER(cr) && bVerbose, mtop);
1667
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 */
1672         int molb = 0;
1673         for (int j = 0; j < g->apm; j++)
1674         {
1675             const t_atom& atom = mtopGetAtomParameters(mtop, g->atomset.globalIndex()[j], &molb);
1676             g->m[j]            = atom.m;
1677             charge += atom.q;
1678         }
1679         /* Total charge of one molecule of this group: */
1680         g->q = charge;
1681     }
1682
1683
1684     /* Need mass-weighted center of split group? */
1685     for (int j = static_cast<int>(SwapGroupSplittingType::Split0);
1686          j <= static_cast<int>(SwapGroupSplittingType::Split1);
1687          j++)
1688     {
1689         g = &(s->group[j]);
1690         if (sc->massw_split[j])
1691         {
1692             /* Save the split group masses if mass-weighting is requested */
1693             snew(g->m, g->atomset.numAtomsGlobal());
1694             int molb = 0;
1695             for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1696             {
1697                 g->m[i] = mtopGetAtomMass(mtop, g->atomset.globalIndex()[i], &molb);
1698             }
1699         }
1700     }
1701
1702     /* Make a t_pbc struct on all nodes so that the molecules
1703      * chosen for an exchange can be made whole. */
1704     snew(s->pbc, 1);
1705
1706     bool restartWithAppending = (startingBehavior == gmx::StartingBehavior::RestartWithAppending);
1707     if (MASTER(cr))
1708     {
1709         if (bVerbose)
1710         {
1711             fprintf(stderr,
1712                     "%s Opening output file %s%s\n",
1713                     SwS.c_str(),
1714                     fn,
1715                     restartWithAppending ? " for appending" : "");
1716         }
1717
1718         s->fpout = gmx_fio_fopen(fn, restartWithAppending ? "a" : "w");
1719
1720         if (!restartWithAppending)
1721         {
1722             xvgr_header(s->fpout, "Molecule counts", "Time (ps)", "counts", exvggtXNY, oenv);
1723
1724             for (int ig = 0; ig < s->ngrp; ig++)
1725             {
1726                 auto enumValue = static_cast<SwapGroupSplittingType>(ig);
1727                 g              = &(s->group[ig]);
1728                 fprintf(s->fpout,
1729                         "# %s group '%s' contains %d atom%s",
1730                         ig < static_cast<int>(SwapGroupSplittingType::Count) ? enumValueToString(enumValue)
1731                                                                              : "Ion",
1732                         g->molname,
1733                         static_cast<int>(g->atomset.numAtomsGlobal()),
1734                         (g->atomset.numAtomsGlobal() > 1) ? "s" : "");
1735                 if (!(SwapGroupSplittingType::Split0 == enumValue
1736                       || SwapGroupSplittingType::Split1 == enumValue))
1737                 {
1738                     fprintf(s->fpout,
1739                             " with %d atom%s in each molecule of charge %g",
1740                             g->apm,
1741                             (g->apm > 1) ? "s" : "",
1742                             g->q);
1743                 }
1744                 fprintf(s->fpout, ".\n");
1745             }
1746
1747             fprintf(s->fpout, "#\n# Initial positions of split groups:\n");
1748         }
1749
1750         for (int j = static_cast<int>(SwapGroupSplittingType::Split0);
1751              j <= static_cast<int>(SwapGroupSplittingType::Split1);
1752              j++)
1753         {
1754             auto enumValue = static_cast<SwapGroupSplittingType>(j);
1755             g              = &(s->group[j]);
1756             for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1757             {
1758                 copy_rvec(globalState->x[sc->grp[j].ind[i]], g->xc[i]);
1759             }
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)
1764             {
1765                 fprintf(s->fpout,
1766                         "# %s group %s-center %5f nm\n",
1767                         enumValueToString(enumValue),
1768                         DimStr[s->swapdim],
1769                         g->center[s->swapdim]);
1770             }
1771         }
1772
1773         if (!restartWithAppending)
1774         {
1775             if ((0 != sc->bulkOffset[Compartment::A]) || (0 != sc->bulkOffset[Compartment::B]))
1776             {
1777                 fprintf(s->fpout, "#\n");
1778                 fprintf(s->fpout,
1779                         "# You provided an offset for the position of the bulk layer(s).\n");
1780                 fprintf(s->fpout,
1781                         "# That means the layers to/from which ions and water molecules are "
1782                         "swapped\n");
1783                 fprintf(s->fpout,
1784                         "# are not midway (= at 0.0) between the compartment-defining layers (at "
1785                         "+/- 1.0).\n");
1786                 fprintf(s->fpout, "# bulk-offsetA = %g\n", sc->bulkOffset[Compartment::A]);
1787                 fprintf(s->fpout, "# bulk-offsetB = %g\n", sc->bulkOffset[Compartment::B]);
1788             }
1789
1790             fprintf(s->fpout, "#\n");
1791             fprintf(s->fpout,
1792                     "# Split0 cylinder radius %f nm, up %f nm, down %f nm\n",
1793                     sc->cyl0r,
1794                     sc->cyl0u,
1795                     sc->cyl0l);
1796             fprintf(s->fpout,
1797                     "# Split1 cylinder radius %f nm, up %f nm, down %f nm\n",
1798                     sc->cyl1r,
1799                     sc->cyl1u,
1800                     sc->cyl1l);
1801
1802             fprintf(s->fpout, "#\n");
1803             if (!mdrunOptions.rerun)
1804             {
1805                 fprintf(s->fpout,
1806                         "# Coupling constant (number of swap attempt steps to average over): %d  "
1807                         "(translates to %f ps).\n",
1808                         sc->nAverage,
1809                         sc->nAverage * sc->nstswap * ir->delta_t);
1810                 fprintf(s->fpout, "# Threshold is %f\n", sc->threshold);
1811                 fprintf(s->fpout, "#\n");
1812                 fprintf(s->fpout,
1813                         "# Remarks about which atoms passed which channel use global atoms numbers "
1814                         "starting at one.\n");
1815             }
1816         }
1817     }
1818     else
1819     {
1820         s->fpout = nullptr;
1821     }
1822
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++)
1825     {
1826         g = &(s->group[ig]);
1827         for (auto ic : keysOf(g->comp))
1828         {
1829             snew(g->comp[ic].nMolPast, sc->nAverage);
1830         }
1831     }
1832
1833     /* Get the initial particle concentrations and let the other nodes know */
1834     if (MASTER(cr))
1835     {
1836         if (startingBehavior != gmx::StartingBehavior::NewSimulation)
1837         {
1838             get_initial_ioncounts_from_cpt(ir, s, swapstate, cr, bVerbose);
1839         }
1840         else
1841         {
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);
1845         }
1846
1847         /* Prepare (further) checkpoint writes ... */
1848         if (startingBehavior != gmx::StartingBehavior::NewSimulation)
1849         {
1850             /* Consistency check */
1851             if (swapstate->nAverage != sc->nAverage)
1852             {
1853                 gmx_fatal(FARGS,
1854                           "%s Ion count averaging steps mismatch! checkpoint: %d, tpr: %d",
1855                           SwS.c_str(),
1856                           swapstate->nAverage,
1857                           sc->nAverage);
1858             }
1859         }
1860         else
1861         {
1862             swapstate->nAverage = sc->nAverage;
1863         }
1864         fprintf(stderr, "%s Setting pointers for checkpoint writing\n", SwS.c_str());
1865         for (auto ic : gmx::EnumerationWrapper<Compartment>{})
1866         {
1867             for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1868             {
1869                 g  = &s->group[ig];
1870                 gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1871
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);
1875             }
1876         }
1877
1878         /* Determine the total charge imbalance */
1879         s->deltaQ = getRequestedChargeImbalance(s);
1880
1881         if (bVerbose)
1882         {
1883             fprintf(stderr, "%s Requested charge imbalance is Q(A) - Q(B) = %g e.\n", SwS.c_str(), s->deltaQ);
1884         }
1885         if (!restartWithAppending)
1886         {
1887             fprintf(s->fpout, "# Requested charge imbalance is Q(A)-Q(B) = %g e.\n", s->deltaQ);
1888         }
1889     }
1890
1891     if (PAR(cr))
1892     {
1893         bc_initial_concentrations(cr, ir->swap, s);
1894     }
1895
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++)
1898     {
1899         g = &s->group[ig];
1900         for (auto ic : keysOf(g->comp))
1901         {
1902             update_time_window(&g->comp[ic], sc->nAverage, -1);
1903         }
1904     }
1905
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);
1908
1909     /* We need to print the legend if we open this file for the first time. */
1910     if (MASTER(cr) && !restartWithAppending)
1911     {
1912         print_ionlist_legend(ir, s, oenv);
1913     }
1914     return s;
1915 }
1916
1917
1918 void finish_swapcoords(t_swap* s)
1919 {
1920     if (s == nullptr)
1921     {
1922         return;
1923     }
1924     if (s->fpout)
1925     {
1926         // Close the swap output file
1927         gmx_fio_fclose(s->fpout);
1928     }
1929 }
1930
1931 /*! \brief Do we need to swap a molecule in any of the ion groups with a water molecule at this step?
1932  *
1933  * From the requested and average molecule counts we determine whether a swap is needed
1934  * at this time step.
1935  */
1936 static gmx_bool need_swap(const t_swapcoords* sc, t_swap* s)
1937 {
1938     for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < sc->ngrp; ig++)
1939     {
1940         t_swapgrp* g = &s->group[ig];
1941
1942         for (auto ic : keysOf(g->comp))
1943         {
1944             if (g->comp[ic].nMolReq - g->comp[ic].nMolAv >= sc->threshold)
1945             {
1946                 return TRUE;
1947             }
1948         }
1949     }
1950     return FALSE;
1951 }
1952
1953
1954 /*! \brief Return the index of an atom or molecule suitable for swapping.
1955  *
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.
1959  *
1960  * \param[in] comp    Structure containing compartment-specific data.
1961  * \param[in] molname Name of the molecule.
1962  *
1963  * \returns Index of the first atom of the molecule chosen for a position exchange.
1964  */
1965 static int get_index_of_distant_atom(t_compartment* comp, const char molname[])
1966 {
1967     int  ibest = -1;
1968     real d     = GMX_REAL_MAX;
1969
1970
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
1974      */
1975     for (int iMol = 0; iMol < comp->nMolBefore; iMol++)
1976     {
1977         if (comp->dist[iMol] < d)
1978         {
1979             ibest = iMol;
1980             d     = comp->dist[ibest];
1981         }
1982     }
1983
1984     if (ibest < 0)
1985     {
1986         gmx_fatal(FARGS,
1987                   "Could not get index of %s atom. Compartment contains %d %s molecules before "
1988                   "swaps.",
1989                   molname,
1990                   comp->nMolBefore,
1991                   molname);
1992     }
1993
1994     /* Set the distance of this index to infinity such that it won't get selected again in
1995      * this time step
1996      */
1997     comp->dist[ibest] = GMX_REAL_MAX;
1998
1999     return comp->ind[ibest];
2000 }
2001
2002
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)
2005 {
2006     int  i;
2007     rvec reference, dx, correctPBCimage;
2008
2009
2010     /* Use the first atom as the reference for PBC */
2011     copy_rvec(x[0], reference);
2012
2013     for (i = 0; i < apm; i++)
2014     {
2015         /* PBC distance between position and reference */
2016         pbc_dx(pbc, x[i], reference, dx);
2017
2018         /* Add PBC distance to reference */
2019         rvec_add(reference, dx, correctPBCimage);
2020
2021         /* Subtract old_com from correct image and add new_com */
2022         rvec_dec(correctPBCimage, old_com);
2023         rvec_inc(correctPBCimage, new_com);
2024
2025         copy_rvec(correctPBCimage, x[i]);
2026     }
2027 }
2028
2029
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[])
2032 {
2033     auto collectiveIndex = g->atomset.collectiveIndex().begin();
2034     for (const auto localIndex : g->atomset.localIndex())
2035     {
2036         /* Copy the possibly modified position */
2037         copy_rvec(g->xc[*collectiveIndex], x[localIndex]);
2038         ++collectiveIndex;
2039     }
2040 }
2041
2042
2043 gmx_bool do_swapcoords(t_commrec*        cr,
2044                        int64_t           step,
2045                        double            t,
2046                        const t_inputrec* ir,
2047                        t_swap*           s,
2048                        gmx_wallcycle*    wcycle,
2049                        rvec              x[],
2050                        matrix            box,
2051                        gmx_bool          bVerbose,
2052                        gmx_bool          bRerun)
2053 {
2054     const t_swapcoords* sc    = ir->swap;
2055     gmx_bool            bSwap = FALSE;
2056     t_swapgrp*          gsol;
2057     int                 isol, iion;
2058     rvec                com_solvent, com_particle; /* solvent and swap molecule's center of mass */
2059
2060
2061     wallcycle_start(wcycle, WallCycleCounter::Swap);
2062
2063     set_pbc(s->pbc, ir->pbcType, box);
2064
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
2068      * any dimension */
2069     for (int ig = static_cast<int>(SwapGroupSplittingType::Split0);
2070          ig <= static_cast<int>(SwapGroupSplittingType::Split1);
2071          ig++)
2072     {
2073         t_swapgrp* g = &(s->group[ig]);
2074         communicate_group_positions(cr,
2075                                     g->xc,
2076                                     g->xc_shifts,
2077                                     g->xc_eshifts,
2078                                     TRUE,
2079                                     x,
2080                                     g->atomset.numAtomsGlobal(),
2081                                     g->atomset.numAtomsLocal(),
2082                                     g->atomset.localIndex().data(),
2083                                     g->atomset.collectiveIndex().data(),
2084                                     g->xc_old,
2085                                     box);
2086
2087         get_center(g->xc, g->m, g->atomset.numAtomsGlobal(), g->center); /* center of split groups == channels */
2088     }
2089
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++)
2094     {
2095         t_swapgrp* g = &(s->group[ig]);
2096         communicate_group_positions(cr,
2097                                     g->xc,
2098                                     nullptr,
2099                                     nullptr,
2100                                     FALSE,
2101                                     x,
2102                                     g->atomset.numAtomsGlobal(),
2103                                     g->atomset.numAtomsLocal(),
2104                                     g->atomset.localIndex().data(),
2105                                     g->atomset.collectiveIndex().data(),
2106                                     nullptr,
2107                                     nullptr);
2108
2109         /* Determine how many ions of this type each compartment contains */
2110         sortMoleculesIntoCompartments(g, cr, sc, s, box, step, s->fpout, bRerun, FALSE);
2111     }
2112
2113     /* Output how many ions are in the compartments */
2114     if (MASTER(cr))
2115     {
2116         print_ionlist(s, t, "");
2117     }
2118
2119     /* If we are doing a rerun, we are finished here, since we cannot perform
2120      * swaps anyway */
2121     if (bRerun)
2122     {
2123         return FALSE;
2124     }
2125
2126     /* Do we have to perform a swap? */
2127     bSwap = need_swap(sc, s);
2128     if (bSwap)
2129     {
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,
2134                                     g->xc,
2135                                     nullptr,
2136                                     nullptr,
2137                                     FALSE,
2138                                     x,
2139                                     g->atomset.numAtomsGlobal(),
2140                                     g->atomset.numAtomsLocal(),
2141                                     g->atomset.localIndex().data(),
2142                                     g->atomset.collectiveIndex().data(),
2143                                     nullptr,
2144                                     nullptr);
2145
2146         /* Determine how many molecules of solvent each compartment contains */
2147         sortMoleculesIntoCompartments(g, cr, sc, s, box, step, s->fpout, bRerun, TRUE);
2148
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;
2152
2153         for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
2154         {
2155             t_swapgrp* g = &(s->group[ig]);
2156
2157             for (auto ic : keysOf(g->comp))
2158             {
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;
2161
2162                 /* Save number of ions per compartment prior to swaps */
2163                 g->comp[ic].nMolBefore = g->comp[ic].nMol;
2164             }
2165         }
2166
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++)
2170         {
2171             int        nswaps = 0;
2172             t_swapgrp* g      = &s->group[ig];
2173             for (auto thisC : gmx::EnumerationWrapper<Compartment>{})
2174             {
2175                 /* Index to the other compartment */
2176                 auto otherC = thisC == Compartment::A ? Compartment::B : Compartment::A;
2177
2178                 while (g->vacancy[thisC] >= sc->threshold)
2179                 {
2180                     /* Swap in an ion */
2181
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);
2184
2185                     /* Get the xc-index of a particle from the other compartment */
2186                     iion = get_index_of_distant_atom(&g->comp[otherC], g->molname);
2187
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);
2190
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);
2195
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++)
2207                     {
2208                         g->comp[thisC].nMolPast[j]++;
2209                         g->comp[otherC].nMolPast[j]--;
2210                     }
2211                     /* Clear ion history */
2212                     if (MASTER(cr))
2213                     {
2214                         int iMol               = iion / g->apm;
2215                         g->channel_label[iMol] = ChannelHistory::None;
2216                         g->comp_from[iMol]     = Domain::Notset;
2217                     }
2218                     /* That was the swap */
2219                     nswaps++;
2220                 }
2221             }
2222
2223             if (nswaps && bVerbose)
2224             {
2225                 fprintf(stderr,
2226                         "%s Performed %d swap%s in step %" PRId64 " for iontype %s.\n",
2227                         SwS.c_str(),
2228                         nswaps,
2229                         nswaps > 1 ? "s" : "",
2230                         step,
2231                         g->molname);
2232             }
2233         }
2234
2235         if (s->fpout != nullptr)
2236         {
2237             print_ionlist(s, t, "  # after swap");
2238         }
2239
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++)
2243         {
2244             t_swapgrp* g = &s->group[ig];
2245             apply_modified_positions(g, x);
2246         }
2247
2248     } /* end of if(bSwap) */
2249
2250     wallcycle_stop(wcycle, WallCycleCounter::Swap);
2251
2252     return bSwap;
2253 }