#include "gmxpre.h"
#include "swapcoords.h"
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <memory>
#include <string>
#include <vector>
#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/domdec/localatomset.h"
#include "gromacs/domdec/localatomsetmanager.h"
#include "gromacs/fileio/confio.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/xvgr.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/groupcoord.h"
#include "gromacs/mdrunutility/handlerestart.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/imdmodule.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/mdrunoptions.h"
#include "gromacs/mdtypes/observableshistory.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/mdtypes/swaphistory.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/timing/wallcycle.h"
#include "gromacs/topology/mtop_lookup.h"
#include "gromacs/topology/topology.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/pleasecite.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/snprintf.h"

Include dependency graph for swapcoords.cpp:

Description

Implements functions in swapcoords.h.

Author: Carsten Kutzner ckutz.nosp@m.ne@g.nosp@m.wdg.d.nosp@m.e

Classes
class	gmx::SwapCoordinates
	Implement Computational Electrophysiology swapping. More...

struct	swap_compartment
	Structure containing compartment-specific data. More...

struct	swap_group
	This structure contains data needed for the groups involved in swapping: split group 0, split group 1, solvent group, ion groups. More...

struct	t_swap
	Main (private) data structure for the position swapping protocol. More...

Typedefs
typedef struct swap_compartment	t_compartment
	Structure containing compartment-specific data. More...

typedef struct swap_group	t_swapgrp
	This structure contains data needed for the groups involved in swapping: split group 0, split group 1, solvent group, ion groups. More...

Enumerations
enum	eChannelHistory { eChHistPassedNone, eChHistPassedCh0, eChHistPassedCh1, eChHistNr }
	Keep track of through which channel the ions have passed.

enum	eDomain { eDomainNotset, eDomainA, eDomainB, eDomainNr }
	Domain identifier. More...

Functions
std::unique_ptr< IMDModule >	gmx::createSwapCoordinatesModule ()
	Creates a module for Computational Electrophysiology swapping.

static gmx_bool	is_in_channel (rvec point, rvec center, real d_up, real d_down, real r_cyl2, t_pbc *pbc, int normal)
	Check whether point is in channel. More...

static void	print_ionlist (t_swap *s, double time, const char comment[])
	Prints output to CompEL output file. More...

static void	get_molecule_center (rvec x[], int nat, const real weights, rvec center, t_pbc pbc)
	Get the center of a group of nat atoms. More...

static gmx_bool	compartment_contains_atom (real w1, real w2, real x, real l, real bulkOffset, real *distance_from_b)
	Return TRUE if position x of ion (or water) is found in the compartment, i.e. between w1 and w2. More...

static void	update_time_window (t_compartment *comp, int values, int replace)
	Updates the time-averaged number of ions in a compartment.

static void	add_to_list (int ci, t_compartment *comp, real distance)
	Add the atom with collective index ci to the atom list in compartment 'comp'. More...

static void	get_compartment_boundaries (int c, t_swap s, const matrix box, real left, real *right)
	Determine the compartment boundaries from the channel centers.

static void	detect_flux_per_channel (t_swapgrp g, int iAtom, int comp, rvec atomPosition, unsigned char comp_now, unsigned char comp_from, unsigned char channel_label, t_swapcoords sc, t_swap s, real cyl0_r2, real cyl1_r2, int64_t step, gmx_bool bRerun, FILE *fpout)
	Determine the per-channel ion flux. More...

static void	sortMoleculesIntoCompartments (t_swapgrp g, t_commrec cr, t_swapcoords sc, t_swap s, const matrix box, int64_t step, FILE *fpout, gmx_bool bRerun, gmx_bool bIsSolvent)
	Determines which ions or solvent molecules are in compartment A and B.

static void	get_initial_ioncounts (const t_inputrec ir, t_swap s, const rvec x[], const matrix box, t_commrec *cr, gmx_bool bRerun)
	Find out how many group atoms are in the compartments initially.

