2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015, 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.
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.
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.
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.
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.
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.
37 * Implements functions in swapcoords.h.
39 * \author Carsten Kutzner <ckutzne@gwdg.de>
40 * \ingroup module_swap
44 #include "swapcoords.h"
51 #include "gromacs/fileio/confio.h"
52 #include "gromacs/fileio/gmxfio.h"
53 #include "gromacs/fileio/xvgr.h"
54 #include "gromacs/legacyheaders/copyrite.h"
55 #include "gromacs/legacyheaders/names.h"
56 #include "gromacs/legacyheaders/network.h"
57 #include "gromacs/legacyheaders/types/inputrec.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/mdlib/groupcoord.h"
60 #include "gromacs/mdlib/mdrun.h"
61 #include "gromacs/mdlib/sim_util.h"
62 #include "gromacs/pbcutil/pbc.h"
63 #include "gromacs/timing/wallcycle.h"
64 #include "gromacs/topology/mtop_util.h"
65 #include "gromacs/utility/cstringutil.h"
66 #include "gromacs/utility/fatalerror.h"
67 #include "gromacs/utility/smalloc.h"
68 #include "gromacs/utility/snprintf.h"
70 static const char *SwS
= {"SWAP:"}; /**< For output that comes from the swap module */
71 static const char *SwSEmpty
= {" "}; /**< Placeholder for multi-line output */
72 static const char* IonString
[eIonNR
] = {"anion", "cation" }; /**< Type of ion, used for verbose output */
73 static const char* IonStr
[eIonNR
] = {"-", "+" }; /**< Type of ion, used for short output */
74 static const char* CompStr
[eCompNR
] = {"A", "B" }; /**< Compartment name */
75 static const char *SwapStr
[eSwapTypesNR
+1] = { "", "X-", "Y-", "Z-", NULL
}; /**< Name for the swap types. */
76 static const char *DimStr
[DIM
+1] = { "X", "Y", "Z", NULL
}; /**< Name for the swap dimension. */
78 /* eGrpSplit0 and eGrpSplit1 _must_ be neighbors in this list because
79 * we sometimes loop from eGrpSplit0 to eGrpSplit1 */
81 eGrpIons
, eGrpSplit0
, eGrpSplit1
, eGrpSolvent
, eGrpNr
82 }; /**< Group identifier */
83 static const char* GrpString
[eGrpNr
] = { "ion", "split0", "split1", "solvent" }; /**< Group name */
85 /** Keep track of through which channel the ions have passed */
86 enum eChannelHistory
{
87 eChHistPassedNone
, eChHistPassedCh0
, eChHistPassedCh1
, eChHistNr
89 static const char* ChannelString
[eChHistNr
] = { "none", "channel0", "channel1" }; /**< Name for the channels */
91 /*! \brief Domain identifier.
93 * Keeps track of from which compartment the ions came before passing the
97 eDomainNotset
, eDomainA
, eDomainB
, eDomainNr
99 static const char* DomainString
[eDomainNr
] = { "not_assigned", "Domain_A", "Domain_B" }; /**< Name for the domains */
104 * Structure containing compartment-specific data.
106 typedef struct swap_compartment
108 int nat
; /**< Number of atoms matching the
109 compartment conditions. */
110 int nat_old
; /**< Number of atoms before swapping. */
111 int nat_req
; /**< Requested number of atoms. */
112 real nat_av
; /**< Time-averaged number of atoms matching
113 the compartment conditions. */
114 int *nat_past
; /**< Past ion counts for time-averaging. */
115 int *ind
; /**< Indices to coll array of atoms. */
116 real
*dist
; /**< Distance of atom to bulk layer, which is
117 normally the center layer of the compartment */
118 int nalloc
; /**< Allocation size for ind array. */
119 int inflow_netto
; /**< Net inflow of ions into this compartment. */
124 * This structure contains data needed for each of the groups involved in swapping: ions, water,
127 typedef struct swap_group
129 int nat
; /**< Number of atoms in the group */
130 int apm
; /**< Number of atoms in each molecule */
131 atom_id
*ind
; /**< Global atom indices of the group */
132 atom_id
*ind_loc
; /**< Local atom indices of the group */
133 int nat_loc
; /**< Number of local group atoms */
134 int nalloc_loc
; /**< Allocation size for ind_loc */
135 rvec
*xc
; /**< Collective array of group atom positions */
136 ivec
*xc_shifts
; /**< Current (collective) shifts */
137 ivec
*xc_eshifts
; /**< Extra shifts since last DD step */
138 rvec
*xc_old
; /**< Old (collective) positions */
139 real
*qc
; /**< Collective array of charges */
140 int *c_ind_loc
; /**< Position of local atoms in the
141 collective array, [0..nat_loc] */
142 real
*m
; /**< Masses (can be omitted) */
143 unsigned char *comp_from
; /**< (Collective) Stores from which compartment this
144 atom has come. This way we keep track of through
145 which channel an ion permeates (only used for
147 unsigned char *comp_now
; /**< In which compartment this ion is now */
148 unsigned char *channel_label
; /**< Which channel was passed at last by this ion? */
149 rvec center
; /**< Center of the group; COM if masses are used */
154 * Main (private) data structure for the position swapping protocol.
156 typedef struct t_swap
158 int swapdim
; /**< One of XX, YY, ZZ */
159 t_pbc
*pbc
; /**< Needed to make molecules whole. */
160 FILE *fpout
; /**< Output file. */
161 t_group group
[eGrpNr
]; /**< Ions, solvent or channels? */
162 t_compartment comp
[eCompNR
][eIonNR
]; /**< Data for a specific compartment and ion type. */
163 t_compartment compsol
[eCompNR
]; /**< Solvent compartments. */
164 int fluxfromAtoB
[eChanNR
][eIonNR
]; /**< Net flux per channels and ion type. */
165 int ncyl0ions
; /**< Number of ions residing in channel 0. */
166 int ncyl1ions
; /**< Same for channel 1. */
167 int cyl0and1
; /**< Ions assigned to cyl0 and cyl1. Not good. */
168 int *fluxleak
; /**< Pointer to a single int value holding the
169 flux not going through any of the channels. */
170 real deltaQ
; /**< The charge imbalance between the compartments. */
175 /*! \brief Check whether point is in channel.
177 * A channel is a cylinder defined by a disc
178 * with radius r around its center c. The thickness of the cylinder is
185 * <---------c--------->
191 * \param[in] point The position (xyz) under consideration.
192 * \param[in] center The center of the cylinder.
193 * \param[in] d_up The upper extension of the cylinder.
194 * \param[in] d_down The lower extension.
195 * \param[in] r_cyl2 Cylinder radius squared.
196 * \param[in] pbc Structure with info about periodic boundary conditions.
197 * \param[in] normal The membrane normal direction is typically 3, i.e. z, but can be x or y also.
199 * \returns Whether the point is inside the defined cylindric channel.
201 static gmx_bool
is_in_channel(
211 int plane1
, plane2
; /* Directions tangential to membrane */
214 plane1
= (normal
+ 1) % 3; /* typically 0, i.e. XX */
215 plane2
= (normal
+ 2) % 3; /* typically 1, i.e. YY */
217 /* Get the distance vector dr between the point and the center of the cylinder */
218 pbc_dx(pbc
, point
, center
, dr
); /* This puts center in the origin */
220 /* Check vertical direction */
221 if ( (dr
[normal
] > d_up
) || (dr
[normal
] < -d_down
) )
226 /* Check radial direction */
227 if ( (dr
[plane1
]*dr
[plane1
] + dr
[plane2
]*dr
[plane2
]) > r_cyl2
)
232 /* All check passed, this point is in the cylinder */
237 /*! \brief Prints to swap output file which ions are in which compartment. */
238 static void print_ionlist(
241 const char comment
[])
243 int itype
, icomp
, i
, j
;
247 fprintf(s
->fpout
, "%12.5e", time
);
248 for (icomp
= 0; icomp
< eCompNR
; icomp
++)
250 for (itype
= 0; itype
< eIonNR
; itype
++)
252 comp
= &(s
->comp
[icomp
][itype
]);
253 fprintf(s
->fpout
, "%7d%7.1f%7d", comp
->nat
, comp
->nat_av
-comp
->nat_req
, comp
->inflow_netto
);
256 fprintf(s
->fpout
, "%12.3e%12.3e",
257 s
->group
[eGrpSplit0
].center
[s
->swapdim
],
258 s
->group
[eGrpSplit1
].center
[s
->swapdim
]);
260 for (i
= 0; i
< eChanNR
; i
++)
262 for (j
= 0; j
< eIonNR
; j
++)
264 fprintf(s
->fpout
, "%12d", s
->fluxfromAtoB
[i
][j
]);
268 /* Also print the number of ions that leaked from A to B: */
269 fprintf(s
->fpout
, "%12d", *s
->fluxleak
);
271 fprintf(s
->fpout
, "%s\n", comment
);
275 /*! \brief Get the center of a group of nat atoms.
277 * Since with PBC an atom group might not be whole, use the first atom as the
278 * reference atom and determine the center with respect to this reference.
280 static void get_molecule_center(
288 rvec weightedPBCimage
;
290 rvec reference
, correctPBCimage
, dx
;
293 /* Use the first atom as the reference and put other atoms near that one */
294 /* This does not work for large molecules that span > half of the box! */
295 copy_rvec(x
[0], reference
);
297 /* Calculate either the weighted center or simply the center of geometry */
300 for (i
= 0; i
< nat
; i
++)
302 /* PBC distance between position and reference */
303 pbc_dx(pbc
, x
[i
], reference
, dx
);
305 /* Add PBC distance to reference */
306 rvec_add(reference
, dx
, correctPBCimage
);
308 /* Take weight into account */
318 svmul(wi
, correctPBCimage
, weightedPBCimage
);
321 rvec_inc(center
, weightedPBCimage
);
325 svmul(1.0/wsum
, center
, center
);
330 /*! \brief Return TRUE if position x of ion (or water) is found in the compartment,
331 * i.e. between w1 and w2.
333 * One can define and additional offset "b" if one wants to exchange ions/water
334 * to or from a plane not directly in the middle of w1 and w2. The offset can be
335 * in ]-1.0, ..., +1.0 [.
336 * A bulkOffset of 0.0 means 'no offset', so the swap-layer is directly in the
337 * middle between w1 and w2. Offsets -1.0 < b < 0.0 will yield swaps nearer to w1,
338 * whereas offsets 0.0 < 0 < +1.0 will yield swaps nearer to w2.
342 * ||--------------+-------------|-------------+------------------------||
343 * w1 ? ? ? ? ? ? ? ? ? ? ? w2
344 * ||--------------+-------------|----b--------+------------------------||
349 * \param[in] w1 Position of 'wall' atom 1.
350 * \param[in] w2 Position of 'wall' atom 2.
351 * \param[in] x Position of the ion or the water molecule under consideration.
352 * \param[in] l Length of the box, from || to || in the sketch.
353 * \param[in] bulkOffset Where is the bulk layer "b" to be found between w1 and w2?
354 * \param[out] distance_from_b Distance of x to the bulk layer "b".
356 * \returns TRUE if x is between w1 and w2.
358 * Also computes the distance of x to the compartment center (the layer that is
359 * normally situated in the middle of w1 and w2 that would be considered as having
360 * the bulk concentration of ions).
362 static gmx_bool
compartment_contains_atom(
368 real
*distance_from_b
)
374 /* First set the origin in the middle of w1 and w2 */
381 /* Now choose the PBC image of x that is closest to the origin: */
392 *distance_from_b
= (real
)fabs(x
- bulkOffset
*0.5*width
);
394 /* Return TRUE if we now are in area "????" */
395 if ( (x
>= w1
) && (x
< w2
) )
406 /*! \brief Updates the time-averaged number of ions in a compartment. */
407 static void update_time_window(t_compartment
*comp
, int values
, int replace
)
413 /* Put in the new value */
416 comp
->nat_past
[replace
] = comp
->nat
;
419 /* Compute the new time-average */
421 for (i
= 0; i
< values
; i
++)
423 average
+= comp
->nat_past
[i
];
426 comp
->nat_av
= average
;
430 /*! \brief Add the atom with collective index ci to the atom list in compartment 'comp'.
432 * \param[in] ci Index of this ion in the collective xc array.
433 * \param[inout] comp Compartment to add this atom to.
434 * \param[in] distance Shortest distance of this atom to the bulk layer,
435 * from which ion/water pairs are selected for swapping.
437 static void add_to_list(
447 if (nr
>= comp
->nalloc
)
449 comp
->nalloc
= over_alloc_dd(nr
+1);
450 srenew(comp
->ind
, comp
->nalloc
);
451 srenew(comp
->dist
, comp
->nalloc
);
454 comp
->dist
[nr
] = distance
;
459 /*! \brief Determine the compartment boundaries from the channel centers. */
460 static void get_compartment_boundaries(
464 real
*left
, real
*right
)
467 real leftpos
, rightpos
, leftpos_orig
;
472 gmx_fatal(FARGS
, "No compartment %c.", c
+'A');
475 pos0
= s
->group
[eGrpSplit0
].center
[s
->swapdim
];
476 pos1
= s
->group
[eGrpSplit1
].center
[s
->swapdim
];
489 /* This gets us the other compartment: */
492 leftpos_orig
= leftpos
;
494 rightpos
= leftpos_orig
+ box
[s
->swapdim
][s
->swapdim
];
502 /*! \brief Determine the per-channel ion flux.
504 * To determine the flux through the individual channels, we
505 * remember the compartment and channel history of each ion. An ion can be
506 * either in channel0 or channel1, or in the remaining volume of compartment
510 * +-----------------+
513 * ||||||||||0|||||||| bilayer with channel 0
518 * |||||1||||||||||||| bilayer with channel 1
521 * +-----------------+
525 static void detect_flux_per_channel(
530 unsigned char *comp_now
,
531 unsigned char *comp_from
,
532 unsigned char *channel_label
,
542 gmx_bool in_cyl0
, in_cyl1
;
549 /* Check whether ion is inside any of the channels */
550 in_cyl0
= is_in_channel(ion_pos
, s
->group
[eGrpSplit0
].center
, sc
->cyl0u
, sc
->cyl0l
, cyl0_r2
, s
->pbc
, sd
);
551 in_cyl1
= is_in_channel(ion_pos
, s
->group
[eGrpSplit1
].center
, sc
->cyl1u
, sc
->cyl1l
, cyl1_r2
, s
->pbc
, sd
);
553 if (in_cyl0
&& in_cyl1
)
555 /* Ion appears to be in both channels. Something is severely wrong! */
557 *comp_now
= eDomainNotset
;
558 *comp_from
= eDomainNotset
;
559 *channel_label
= eChHistPassedNone
;
563 /* Ion is in channel 0 now */
564 *channel_label
= eChHistPassedCh0
;
565 *comp_now
= eDomainNotset
;
570 /* Ion is in channel 1 now */
571 *channel_label
= eChHistPassedCh1
;
572 *comp_now
= eDomainNotset
;
577 /* Ion is not in any of the channels, so it must be in domain A or B */
580 *comp_now
= eDomainA
;
584 *comp_now
= eDomainB
;
588 /* Only take action, if ion is now in domain A or B, and was before
589 * in the other domain!
591 if (eDomainNotset
== *comp_from
)
593 /* Maybe we can set the domain now */
594 *comp_from
= *comp_now
; /* Could still be eDomainNotset, though */
596 else if ( (*comp_now
!= eDomainNotset
) /* if in channel */
597 && (*comp_from
!= *comp_now
) )
599 /* Obviously the ion changed its domain.
600 * Count this for the channel through which it has passed. */
601 switch (*channel_label
)
603 case eChHistPassedNone
:
604 *s
->fluxleak
= *s
->fluxleak
+ 1;
606 fprintf(stderr
, " %s Warning! Step %s, ion %d (%s) moved from %s to %s\n",
607 SwS
, gmx_step_str(step
, buf
), iion
, IonStr
[iontype
], DomainString
[*comp_from
], DomainString
[*comp_now
]);
610 fprintf(stderr
, ", possibly due to a swap in the original simulation.\n");
614 fprintf(stderr
, "but did not pass cyl0 or cyl1 as defined in the .mdp file.\n"
615 "Do you have an ion somewhere within the membrane?\n");
616 /* Write this info to the CompEL output file: */
617 fprintf(s
->fpout
, " # Warning: step %s, ion %d (%s) moved from %s to %s (probably through the membrane)\n",
618 gmx_step_str(step
, buf
), iion
, IonStr
[iontype
],
619 DomainString
[*comp_from
], DomainString
[*comp_now
]);
623 case eChHistPassedCh0
:
624 case eChHistPassedCh1
:
625 if (*channel_label
== eChHistPassedCh0
)
634 if (eDomainA
== *comp_from
)
636 s
->fluxfromAtoB
[chan_nr
][iontype
]++;
640 s
->fluxfromAtoB
[chan_nr
][iontype
]--;
642 fprintf(fpout
, "# Atom nr. %d finished passing %s.\n", iion
, ChannelString
[*channel_label
]);
645 gmx_fatal(FARGS
, "%s Unknown channel history entry!\n", SwS
);
649 /* This ion has moved to the _other_ compartment ... */
650 *comp_from
= *comp_now
;
651 /* ... and it did not pass any channel yet */
652 *channel_label
= eChHistPassedNone
;
657 /*! \brief Get the lists of ions for the two compartments */
658 static void compartmentalize_ions(
671 real cyl0_r2
, cyl1_r2
;
673 int sum
, not_in_comp
[eCompNR
]; /* consistency check */
678 iong
= &s
->group
[eGrpIons
];
681 cyl0_r2
= sc
->cyl0r
* sc
->cyl0r
;
682 cyl1_r2
= sc
->cyl1r
* sc
->cyl1r
;
685 /* Get us a counter that cycles in the range of [0 ... sc->nAverage[ */
686 replace
= (step
/sc
->nstswap
) % sc
->nAverage
;
688 for (comp
= eCompA
; comp
<= eCompB
; comp
++)
690 /* Get lists of atoms that match criteria for this compartment */
691 get_compartment_boundaries(comp
, sc
->si_priv
, box
, &left
, &right
);
693 /* First clear the ion lists */
694 s
->comp
[comp
][eIonNEG
].nat
= 0;
695 s
->comp
[comp
][eIonPOS
].nat
= 0;
696 not_in_comp
[comp
] = 0; /* consistency check */
698 /* Loop over the IONS */
699 for (i
= 0; i
< iong
->nat
; i
++)
701 /* Anion or cation? */
702 type
= iong
->qc
[i
] < 0 ? eIonNEG
: eIonPOS
;
704 /* Is this ion in the compartment that we look at? */
705 if (compartment_contains_atom(left
, right
, iong
->xc
[i
][sd
], box
[sd
][sd
], sc
->bulkOffset
[comp
], &dist
) )
707 /* Now put it into the list containing only ions of its type */
708 add_to_list(i
, &s
->comp
[comp
][type
], dist
);
710 /* Master also checks through which channel each ion has passed */
711 if (MASTER(cr
) && (iong
->comp_now
!= NULL
))
713 ion_nr_global
= iong
->ind
[i
] + 1; /* PDB index starts at 1 ... */
714 detect_flux_per_channel(ion_nr_global
, comp
, type
, iong
->xc
[i
],
715 &iong
->comp_now
[i
], &iong
->comp_from
[i
], &iong
->channel_label
[i
],
716 sc
, cyl0_r2
, cyl1_r2
, step
, bRerun
, fpout
);
721 not_in_comp
[comp
] += 1;
724 /* Correct the time-averaged number of ions in both compartments */
725 update_time_window(&s
->comp
[comp
][eIonNEG
], sc
->nAverage
, replace
);
726 update_time_window(&s
->comp
[comp
][eIonPOS
], sc
->nAverage
, replace
);
729 /* Flux detection warnings */
735 "%s Warning: %d atoms were detected as being in both channels! Probably your split\n"
736 "%s cylinder is way too large, or one compartment has collapsed (step %" GMX_PRId64
")\n",
737 SwS
, s
->cyl0and1
, SwS
, step
);
739 fprintf(s
->fpout
, "Warning: %d atoms were assigned to both channels!\n", s
->cyl0and1
);
746 /* Consistency checks */
747 if (not_in_comp
[eCompA
] + not_in_comp
[eCompB
] != iong
->nat
)
751 fprintf(fpout
, "# Warning: Inconsistency during ion compartmentalization. !inA: %d, !inB: %d, total ions %d\n",
752 not_in_comp
[eCompA
], not_in_comp
[eCompB
], iong
->nat
);
757 fprintf(stderr
, "%s rank %d: Inconsistency during ion compartmentalization. !inA: %d, !inB: %d, total ions %d\n",
758 SwS
, cr
->nodeid
, not_in_comp
[eCompA
], not_in_comp
[eCompB
], iong
->nat
);
762 sum
= s
->comp
[eCompA
][eIonNEG
].nat
+ s
->comp
[eCompA
][eIonPOS
].nat
763 + s
->comp
[eCompB
][eIonNEG
].nat
+ s
->comp
[eCompB
][eIonPOS
].nat
;
764 if (sum
!= iong
->nat
)
768 fprintf(fpout
, "# Warning: %d atoms are in the ion group, but altogether %d have been assigned to the compartments.\n",
774 fprintf(stderr
, "%s rank %d: %d atoms are in the ion group, but altogether %d have been assigned to the compartments.\n",
775 SwS
, cr
->nodeid
, iong
->nat
, sum
);
783 /*! \brief Set up the compartments and get lists of solvent atoms in each compartment */
784 static void compartmentalize_solvent(
796 int sum
, not_in_comp
[eCompNR
]; /* consistency check */
800 solg
= &s
->group
[eGrpSolvent
];
804 for (comp
= eCompA
; comp
<= eCompB
; comp
++)
806 /* Get lists of atoms that match criteria for this compartment */
807 get_compartment_boundaries(comp
, sc
->si_priv
, box
, &left
, &right
);
809 /* First clear the solvent molecule lists */
810 s
->compsol
[comp
].nat
= 0;
811 not_in_comp
[comp
] = 0; /* consistency check */
813 /* Loop over the solvent MOLECULES */
814 for (i
= 0; i
< sc
->nat_sol
; i
+= apm
)
816 if (compartment_contains_atom(left
, right
, solg
->xc
[i
][sd
], box
[sd
][sd
], sc
->bulkOffset
[comp
], &dist
))
818 /* Add the whole molecule to the list */
819 for (j
= 0; j
< apm
; j
++)
821 add_to_list(i
+j
, &s
->compsol
[comp
], dist
);
826 not_in_comp
[comp
] += apm
;
833 fprintf(fpout
, "# Solv. molecules in comp.%s: %d comp.%s: %d\n",
834 CompStr
[eCompA
], s
->compsol
[eCompA
].nat
/apm
,
835 CompStr
[eCompB
], s
->compsol
[eCompB
].nat
/apm
);
838 /* Consistency checks */
839 if (not_in_comp
[eCompA
] + not_in_comp
[eCompB
] != solg
->nat
)
843 fprintf(fpout
, "# Warning: Inconsistency during solvent compartmentalization. !inA: %d, !inB: %d, solvent atoms %d\n",
844 not_in_comp
[eCompA
], not_in_comp
[eCompB
], solg
->nat
);
849 fprintf(stderr
, "%s rank %d: Inconsistency during solvent compartmentalization. !inA: %d, !inB: %d, solvent atoms %d\n",
850 SwS
, cr
->nodeid
, not_in_comp
[eCompA
], not_in_comp
[eCompB
], solg
->nat
);
853 sum
= s
->compsol
[eCompA
].nat
+ s
->compsol
[eCompB
].nat
;
854 if (sum
!= solg
->nat
)
858 fprintf(fpout
, "# Warning: %d atoms in solvent group, but %d have been assigned to the compartments.\n",
864 fprintf(stderr
, "%s rank %d: %d atoms in solvent group, but %d have been assigned to the compartments.\n",
865 SwS
, cr
->nodeid
, solg
->nat
, sum
);
871 /*! \brief Find out how many group atoms are in the compartments initially */
872 static void get_initial_ioncounts(
874 rvec x
[], /* the initial positions */
888 /* Copy the initial swap group positions to the collective array so
889 * that we can compartmentalize */
890 for (i
= 0; i
< sc
->nat
; i
++)
893 copy_rvec(x
[ind
], s
->group
[eGrpIons
].xc
[i
]);
896 /* Set up the compartments and get lists of atoms in each compartment */
897 compartmentalize_ions(cr
, sc
, box
, 0, s
->fpout
, bRerun
);
899 /* Set initial concentrations if requested */
900 for (ic
= 0; ic
< eCompNR
; ic
++)
902 s
->comp
[ic
][eIonPOS
].nat_req
= sc
->ncations
[ic
];
903 s
->comp
[ic
][eIonNEG
].nat_req
= sc
->nanions
[ic
];
905 for (ic
= 0; ic
< eCompNR
; ic
++)
907 for (ii
= 0; ii
< eIonNR
; ii
++)
909 if (s
->comp
[ic
][ii
].nat_req
< 0)
911 s
->comp
[ic
][ii
].nat_req
= s
->comp
[ic
][ii
].nat
;
916 /* Check whether the number of requested ions adds up to the total number of ions */
917 for (ii
= 0; ii
< eIonNR
; ii
++)
919 req
[ii
] = s
->comp
[eCompA
][ii
].nat_req
+ s
->comp
[eCompB
][ii
].nat_req
;
920 tot
[ii
] = s
->comp
[eCompA
][ii
].nat
+ s
->comp
[eCompB
][ii
].nat
;
922 if ( (req
[eCompA
] != tot
[eCompA
]) || (req
[eCompB
] != tot
[eCompB
]) )
924 gmx_fatal(FARGS
, "Mismatch of the number of ions summed over both compartments.\n"
925 "You requested a total of %d anions and %d cations,\n"
926 "but there are a total of %d anions and %d cations in the system.\n",
927 req
[eIonNEG
], req
[eIonPOS
],
928 tot
[eIonNEG
], tot
[eIonPOS
]);
931 /* Initialize time-averaging:
932 * Write initial concentrations to all time bins to start with */
933 for (ic
= 0; ic
< eCompNR
; ic
++)
935 for (ii
= 0; ii
< eIonNR
; ii
++)
937 s
->comp
[ic
][ii
].nat_av
= s
->comp
[ic
][ii
].nat
;
938 for (i
= 0; i
< sc
->nAverage
; i
++)
940 s
->comp
[ic
][ii
].nat_past
[i
] = s
->comp
[ic
][ii
].nat
;
947 /*! \brief Copy history of ion counts from checkpoint file.
949 * When called, the checkpoint file has already been read in. Here we copy
950 * over the values from .cpt file to the swap data structure.
952 static void get_initial_ioncounts_from_cpt(
953 t_inputrec
*ir
, swapstate_t
*swapstate
,
954 t_commrec
*cr
, gmx_bool bVerbose
)
965 /* Copy the past values from the checkpoint values that have been read in already */
968 fprintf(stderr
, "%s Copying values from checkpoint\n", SwS
);
971 for (ic
= 0; ic
< eCompNR
; ic
++)
973 for (ii
= 0; ii
< eIonNR
; ii
++)
975 s
->comp
[ic
][ii
].nat_req
= swapstate
->nat_req
[ic
][ii
];
976 s
->comp
[ic
][ii
].inflow_netto
= swapstate
->inflow_netto
[ic
][ii
];
980 fprintf(stderr
, "%s ... Influx netto: %d Requested: %d Past values: ", SwS
,
981 s
->comp
[ic
][ii
].inflow_netto
, s
->comp
[ic
][ii
].nat_req
);
984 for (j
= 0; j
< sc
->nAverage
; j
++)
986 s
->comp
[ic
][ii
].nat_past
[j
] = swapstate
->nat_past
[ic
][ii
][j
];
989 fprintf(stderr
, "%d ", s
->comp
[ic
][ii
].nat_past
[j
]);
994 fprintf(stderr
, "\n");
1002 /*! \brief The master lets all others know about the initial ion counts. */
1003 static void bc_initial_concentrations(
1012 for (ic
= 0; ic
< eCompNR
; ic
++)
1014 for (ii
= 0; ii
< eIonNR
; ii
++)
1016 gmx_bcast(sizeof(s
->comp
[ic
][ii
].nat_req
), &(s
->comp
[ic
][ii
].nat_req
), cr
);
1017 gmx_bcast(sizeof(s
->comp
[ic
][ii
].nat
), &(s
->comp
[ic
][ii
].nat
), cr
);
1018 gmx_bcast( swap
->nAverage
* sizeof(s
->comp
[ic
][ii
].nat_past
[0]), s
->comp
[ic
][ii
].nat_past
, cr
);
1024 /*! \brief Ensure that each atom belongs to at most one of the swap groups. */
1025 static void check_swap_groups(t_swap
*s
, int nat
, gmx_bool bVerbose
)
1029 atom_id
*nGroup
= NULL
; /* This array counts for each atom in the MD system to
1030 how many swap groups it belongs (should be 0 or 1!) */
1032 int nMultiple
= 0; /* Number of atoms belonging to multiple groups */
1037 fprintf(stderr
, "%s Making sure each atom belongs to at most one of the swap groups.\n", SwS
);
1040 /* Add one to the group count of atoms belonging to a swap group: */
1042 for (i
= 0; i
< eGrpNr
; i
++)
1045 for (j
= 0; j
< g
->nat
; j
++)
1047 /* Get the global index of this atom of this group: */
1052 /* Make sure each atom belongs to at most one swap group: */
1053 for (j
= 0; j
< g
->nat
; j
++)
1064 gmx_fatal(FARGS
, "%s Cannot perform swapping since %d atom%s allocated to more than one swap index group.\n"
1065 "%s Each atom must be allocated to at most one of the split groups, the swap group, or the solvent.\n"
1066 "%s Check the .mdp file settings regarding the swap index groups or the index groups themselves.\n",
1067 SwS
, nMultiple
, (1 == nMultiple
) ? " is" : "s are", SwSEmpty
, SwSEmpty
);
1072 /*! \brief Get the number of atoms per molecule for this group.
1074 * Also ensure that all the molecules in this group have this number of atoms.
1076 static int get_group_apm_check(
1080 const gmx_mtop_atomlookup_t alook
,
1085 int molb
, molnr
, atnr_mol
;
1088 ind
= s
->group
[group
].ind
;
1089 nat
= s
->group
[group
].nat
;
1091 /* Determine the number of solvent atoms per solvent molecule from the
1092 * first solvent atom: */
1094 gmx_mtop_atomnr_to_molblock_ind(alook
, ind
[i
], &molb
, &molnr
, &atnr_mol
);
1095 apm
= mtop
->molblock
[molb
].natoms_mol
;
1099 fprintf(stderr
, "%s Checking whether all %s molecules consist of %d atom%s\n",
1100 SwS
, GrpString
[group
], apm
, apm
> 1 ? "s" : "");
1103 /* Check whether this is also true for all other solvent atoms */
1104 for (i
= 1; i
< nat
; i
++)
1106 gmx_mtop_atomnr_to_molblock_ind(alook
, ind
[i
], &molb
, &molnr
, &atnr_mol
);
1107 if (apm
!= mtop
->molblock
[molb
].natoms_mol
)
1109 gmx_fatal(FARGS
, "Not all %s group molecules consist of %d atoms.",
1110 GrpString
[group
], apm
);
1118 /*! \brief Print the legend to the swap output file.
1120 * Also print the initial ion counts
1122 static void print_ionlist_legend(t_inputrec
*ir
,
1123 const gmx_output_env_t
*oenv
)
1125 const char **legend
;
1131 s
= ir
->swap
->si_priv
;
1133 snew(legend
, eCompNR
*eIonNR
*3 + 2 + eChanNR
*eIonNR
+ 1);
1134 for (ic
= count
= 0; ic
< eCompNR
; ic
++)
1136 for (ii
= 0; ii
< eIonNR
; ii
++)
1138 snprintf(buf
, STRLEN
, "%s %ss", CompStr
[ic
], IonString
[ii
]);
1139 legend
[count
++] = gmx_strdup(buf
);
1140 snprintf(buf
, STRLEN
, "%s av. mismatch to %d%s",
1141 CompStr
[ic
], s
->comp
[ic
][ii
].nat_req
, IonStr
[ii
]);
1142 legend
[count
++] = gmx_strdup(buf
);
1143 snprintf(buf
, STRLEN
, "%s netto %s influx", CompStr
[ic
], IonString
[ii
]);
1144 legend
[count
++] = gmx_strdup(buf
);
1147 snprintf(buf
, STRLEN
, "%scenter of %s of split group 0", SwapStr
[ir
->eSwapCoords
], (NULL
!= s
->group
[eGrpSplit0
].m
) ? "mass" : "geometry");
1148 legend
[count
++] = gmx_strdup(buf
);
1149 snprintf(buf
, STRLEN
, "%scenter of %s of split group 1", SwapStr
[ir
->eSwapCoords
], (NULL
!= s
->group
[eGrpSplit1
].m
) ? "mass" : "geometry");
1150 legend
[count
++] = gmx_strdup(buf
);
1152 for (ic
= 0; ic
< eChanNR
; ic
++)
1154 for (ii
= 0; ii
< eIonNR
; ii
++)
1156 snprintf(buf
, STRLEN
, "A->ch%d->B %s permeations", ic
, IonString
[ii
]);
1157 legend
[count
++] = gmx_strdup(buf
);
1161 snprintf(buf
, STRLEN
, "leakage");
1162 legend
[count
++] = gmx_strdup(buf
);
1164 xvgr_legend(s
->fpout
, count
, legend
, oenv
);
1166 fprintf(s
->fpout
, "# Instantaneous ion counts and time-averaged differences to requested numbers\n");
1167 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 ");
1168 fprintf(s
->fpout
, " %s-Split0 %s-Split1", DimStr
[s
->swapdim
], DimStr
[s
->swapdim
]);
1169 fprintf(s
->fpout
, " A-ch0-B_an A-ch0-B_cat A-ch1-B_an A-ch1-B_cat ion_leakage\n");
1174 /*! \brief Initialize channel ion flux detection routine.
1176 * Initialize arrays that keep track of where the ions come from and where
1179 static void detect_flux_per_channel_init(
1182 swapstate_t
*swapstate
,
1183 gmx_bool bStartFromCpt
)
1189 g
= &(s
->group
[eGrpIons
]);
1191 /* All these flux detection routines run on the master only */
1195 g
->comp_from
= NULL
;
1196 g
->channel_label
= NULL
;
1201 /******************************************************/
1202 /* Channel and domain history for the individual ions */
1203 /******************************************************/
1204 if (bStartFromCpt
) /* set the pointers right */
1206 g
->comp_from
= swapstate
->comp_from
;
1207 g
->channel_label
= swapstate
->channel_label
;
1209 else /* allocate memory */
1211 snew(g
->comp_from
, g
->nat
);
1212 swapstate
->comp_from
= g
->comp_from
;
1213 snew(g
->channel_label
, g
->nat
);
1214 swapstate
->channel_label
= g
->channel_label
;
1216 snew(g
->comp_now
, g
->nat
);
1218 /* Initialize the channel and domain history counters */
1219 for (i
= 0; i
< g
->nat
; i
++)
1221 g
->comp_now
[i
] = eDomainNotset
;
1224 g
->comp_from
[i
] = eDomainNotset
;
1225 g
->channel_label
[i
] = eChHistPassedNone
;
1229 /************************************/
1230 /* Channel fluxes for both channels */
1231 /************************************/
1238 fprintf(stderr
, "%s Copying channel fluxes from checkpoint file data\n", SwS
);
1241 for (ic
= 0; ic
< eChanNR
; ic
++)
1243 fprintf(stderr
, "%s Channel %d flux history: ", SwS
, ic
);
1244 for (ii
= 0; ii
< eIonNR
; ii
++)
1248 s
->fluxfromAtoB
[ic
][ii
] = swapstate
->fluxfromAtoB
[ic
][ii
];
1252 s
->fluxfromAtoB
[ic
][ii
] = 0;
1255 fprintf(stderr
, "%d %s%s ", s
->fluxfromAtoB
[ic
][ii
], IonString
[ii
], s
->fluxfromAtoB
[ic
][ii
] == 1 ? "" : "s");
1257 fprintf(stderr
, "\n");
1261 s
->fluxleak
= swapstate
->fluxleak
;
1265 snew(s
->fluxleak
, 1);
1267 /* Set pointer for checkpoint writing */
1268 swapstate
->fluxleak
= s
->fluxleak
;
1271 /* Set pointers for checkpoint writing */
1272 for (ic
= 0; ic
< eChanNR
; ic
++)
1274 for (ii
= 0; ii
< eIonNR
; ii
++)
1276 swapstate
->fluxfromAtoB_p
[ic
][ii
] = &(s
->fluxfromAtoB
[ic
][ii
]);
1282 /*! \brief Outputs the initial structure to PDB file for debugging reasons.
1284 * Output the starting structure so that in case of multimeric channels
1285 * the user can check whether we have the correct PBC image for all atoms.
1286 * If this is not correct, the ion counts per channel will be very likely
1289 static void outputStartStructureIfWanted(gmx_mtop_t
*mtop
, rvec
*x
, int ePBC
, matrix box
)
1291 char *env
= getenv("GMX_COMPELDUMP");
1295 fprintf(stderr
, "\n%s Found env.var. GMX_COMPELDUMP, will output CompEL starting structure made whole.\n"
1296 "%s In case of multimeric channels, please check whether they have the correct PBC representation.\n",
1299 write_sto_conf_mtop("CompELAssumedWholeConfiguration.pdb", *mtop
->name
, mtop
, x
, NULL
, ePBC
, box
);
1304 /*! \brief Initialize the swapstate structure, used for checkpoint writing.
1306 * The swapstate struct stores the information we need to make the channels
1307 * whole again after restarts from a checkpoint file. Here we do the following:\n
1308 * a) If we did not start from .cpt, we prepare the struct for proper .cpt writing,\n
1309 * b) if we did start from .cpt, we copy over the last whole structures from .cpt,\n
1310 * c) in any case, for subsequent checkpoint writing, we set the pointers in\n
1311 * swapstate to the x_old arrays, which contain the correct PBC representation of
1312 * multimeric channels at the last time step.
1314 static void init_swapstate(
1315 swapstate_t
*swapstate
,
1318 rvec x
[], /* the initial positions */
1323 rvec
*x_pbc
= NULL
; /* positions of the whole MD system with molecules made whole */
1330 /* We always need the last whole positions such that
1331 * in the next time step we can make the channels whole again in PBC */
1332 if (swapstate
->bFromCpt
)
1334 /* Copy the last whole positions of each channel from .cpt */
1335 g
= &(s
->group
[eGrpSplit0
]);
1336 for (i
= 0; i
< g
->nat
; i
++)
1338 copy_rvec(swapstate
->xc_old_whole
[eChan0
][i
], g
->xc_old
[i
]);
1340 g
= &(s
->group
[eGrpSplit1
]);
1341 for (i
= 0; i
< g
->nat
; i
++)
1343 copy_rvec(swapstate
->xc_old_whole
[eChan1
][i
], g
->xc_old
[i
]);
1348 /* Extract the initial split group positions. */
1350 /* Remove pbc, make molecule whole. */
1351 snew(x_pbc
, mtop
->natoms
);
1352 m_rveccopy(mtop
->natoms
, x
, x_pbc
);
1354 /* This can only make individual molecules whole, not multimers */
1355 do_pbc_mtop(NULL
, ePBC
, box
, mtop
, x_pbc
);
1357 /* Output the starting structure? */
1358 outputStartStructureIfWanted(mtop
, x_pbc
, ePBC
, box
);
1360 /* If this is the first run (i.e. no checkpoint present) we assume
1361 * that the starting positions give us the correct PBC representation */
1362 for (ig
= eGrpSplit0
; ig
<= eGrpSplit1
; ig
++)
1364 g
= &(s
->group
[ig
]);
1365 for (i
= 0; i
< g
->nat
; i
++)
1367 copy_rvec(x_pbc
[g
->ind
[i
]], g
->xc_old
[i
]);
1372 /* Prepare swapstate arrays for later checkpoint writing */
1373 swapstate
->nat
[eChan0
] = s
->group
[eGrpSplit0
].nat
;
1374 swapstate
->nat
[eChan1
] = s
->group
[eGrpSplit1
].nat
;
1377 /* For subsequent checkpoint writing, set the swapstate pointers to the xc_old
1378 * arrays that get updated at every swapping step */
1379 swapstate
->xc_old_whole_p
[eChan0
] = &s
->group
[eGrpSplit0
].xc_old
;
1380 swapstate
->xc_old_whole_p
[eChan1
] = &s
->group
[eGrpSplit1
].xc_old
;
1384 void init_swapcoords(FILE *fplog
,
1391 swapstate_t
*swapstate
,
1393 const gmx_output_env_t
*oenv
,
1394 unsigned long Flags
)
1396 int i
, ic
, ig
, ii
, j
;
1401 gmx_bool bAppend
, bStartFromCpt
, bRerun
;
1402 gmx_mtop_atomlookup_t alook
= NULL
;
1406 alook
= gmx_mtop_atomlookup_init(mtop
);
1408 if ( (PAR(cr
)) && !DOMAINDECOMP(cr
) )
1410 gmx_fatal(FARGS
, "Position swapping is only implemented for domain decomposition!");
1413 bAppend
= Flags
& MD_APPENDFILES
;
1414 bStartFromCpt
= Flags
& MD_STARTFROMCPT
;
1415 bRerun
= Flags
& MD_RERUN
;
1418 snew(sc
->si_priv
, 1);
1425 gmx_fatal(FARGS
, "%s This module does not support reruns in parallel\nPlease request a serial run with -nt 1 / -np 1\n", SwS
);
1428 fprintf(stderr
, "%s Rerun - using every available frame\n", SwS
);
1430 sc
->nAverage
= 1; /* averaging makes no sense for reruns */
1433 if (MASTER(cr
) && !bAppend
)
1435 fprintf(fplog
, "\nInitializing ion/water position exchanges\n");
1436 please_cite(fplog
, "Kutzner2011b");
1439 switch (ir
->eSwapCoords
)
1455 /* Copy some data to the group structures for convenience */
1456 /* Number of atoms in the group */
1457 s
->group
[eGrpIons
].nat
= sc
->nat
;
1458 s
->group
[eGrpSplit0
].nat
= sc
->nat_split
[0];
1459 s
->group
[eGrpSplit1
].nat
= sc
->nat_split
[1];
1460 s
->group
[eGrpSolvent
].nat
= sc
->nat_sol
;
1461 /* Pointer to the indices */
1462 s
->group
[eGrpIons
].ind
= sc
->ind
;
1463 s
->group
[eGrpSplit0
].ind
= sc
->ind_split
[0];
1464 s
->group
[eGrpSplit1
].ind
= sc
->ind_split
[1];
1465 s
->group
[eGrpSolvent
].ind
= sc
->ind_sol
;
1467 check_swap_groups(s
, mtop
->natoms
, bVerbose
&& MASTER(cr
));
1469 /* Allocate space for the collective arrays for all groups */
1470 for (ig
= 0; ig
< eGrpNr
; ig
++)
1472 g
= &(s
->group
[ig
]);
1473 snew(g
->xc
, g
->nat
);
1474 snew(g
->c_ind_loc
, g
->nat
);
1475 /* For the split groups (the channels) we need some extra memory to
1476 * be able to make the molecules whole even if they span more than
1477 * half of the box size. */
1478 if (eGrpSplit0
== ig
|| eGrpSplit1
== ig
)
1480 snew(g
->xc_shifts
, g
->nat
);
1481 snew(g
->xc_eshifts
, g
->nat
);
1482 snew(g
->xc_old
, g
->nat
);
1488 init_swapstate(swapstate
, sc
, mtop
, x
, box
, ir
->ePBC
);
1491 /* After init_swapstate we have a set of (old) whole positions for our
1492 * channels. Now transfer that to all nodes */
1495 for (ig
= eGrpSplit0
; ig
<= eGrpSplit1
; ig
++)
1497 g
= &(s
->group
[ig
]);
1498 gmx_bcast((g
->nat
)*sizeof((g
->xc_old
)[0]), g
->xc_old
, (cr
));
1502 /* Make sure that all molecules in the ion and solvent groups contain the
1503 * same number of atoms each */
1504 s
->group
[eGrpIons
].apm
= get_group_apm_check(eGrpIons
, s
, MASTER(cr
) && bVerbose
, alook
, mtop
);
1505 s
->group
[eGrpSolvent
].apm
= get_group_apm_check(eGrpSolvent
, s
, MASTER(cr
) && bVerbose
, alook
, mtop
);
1507 /* Save masses where needed */
1508 s
->group
[eGrpIons
].m
= NULL
;
1509 /* We only need enough space to determine a single solvent molecule's
1510 * center at at time */
1511 g
= &(s
->group
[eGrpSolvent
]);
1514 /* Need mass-weighted center of split group? */
1515 for (j
= 0, ig
= eGrpSplit0
; j
< eChanNR
; ig
++, j
++)
1517 g
= &(s
->group
[ig
]);
1518 if (TRUE
== sc
->massw_split
[j
])
1520 /* Save the split group charges if mass-weighting is requested */
1522 for (i
= 0; i
< g
->nat
; i
++)
1524 gmx_mtop_atomnr_to_atom(alook
, g
->ind
[i
], &atom
);
1534 /* Save the ionic charges */
1535 g
= &(s
->group
[eGrpIons
]);
1536 snew(g
->qc
, g
->nat
);
1537 for (i
= 0; i
< g
->nat
; i
++)
1539 gmx_mtop_atomnr_to_atom(alook
, g
->ind
[i
], &atom
);
1543 /* Make a t_pbc struct on all nodes so that the molecules
1544 * chosen for an exchange can be made whole. */
1546 /* Every node needs to call set_pbc() and therefore every node needs
1547 * to know the box dimensions */
1548 copy_mat(box
, boxCopy
);
1551 gmx_bcast(sizeof(boxCopy
), boxCopy
, cr
);
1553 set_pbc(s
->pbc
, -1, boxCopy
);
1559 fprintf(stderr
, "%s Opening output file %s%s\n", SwS
, fn
, bAppend
? " for appending" : "");
1562 s
->fpout
= gmx_fio_fopen(fn
, bAppend
? "a" : "w" );
1566 xvgr_header(s
->fpout
, "Ion counts", "Time (ps)", "counts", exvggtXNY
, oenv
);
1568 for (ig
= 0; ig
< eGrpNr
; ig
++)
1570 g
= &(s
->group
[ig
]);
1571 fprintf(s
->fpout
, "# %s group contains %d atom%s", GrpString
[ig
], g
->nat
, (g
->nat
> 1) ? "s" : "");
1572 if (eGrpSolvent
== ig
|| eGrpIons
== ig
)
1574 fprintf(s
->fpout
, " with %d atom%s in each molecule", g
->apm
, (g
->apm
> 1) ? "s" : "");
1576 fprintf(s
->fpout
, ".\n");
1579 fprintf(s
->fpout
, "#\n# Initial positions of split groups:\n");
1582 for (j
= 0, ig
= eGrpSplit0
; j
< eChanNR
; j
++, ig
++)
1584 g
= &(s
->group
[ig
]);
1585 for (i
= 0; i
< g
->nat
; i
++)
1587 copy_rvec(x
[sc
->ind_split
[j
][i
]], g
->xc
[i
]);
1589 if (eGrpSplit0
== ig
|| eGrpSplit1
== ig
)
1591 /* xc has the correct PBC representation for the two channels, so we do
1592 * not need to correct for that */
1593 get_center(g
->xc
, g
->m
, g
->nat
, g
->center
);
1597 /* For the water molecules, we need to make the molecules whole */
1598 get_molecule_center(g
->xc
, g
->nat
, g
->m
, g
->center
, s
->pbc
);
1602 fprintf(s
->fpout
, "# %s group %s-center %5f nm\n", GrpString
[ig
],
1603 DimStr
[s
->swapdim
], g
->center
[s
->swapdim
]);
1609 if ( (0 != sc
->bulkOffset
[eCompA
]) || (0 != sc
->bulkOffset
[eCompB
]) )
1611 fprintf(s
->fpout
, "#\n");
1612 fprintf(s
->fpout
, "# You provided an offset for the position of the bulk layer(s).\n");
1613 fprintf(s
->fpout
, "# That means the layers to/from which ions and water molecules are swapped\n");
1614 fprintf(s
->fpout
, "# are not midway (= at 0.0) between the compartment-defining layers (at +/- 1.0).\n");
1615 fprintf(s
->fpout
, "# bulk-offsetA = %g\n", sc
->bulkOffset
[eCompA
]);
1616 fprintf(s
->fpout
, "# bulk-offsetB = %g\n", sc
->bulkOffset
[eCompB
]);
1619 fprintf(s
->fpout
, "#\n");
1620 fprintf(s
->fpout
, "# split0 cylinder radius %f nm, up %f nm, down %f nm\n",
1621 sc
->cyl0r
, sc
->cyl0u
, sc
->cyl0l
);
1622 fprintf(s
->fpout
, "# split1 cylinder radius %f nm, up %f nm, down %f nm\n",
1623 sc
->cyl1r
, sc
->cyl1u
, sc
->cyl1l
);
1625 fprintf(s
->fpout
, "#\n");
1628 fprintf(s
->fpout
, "# Coupling constant (number of swap attempt steps to average over): %d (translates to %f ps).\n",
1629 sc
->nAverage
, sc
->nAverage
*sc
->nstswap
*ir
->delta_t
);
1630 fprintf(s
->fpout
, "# Threshold is %f\n", sc
->threshold
);
1631 fprintf(s
->fpout
, "#\n");
1632 fprintf(s
->fpout
, "# Remarks about which atoms passed which channel use global atoms numbers starting at one.\n");
1641 /* Prepare for parallel or serial run */
1644 for (ig
= 0; ig
< eGrpNr
; ig
++)
1646 g
= &(s
->group
[ig
]);
1654 for (ig
= 0; ig
< eGrpNr
; ig
++)
1656 g
= &(s
->group
[ig
]);
1657 g
->nat_loc
= g
->nat
;
1658 g
->ind_loc
= g
->ind
;
1659 /* c_ind_loc needs to be set to identity in the serial case */
1660 for (i
= 0; i
< g
->nat
; i
++)
1662 g
->c_ind_loc
[i
] = i
;
1667 /* Allocate memory for the ion counts time window */
1668 for (ic
= 0; ic
< eCompNR
; ic
++)
1670 for (ii
= 0; ii
< eIonNR
; ii
++)
1672 snew(s
->comp
[ic
][ii
].nat_past
, sc
->nAverage
);
1676 /* Get the initial ion concentrations and let the other nodes know */
1679 swapstate
->nions
= s
->group
[eGrpIons
].nat
;
1683 get_initial_ioncounts_from_cpt(ir
, swapstate
, cr
, bVerbose
);
1687 fprintf(stderr
, "%s Determining initial ion counts.\n", SwS
);
1688 get_initial_ioncounts(ir
, x
, box
, cr
, bRerun
);
1691 /* Prepare (further) checkpoint writes ... */
1694 /* Consistency check */
1695 if (swapstate
->nAverage
!= sc
->nAverage
)
1697 gmx_fatal(FARGS
, "%s Ion count averaging steps mismatch! checkpoint: %d, tpr: %d",
1698 SwS
, swapstate
->nAverage
, sc
->nAverage
);
1703 swapstate
->nAverage
= sc
->nAverage
;
1705 fprintf(stderr
, "%s Setting pointers for checkpoint writing\n", SwS
);
1706 for (ic
= 0; ic
< eCompNR
; ic
++)
1708 for (ii
= 0; ii
< eIonNR
; ii
++)
1710 swapstate
->nat_req_p
[ic
][ii
] = &(s
->comp
[ic
][ii
].nat_req
);
1711 swapstate
->nat_past_p
[ic
][ii
] = &(s
->comp
[ic
][ii
].nat_past
[0]);
1712 swapstate
->inflow_netto_p
[ic
][ii
] = &(s
->comp
[ic
][ii
].inflow_netto
);
1716 /* Determine the total charge imbalance */
1717 s
->deltaQ
= ( (-1) * s
->comp
[eCompA
][eIonNEG
].nat_req
+ s
->comp
[eCompA
][eIonPOS
].nat_req
)
1718 - ( (-1) * s
->comp
[eCompB
][eIonNEG
].nat_req
+ s
->comp
[eCompB
][eIonPOS
].nat_req
);
1722 fprintf(stderr
, "%s Requested charge imbalance is Q(A) - Q(B) = %g e.\n", SwS
, s
->deltaQ
);
1726 fprintf(s
->fpout
, "# Requested charge imbalance is Q(A)-Q(B) = %g e.\n", s
->deltaQ
);
1732 bc_initial_concentrations(cr
, ir
->swap
);
1735 /* Put the time-averaged number of ions for all compartments */
1736 for (ic
= 0; ic
< eCompNR
; ic
++)
1738 for (ii
= 0; ii
< eIonNR
; ii
++)
1740 update_time_window(&(s
->comp
[ic
][ii
]), sc
->nAverage
, -1);
1744 /* Initialize arrays that keep track of through which channel the ions go */
1745 detect_flux_per_channel_init(cr
, s
, swapstate
, bStartFromCpt
);
1747 /* We need to print the legend if we open this file for the first time. */
1748 if (MASTER(cr
) && !bAppend
)
1750 print_ionlist_legend(ir
, oenv
);
1755 void dd_make_local_swap_groups(gmx_domdec_t
*dd
, t_swapcoords
*sc
)
1761 /* Make ion group, split groups and solvent group */
1762 for (ig
= 0; ig
< eGrpNr
; ig
++)
1764 g
= &(sc
->si_priv
->group
[ig
]);
1765 dd_make_local_group_indices(dd
->ga2la
, g
->nat
, g
->ind
,
1766 &(g
->nat_loc
), &(g
->ind_loc
), &(g
->nalloc_loc
), g
->c_ind_loc
);
1771 /*! \brief Do we need to swap ions with water molecules at this step?
1773 * From the requested and average ion counts we determine whether a swap is needed
1774 * at this time step.
1776 static gmx_bool
need_swap(t_swapcoords
*sc
)
1783 for (ic
= 0; ic
< eCompNR
; ic
++)
1785 for (ii
= 0; ii
< eIonNR
; ii
++)
1787 if (s
->comp
[ic
][ii
].nat_req
- s
->comp
[ic
][ii
].nat_av
>= sc
->threshold
)
1797 /*! \brief Return the index of an atom or molecule suitable for swapping.
1799 * Returns the index of an atom that is far off the compartment boundaries,
1800 * that is near to the bulk layer to/from which the swaps take place.
1801 * Other atoms of the molecule (if any) will directly follow the returned index.
1803 * \param[in] comp Structure containing compartment-specific data.
1804 * \param[in] apm Atoms per molecule - just return the first atom index of a molecule
1806 * \returns Index of the first atom of the molecule chosen for a position exchange.
1808 static int get_index_of_distant_atom(
1809 t_compartment
*comp
,
1813 real d
= GMX_REAL_MAX
;
1816 /* comp->nat contains the original number of atoms in this compartment
1817 * prior to doing any swaps. Some of these atoms may already have been
1818 * swapped out, but then they are marked with a distance of GMX_REAL_MAX
1820 for (i
= 0; i
< comp
->nat_old
; i
+= apm
)
1822 if (comp
->dist
[i
] < d
)
1825 d
= comp
->dist
[ibest
];
1831 gmx_fatal(FARGS
, "Could not get index of swap atom. Compartment atoms %d before swaps, atoms per molecule %d.",
1832 comp
->nat_old
, apm
);
1835 /* Set the distance of this index to infinity such that it won't get selected again in
1838 comp
->dist
[ibest
] = GMX_REAL_MAX
;
1840 return comp
->ind
[ibest
];
1844 /*! \brief Swaps centers of mass and makes molecule whole if broken */
1845 static void translate_positions(
1853 rvec reference
, dx
, correctPBCimage
;
1856 /* Use the first atom as the reference for PBC */
1857 copy_rvec(x
[0], reference
);
1859 for (i
= 0; i
< apm
; i
++)
1861 /* PBC distance between position and reference */
1862 pbc_dx(pbc
, x
[i
], reference
, dx
);
1864 /* Add PBC distance to reference */
1865 rvec_add(reference
, dx
, correctPBCimage
);
1867 /* Subtract old_com from correct image and add new_com */
1868 rvec_dec(correctPBCimage
, old_com
);
1869 rvec_inc(correctPBCimage
, new_com
);
1871 copy_rvec(correctPBCimage
, x
[i
]);
1876 /*! \brief Write back the the modified local positions from the collective array to the official positions. */
1877 static void apply_modified_positions(
1884 for (l
= 0; l
< g
->nat_loc
; l
++)
1886 /* Get the right local index to write to */
1888 /* Where is the local atom in the collective array? */
1889 cind
= g
->c_ind_loc
[l
];
1891 /* Copy the possibly modified position */
1892 copy_rvec(g
->xc
[cind
], x
[ii
]);
1897 gmx_bool
do_swapcoords(t_commrec
*cr
,
1901 gmx_wallcycle_t wcycle
,
1910 int j
, ii
, ic
, ig
, im
, nswaps
;
1911 gmx_bool bSwap
= FALSE
;
1913 real vacancy
[eCompNR
][eIonNR
];
1915 rvec solvent_center
, ion_center
;
1917 gmx_mtop_atomlookup_t alook
= NULL
;
1920 wallcycle_start(wcycle
, ewcSWAP
);
1925 /* Assemble the positions of the ions. Since we do not need to make these
1926 * whole, we pass NULL as third argument. */
1927 g
= &(s
->group
[eGrpIons
]);
1928 communicate_group_positions(cr
, g
->xc
, NULL
, NULL
, FALSE
,
1929 x
, g
->nat
, g
->nat_loc
, g
->ind_loc
, g
->c_ind_loc
, NULL
, NULL
);
1931 /* Assemble the positions of the split groups, i.e. the channels.
1932 * Here we also pass a shifts array to communicate_group_positions(), so that it can make
1933 * the molecules whole even in cases where they span more than half of the box in
1935 for (ig
= eGrpSplit0
; ig
<= eGrpSplit1
; ig
++)
1937 g
= &(s
->group
[ig
]);
1938 communicate_group_positions(cr
, g
->xc
, g
->xc_shifts
, g
->xc_eshifts
, TRUE
,
1939 x
, g
->nat
, g
->nat_loc
, g
->ind_loc
, g
->c_ind_loc
, g
->xc_old
, box
);
1941 get_center(g
->xc
, g
->m
, g
->nat
, g
->center
); /* center of split groups == channels */
1944 /* Set up the compartments and get lists of atoms in each compartment,
1945 * determine how many ions each compartment contains */
1946 compartmentalize_ions(cr
, sc
, box
, step
, s
->fpout
, bRerun
);
1948 /* Output how many ions are in the compartments */
1951 print_ionlist(s
, t
, "");
1954 /* If we are doing a rerun, we are finished here, since we cannot perform
1961 /* Do we have to perform a swap? */
1962 bSwap
= need_swap(sc
);
1965 /* Since we here know that we have to perform ion/water position exchanges,
1966 * we now assemble the solvent positions */
1967 g
= &(s
->group
[eGrpSolvent
]);
1968 communicate_group_positions(cr
, g
->xc
, NULL
, NULL
, FALSE
,
1969 x
, g
->nat
, g
->nat_loc
, g
->ind_loc
, g
->c_ind_loc
, NULL
, NULL
);
1971 compartmentalize_solvent(cr
, sc
, box
, s
->fpout
);
1973 /* Determine where ions are missing and where ions are too many */
1974 for (ic
= 0; ic
< eCompNR
; ic
++)
1976 for (ii
= 0; ii
< eIonNR
; ii
++)
1978 vacancy
[ic
][ii
] = s
->comp
[ic
][ii
].nat_req
- s
->comp
[ic
][ii
].nat_av
;
1982 /* Remember the original number of ions per compartment */
1983 for (ic
= 0; ic
< eCompNR
; ic
++)
1985 s
->compsol
[ic
].nat_old
= s
->compsol
[ic
].nat
;
1986 for (ii
= 0; ii
< eIonNR
; ii
++)
1988 s
->comp
[ic
][ii
].nat_old
= s
->comp
[ic
][ii
].nat
;
1992 /* Now actually correct the number of ions */
1994 alook
= gmx_mtop_atomlookup_init(mtop
);
1995 for (ic
= 0; ic
< eCompNR
; ic
++)
1997 for (ii
= 0; ii
< eIonNR
; ii
++)
1999 while (vacancy
[ic
][ii
] >= sc
->threshold
)
2001 /* Swap in an ion */
2003 /* Get the xc-index of the first atom of a solvent molecule of this compartment */
2004 isol
= get_index_of_distant_atom(&(s
->compsol
[ic
]), s
->group
[eGrpSolvent
].apm
);
2006 /* Get the xc-index of an ion from the other compartment */
2007 iion
= get_index_of_distant_atom(&(s
->comp
[(ic
+1)%eCompNR
][ii
]), s
->group
[eGrpIons
].apm
);
2009 /* Get the solvent molecule's center of mass */
2010 for (im
= 0; im
< s
->group
[eGrpSolvent
].apm
; im
++)
2012 gmx_mtop_atomnr_to_atom(alook
, s
->group
[eGrpSolvent
].ind
[isol
+im
], &atom
);
2013 s
->group
[eGrpSolvent
].m
[im
] = atom
->m
;
2015 get_molecule_center(&(s
->group
[eGrpSolvent
].xc
[isol
]), s
->group
[eGrpSolvent
].apm
, s
->group
[eGrpSolvent
].m
, solvent_center
, s
->pbc
);
2016 get_molecule_center(&(s
->group
[eGrpIons
].xc
[iion
]), s
->group
[eGrpIons
].apm
, NULL
, ion_center
, s
->pbc
);
2018 /* subtract com_solvent and add com_ion */
2019 translate_positions(&(s
->group
[eGrpSolvent
].xc
[isol
]), s
->group
[eGrpSolvent
].apm
, solvent_center
, ion_center
, s
->pbc
);
2020 /* For the ion, subtract com_ion and add com_solvent */
2021 translate_positions(&(s
->group
[eGrpIons
].xc
[iion
]), s
->group
[eGrpIons
].apm
, ion_center
, solvent_center
, s
->pbc
);
2024 vacancy
[(ic
+1) % eCompNR
][ii
]++;
2026 /* Keep track of the changes */
2027 s
->comp
[ic
][ii
].nat
++;
2028 s
->comp
[(ic
+1) % eCompNR
][ii
].nat
--;
2029 s
->comp
[ic
][ii
].inflow_netto
++;
2030 s
->comp
[(ic
+1) % eCompNR
][ii
].inflow_netto
--;
2031 /* Correct the past time window to still get the right averages from now on */
2032 s
->comp
[ic
][ii
].nat_av
++;
2033 s
->comp
[(ic
+1) % eCompNR
][ii
].nat_av
--;
2034 for (j
= 0; j
< sc
->nAverage
; j
++)
2036 s
->comp
[ic
][ii
].nat_past
[j
]++;
2037 s
->comp
[(ic
+1) % eCompNR
][ii
].nat_past
[j
]--;
2039 /* Clear ion history */
2042 s
->group
[eGrpIons
].channel_label
[iion
] = eChHistPassedNone
;
2043 s
->group
[eGrpIons
].comp_from
[iion
] = eDomainNotset
;
2045 /* That was the swap */
2050 gmx_mtop_atomlookup_destroy(alook
);
2054 fprintf(stderr
, "%s Performed %d swap%s in step %" GMX_PRId64
".\n", SwS
, nswaps
, nswaps
> 1 ? "s" : "", step
);
2056 if (s
->fpout
!= NULL
)
2058 print_ionlist(s
, t
, " # after swap");
2061 /* Write back the the modified local positions from the collective array to the official coordinates */
2062 apply_modified_positions(&(s
->group
[eGrpIons
]), x
);
2063 apply_modified_positions(&(s
->group
[eGrpSolvent
]), x
);
2064 } /* end of if(bSwap) */
2066 wallcycle_stop(wcycle
, ewcSWAP
);