Gromacs  2021.2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Classes | Typedefs | Enumerations | Functions | Variables
#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, PbcType pbcType, 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_swapinit_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

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]ciIndex of this ion in the collective xc array.
[in,out]compCompartment to add this atom to.
[in]distanceShortest 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]w1Position of 'wall' atom 1.
[in]w2Position of 'wall' atom 2.
[in]xPosition of the ion or the water molecule under consideration.
[in]lLength of the box, from || to || in the sketch.
[in]bulkOffsetWhere is the bulk layer "b" to be found between w1 and w2?
[out]distance_from_bDistance 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]crPointer to MPI communication data.
[in]stepThe number of the MD time step.
[in]tThe time.
[in]irStructure containing MD input parameters
[in,out]sThe structure needed for position swapping.
[in]wcycleCount wallcycles of swap routines for diagnostic output.
[in]xPositions of home particles this node owns.
[in]boxThe simulation box.
[in]bVerboseShould we be quiet or verbose?
[in]bRerunAre 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]sPointer 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]compStructure containing compartment-specific data.
[in]molnameName 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]fplogGeneral output file, normally md.log.
[in]irStructure containing MD input parameters, among those also the structure needed for position swapping.
[in]fnOutput file name for swap data.
[in]mtopMolecular topology.
[in]globalStateThe global state, only used on the master rank.
[in]ohContains struct with swap data that is read from or written to checkpoint.
[in]crPointer to MPI communication data.
[in]atomSetsManager tending to swap atom indices.
[in]oenvNeeded to open the swap output XVGR file.
[in]mdrunOptionsOptions for mdrun.
[in]startingBehaviorDescribes 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]pointThe position (xyz) under consideration.
[in]centerThe center of the cylinder.
[in]d_upThe upper extension of the cylinder.
[in]d_downThe lower extension.
[in]r_cyl2Cylinder radius squared.
[in]pbcStructure with info about periodic boundary conditions.
[in]normalThe 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,
PbcType  pbcType,
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.