static void	get_initial_ioncounts_from_cpt (const t_inputrec ir, t_swap s, swaphistory_t swapstate, t_commrec cr, gmx_bool bVerbose)
	Copy history of ion counts from checkpoint file. More...

static void	bc_initial_concentrations (t_commrec cr, t_swapcoords swap, t_swap *s)
	The master lets all others know about the initial ion counts.

static void	check_swap_groups (t_swap *s, int nat, gmx_bool bVerbose)
	Ensure that each atom belongs to at most one of the swap groups.

static int	get_group_apm_check (int igroup, t_swap s, gmx_bool bVerbose, gmx_mtop_t mtop)
	Get the number of atoms per molecule for this group. More...

static void	print_ionlist_legend (const t_inputrec ir, t_swap s, const gmx_output_env_t *oenv)
	Print the legend to the swap output file. More...

static void	detect_flux_per_channel_init (t_swap s, swaphistory_t swapstate, const bool isRestart)
	Initialize channel ion flux detection routine. More...

static void	outputStartStructureIfWanted (gmx_mtop_t mtop, rvec x, int ePBC, const matrix box)
	Outputs the initial structure to PDB file for debugging reasons. More...

static void	init_swapstate (swaphistory_t swapstate, t_swapcoords sc, t_swap s, gmx_mtop_t mtop, const rvec x, const matrix box, const t_inputrec ir)
	Initialize the swapstate structure, used for checkpoint writing. More...

static real	getRequestedChargeImbalance (t_swap *s)
	Determine the total charge imbalance resulting from the swap groups.

static void	copyIndicesToGroup (const int indIons, int nIons, t_swapGroup g, t_commrec *cr)
	Sorts anions and cations into two separate groups. More...

static void	convertOldToNewGroupFormat (t_swapcoords sc, gmx_mtop_t mtop, gmx_bool bVerbose, t_commrec *cr)
	Converts old .tpr file CompEL contents to new data layout. More...

static gmx_bool	bConvertFromOldTpr (t_swapcoords *sc)
	Returns TRUE if we started from an old .tpr. More...

t_swap *	init_swapcoords (FILE fplog, const t_inputrec ir, const char fn, gmx_mtop_t mtop, const t_state globalState, ObservablesHistory oh, t_commrec cr, gmx::LocalAtomSetManager atomSets, const gmx_output_env_t *oenv, const gmx::MdrunOptions &mdrunOptions, const gmx::StartingBehavior startingBehavior)
	Initialize ion / water position swapping ("Computational Electrophysiology"). More...

void	finish_swapcoords (t_swap *s)
	Finalizes ion / water position swapping, if it was active. More...

static gmx_bool	need_swap (t_swapcoords sc, t_swap s)
	Do we need to swap a molecule in any of the ion groups with a water molecule at this step? More...

static int	get_index_of_distant_atom (t_compartment *comp, const char molname[])
	Return the index of an atom or molecule suitable for swapping. More...

static void	translate_positions (rvec x, int apm, rvec old_com, rvec new_com, t_pbc pbc)
	Swaps centers of mass and makes molecule whole if broken.

static void	apply_modified_positions (swap_group *g, rvec x[])
	Write back the modified local positions from the collective array to the official positions.

gmx_bool	do_swapcoords (t_commrec cr, int64_t step, double t, t_inputrec ir, t_swap s, gmx_wallcycle wcycle, rvec x[], matrix box, gmx_bool bVerbose, gmx_bool bRerun)
	"Computational Electrophysiology" main routine within MD loop. More...

Variables
static const char *	SwS = { "SWAP:" }
	For output that comes from the swap module.

static const char *	SwSEmpty = { " " }
	Placeholder for multi-line output.

static const char *	CompStr [eCompNR] = { "A", "B" }
	Compartment name.

static const char *	SwapStr [eSwapTypesNR+1]
	Name for the swap types. More...

static const char *	DimStr [3+1] = { "X", "Y", "Z", nullptr }
	Name for the swap dimension. More...

static const char *	ChannelString [eChHistNr] = { "none", "channel0", "channel1" }
	Name for the channels.

