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