Bug Summary

File:gromacs/swap/swapcoords.c
Location:line 1935, column 9
Description:Value stored to 'g' is never read

Annotated Source Code

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