src/gromacs/swap/swapcoords.c

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