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