src/gromacs/swap/swapcoords.cpp

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