static const char *	DomainString [eDomainNr]
	Name for the domains. More...

Typedef Documentation

typedef struct swap_compartment t_compartment

Structure containing compartment-specific data.

typedef struct swap_group t_swapgrp

This structure contains data needed for the groups involved in swapping: split group 0, split group 1, solvent group, ion groups.

Enumeration Type Documentation

enum eDomain

Domain identifier.

Keeps track of from which compartment the ions came before passing the channel.

Function Documentation

static void add_to_list	(	int	ci,
		t_compartment *	comp,
		real	distance
	)

static

Add the atom with collective index ci to the atom list in compartment 'comp'.

Parameters

[in]	ci	Index of this ion in the collective xc array.
[in,out]	comp	Compartment to add this atom to.
[in]	distance	Shortest distance of this atom to the bulk layer, from which ion/water pairs are selected for swapping.

static gmx_bool bConvertFromOldTpr ( t_swapcoords * sc )

static

Returns TRUE if we started from an old .tpr.

Then we need to re-sort anions and cations into separate groups

static gmx_bool compartment_contains_atom	(	real	w1,
		real	w2,
		real	x,
		real	l,
		real	bulkOffset,
		real *	distance_from_b
	)

static

Return TRUE if position x of ion (or water) is found in the compartment, i.e. between w1 and w2.

