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 char* SwS      = { "SWAP:" }; /**< For output that comes from the swap module */
  85 static const char* SwSEmpty = { "     " }; /**< Placeholder for multi-line output */
  86 static const char* CompStr[eCompNR] = { "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* 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 eChannelHistory
  92 {
  93     eChHistPassedNone,
  94     eChHistPassedCh0,
  95     eChHistPassedCh1,
  96     eChHistNr
  97 };
  98 static const char* ChannelString[eChHistNr] = { "none", "channel0", "channel1" }; /**< Name for the channels */
  99
 100 /*! \brief Domain identifier.
 101  *
 102  * Keeps track of from which compartment the ions came before passing the
 103  * channel.
 104  */
 105 enum eDomain
 106 {
 107     eDomainNotset,
 108     eDomainA,
 109     eDomainB,
 110     eDomainNr
 111 };
 112 static const char* DomainString[eDomainNr] = { "not_assigned", "Domain_A", "Domain_B" }; /**< Name for the domains */
 113
 114 namespace gmx
 115 {
 116
 117 extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const int> globalAtomIndex);
 118
 119 /*! \internal
 120  * \brief Implement Computational Electrophysiology swapping.
 121  */
 122 class SwapCoordinates final : public IMDModule
 123 {
 124     // From IMDModule
 125     IMdpOptionProvider* mdpOptionProvider() override { return nullptr; }
 126     IMDOutputProvider*  outputProvider() override { return nullptr; }
 127     void                initForceProviders(ForceProviders* /* forceProviders */) override {}
 128     void subscribeToSimulationSetupNotifications(MdModulesNotifier* /* notifier */) override {}
 129     void subscribeToPreProcessingNotifications(MdModulesNotifier* /* notifier */) override {}
 130 };
 131
 132 std::unique_ptr<IMDModule> createSwapCoordinatesModule()
 133 {
 134     return std::make_unique<SwapCoordinates>();
 135 }
 136
 137
 138 } // namespace gmx
 139
 140 /*! \internal \brief
 141  * Structure containing compartment-specific data.
 142  */
 143 typedef struct swap_compartment
 144 {
 145     int nMol;        /**< Number of ion or water molecules detected
 146                           in this compartment.                          */
 147     int  nMolBefore; /**< Number of molecules before swapping.          */
 148     int  nMolReq;    /**< Requested number of molecules in compartment. */
 149     real nMolAv;     /**< Time-averaged number of molecules matching
 150                           the compartment conditions.                   */
 151     int*  nMolPast;  /**< Past molecule counts for time-averaging.      */
 152     int*  ind;       /**< Indices to collective array of atoms.         */
 153     real* dist;      /**< Distance of atom to bulk layer, which is
 154                           normally the center layer of the compartment  */
 155     int nalloc;      /**< Allocation size for ind array.                */
 156     int inflow_net;  /**< Net inflow of ions into this compartment.     */
 157 } t_compartment;
 158
 159
 160 /*! \internal \brief
 161  * This structure contains data needed for the groups involved in swapping:
 162  * split group 0, split group 1, solvent group, ion groups.
 163  */
 164 typedef struct swap_group
 165 {
 166     /*!\brief Construct a swap group given the managed swap atoms.
 167      *
 168      * \param[in] atomset Managed indices of atoms that are part of the swap group.
 169      */
 170     swap_group(const gmx::LocalAtomSet& atomset);
 171     char* molname = nullptr;        /**< Name of the group or ion type                         */
 172     int   apm     = 0;              /**< Number of atoms in each molecule                      */
 173     gmx::LocalAtomSet atomset;      /**< The atom indices in the swap group                    */
 174     rvec*             xc = nullptr; /**< Collective array of group atom positions (size nat)   */
 175     ivec* xc_shifts      = nullptr; /**< Current (collective) shifts (size nat)                */
 176     ivec* xc_eshifts     = nullptr; /**< Extra shifts since last DD step (size nat)            */
 177     rvec* xc_old         = nullptr; /**< Old (collective) positions (size nat)                 */
 178     real  q              = 0.;      /**< Total charge of one molecule of this group            */
 179     real* m              = nullptr; /**< Masses (can be omitted, size apm)                     */
 180     unsigned char* comp_from = nullptr; /**< (Collective) Stores from which compartment this
 181                                            molecule has come. This way we keep track of
 182                                            through which channel an ion permeates
 183                                            (size nMol = nat/apm)                                 */
 184     unsigned char* comp_now = nullptr; /**< In which compartment this ion is now (size nMol)      */
 185     unsigned char* channel_label = nullptr; /**< Which channel was passed at last by this ion?
 186                                                (size nMol) */
 187     rvec          center;        /**< Center of the group; COM if masses are used           */
 188     t_compartment comp[eCompNR]; /**< Distribution of particles of this group across
 189                                        the two compartments                                 */
 190     real vacancy[eCompNR];       /**< How many molecules need to be swapped in?             */
 191     int  fluxfromAtoB[eChanNR];  /**< Net flux of ions per channel                          */
 192     int  nCyl[eChanNR];          /**< Number of ions residing in a channel                  */
 193     int  nCylBoth = 0;           /**< Ions assigned to cyl0 and cyl1. Not good.             */
 194 } t_swapgrp;
 195
 196 t_swapgrp::swap_group(const gmx::LocalAtomSet& atomset) : atomset{ atomset }
 197 {
 198     center[0] = 0;
 199     center[1] = 0;
 200     center[2] = 0;
 201     for (int compartment = eCompA; compartment < eCompNR; ++compartment)
 202     {
 203         comp[compartment]    = {};
 204         vacancy[compartment] = 0;
 205     }
 206     for (int channel = eChan0; channel < eChanNR; ++channel)
 207     {
 208         fluxfromAtoB[channel] = 0;
 209         nCyl[channel]         = 0;
 210     }
 211 };
 212
 213 /*! \internal \brief
 214  * Main (private) data structure for the position swapping protocol.
 215  */
 216 struct t_swap
 217 {
 218     int                    swapdim;  /**< One of XX, YY, ZZ                               */
 219     t_pbc*                 pbc;      /**< Needed to make molecules whole.                 */
 220     FILE*                  fpout;    /**< Output file.                                    */
 221     int                    ngrp;     /**< Number of t_swapgrp groups                      */
 222     std::vector<t_swapgrp> group;    /**< Separate groups for channels, solvent, ions     */
 223     int                    fluxleak; /**< Flux not going through any of the channels.     */
 224     real                   deltaQ;   /**< The charge imbalance between the compartments.  */
 225 };
 226
 227
 228 /*! \brief Check whether point is in channel.
 229  *
 230  * A channel is a cylinder defined by a disc
 231  * with radius r around its center c. The thickness of the cylinder is
 232  * d_up - d_down.
 233  *
 234  * \code
 235  *               ^  d_up
 236  *               |
 237  *     r         |
 238  *     <---------c--------->
 239  *               |
 240  *               v  d_down
 241  *
 242  * \endcode
 243  *
 244  * \param[in] point    The position (xyz) under consideration.
 245  * \param[in] center   The center of the cylinder.
 246  * \param[in] d_up     The upper extension of the cylinder.
 247  * \param[in] d_down   The lower extension.
 248  * \param[in] r_cyl2   Cylinder radius squared.
 249  * \param[in] pbc      Structure with info about periodic boundary conditions.
 250  * \param[in] normal   The membrane normal direction is typically 3, i.e. z, but can be x or y also.
 251  *
 252  * \returns   Whether the point is inside the defined cylindric channel.
 253  */
 254 static gmx_bool is_in_channel(rvec point, rvec center, real d_up, real d_down, real r_cyl2, t_pbc* pbc, int normal)
 255 {
 256     rvec dr;
 257     int  plane1, plane2; /* Directions tangential to membrane */
 258
 259
 260     plane1 = (normal + 1) % 3; /* typically 0, i.e. XX */
 261     plane2 = (normal + 2) % 3; /* typically 1, i.e. YY */
 262
 263     /* Get the distance vector dr between the point and the center of the cylinder */
 264     pbc_dx(pbc, point, center, dr); /* This puts center in the origin */
 265
 266     /* Check vertical direction */
 267     if ((dr[normal] > d_up) || (dr[normal] < -d_down))
 268     {
 269         return FALSE;
 270     }
 271
 272     /* Check radial direction */
 273     if ((dr[plane1] * dr[plane1] + dr[plane2] * dr[plane2]) > r_cyl2)
 274     {
 275         return FALSE;
 276     }
 277
 278     /* All check passed, this point is in the cylinder */
 279     return TRUE;
 280 }
 281
 282
 283 /*! \brief Prints output to CompEL output file.
 284  *
 285  * Prints to swap output file how many ions are in each compartment,
 286  * where the centers of the split groups are, and how many ions of each type
 287  * passed the channels.
 288  */
 289 static void print_ionlist(t_swap* s, double time, const char comment[])
 290 {
 291     // Output time
 292     fprintf(s->fpout, "%12.5e", time);
 293
 294     // Output number of molecules and difference to reference counts for each
 295     // compartment and ion type
 296     for (int iComp = 0; iComp < eCompNR; iComp++)
 297     {
 298         for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
 299         {
 300             t_compartment* comp = &s->group[ig].comp[iComp];
 301
 302             fprintf(s->fpout, "%10d%10.1f%10d", comp->nMol, comp->nMolAv - comp->nMolReq, comp->inflow_net);
 303         }
 304     }
 305
 306     // Output center of split groups
 307     fprintf(s->fpout,
 308             "%10g%10g",
 309             s->group[static_cast<int>(SwapGroupSplittingType::Split0)].center[s->swapdim],
 310             s->group[static_cast<int>(SwapGroupSplittingType::Split1)].center[s->swapdim]);
 311
 312     // Output ion flux for each channel and ion type
 313     for (int iChan = 0; iChan < eChanNR; iChan++)
 314     {
 315         for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
 316         {
 317             t_swapgrp* g = &s->group[ig];
 318             fprintf(s->fpout, "%10d", g->fluxfromAtoB[iChan]);
 319         }
 320     }
 321
 322     /* Output the number of molecules that leaked from A to B */
 323     fprintf(s->fpout, "%10d", s->fluxleak);
 324
 325     fprintf(s->fpout, "%s\n", comment);
 326 }
 327
 328
 329 /*! \brief Get the center of a group of nat atoms.
 330  *
 331  * Since with PBC an atom group might not be whole, use the first atom as the
 332  * reference atom and determine the center with respect to this reference.
 333  */
 334 static void get_molecule_center(rvec x[], int nat, const real* weights, rvec center, t_pbc* pbc)
 335 {
 336     int  i;
 337     rvec weightedPBCimage;
 338     real wi, wsum;
 339     rvec reference, correctPBCimage, dx;
 340
 341
 342     /* Use the first atom as the reference and put other atoms near that one */
 343     /* This does not work for large molecules that span > half of the box! */
 344     copy_rvec(x[0], reference);
 345
 346     /* Calculate either the weighted center or simply the center of geometry */
 347     wsum = 0.0;
 348     clear_rvec(center);
 349     for (i = 0; i < nat; i++)
 350     {
 351         /* PBC distance between position and reference */
 352         pbc_dx(pbc, x[i], reference, dx);
 353
 354         /* Add PBC distance to reference */
 355         rvec_add(reference, dx, correctPBCimage);
 356
 357         /* Take weight into account */
 358         if (nullptr == weights)
 359         {
 360             wi = 1.0;
 361         }
 362         else
 363         {
 364             wi = weights[i];
 365         }
 366         wsum += wi;
 367         svmul(wi, correctPBCimage, weightedPBCimage);
 368
 369         /* Add to center */
 370         rvec_inc(center, weightedPBCimage);
 371     }
 372
 373     /* Normalize */
 374     svmul(1.0 / wsum, center, center);
 375 }
 376
 377
 378 /*! \brief Return TRUE if position x of ion (or water) is found in the compartment,
 379  * i.e. between w1 and w2.
 380  *
 381  * One can define and additional offset "b" if one wants to exchange ions/water
 382  * to or from a plane not directly in the middle of w1 and w2. The offset can be
 383  * in  ]-1.0, ..., +1.0 [.
 384  * A bulkOffset of 0.0 means 'no offset', so the swap-layer is directly in the
 385  * middle between w1 and w2. Offsets -1.0 < b <  0.0 will yield swaps nearer to w1,
 386  * whereas offsets  0.0 < 0 < +1.0 will yield swaps nearer to w2.
 387  *
 388  * \code
 389  *
 390  * ||--------------+-------------|-------------+------------------------||
 391  *                w1  ? ? ? ? ? ? ? ? ? ? ?   w2
 392  * ||--------------+-------------|----b--------+------------------------||
 393  *                -1            0.0           +1
 394  *
 395  * \endcode
 396  *
 397  * \param[in]  w1               Position of 'wall' atom 1.
 398  * \param[in]  w2               Position of 'wall' atom 2.
 399  * \param[in]  x                Position of the ion or the water molecule under consideration.
 400  * \param[in]  l                Length of the box, from || to || in the sketch.
 401  * \param[in]  bulkOffset       Where is the bulk layer "b" to be found between w1 and w2?
 402  * \param[out] distance_from_b  Distance of x to the bulk layer "b".
 403  *
 404  * \returns TRUE if x is between w1 and w2.
 405  *
 406  * Also computes the distance of x to the compartment center (the layer that is
 407  * normally situated in the middle of w1 and w2 that would be considered as having
 408  * the bulk concentration of ions).
 409  */
 410 static gmx_bool compartment_contains_atom(real w1, real w2, real x, real l, real bulkOffset, real* distance_from_b)
 411 {
 412     real m, l_2;
 413     real width;
 414
 415
 416     /* First set the origin in the middle of w1 and w2 */
 417     m = 0.5 * (w1 + w2);
 418     w1 -= m;
 419     w2 -= m;
 420     x -= m;
 421     width = w2 - w1;
 422
 423     /* Now choose the PBC image of x that is closest to the origin: */
 424     l_2 = 0.5 * l;
 425     while (x > l_2)
 426     {
 427         x -= l;
 428     }
 429     while (x <= -l_2)
 430     {
 431         x += l;
 432     }
 433
 434     *distance_from_b = static_cast<real>(fabs(x - bulkOffset * 0.5 * width));
 435
 436     /* Return TRUE if we now are in area "????" */
 437     return (x >= w1) && (x < w2);
 438 }
 439
 440
 441 /*! \brief Updates the time-averaged number of ions in a compartment. */
 442 static void update_time_window(t_compartment* comp, int values, int replace)
 443 {
 444     real average;
 445     int  i;
 446
 447
 448     /* Put in the new value */
 449     if (replace >= 0)
 450     {
 451         comp->nMolPast[replace] = comp->nMol;
 452     }
 453
 454     /* Compute the new time-average */
 455     average = 0.0;
 456     for (i = 0; i < values; i++)
 457     {
 458         average += comp->nMolPast[i];
 459     }
 460     average /= values;
 461     comp->nMolAv = average;
 462 }
 463
 464
 465 /*! \brief Add the atom with collective index ci to the atom list in compartment 'comp'.
 466  *
 467  * \param[in]     ci        Index of this ion in the collective xc array.
 468  * \param[inout]  comp      Compartment to add this atom to.
 469  * \param[in]     distance  Shortest distance of this atom to the bulk layer,
 470  *                          from which ion/water pairs are selected for swapping.
 471  */
 472 static void add_to_list(int ci, t_compartment* comp, real distance)
 473 {
 474     int nr = comp->nMol;
 475
 476     if (nr >= comp->nalloc)
 477     {
 478         comp->nalloc = over_alloc_dd(nr + 1);
 479         srenew(comp->ind, comp->nalloc);
 480         srenew(comp->dist, comp->nalloc);
 481     }
 482     comp->ind[nr]  = ci;
 483     comp->dist[nr] = distance;
 484     comp->nMol++;
 485 }
 486
 487
 488 /*! \brief Determine the compartment boundaries from the channel centers. */
 489 static void get_compartment_boundaries(int c, t_swap* s, const matrix box, real* left, real* right)
 490 {
 491     real pos0, pos1;
 492     real leftpos, rightpos, leftpos_orig;
 493
 494
 495     if (c >= eCompNR)
 496     {
 497         gmx_fatal(FARGS, "No compartment %c.", c + 'A');
 498     }
 499
 500     pos0 = s->group[static_cast<int>(SwapGroupSplittingType::Split0)].center[s->swapdim];
 501     pos1 = s->group[static_cast<int>(SwapGroupSplittingType::Split1)].center[s->swapdim];
 502
 503     if (pos0 < pos1)
 504     {
 505         leftpos  = pos0;
 506         rightpos = pos1;
 507     }
 508     else
 509     {
 510         leftpos  = pos1;
 511         rightpos = pos0;
 512     }
 513
 514     /* This gets us the other compartment: */
 515     if (c == eCompB)
 516     {
 517         leftpos_orig = leftpos;
 518         leftpos      = rightpos;
 519         rightpos     = leftpos_orig + box[s->swapdim][s->swapdim];
 520     }
 521
 522     *left  = leftpos;
 523     *right = rightpos;
 524 }
 525
 526
 527 /*! \brief Determine the per-channel ion flux.
 528  *
 529  * To determine the flux through the individual channels, we
 530  * remember the compartment and channel history of each ion. An ion can be
 531  * either in channel0 or channel1, or in the remaining volume of compartment
 532  * A or B.
 533  *
 534  * \code
 535  *    +-----------------+
 536  *    |                 | B
 537  *    |                 | B compartment
 538  *    ||||||||||0|||||||| bilayer with channel 0
 539  *    |                 | A
 540  *    |                 | A
 541  *    |                 | A compartment
 542  *    |                 | A
 543  *    |||||1||||||||||||| bilayer with channel 1
 544  *    |                 | B
 545  *    |                 | B compartment
 546  *    +-----------------+
 547  *
 548  * \endcode
 549  */
 550 static void detect_flux_per_channel(t_swapgrp*          g,
 551                                     int                 iAtom,
 552                                     int                 comp,
 553                                     rvec                atomPosition,
 554                                     unsigned char*      comp_now,
 555                                     unsigned char*      comp_from,
 556                                     unsigned char*      channel_label,
 557                                     const t_swapcoords* sc,
 558                                     t_swap*             s,
 559                                     real                cyl0_r2,
 560                                     real                cyl1_r2,
 561                                     int64_t             step,
 562                                     gmx_bool            bRerun,
 563                                     FILE*               fpout)
 564 {
 565     int      sd, chan_nr;
 566     gmx_bool in_cyl0, in_cyl1;
 567     char     buf[STRLEN];
 568
 569
 570     sd = s->swapdim;
 571
 572     /* Check whether ion is inside any of the channels */
 573     in_cyl0 = is_in_channel(atomPosition,
 574                             s->group[static_cast<int>(SwapGroupSplittingType::Split0)].center,
 575                             sc->cyl0u,
 576                             sc->cyl0l,
 577                             cyl0_r2,
 578                             s->pbc,
 579                             sd);
 580     in_cyl1 = is_in_channel(atomPosition,
 581                             s->group[static_cast<int>(SwapGroupSplittingType::Split1)].center,
 582                             sc->cyl1u,
 583                             sc->cyl1l,
 584                             cyl1_r2,
 585                             s->pbc,
 586                             sd);
 587
 588     if (in_cyl0 && in_cyl1)
 589     {
 590         /* Ion appears to be in both channels. Something is severely wrong! */
 591         g->nCylBoth++;
 592         *comp_now      = eDomainNotset;
 593         *comp_from     = eDomainNotset;
 594         *channel_label = eChHistPassedNone;
 595     }
 596     else if (in_cyl0)
 597     {
 598         /* Ion is in channel 0 now */
 599         *channel_label = eChHistPassedCh0;
 600         *comp_now      = eDomainNotset;
 601         g->nCyl[eChan0]++;
 602     }
 603     else if (in_cyl1)
 604     {
 605         /* Ion is in channel 1 now */
 606         *channel_label = eChHistPassedCh1;
 607         *comp_now      = eDomainNotset;
 608         g->nCyl[eChan1]++;
 609     }
 610     else
 611     {
 612         /* Ion is not in any of the channels, so it must be in domain A or B */
 613         if (eCompA == comp)
 614         {
 615             *comp_now = eDomainA;
 616         }
 617         else
 618         {
 619             *comp_now = eDomainB;
 620         }
 621     }
 622
 623     /* Only take action, if ion is now in domain A or B, and was before
 624      * in the other domain!
 625      */
 626     if (eDomainNotset == *comp_from)
 627     {
 628         /* Maybe we can set the domain now */
 629         *comp_from = *comp_now; /* Could still be eDomainNotset, though */
 630     }
 631     else if ((*comp_now != eDomainNotset) /* if in channel */
 632              && (*comp_from != *comp_now))
 633     {
 634         /* Obviously the ion changed its domain.
 635          * Count this for the channel through which it has passed. */
 636         switch (*channel_label)
 637         {
 638             case eChHistPassedNone:
 639                 ++s->fluxleak;
 640
 641                 fprintf(stderr,
 642                         " %s Warning! Step %s, ion %d moved from %s to %s\n",
 643                         SwS,
 644                         gmx_step_str(step, buf),
 645                         iAtom,
 646                         DomainString[*comp_from],
 647                         DomainString[*comp_now]);
 648                 if (bRerun)
 649                 {
 650                     fprintf(stderr, ", possibly due to a swap in the original simulation.\n");
 651                 }
 652                 else
 653                 {
 654                     fprintf(stderr,
 655                             "but did not pass cyl0 or cyl1 as defined in the .mdp file.\n"
 656                             "Do you have an ion somewhere within the membrane?\n");
 657                     /* Write this info to the CompEL output file: */
 658                     fprintf(s->fpout,
 659                             " # Warning: step %s, ion %d moved from %s to %s (probably through the "
 660                             "membrane)\n",
 661                             gmx_step_str(step, buf),
 662                             iAtom,
 663                             DomainString[*comp_from],
 664                             DomainString[*comp_now]);
 665                 }
 666                 break;
 667             case eChHistPassedCh0:
 668             case eChHistPassedCh1:
 669                 if (*channel_label == eChHistPassedCh0)
 670                 {
 671                     chan_nr = 0;
 672                 }
 673                 else
 674                 {
 675                     chan_nr = 1;
 676                 }
 677
 678                 if (eDomainA == *comp_from)
 679                 {
 680                     g->fluxfromAtoB[chan_nr]++;
 681                 }
 682                 else
 683                 {
 684                     g->fluxfromAtoB[chan_nr]--;
 685                 }
 686                 fprintf(fpout, "# Atom nr. %d finished passing %s.\n", iAtom, ChannelString[*channel_label]);
 687                 break;
 688             default:
 689                 gmx_fatal(FARGS, "%s Unknown channel history entry for ion type '%s'\n", SwS, g->molname);
 690         }
 691
 692         /* This ion has moved to the _other_ compartment ... */
 693         *comp_from = *comp_now;
 694         /* ... and it did not pass any channel yet */
 695         *channel_label = eChHistPassedNone;
 696     }
 697 }
 698
 699
 700 /*! \brief Determines which ions or solvent molecules are in compartment A and B */
 701 static void sortMoleculesIntoCompartments(t_swapgrp*          g,
 702                                           t_commrec*          cr,
 703                                           const t_swapcoords* sc,
 704                                           t_swap*             s,
 705                                           const matrix        box,
 706                                           int64_t             step,
 707                                           FILE*               fpout,
 708                                           gmx_bool            bRerun,
 709                                           gmx_bool            bIsSolvent)
 710 {
 711     int  nMolNotInComp[eCompNR]; /* consistency check */
 712     real cyl0_r2 = sc->cyl0r * sc->cyl0r;
 713     real cyl1_r2 = sc->cyl1r * sc->cyl1r;
 714
 715     /* Get us a counter that cycles in the range of [0 ... sc->nAverage[ */
 716     int replace = (step / sc->nstswap) % sc->nAverage;
 717
 718     for (int comp = eCompA; comp <= eCompB; comp++)
 719     {
 720         real left, right;
 721
 722         /* Get lists of atoms that match criteria for this compartment */
 723         get_compartment_boundaries(comp, s, box, &left, &right);
 724
 725         /* First clear the ion molecule lists */
 726         g->comp[comp].nMol  = 0;
 727         nMolNotInComp[comp] = 0; /* consistency check */
 728
 729         /* Loop over the molecules and atoms of this group */
 730         for (int iMol = 0, iAtom = 0; iAtom < static_cast<int>(g->atomset.numAtomsGlobal());
 731              iAtom += g->apm, iMol++)
 732         {
 733             real dist;
 734             int  sd = s->swapdim;
 735
 736             /* Is this first atom of the molecule in the compartment that we look at? */
 737             if (compartment_contains_atom(
 738                         left, right, g->xc[iAtom][sd], box[sd][sd], sc->bulkOffset[comp], &dist))
 739             {
 740                 /* Add the first atom of this molecule to the list of molecules in this compartment */
 741                 add_to_list(iAtom, &g->comp[comp], dist);
 742
 743                 /* Master also checks for ion groups through which channel each ion has passed */
 744                 if (MASTER(cr) && (g->comp_now != nullptr) && !bIsSolvent)
 745                 {
 746                     int globalAtomNr = g->atomset.globalIndex()[iAtom] + 1; /* PDB index starts at 1 ... */
 747                     detect_flux_per_channel(g,
 748                                             globalAtomNr,
 749                                             comp,
 750                                             g->xc[iAtom],
 751                                             &g->comp_now[iMol],
 752                                             &g->comp_from[iMol],
 753                                             &g->channel_label[iMol],
 754                                             sc,
 755                                             s,
 756                                             cyl0_r2,
 757                                             cyl1_r2,
 758                                             step,
 759                                             bRerun,
 760                                             fpout);
 761                 }
 762             }
 763             else
 764             {
 765                 nMolNotInComp[comp]++;
 766             }
 767         }
 768         /* Correct the time-averaged number of ions in the compartment */
 769         if (!bIsSolvent)
 770         {
 771             update_time_window(&g->comp[comp], sc->nAverage, replace);
 772         }
 773     }
 774
 775     /* Flux detection warnings */
 776     if (MASTER(cr) && !bIsSolvent)
 777     {
 778         if (g->nCylBoth > 0)
 779         {
 780             fprintf(stderr,
 781                     "\n"
 782                     "%s Warning: %d atoms were detected as being in both channels! Probably your "
 783                     "split\n"
 784                     "%s          cylinder is way too large, or one compartment has collapsed (step "
 785                     "%" PRId64 ")\n",
 786                     SwS,
 787                     g->nCylBoth,
 788                     SwS,
 789                     step);
 790
 791             fprintf(s->fpout, "Warning: %d atoms were assigned to both channels!\n", g->nCylBoth);
 792
 793             g->nCylBoth = 0;
 794         }
 795     }
 796
 797     if (bIsSolvent && nullptr != fpout)
 798     {
 799         fprintf(fpout,
 800                 "# Solv. molecules in comp.%s: %d   comp.%s: %d\n",
 801                 CompStr[eCompA],
 802                 g->comp[eCompA].nMol,
 803                 CompStr[eCompB],
 804                 g->comp[eCompB].nMol);
 805     }
 806
 807     /* Consistency checks */
 808     const auto numMolecules = static_cast<int>(g->atomset.numAtomsGlobal() / g->apm);
 809     if (nMolNotInComp[eCompA] + nMolNotInComp[eCompB] != numMolecules)
 810     {
 811         fprintf(stderr,
 812                 "%s Warning: Inconsistency while assigning '%s' molecules to compartments. !inA: "
 813                 "%d, !inB: %d, total molecules %d\n",
 814                 SwS,
 815                 g->molname,
 816                 nMolNotInComp[eCompA],
 817                 nMolNotInComp[eCompB],
 818                 numMolecules);
 819     }
 820
 821     int sum = g->comp[eCompA].nMol + g->comp[eCompB].nMol;
 822     if (sum != numMolecules)
 823     {
 824         fprintf(stderr,
 825                 "%s Warning: %d molecules are in group '%s', but altogether %d have been assigned "
 826                 "to the compartments.\n",
 827                 SwS,
 828                 numMolecules,
 829                 g->molname,
 830                 sum);
 831     }
 832 }
 833
 834
 835 /*! \brief Find out how many group atoms are in the compartments initially */
 836 static void get_initial_ioncounts(const t_inputrec* ir,
 837                                   t_swap*           s,
 838                                   const rvec        x[], /* the initial positions */
 839                                   const matrix      box,
 840                                   t_commrec*        cr,
 841                                   gmx_bool          bRerun)
 842 {
 843     t_swapcoords* sc;
 844     t_swapgrp*    g;
 845
 846     sc = ir->swap;
 847
 848     /* Loop over the user-defined (ion) groups */
 849     for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
 850     {
 851         g = &s->group[ig];
 852
 853         /* Copy the initial positions of the atoms in the group
 854          * to the collective array so that we can compartmentalize */
 855         for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
 856         {
 857             int ind = g->atomset.globalIndex()[i];
 858             copy_rvec(x[ind], g->xc[i]);
 859         }
 860
 861         /* Set up the compartments and get lists of atoms in each compartment */
 862         sortMoleculesIntoCompartments(g, cr, sc, s, box, 0, s->fpout, bRerun, FALSE);
 863
 864         /* Set initial molecule counts if requested (as signaled by "-1" value) */
 865         for (int ic = 0; ic < eCompNR; ic++)
 866         {
 867             int requested = sc->grp[ig].nmolReq[ic];
 868             if (requested < 0)
 869             {
 870                 g->comp[ic].nMolReq = g->comp[ic].nMol;
 871             }
 872             else
 873             {
 874                 g->comp[ic].nMolReq = requested;
 875             }
 876         }
 877
 878         /* Check whether the number of requested molecules adds up to the total number */
 879         int req = g->comp[eCompA].nMolReq + g->comp[eCompB].nMolReq;
 880         int tot = g->comp[eCompA].nMol + g->comp[eCompB].nMol;
 881
 882         if ((req != tot))
 883         {
 884             gmx_fatal(FARGS,
 885                       "Mismatch of the number of %s ions summed over both compartments.\n"
 886                       "You requested a total of %d ions (%d in A and %d in B),\n"
 887                       "but there are a total of %d ions of this type in the system.\n",
 888                       g->molname,
 889                       req,
 890                       g->comp[eCompA].nMolReq,
 891                       g->comp[eCompB].nMolReq,
 892                       tot);
 893         }
 894
 895         /* Initialize time-averaging:
 896          * Write initial concentrations to all time bins to start with */
 897         for (int ic = 0; ic < eCompNR; ic++)
 898         {
 899             g->comp[ic].nMolAv = g->comp[ic].nMol;
 900             for (int i = 0; i < sc->nAverage; i++)
 901             {
 902                 g->comp[ic].nMolPast[i] = g->comp[ic].nMol;
 903             }
 904         }
 905     }
 906 }
 907
 908
 909 /*! \brief Copy history of ion counts from checkpoint file.
 910  *
 911  * When called, the checkpoint file has already been read in. Here we copy
 912  * over the values from .cpt file to the swap data structure.
 913  */
 914 static void get_initial_ioncounts_from_cpt(const t_inputrec* ir,
 915                                            t_swap*           s,
 916                                            swaphistory_t*    swapstate,
 917                                            t_commrec*        cr,
 918                                            gmx_bool          bVerbose)
 919 {
 920     t_swapcoords*    sc;
 921     t_swapgrp*       g;
 922     swapstateIons_t* gs;
 923
 924     sc = ir->swap;
 925
 926     if (MASTER(cr))
 927     {
 928         /* Copy the past values from the checkpoint values that have been read in already */
 929         if (bVerbose)
 930         {
 931             fprintf(stderr, "%s Copying values from checkpoint\n", SwS);
 932         }
 933
 934         for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
 935         {
 936             g  = &s->group[ig];
 937             gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
 938
 939             for (int ic = 0; ic < eCompNR; ic++)
 940             {
 941                 g->comp[ic].nMolReq    = gs->nMolReq[ic];
 942                 g->comp[ic].inflow_net = gs->inflow_net[ic];
 943
 944                 if (bVerbose)
 945                 {
 946                     fprintf(stderr,
 947                             "%s ... Influx netto: %d   Requested: %d   Past values: ",
 948                             SwS,
 949                             g->comp[ic].inflow_net,
 950                             g->comp[ic].nMolReq);
 951                 }
 952
 953                 for (int j = 0; j < sc->nAverage; j++)
 954                 {
 955                     g->comp[ic].nMolPast[j] = gs->nMolPast[ic][j];
 956                     if (bVerbose)
 957                     {
 958                         fprintf(stderr, "%d ", g->comp[ic].nMolPast[j]);
 959                     }
 960                 }
 961                 if (bVerbose)
 962                 {
 963                     fprintf(stderr, "\n");
 964                 }
 965             }
 966         }
 967     }
 968 }
 969
 970
 971 /*! \brief The master lets all others know about the initial ion counts. */
 972 static void bc_initial_concentrations(t_commrec* cr, t_swapcoords* swap, t_swap* s)
 973 {
 974     int        ic, ig;
 975     t_swapgrp* g;
 976
 977
 978     for (ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
 979     {
 980         g = &s->group[ig];
 981
 982         for (ic = 0; ic < eCompNR; ic++)
 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);
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,
1038                   nMultiple,
1039                   (1 == nMultiple) ? " is" : "s are",
1040                   SwSEmpty,
1041                   SwSEmpty);
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,
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, eCompNR * nIonTypes * 3 + 2 + eChanNR * nIonTypes + 1);
1099
1100     // Number of molecules and difference to reference counts for each
1101     // compartment and ion type
1102     for (int ic = count = 0; ic < eCompNR; ic++)
1103     {
1104         for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1105         {
1106             t_swapGroup* g = &ir->swap->grp[ig];
1107             real         q = s->group[ig].q;
1108
1109             snprintf(buf, STRLEN, "%s %s ions (charge %s%g)", CompStr[ic], g->molname, q > 0 ? "+" : "", q);
1110             legend[count++] = gmx_strdup(buf);
1111
1112             snprintf(buf,
1113                      STRLEN,
1114                      "%s av. mismatch to %d %s ions",
1115                      CompStr[ic],
1116                      s->group[ig].comp[ic].nMolReq,
1117                      g->molname);
1118             legend[count++] = gmx_strdup(buf);
1119
1120             snprintf(buf, STRLEN, "%s net %s ion influx", CompStr[ic], g->molname);
1121             legend[count++] = gmx_strdup(buf);
1122         }
1123     }
1124
1125     // Center of split groups
1126     snprintf(buf,
1127              STRLEN,
1128              "%scenter of %s of split group 0",
1129              SwapStr[ir->eSwapCoords],
1130              (nullptr != s->group[static_cast<int>(SwapGroupSplittingType::Split0)].m) ? "mass" : "geometry");
1131     legend[count++] = gmx_strdup(buf);
1132     snprintf(buf,
1133              STRLEN,
1134              "%scenter of %s of split group 1",
1135              SwapStr[ir->eSwapCoords],
1136              (nullptr != s->group[static_cast<int>(SwapGroupSplittingType::Split1)].m) ? "mass" : "geometry");
1137     legend[count++] = gmx_strdup(buf);
1138
1139     // Ion flux for each channel and ion type
1140     for (int ic = 0; ic < eChanNR; ic++)
1141     {
1142         for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1143         {
1144             t_swapGroup* g = &ir->swap->grp[ig];
1145             snprintf(buf, STRLEN, "A->ch%d->B %s permeations", ic, g->molname);
1146             legend[count++] = gmx_strdup(buf);
1147         }
1148     }
1149
1150     // Number of molecules that leaked from A to B
1151     snprintf(buf, STRLEN, "leakage");
1152     legend[count++] = gmx_strdup(buf);
1153
1154     xvgr_legend(s->fpout, count, legend, oenv);
1155
1156     fprintf(s->fpout,
1157             "# Instantaneous ion counts and time-averaged differences to requested numbers\n");
1158
1159     // We add a simple text legend helping to identify the columns with xvgr legend strings
1160     fprintf(s->fpout, "#  time (ps)");
1161     for (int i = 0; i < count; i++)
1162     {
1163         snprintf(buf, STRLEN, "s%d", i);
1164         fprintf(s->fpout, "%10s", buf);
1165     }
1166     fprintf(s->fpout, "\n");
1167     fflush(s->fpout);
1168 }
1169
1170
1171 /*! \brief Initialize channel ion flux detection routine.
1172  *
1173  * Initialize arrays that keep track of where the ions come from and where
1174  * they go.
1175  */
1176 static void detect_flux_per_channel_init(t_swap* s, swaphistory_t* swapstate, const bool isRestart)
1177 {
1178     t_swapgrp*       g;
1179     swapstateIons_t* gs;
1180
1181     /* All these flux detection routines run on the master only */
1182     if (swapstate == nullptr)
1183     {
1184         return;
1185     }
1186
1187     for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1188     {
1189         g  = &s->group[ig];
1190         gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1191
1192         /******************************************************/
1193         /* Channel and domain history for the individual ions */
1194         /******************************************************/
1195         if (isRestart) /* set the pointers right */
1196         {
1197             g->comp_from     = gs->comp_from;
1198             g->channel_label = gs->channel_label;
1199         }
1200         else /* allocate memory for molecule counts */
1201         {
1202             snew(g->comp_from, g->atomset.numAtomsGlobal() / g->apm);
1203             gs->comp_from = g->comp_from;
1204             snew(g->channel_label, g->atomset.numAtomsGlobal() / g->apm);
1205             gs->channel_label = g->channel_label;
1206         }
1207         snew(g->comp_now, g->atomset.numAtomsGlobal() / g->apm);
1208
1209         /* Initialize the channel and domain history counters */
1210         for (size_t i = 0; i < g->atomset.numAtomsGlobal() / g->apm; i++)
1211         {
1212             g->comp_now[i] = eDomainNotset;
1213             if (!isRestart)
1214             {
1215                 g->comp_from[i]     = eDomainNotset;
1216                 g->channel_label[i] = eChHistPassedNone;
1217             }
1218         }
1219
1220         /************************************/
1221         /* Channel fluxes for both channels */
1222         /************************************/
1223         g->nCyl[eChan0] = 0;
1224         g->nCyl[eChan1] = 0;
1225         g->nCylBoth     = 0;
1226     }
1227
1228     if (isRestart)
1229     {
1230         fprintf(stderr, "%s Copying channel fluxes from checkpoint file data\n", SwS);
1231     }
1232
1233
1234     // Loop over ion types (and both channels)
1235     for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1236     {
1237         g  = &s->group[ig];
1238         gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1239
1240         for (int ic = 0; ic < eChanNR; ic++)
1241         {
1242             fprintf(stderr, "%s Channel %d flux history for ion type %s (charge %g): ", SwS, ic, g->molname, g->q);
1243             if (isRestart)
1244             {
1245                 g->fluxfromAtoB[ic] = gs->fluxfromAtoB[ic];
1246             }
1247             else
1248             {
1249                 g->fluxfromAtoB[ic] = 0;
1250             }
1251
1252             fprintf(stderr, "%d molecule%s", g->fluxfromAtoB[ic], g->fluxfromAtoB[ic] == 1 ? "" : "s");
1253             fprintf(stderr, "\n");
1254         }
1255     }
1256
1257     /* Set pointers for checkpoint writing */
1258     swapstate->fluxleak_p = &s->fluxleak;
1259     for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1260     {
1261         g  = &s->group[ig];
1262         gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1263
1264         for (int ic = 0; ic < eChanNR; ic++)
1265         {
1266             gs->fluxfromAtoB_p[ic] = &g->fluxfromAtoB[ic];
1267         }
1268     }
1269 }
1270
1271
1272 /*! \brief Outputs the initial structure to PDB file for debugging reasons.
1273  *
1274  * Output the starting structure so that in case of multimeric channels
1275  * the user can check whether we have the correct PBC image for all atoms.
1276  * If this is not correct, the ion counts per channel will be very likely
1277  * wrong.
1278  */
1279 static void outputStartStructureIfWanted(const gmx_mtop_t& mtop, rvec* x, PbcType pbcType, const matrix box)
1280 {
1281     char* env = getenv("GMX_COMPELDUMP");
1282
1283     if (env != nullptr)
1284     {
1285         fprintf(stderr,
1286                 "\n%s Found env.var. GMX_COMPELDUMP, will output CompEL starting structure made "
1287                 "whole.\n"
1288                 "%s In case of multimeric channels, please check whether they have the correct PBC "
1289                 "representation.\n",
1290                 SwS,
1291                 SwSEmpty);
1292
1293         write_sto_conf_mtop(
1294                 "CompELAssumedWholeConfiguration.pdb", *mtop.name, mtop, x, nullptr, pbcType, box);
1295     }
1296 }
1297
1298
1299 /*! \brief Initialize the swapstate structure, used for checkpoint writing.
1300  *
1301  * The swapstate struct stores the information we need to make the channels
1302  * whole again after restarts from a checkpoint file. Here we do the following:
1303  * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,
1304  * b) if we did start from .cpt, we copy over the last whole structures from .cpt,
1305  * c) in any case, for subsequent checkpoint writing, we set the pointers in
1306  * swapstate to the x_old arrays, which contain the correct PBC representation of
1307  * multimeric channels at the last time step.
1308  */
1309 static void init_swapstate(swaphistory_t*    swapstate,
1310                            t_swapcoords*     sc,
1311                            t_swap*           s,
1312                            const gmx_mtop_t& mtop,
1313                            const rvec*       x, /* the initial positions */
1314                            const matrix      box,
1315                            const t_inputrec* ir)
1316 {
1317     rvec*      x_pbc = nullptr; /* positions of the whole MD system with molecules made whole */
1318     t_swapgrp* g;
1319
1320
1321     /* We always need the last whole positions such that
1322      * in the next time step we can make the channels whole again in PBC */
1323     if (swapstate->bFromCpt)
1324     {
1325         /* Copy the last whole positions of each channel from .cpt */
1326         g = &(s->group[static_cast<int>(SwapGroupSplittingType::Split0)]);
1327         for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1328         {
1329             copy_rvec(swapstate->xc_old_whole[eChan0][i], g->xc_old[i]);
1330         }
1331         g = &(s->group[static_cast<int>(SwapGroupSplittingType::Split1)]);
1332         for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1333         {
1334             copy_rvec(swapstate->xc_old_whole[eChan1][i], g->xc_old[i]);
1335         }
1336     }
1337     else
1338     {
1339         swapstate->eSwapCoords = ir->eSwapCoords;
1340
1341         /* Set the number of ion types and allocate memory for checkpointing */
1342         swapstate->nIonTypes = s->ngrp - static_cast<int>(SwapGroupSplittingType::Count);
1343         snew(swapstate->ionType, swapstate->nIonTypes);
1344
1345         /* Store the total number of ions of each type in the swapstateIons
1346          * structure that is accessible during checkpoint writing */
1347         for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1348         {
1349             swapstateIons_t* gs = &swapstate->ionType[ii];
1350             gs->nMol            = sc->grp[ii + static_cast<int>(SwapGroupSplittingType::Count)].nat;
1351         }
1352
1353         /* Extract the initial split group positions. */
1354
1355         /* Remove pbc, make molecule whole. */
1356         snew(x_pbc, mtop.natoms);
1357         copy_rvecn(x, x_pbc, 0, mtop.natoms);
1358
1359         /* This can only make individual molecules whole, not multimers */
1360         do_pbc_mtop(ir->pbcType, box, &mtop, x_pbc);
1361
1362         /* Output the starting structure? */
1363         outputStartStructureIfWanted(mtop, x_pbc, ir->pbcType, box);
1364
1365         /* If this is the first run (i.e. no checkpoint present) we assume
1366          * that the starting positions give us the correct PBC representation */
1367         for (int ig = static_cast<int>(SwapGroupSplittingType::Split0);
1368              ig <= static_cast<int>(SwapGroupSplittingType::Split1);
1369              ig++)
1370         {
1371             g = &(s->group[ig]);
1372             for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1373             {
1374                 copy_rvec(x_pbc[g->atomset.globalIndex()[i]], g->xc_old[i]);
1375             }
1376         }
1377         sfree(x_pbc);
1378
1379         /* Prepare swapstate arrays for later checkpoint writing */
1380         swapstate->nat[eChan0] =
1381                 s->group[static_cast<int>(SwapGroupSplittingType::Split0)].atomset.numAtomsGlobal();
1382         swapstate->nat[eChan1] =
1383                 s->group[static_cast<int>(SwapGroupSplittingType::Split1)].atomset.numAtomsGlobal();
1384     }
1385
1386     /* For subsequent checkpoint writing, set the swapstate pointers to the xc_old
1387      * arrays that get updated at every swapping step */
1388     swapstate->xc_old_whole_p[eChan0] = &s->group[static_cast<int>(SwapGroupSplittingType::Split0)].xc_old;
1389     swapstate->xc_old_whole_p[eChan1] = &s->group[static_cast<int>(SwapGroupSplittingType::Split1)].xc_old;
1390 }
1391
1392 /*! \brief Determine the total charge imbalance resulting from the swap groups */
1393 static real getRequestedChargeImbalance(t_swap* s)
1394 {
1395     int        ig;
1396     real       DeltaQ = 0.0;
1397     t_swapgrp* g;
1398     real       particle_charge;
1399     real       particle_number[eCompNR];
1400
1401     //        s->deltaQ =  ( (-1) * s->comp[eCompA][eIonNEG].nat_req + s->comp[eCompA][eIonPOS].nat_req )
1402     //                   - ( (-1) * s->comp[eCompB][eIonNEG].nat_req + s->comp[eCompB][eIonPOS].nat_req );
1403
1404     for (ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1405     {
1406         g = &s->group[ig];
1407
1408         particle_charge         = g->q;
1409         particle_number[eCompA] = g->comp[eCompA].nMolReq;
1410         particle_number[eCompB] = g->comp[eCompB].nMolReq;
1411
1412         DeltaQ += particle_charge * (particle_number[eCompA] - particle_number[eCompB]);
1413     }
1414
1415     return DeltaQ;
1416 }
1417
1418
1419 /*! \brief Sorts anions and cations into two separate groups
1420  *
1421  * This routine should be called for the 'anions' and 'cations' group,
1422  * of which the indices were lumped together in the older version of the code.
1423  */
1424 static void copyIndicesToGroup(const int* indIons, int nIons, t_swapGroup* g, t_commrec* cr)
1425 {
1426     g->nat = nIons;
1427
1428     /* If explicit ion counts were requested in the .mdp file
1429      * (by setting positive values for the number of ions),
1430      * we can make an additional consistency check here */
1431     if ((g->nmolReq[eCompA] < 0) && (g->nmolReq[eCompB] < 0))
1432     {
1433         if (g->nat != (g->nmolReq[eCompA] + g->nmolReq[eCompB]))
1434         {
1435             gmx_fatal_collective(FARGS,
1436                                  cr->mpi_comm_mysim,
1437                                  MASTER(cr),
1438                                  "%s Inconsistency while importing swap-related data from an old "
1439                                  "input file version.\n"
1440                                  "%s The requested ion counts in compartments A (%d) and B (%d)\n"
1441                                  "%s do not add up to the number of ions (%d) of this type for the "
1442                                  "group '%s'.\n",
1443                                  SwS,
1444                                  SwSEmpty,
1445                                  g->nmolReq[eCompA],
1446                                  g->nmolReq[eCompB],
1447                                  SwSEmpty,
1448                                  g->nat,
1449                                  g->molname);
1450         }
1451     }
1452
1453     srenew(g->ind, g->nat);
1454     for (int i = 0; i < g->nat; i++)
1455     {
1456         g->ind[i] = indIons[i];
1457     }
1458 }
1459
1460
1461 /*! \brief Converts old .tpr file CompEL contents to new data layout.
1462  *
1463  *  If we have read an old .tpr file (tpxv <= tpxv_CompElPolyatomicIonsAndMultipleIonTypes),
1464  * anions and cations are stored together in group #3. In the new
1465  * format we store each ion type in a separate group.
1466  * The 'classic' groups are:
1467  * #0 split group 0  - OK
1468  * #1 split group 1  - OK
1469  * #2 solvent        - OK
1470  * #3 anions         - contains also cations, needs to be converted
1471  * #4 cations        - empty before conversion
1472  *
1473  */
1474 static void convertOldToNewGroupFormat(t_swapcoords* sc, const gmx_mtop_t& mtop, gmx_bool bVerbose, t_commrec* cr)
1475 {
1476     t_swapGroup* g = &sc->grp[3];
1477
1478     /* Loop through the atom indices of group #3 (anions) and put all indices
1479      * that belong to cations into the cation group.
1480      */
1481     int  nAnions    = 0;
1482     int  nCations   = 0;
1483     int* indAnions  = nullptr;
1484     int* indCations = nullptr;
1485     snew(indAnions, g->nat);
1486     snew(indCations, g->nat);
1487
1488     int molb = 0;
1489     for (int i = 0; i < g->nat; i++)
1490     {
1491         const t_atom& atom = mtopGetAtomParameters(mtop, g->ind[i], &molb);
1492         if (atom.q < 0)
1493         {
1494             // This is an anion, add it to the list of anions
1495             indAnions[nAnions++] = g->ind[i];
1496         }
1497         else
1498         {
1499             // This is a cation, add it to the list of cations
1500             indCations[nCations++] = g->ind[i];
1501         }
1502     }
1503
1504     if (bVerbose)
1505     {
1506         fprintf(stdout,
1507                 "%s Sorted %d ions into separate groups of %d anions and %d cations.\n",
1508                 SwS,
1509                 g->nat,
1510                 nAnions,
1511                 nCations);
1512     }
1513
1514
1515     /* Now we have the correct lists of anions and cations.
1516      * Copy it to the right groups.
1517      */
1518     copyIndicesToGroup(indAnions, nAnions, g, cr);
1519     g = &sc->grp[4];
1520     copyIndicesToGroup(indCations, nCations, g, cr);
1521     sfree(indAnions);
1522     sfree(indCations);
1523 }
1524
1525
1526 /*! \brief Returns TRUE if we started from an old .tpr
1527  *
1528  * Then we need to re-sort anions and cations into separate groups */
1529 static gmx_bool bConvertFromOldTpr(t_swapcoords* sc)
1530 {
1531     // If the last group has no atoms it means we need to convert!
1532     return (sc->ngrp >= 5) && (0 == sc->grp[4].nat);
1533 }
1534
1535
1536 t_swap* init_swapcoords(FILE*                       fplog,
1537                         const t_inputrec*           ir,
1538                         const char*                 fn,
1539                         const gmx_mtop_t&           mtop,
1540                         const t_state*              globalState,
1541                         ObservablesHistory*         oh,
1542                         t_commrec*                  cr,
1543                         gmx::LocalAtomSetManager*   atomSets,
1544                         const gmx_output_env_t*     oenv,
1545                         const gmx::MdrunOptions&    mdrunOptions,
1546                         const gmx::StartingBehavior startingBehavior)
1547 {
1548     t_swapgrp*       g;
1549     swapstateIons_t* gs;
1550     swaphistory_t*   swapstate = nullptr;
1551
1552     if ((PAR(cr)) && !DOMAINDECOMP(cr))
1553     {
1554         gmx_fatal(FARGS, "Position swapping is only implemented for domain decomposition!");
1555     }
1556
1557     auto sc = ir->swap;
1558     auto s  = new t_swap();
1559
1560     if (mdrunOptions.rerun)
1561     {
1562         if (PAR(cr))
1563         {
1564             gmx_fatal(FARGS,
1565                       "%s This module does not support reruns in parallel\nPlease request a serial "
1566                       "run with -nt 1 / -np 1\n",
1567                       SwS);
1568         }
1569
1570         fprintf(stderr, "%s Rerun - using every available frame\n", SwS);
1571         sc->nstswap  = 1;
1572         sc->nAverage = 1; /* averaging makes no sense for reruns */
1573     }
1574
1575     if (MASTER(cr) && startingBehavior == gmx::StartingBehavior::NewSimulation)
1576     {
1577         fprintf(fplog, "\nInitializing ion/water position exchanges\n");
1578         please_cite(fplog, "Kutzner2011b");
1579     }
1580
1581     switch (ir->eSwapCoords)
1582     {
1583         case SwapType::X: s->swapdim = XX; break;
1584         case SwapType::Y: s->swapdim = YY; break;
1585         case SwapType::Z: s->swapdim = ZZ; break;
1586         default: s->swapdim = -1; break;
1587     }
1588
1589     const gmx_bool bVerbose = mdrunOptions.verbose;
1590
1591     // For compatibility with old .tpr files
1592     if (bConvertFromOldTpr(sc))
1593     {
1594         convertOldToNewGroupFormat(sc, mtop, bVerbose && MASTER(cr), cr);
1595     }
1596
1597     /* Copy some data and pointers to the group structures for convenience */
1598     /* Number of atoms in the group */
1599     s->ngrp = sc->ngrp;
1600     for (int i = 0; i < s->ngrp; i++)
1601     {
1602         s->group.emplace_back(atomSets->add(
1603                 gmx::ArrayRef<const int>(sc->grp[i].ind, sc->grp[i].ind + sc->grp[i].nat)));
1604         s->group[i].molname = sc->grp[i].molname;
1605     }
1606
1607     /* Check for overlapping atoms */
1608     check_swap_groups(s, mtop.natoms, bVerbose && MASTER(cr));
1609
1610     /* Allocate space for the collective arrays for all groups */
1611     /* For the collective position array */
1612     for (int i = 0; i < s->ngrp; i++)
1613     {
1614         g = &s->group[i];
1615         snew(g->xc, g->atomset.numAtomsGlobal());
1616
1617         /* For the split groups (the channels) we need some extra memory to
1618          * be able to make the molecules whole even if they span more than
1619          * half of the box size. */
1620         if ((i == static_cast<int>(SwapGroupSplittingType::Split0))
1621             || (i == static_cast<int>(SwapGroupSplittingType::Split1)))
1622         {
1623             snew(g->xc_shifts, g->atomset.numAtomsGlobal());
1624             snew(g->xc_eshifts, g->atomset.numAtomsGlobal());
1625             snew(g->xc_old, g->atomset.numAtomsGlobal());
1626         }
1627     }
1628
1629     if (MASTER(cr))
1630     {
1631         if (oh->swapHistory == nullptr)
1632         {
1633             oh->swapHistory = std::make_unique<swaphistory_t>(swaphistory_t{});
1634         }
1635         swapstate = oh->swapHistory.get();
1636
1637         init_swapstate(swapstate, sc, s, mtop, globalState->x.rvec_array(), globalState->box, ir);
1638     }
1639
1640     /* After init_swapstate we have a set of (old) whole positions for our
1641      * channels. Now transfer that to all nodes */
1642     if (PAR(cr))
1643     {
1644         for (int ig = static_cast<int>(SwapGroupSplittingType::Split0);
1645              ig <= static_cast<int>(SwapGroupSplittingType::Split1);
1646              ig++)
1647         {
1648             g = &(s->group[ig]);
1649             gmx_bcast((g->atomset.numAtomsGlobal()) * sizeof((g->xc_old)[0]), g->xc_old, cr->mpi_comm_mygroup);
1650         }
1651     }
1652
1653     /* Make sure that all molecules in the solvent and ion groups contain the
1654      * same number of atoms each */
1655     for (int ig = static_cast<int>(SwapGroupSplittingType::Solvent); ig < s->ngrp; ig++)
1656     {
1657         real charge;
1658
1659         g      = &(s->group[ig]);
1660         g->apm = get_group_apm_check(ig, s, MASTER(cr) && bVerbose, mtop);
1661
1662         /* Since all molecules of a group are equal, we only need enough space
1663          * to determine properties of a single molecule at at time */
1664         snew(g->m, g->apm); /* For the center of mass */
1665         charge   = 0;       /* To determine the charge imbalance */
1666         int molb = 0;
1667         for (int j = 0; j < g->apm; j++)
1668         {
1669             const t_atom& atom = mtopGetAtomParameters(mtop, g->atomset.globalIndex()[j], &molb);
1670             g->m[j]            = atom.m;
1671             charge += atom.q;
1672         }
1673         /* Total charge of one molecule of this group: */
1674         g->q = charge;
1675     }
1676
1677
1678     /* Need mass-weighted center of split group? */
1679     for (int j = static_cast<int>(SwapGroupSplittingType::Split0);
1680          j <= static_cast<int>(SwapGroupSplittingType::Split1);
1681          j++)
1682     {
1683         g = &(s->group[j]);
1684         if (sc->massw_split[j])
1685         {
1686             /* Save the split group masses if mass-weighting is requested */
1687             snew(g->m, g->atomset.numAtomsGlobal());
1688             int molb = 0;
1689             for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1690             {
1691                 g->m[i] = mtopGetAtomMass(mtop, g->atomset.globalIndex()[i], &molb);
1692             }
1693         }
1694     }
1695
1696     /* Make a t_pbc struct on all nodes so that the molecules
1697      * chosen for an exchange can be made whole. */
1698     snew(s->pbc, 1);
1699
1700     bool restartWithAppending = (startingBehavior == gmx::StartingBehavior::RestartWithAppending);
1701     if (MASTER(cr))
1702     {
1703         if (bVerbose)
1704         {
1705             fprintf(stderr, "%s Opening output file %s%s\n", SwS, fn, restartWithAppending ? " for appending" : "");
1706         }
1707
1708         s->fpout = gmx_fio_fopen(fn, restartWithAppending ? "a" : "w");
1709
1710         if (!restartWithAppending)
1711         {
1712             xvgr_header(s->fpout, "Molecule counts", "Time (ps)", "counts", exvggtXNY, oenv);
1713
1714             for (int ig = 0; ig < s->ngrp; ig++)
1715             {
1716                 auto enumValue = static_cast<SwapGroupSplittingType>(ig);
1717                 g              = &(s->group[ig]);
1718                 fprintf(s->fpout,
1719                         "# %s group '%s' contains %d atom%s",
1720                         ig < static_cast<int>(SwapGroupSplittingType::Count) ? enumValueToString(enumValue)
1721                                                                              : "Ion",
1722                         g->molname,
1723                         static_cast<int>(g->atomset.numAtomsGlobal()),
1724                         (g->atomset.numAtomsGlobal() > 1) ? "s" : "");
1725                 if (!(SwapGroupSplittingType::Split0 == enumValue
1726                       || SwapGroupSplittingType::Split1 == enumValue))
1727                 {
1728                     fprintf(s->fpout,
1729                             " with %d atom%s in each molecule of charge %g",
1730                             g->apm,
1731                             (g->apm > 1) ? "s" : "",
1732                             g->q);
1733                 }
1734                 fprintf(s->fpout, ".\n");
1735             }
1736
1737             fprintf(s->fpout, "#\n# Initial positions of split groups:\n");
1738         }
1739
1740         for (int j = static_cast<int>(SwapGroupSplittingType::Split0);
1741              j <= static_cast<int>(SwapGroupSplittingType::Split1);
1742              j++)
1743         {
1744             auto enumValue = static_cast<SwapGroupSplittingType>(j);
1745             g              = &(s->group[j]);
1746             for (size_t i = 0; i < g->atomset.numAtomsGlobal(); i++)
1747             {
1748                 copy_rvec(globalState->x[sc->grp[j].ind[i]], g->xc[i]);
1749             }
1750             /* xc has the correct PBC representation for the two channels, so we do
1751              * not need to correct for that */
1752             get_center(g->xc, g->m, g->atomset.numAtomsGlobal(), g->center);
1753             if (!restartWithAppending)
1754             {
1755                 fprintf(s->fpout,
1756                         "# %s group %s-center %5f nm\n",
1757                         enumValueToString(enumValue),
1758                         DimStr[s->swapdim],
1759                         g->center[s->swapdim]);
1760             }
1761         }
1762
1763         if (!restartWithAppending)
1764         {
1765             if ((0 != sc->bulkOffset[eCompA]) || (0 != sc->bulkOffset[eCompB]))
1766             {
1767                 fprintf(s->fpout, "#\n");
1768                 fprintf(s->fpout,
1769                         "# You provided an offset for the position of the bulk layer(s).\n");
1770                 fprintf(s->fpout,
1771                         "# That means the layers to/from which ions and water molecules are "
1772                         "swapped\n");
1773                 fprintf(s->fpout,
1774                         "# are not midway (= at 0.0) between the compartment-defining layers (at "
1775                         "+/- 1.0).\n");
1776                 fprintf(s->fpout, "# bulk-offsetA = %g\n", sc->bulkOffset[eCompA]);
1777                 fprintf(s->fpout, "# bulk-offsetB = %g\n", sc->bulkOffset[eCompB]);
1778             }
1779
1780             fprintf(s->fpout, "#\n");
1781             fprintf(s->fpout,
1782                     "# Split0 cylinder radius %f nm, up %f nm, down %f nm\n",
1783                     sc->cyl0r,
1784                     sc->cyl0u,
1785                     sc->cyl0l);
1786             fprintf(s->fpout,
1787                     "# Split1 cylinder radius %f nm, up %f nm, down %f nm\n",
1788                     sc->cyl1r,
1789                     sc->cyl1u,
1790                     sc->cyl1l);
1791
1792             fprintf(s->fpout, "#\n");
1793             if (!mdrunOptions.rerun)
1794             {
1795                 fprintf(s->fpout,
1796                         "# Coupling constant (number of swap attempt steps to average over): %d  "
1797                         "(translates to %f ps).\n",
1798                         sc->nAverage,
1799                         sc->nAverage * sc->nstswap * ir->delta_t);
1800                 fprintf(s->fpout, "# Threshold is %f\n", sc->threshold);
1801                 fprintf(s->fpout, "#\n");
1802                 fprintf(s->fpout,
1803                         "# Remarks about which atoms passed which channel use global atoms numbers "
1804                         "starting at one.\n");
1805             }
1806         }
1807     }
1808     else
1809     {
1810         s->fpout = nullptr;
1811     }
1812
1813     /* Allocate memory to remember the past particle counts for time averaging */
1814     for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1815     {
1816         g = &(s->group[ig]);
1817         for (int ic = 0; ic < eCompNR; ic++)
1818         {
1819             snew(g->comp[ic].nMolPast, sc->nAverage);
1820         }
1821     }
1822
1823     /* Get the initial particle concentrations and let the other nodes know */
1824     if (MASTER(cr))
1825     {
1826         if (startingBehavior != gmx::StartingBehavior::NewSimulation)
1827         {
1828             get_initial_ioncounts_from_cpt(ir, s, swapstate, cr, bVerbose);
1829         }
1830         else
1831         {
1832             fprintf(stderr, "%s Determining initial numbers of ions per compartment.\n", SwS);
1833             get_initial_ioncounts(
1834                     ir, s, globalState->x.rvec_array(), globalState->box, cr, mdrunOptions.rerun);
1835         }
1836
1837         /* Prepare (further) checkpoint writes ... */
1838         if (startingBehavior != gmx::StartingBehavior::NewSimulation)
1839         {
1840             /* Consistency check */
1841             if (swapstate->nAverage != sc->nAverage)
1842             {
1843                 gmx_fatal(FARGS,
1844                           "%s Ion count averaging steps mismatch! checkpoint: %d, tpr: %d",
1845                           SwS,
1846                           swapstate->nAverage,
1847                           sc->nAverage);
1848             }
1849         }
1850         else
1851         {
1852             swapstate->nAverage = sc->nAverage;
1853         }
1854         fprintf(stderr, "%s Setting pointers for checkpoint writing\n", SwS);
1855         for (int ic = 0; ic < eCompNR; ic++)
1856         {
1857             for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
1858             {
1859                 g  = &s->group[ig];
1860                 gs = &swapstate->ionType[ig - static_cast<int>(SwapGroupSplittingType::Count)];
1861
1862                 gs->nMolReq_p[ic]    = &(g->comp[ic].nMolReq);
1863                 gs->nMolPast_p[ic]   = &(g->comp[ic].nMolPast[0]);
1864                 gs->inflow_net_p[ic] = &(g->comp[ic].inflow_net);
1865             }
1866         }
1867
1868         /* Determine the total charge imbalance */
1869         s->deltaQ = getRequestedChargeImbalance(s);
1870
1871         if (bVerbose)
1872         {
1873             fprintf(stderr, "%s Requested charge imbalance is Q(A) - Q(B) = %g e.\n", SwS, s->deltaQ);
1874         }
1875         if (!restartWithAppending)
1876         {
1877             fprintf(s->fpout, "# Requested charge imbalance is Q(A)-Q(B) = %g e.\n", s->deltaQ);
1878         }
1879     }
1880
1881     if (PAR(cr))
1882     {
1883         bc_initial_concentrations(cr, ir->swap, s);
1884     }
1885
1886     /* Update the time-averaged number of molecules for all groups and compartments */
1887     for (int ig = static_cast<int>(SwapGroupSplittingType::Count); ig < sc->ngrp; ig++)
1888     {
1889         g = &s->group[ig];
1890         for (int ic = 0; ic < eCompNR; ic++)
1891         {
1892             update_time_window(&g->comp[ic], sc->nAverage, -1);
1893         }
1894     }
1895
1896     /* Initialize arrays that keep track of through which channel the ions go */
1897     detect_flux_per_channel_init(s, swapstate, startingBehavior != gmx::StartingBehavior::NewSimulation);
1898
1899     /* We need to print the legend if we open this file for the first time. */
1900     if (MASTER(cr) && !restartWithAppending)
1901     {
1902         print_ionlist_legend(ir, s, oenv);
1903     }
1904     return s;
1905 }
1906
1907
1908 void finish_swapcoords(t_swap* s)
1909 {
1910     if (s == nullptr)
1911     {
1912         return;
1913     }
1914     if (s->fpout)
1915     {
1916         // Close the swap output file
1917         gmx_fio_fclose(s->fpout);
1918     }
1919 }
1920
1921 /*! \brief Do we need to swap a molecule in any of the ion groups with a water molecule at this step?
1922  *
1923  * From the requested and average molecule counts we determine whether a swap is needed
1924  * at this time step.
1925  */
1926 static gmx_bool need_swap(const t_swapcoords* sc, t_swap* s)
1927 {
1928     int        ic, ig;
1929     t_swapgrp* g;
1930
1931     for (ig = static_cast<int>(SwapGroupSplittingType::Count); ig < sc->ngrp; ig++)
1932     {
1933         g = &s->group[ig];
1934
1935         for (ic = 0; ic < eCompNR; ic++)
1936         {
1937             if (g->comp[ic].nMolReq - g->comp[ic].nMolAv >= sc->threshold)
1938             {
1939                 return TRUE;
1940             }
1941         }
1942     }
1943     return FALSE;
1944 }
1945
1946
1947 /*! \brief Return the index of an atom or molecule suitable for swapping.
1948  *
1949  * Returns the index of an atom that is far off the compartment boundaries,
1950  * that is near to the bulk layer to/from which the swaps take place.
1951  * Other atoms of the molecule (if any) will directly follow the returned index.
1952  *
1953  * \param[in] comp    Structure containing compartment-specific data.
1954  * \param[in] molname Name of the molecule.
1955  *
1956  * \returns Index of the first atom of the molecule chosen for a position exchange.
1957  */
1958 static int get_index_of_distant_atom(t_compartment* comp, const char molname[])
1959 {
1960     int  ibest = -1;
1961     real d     = GMX_REAL_MAX;
1962
1963
1964     /* comp->nat contains the original number of atoms in this compartment
1965      * prior to doing any swaps. Some of these atoms may already have been
1966      * swapped out, but then they are marked with a distance of GMX_REAL_MAX
1967      */
1968     for (int iMol = 0; iMol < comp->nMolBefore; iMol++)
1969     {
1970         if (comp->dist[iMol] < d)
1971         {
1972             ibest = iMol;
1973             d     = comp->dist[ibest];
1974         }
1975     }
1976
1977     if (ibest < 0)
1978     {
1979         gmx_fatal(FARGS,
1980                   "Could not get index of %s atom. Compartment contains %d %s molecules before "
1981                   "swaps.",
1982                   molname,
1983                   comp->nMolBefore,
1984                   molname);
1985     }
1986
1987     /* Set the distance of this index to infinity such that it won't get selected again in
1988      * this time step
1989      */
1990     comp->dist[ibest] = GMX_REAL_MAX;
1991
1992     return comp->ind[ibest];
1993 }
1994
1995
1996 /*! \brief Swaps centers of mass and makes molecule whole if broken */
1997 static void translate_positions(rvec* x, int apm, rvec old_com, rvec new_com, t_pbc* pbc)
1998 {
1999     int  i;
2000     rvec reference, dx, correctPBCimage;
2001
2002
2003     /* Use the first atom as the reference for PBC */
2004     copy_rvec(x[0], reference);
2005
2006     for (i = 0; i < apm; i++)
2007     {
2008         /* PBC distance between position and reference */
2009         pbc_dx(pbc, x[i], reference, dx);
2010
2011         /* Add PBC distance to reference */
2012         rvec_add(reference, dx, correctPBCimage);
2013
2014         /* Subtract old_com from correct image and add new_com */
2015         rvec_dec(correctPBCimage, old_com);
2016         rvec_inc(correctPBCimage, new_com);
2017
2018         copy_rvec(correctPBCimage, x[i]);
2019     }
2020 }
2021
2022
2023 /*! \brief Write back the modified local positions from the collective array to the official positions. */
2024 static void apply_modified_positions(swap_group* g, rvec x[])
2025 {
2026     auto collectiveIndex = g->atomset.collectiveIndex().begin();
2027     for (const auto localIndex : g->atomset.localIndex())
2028     {
2029         /* Copy the possibly modified position */
2030         copy_rvec(g->xc[*collectiveIndex], x[localIndex]);
2031         ++collectiveIndex;
2032     }
2033 }
2034
2035
2036 gmx_bool do_swapcoords(t_commrec*        cr,
2037                        int64_t           step,
2038                        double            t,
2039                        const t_inputrec* ir,
2040                        t_swap*           s,
2041                        gmx_wallcycle*    wcycle,
2042                        rvec              x[],
2043                        matrix            box,
2044                        gmx_bool          bVerbose,
2045                        gmx_bool          bRerun)
2046 {
2047     const t_swapcoords* sc = ir->swap;
2048     int                 j, ic, ig, nswaps;
2049     int                 thisC, otherC; /* Index into this compartment and the other one */
2050     gmx_bool            bSwap = FALSE;
2051     t_swapgrp *         g, *gsol;
2052     int                 isol, iion;
2053     rvec                com_solvent, com_particle; /* solvent and swap molecule's center of mass */
2054
2055
2056     wallcycle_start(wcycle, ewcSWAP);
2057
2058     set_pbc(s->pbc, ir->pbcType, box);
2059
2060     /* Assemble the positions of the split groups, i.e. the channels.
2061      * Here we also pass a shifts array to communicate_group_positions(), so that it can make
2062      * the molecules whole even in cases where they span more than half of the box in
2063      * any dimension */
2064     for (ig = static_cast<int>(SwapGroupSplittingType::Split0);
2065          ig <= static_cast<int>(SwapGroupSplittingType::Split1);
2066          ig++)
2067     {
2068         g = &(s->group[ig]);
2069         communicate_group_positions(cr,
2070                                     g->xc,
2071                                     g->xc_shifts,
2072                                     g->xc_eshifts,
2073                                     TRUE,
2074                                     x,
2075                                     g->atomset.numAtomsGlobal(),
2076                                     g->atomset.numAtomsLocal(),
2077                                     g->atomset.localIndex().data(),
2078                                     g->atomset.collectiveIndex().data(),
2079                                     g->xc_old,
2080                                     box);
2081
2082         get_center(g->xc, g->m, g->atomset.numAtomsGlobal(), g->center); /* center of split groups == channels */
2083     }
2084
2085     /* Assemble the positions of the ions (ig = 3, 4, ...). These molecules should
2086      * be small and we can always make them whole with a simple distance check.
2087      * Therefore we pass NULL as third argument. */
2088     for (ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
2089     {
2090         g = &(s->group[ig]);
2091         communicate_group_positions(cr,
2092                                     g->xc,
2093                                     nullptr,
2094                                     nullptr,
2095                                     FALSE,
2096                                     x,
2097                                     g->atomset.numAtomsGlobal(),
2098                                     g->atomset.numAtomsLocal(),
2099                                     g->atomset.localIndex().data(),
2100                                     g->atomset.collectiveIndex().data(),
2101                                     nullptr,
2102                                     nullptr);
2103
2104         /* Determine how many ions of this type each compartment contains */
2105         sortMoleculesIntoCompartments(g, cr, sc, s, box, step, s->fpout, bRerun, FALSE);
2106     }
2107
2108     /* Output how many ions are in the compartments */
2109     if (MASTER(cr))
2110     {
2111         print_ionlist(s, t, "");
2112     }
2113
2114     /* If we are doing a rerun, we are finished here, since we cannot perform
2115      * swaps anyway */
2116     if (bRerun)
2117     {
2118         return FALSE;
2119     }
2120
2121     /* Do we have to perform a swap? */
2122     bSwap = need_swap(sc, s);
2123     if (bSwap)
2124     {
2125         /* Since we here know that we have to perform ion/water position exchanges,
2126          * we now assemble the solvent positions */
2127         g = &(s->group[static_cast<int>(SwapGroupSplittingType::Solvent)]);
2128         communicate_group_positions(cr,
2129                                     g->xc,
2130                                     nullptr,
2131                                     nullptr,
2132                                     FALSE,
2133                                     x,
2134                                     g->atomset.numAtomsGlobal(),
2135                                     g->atomset.numAtomsLocal(),
2136                                     g->atomset.localIndex().data(),
2137                                     g->atomset.collectiveIndex().data(),
2138                                     nullptr,
2139                                     nullptr);
2140
2141         /* Determine how many molecules of solvent each compartment contains */
2142         sortMoleculesIntoCompartments(g, cr, sc, s, box, step, s->fpout, bRerun, TRUE);
2143
2144         /* Save number of solvent molecules per compartment prior to any swaps */
2145         g->comp[eCompA].nMolBefore = g->comp[eCompA].nMol;
2146         g->comp[eCompB].nMolBefore = g->comp[eCompB].nMol;
2147
2148         for (ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
2149         {
2150             g = &(s->group[ig]);
2151
2152             for (ic = 0; ic < eCompNR; ic++)
2153             {
2154                 /* Determine in which compartment ions are missing and where they are too many */
2155                 g->vacancy[ic] = g->comp[ic].nMolReq - g->comp[ic].nMolAv;
2156
2157                 /* Save number of ions per compartment prior to swaps */
2158                 g->comp[ic].nMolBefore = g->comp[ic].nMol;
2159             }
2160         }
2161
2162         /* Now actually perform the particle exchanges, one swap group after another */
2163         gsol = &s->group[static_cast<int>(SwapGroupSplittingType::Solvent)];
2164         for (ig = static_cast<int>(SwapGroupSplittingType::Count); ig < s->ngrp; ig++)
2165         {
2166             nswaps = 0;
2167             g      = &s->group[ig];
2168             for (thisC = 0; thisC < eCompNR; thisC++)
2169             {
2170                 /* Index to the other compartment */
2171                 otherC = (thisC + 1) % eCompNR;
2172
2173                 while (g->vacancy[thisC] >= sc->threshold)
2174                 {
2175                     /* Swap in an ion */
2176
2177                     /* Get the xc-index of the first atom of a solvent molecule of this compartment */
2178                     isol = get_index_of_distant_atom(&gsol->comp[thisC], gsol->molname);
2179
2180                     /* Get the xc-index of a particle from the other compartment */
2181                     iion = get_index_of_distant_atom(&g->comp[otherC], g->molname);
2182
2183                     get_molecule_center(&gsol->xc[isol], gsol->apm, gsol->m, com_solvent, s->pbc);
2184                     get_molecule_center(&g->xc[iion], g->apm, g->m, com_particle, s->pbc);
2185
2186                     /* Subtract solvent molecule's center of mass and add swap particle's center of mass */
2187                     translate_positions(&gsol->xc[isol], gsol->apm, com_solvent, com_particle, s->pbc);
2188                     /* Similarly for the swap particle, subtract com_particle and add com_solvent */
2189                     translate_positions(&g->xc[iion], g->apm, com_particle, com_solvent, s->pbc);
2190
2191                     /* Keep track of the changes */
2192                     g->vacancy[thisC]--;
2193                     g->vacancy[otherC]++;
2194                     g->comp[thisC].nMol++;
2195                     g->comp[otherC].nMol--;
2196                     g->comp[thisC].inflow_net++;
2197                     g->comp[otherC].inflow_net--;
2198                     /* Correct the past time window to still get the right averages from now on */
2199                     g->comp[thisC].nMolAv++;
2200                     g->comp[otherC].nMolAv--;
2201                     for (j = 0; j < sc->nAverage; j++)
2202                     {
2203                         g->comp[thisC].nMolPast[j]++;
2204                         g->comp[otherC].nMolPast[j]--;
2205                     }
2206                     /* Clear ion history */
2207                     if (MASTER(cr))
2208                     {
2209                         int iMol               = iion / g->apm;
2210                         g->channel_label[iMol] = eChHistPassedNone;
2211                         g->comp_from[iMol]     = eDomainNotset;
2212                     }
2213                     /* That was the swap */
2214                     nswaps++;
2215                 }
2216             }
2217
2218             if (nswaps && bVerbose)
2219             {
2220                 fprintf(stderr,
2221                         "%s Performed %d swap%s in step %" PRId64 " for iontype %s.\n",
2222                         SwS,
2223                         nswaps,
2224                         nswaps > 1 ? "s" : "",
2225                         step,
2226                         g->molname);
2227             }
2228         }
2229
2230         if (s->fpout != nullptr)
2231         {
2232             print_ionlist(s, t, "  # after swap");
2233         }
2234
2235         /* For the solvent and user-defined swap groups, each rank writes back its
2236          * (possibly modified) local positions to the official position array. */
2237         for (ig = static_cast<int>(SwapGroupSplittingType::Solvent); ig < s->ngrp; ig++)
2238         {
2239             g = &s->group[ig];
2240             apply_modified_positions(g, x);
2241         }
2242
2243     } /* end of if(bSwap) */
2244
2245     wallcycle_stop(wcycle, ewcSWAP);
2246
2247     return bSwap;
2248 }