One can define and additional offset "b" if one wants to exchange ions/water to or from a plane not directly in the middle of w1 and w2. The offset can be in ]-1.0, ..., +1.0 [. A bulkOffset of 0.0 means 'no offset', so the swap-layer is directly in the middle between w1 and w2. Offsets -1.0 < b < 0.0 will yield swaps nearer to w1, whereas offsets 0.0 < 0 < +1.0 will yield swaps nearer to w2.

*
* ||--------------+-------------|-------------+------------------------||
*                w1  ? ? ? ? ? ? ? ? ? ? ?   w2
* ||--------------+-------------|----b--------+------------------------||
*                -1            0.0           +1
*
* 

Parameters

[in]	w1	Position of 'wall' atom 1.
[in]	w2	Position of 'wall' atom 2.
[in]	x	Position of the ion or the water molecule under consideration.
[in]	l	Length of the box, from \|\| to \|\| in the sketch.
[in]	bulkOffset	Where is the bulk layer "b" to be found between w1 and w2?
[out]	distance_from_b	Distance of x to the bulk layer "b".

Returns: TRUE if x is between w1 and w2.

Also computes the distance of x to the compartment center (the layer that is normally situated in the middle of w1 and w2 that would be considered as having the bulk concentration of ions).

static void convertOldToNewGroupFormat	(	t_swapcoords *	sc,
		gmx_mtop_t *	mtop,
		gmx_bool	bVerbose,
		t_commrec *	cr
	)

static

Converts old .tpr file CompEL contents to new data layout.

If we have read an old .tpr file (tpxv <= tpxv_CompElPolyatomicIonsAndMultipleIonTypes), anions and cations are stored together in group #3. In the new format we store each ion type in a separate group. The 'classic' groups are: #0 split group 0 - OK #1 split group 1 - OK #2 solvent - OK #3 anions - contains also cations, needs to be converted #4 cations - empty before conversion

static void copyIndicesToGroup	(	const int *	indIons,
		int	nIons,
		t_swapGroup *	g,
		t_commrec *	cr
	)

static

Sorts anions and cations into two separate groups.

This routine should be called for the 'anions' and 'cations' group, of which the indices were lumped together in the older version of the code.

static void detect_flux_per_channel	(	t_swapgrp *	g,
		int	iAtom,
		int	comp,
		rvec	atomPosition,
		unsigned char *	comp_now,
		unsigned char *	comp_from,
		unsigned char *	channel_label,
		t_swapcoords *	sc,
		t_swap *	s,
		real	cyl0_r2,
		real	cyl1_r2,
		int64_t	step,
		gmx_bool	bRerun,
		FILE *	fpout
	)

static

Determine the per-channel ion flux.

To determine the flux through the individual channels, we remember the compartment and channel history of each ion. An ion can be either in channel0 or channel1, or in the remaining volume of compartment A or B.

*    +-----------------+
*    |                 | B
*    |                 | B compartment
*    ||||||||||0|||||||| bilayer with channel 0
*    |                 | A
*    |                 | A
*    |                 | A compartment
*    |                 | A
*    |||||1||||||||||||| bilayer with channel 1
*    |                 | B
*    |                 | B compartment
*    +-----------------+
*
* 

static void detect_flux_per_channel_init	(	t_swap *	s,
		swaphistory_t *	swapstate,
		const bool	isRestart
	)

static

Initialize channel ion flux detection routine.

Initialize arrays that keep track of where the ions come from and where they go.

gmx_bool do_swapcoords	(	t_commrec *	cr,
		int64_t	step,
		double	t,
		t_inputrec *	ir,
		t_swap *	s,
		gmx_wallcycle *	wcycle,
		rvec	x[],
		matrix	box,
		gmx_bool	bVerbose,
		gmx_bool	bRerun
	)

"Computational Electrophysiology" main routine within MD loop.

Parameters

[in]	cr	Pointer to MPI communication data.
[in]	step	The number of the MD time step.
[in]	t	The time.
[in]	ir	Structure containing MD input parameters
[in,out]	s	The structure needed for position swapping.
[in]	wcycle	Count wallcycles of swap routines for diagnostic output.
[in]	x	Positions of home particles this node owns.
[in]	box	The simulation box.
[in]	bVerbose	Should we be quiet or verbose?
[in]	bRerun	Are we doing a rerun?

Returns: Whether at least one pair of molecules was swapped.

void finish_swapcoords ( t_swap * s )

Finalizes ion / water position swapping, if it was active.

Parameters

[in] s Pointer to swap data.

static int get_group_apm_check	(	int	igroup,
		t_swap *	s,
		gmx_bool	bVerbose,
		gmx_mtop_t *	mtop
	)

static

Get the number of atoms per molecule for this group.

Also ensure that all the molecules in this group have this number of atoms.

static int get_index_of_distant_atom	(	t_compartment *	comp,
		const char	molname[]
	)

static

Return the index of an atom or molecule suitable for swapping.

Returns the index of an atom that is far off the compartment boundaries, that is near to the bulk layer to/from which the swaps take place. Other atoms of the molecule (if any) will directly follow the returned index.

Parameters

[in]	comp	Structure containing compartment-specific data.
[in]	molname	Name of the molecule.

Returns: Index of the first atom of the molecule chosen for a position exchange.

static void get_initial_ioncounts_from_cpt	(	const t_inputrec *	ir,
		t_swap *	s,
		swaphistory_t *	swapstate,
		t_commrec *	cr,
		gmx_bool	bVerbose
	)

static

Copy history of ion counts from checkpoint file.

When called, the checkpoint file has already been read in. Here we copy over the values from .cpt file to the swap data structure.

static void get_molecule_center	(	rvec	x[],
		int	nat,
		const real *	weights,
		rvec	center,
		t_pbc *	pbc
	)

static

Get the center of a group of nat atoms.

Since with PBC an atom group might not be whole, use the first atom as the reference atom and determine the center with respect to this reference.

t_swap* init_swapcoords	(	FILE *	fplog,
		const t_inputrec *	ir,
		const char *	fn,
		gmx_mtop_t *	mtop,
		const t_state *	globalState,
		ObservablesHistory *	oh,
		t_commrec *	cr,
		gmx::LocalAtomSetManager *	atomSets,
		const gmx_output_env_t *	oenv,
		const gmx::MdrunOptions &	mdrunOptions,
		gmx::StartingBehavior	startingBehavior
	)

Initialize ion / water position swapping ("Computational Electrophysiology").

This routine does the memory allocation for various helper arrays, opens the output file, sets up swap data checkpoint writing, etc. and returns it.

Parameters

[in]	fplog	General output file, normally md.log.
[in]	ir	Structure containing MD input parameters, among those also the structure needed for position swapping.
[in]	fn	Output file name for swap data.
[in]	mtop	Molecular topology.
[in]	globalState	The global state, only used on the master rank.
[in]	oh	Contains struct with swap data that is read from or written to checkpoint.
[in]	cr	Pointer to MPI communication data.
[in]	atomSets	Manager tending to swap atom indices.
[in]	oenv	Needed to open the swap output XVGR file.
[in]	mdrunOptions	Options for mdrun.
[in]	startingBehavior	Describes whether this is a restart appending to output files

static void init_swapstate	(	swaphistory_t *	swapstate,
		t_swapcoords *	sc,
		t_swap *	s,
		gmx_mtop_t *	mtop,
		const rvec *	x,
		const matrix	box,
		const t_inputrec *	ir
	)

static

Initialize the swapstate structure, used for checkpoint writing.

The swapstate struct stores the information we need to make the channels whole again after restarts from a checkpoint file. Here we do the following: a) If we did not start from .cpt, we prepare the struct for proper .cpt writing, b) if we did start from .cpt, we copy over the last whole structures from .cpt, c) in any case, for subsequent checkpoint writing, we set the pointers in swapstate to the x_old arrays, which contain the correct PBC representation of multimeric channels at the last time step.

static gmx_bool is_in_channel	(	rvec	point,
		rvec	center,
		real	d_up,
		real	d_down,
		real	r_cyl2,
		t_pbc *	pbc,
		int	normal
	)

static

Check whether point is in channel.

A channel is a cylinder defined by a disc with radius r around its center c. The thickness of the cylinder is d_up - d_down.

*               ^  d_up
*               |
*     r         |
*     <---------c--------->
*               |
*               v  d_down
*
* 

Parameters

[in]	point	The position (xyz) under consideration.
[in]	center	The center of the cylinder.
[in]	d_up	The upper extension of the cylinder.
[in]	d_down	The lower extension.
[in]	r_cyl2	Cylinder radius squared.
[in]	pbc	Structure with info about periodic boundary conditions.
[in]	normal	The membrane normal direction is typically 3, i.e. z, but can be x or y also.

Returns: Whether the point is inside the defined cylindric channel.

static gmx_bool need_swap	(	t_swapcoords *	sc,
		t_swap *	s
	)

static

Do we need to swap a molecule in any of the ion groups with a water molecule at this step?

From the requested and average molecule counts we determine whether a swap is needed at this time step.

static void outputStartStructureIfWanted	(	gmx_mtop_t *	mtop,
		rvec *	x,
		int	ePBC,
		const matrix	box
	)

static

Outputs the initial structure to PDB file for debugging reasons.

Output the starting structure so that in case of multimeric channels the user can check whether we have the correct PBC image for all atoms. If this is not correct, the ion counts per channel will be very likely wrong.

static void print_ionlist	(	t_swap *	s,
		double	time,
		const char	comment[]
	)

static

Prints output to CompEL output file.

Prints to swap output file how many ions are in each compartment, where the centers of the split groups are, and how many ions of each type passed the channels.

static void print_ionlist_legend	(	const t_inputrec *	ir,
		t_swap *	s,
		const gmx_output_env_t *	oenv
	)

static

Print the legend to the swap output file.

Also print the initial values of ion counts and position of split groups.

Variable Documentation

const char* DimStr[3+1] = { "X", "Y", "Z", nullptr }

static

Name for the swap dimension.

const char* DomainString[eDomainNr]

static

Initial value:

= { "not_assigned", "Domain_A",

"Domain_B" }

Name for the domains.

const char* SwapStr[eSwapTypesNR+1]

static

Initial value:

= { "", "X-", "Y-", "Z-",

nullptr }

Name for the swap types.

Description

Classes

Typedefs

Enumerations

Functions

Variables

Typedef Documentation

Enumeration Type Documentation

Function Documentation

Variable Documentation