#include "gmxpre.h"
#include "gromacs/pbcutil/pbc.h"
#include <cmath>
#include <algorithm>
#include "gromacs/math/functions.h"
#include "gromacs/math/units.h"
#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vecdump.h"
#include "gromacs/pbcutil/ishift.h"
#include "gromacs/pbcutil/mshift.h"
#include "gromacs/topology/topology.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"

Include dependency graph for pbc.cpp:

Description

Implements routines in pbc.h.

Utility functions for handling periodic boundary conditions. Mainly used in analysis tools.

Enumerations
enum	{ epbcdxRECTANGULAR = 1, epbcdxTRICLINIC, epbcdx2D_RECT, epbcdx2D_TRIC, epbcdx1D_RECT, epbcdx1D_TRIC, epbcdxSCREW_RECT, epbcdxSCREW_TRIC, epbcdxNOPBC, epbcdxUNSUPPORTED }

Functions
int	numPbcDimensions (PbcType pbcType)
	Returns the number of dimensions that use pbc. More...

void	dump_pbc (FILE fp, t_pbc pbc)
	Dump the contents of the pbc structure to the file. More...

const char *	check_box (PbcType pbcType, const matrix box)
	Check the box for consistency. More...

void	matrix_convert (matrix box, const rvec vec, const rvec angleInDegrees)
	Creates box matrix from edge lengths and angles. More...

real	max_cutoff2 (PbcType pbcType, const matrix box)
	Compute the maximum cutoff for the box. More...

PbcType	guessPbcType (const matrix box)
	Guess PBC type. More...

static int	correct_box_elem (FILE *fplog, const int64_t step, matrix box, const int v, const int d)
	Check if the box still obeys the restrictions, if not, correct it.

gmx_bool	correct_box (FILE *fplog, const int64_t step, matrix box)
	Corrects the box if necessary. More...

static void	low_set_pbc (t_pbc *pbc, PbcType pbcType, const ivec dd_pbc, const matrix box)
	Do the real arithmetic for filling the pbc struct.

void	set_pbc (t_pbc *pbc, PbcType pbcType, const matrix box)
	Initiate the periodic boundary condition algorithms. More...

t_pbc *	set_pbc_dd (t_pbc *pbc, PbcType pbcType, const ivec domdecCells, gmx_bool bSingleDir, const matrix box)
	Initiate the periodic boundary condition algorithms. More...

void	pbc_dx (const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx)
	Compute distance with PBC. More...

int	pbc_dx_aiuc (const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx)
	Compute distance vector for simple PBC types. More...

void	pbc_dx_d (const t_pbc *pbc, const dvec x1, const dvec x2, dvec dx)
	Compute distance vector in double precision. More...

void	calc_shifts (const matrix box, gmx::ArrayRef< gmx::RVec > shift_vec)
	Computes shift vectors. More...

void	calc_box_center (int ecenter, const matrix box, rvec box_center)
	Calculates the center of the box. More...

void	calc_triclinic_images (const matrix box, rvec img[])
	Calculates the NTRICIMG box images. More...

void	calc_compact_unitcell_vertices (int ecenter, const matrix box, rvec vert[])
	Calculates the NCUCVERT vertices of a compact unitcell. More...

int *	compact_unitcell_edges ()
	Compute unitcell edges. More...

template<bool haveBoxDeformation>
static void	putAtomsInBoxTemplated (PbcType pbcType, const matrix box, const matrix boxDeformation, gmx::ArrayRef< gmx::RVec > x, gmx::ArrayRef< gmx::RVec > v)

void	put_atoms_in_box (PbcType pbcType, const matrix box, gmx::ArrayRef< gmx::RVec > x)
	Put atoms inside the simulations box. More...

void	put_atoms_in_box_omp (PbcType pbcType, const matrix box, const bool haveBoxDeformation, const matrix boxDeformation, gmx::ArrayRef< gmx::RVec > x, gmx::ArrayRef< gmx::RVec > v, gmx_unused int nth)
	Parallellizes put_atoms_in_box() More...

void	put_atoms_in_triclinic_unitcell (int ecenter, const matrix box, gmx::ArrayRef< gmx::RVec > x)
	Put atoms inside triclinic box. More...

void	put_atoms_in_compact_unitcell (PbcType pbcType, int ecenter, const matrix box, gmx::ArrayRef< gmx::RVec > x)
	Put atoms inside the unitcell. More...

static void	low_do_pbc_mtop (FILE fplog, PbcType pbcType, const bool correctVelocitiesForBoxDeformation, const matrix boxDeformation, const matrix box, const gmx_mtop_t mtop, gmx::ArrayRef< gmx::RVec > x, gmx::ArrayRef< gmx::RVec > v, gmx_bool bFirst)
	Make molecules whole by shifting positions. More...

void	do_pbc_first_mtop (FILE fplog, PbcType pbcType, const bool correctVelocitiesForBoxDeformation, const matrix boxDeformation, const matrix box, const gmx_mtop_t mtop, gmx::ArrayRef< gmx::RVec > x, gmx::ArrayRef< gmx::RVec > v)
	Make all molecules whole by shifting positions. More...

void	do_pbc_mtop (PbcType pbcType, const matrix box, const gmx_mtop_t *mtop, rvec x[])
	Make molecules consisting of multiple charge groups whole by shifting positions. More...

void	setBoxDeformationRate (const matrix boxDeformation, const matrix box, matrix boxDeformationRate)
	Sets the box deformation rate. More...

Variables
const gmx::EnumerationArray < PbcType, std::string >	c_pbcTypeNames
	Names for all values in PBC types enumeration. More...

static constexpr real	sc_skewnessMargin = 1.001
	Margin for correction when the box is too skewed.

static constexpr real	sc_boxSkewnessMarginForWarning = sc_skewnessMargin + 0.004
	Margin factor for warning/error message. More...

static gmx_bool	bWarnedGuess = FALSE
	Set to true if warning has been printed.

Function Documentation

void calc_box_center	(	int	ecenter,
		const matrix	box,
		rvec	box_center
	)

Calculates the center of the box.

See the description for the enum ecenter above.

Parameters

[in]	ecenter	Description of center type
[in]	box	The simulation box
[out]	box_center	The center of the box

void calc_compact_unitcell_vertices	(	int	ecenter,
		const matrix	box,
		rvec	vert[]
	)

Calculates the NCUCVERT vertices of a compact unitcell.

Parameters

[in]	ecenter	The center type
[in]	box	The simulation box
[out]	vert	The vertices

void calc_shifts	(	const matrix	box,
		gmx::ArrayRef< gmx::RVec >	shift_vec
	)

Computes shift vectors.

This routine calculates ths shift vectors necessary to use the neighbor searching routine.

Parameters

[in]	box	The simulation box
[out]	shift_vec	The shifting vectors

void calc_triclinic_images	(	const matrix	box,
		rvec	img[]
	)

Calculates the NTRICIMG box images.

Parameters

[in]	box	The simulation box
[out]	img	The triclinic box images

const char* check_box	(	PbcType	pbcType,
		const matrix	box
	)

Check the box for consistency.

When pbcType=PbcTypes::Unset, the type of pbc is guessed from the box matrix.

Parameters

[in]	pbcType	The pbc identifier
[in]	box	The box matrix

Returns: NULL if the box is supported by Gromacs. Otherwise returns a string with the problem.

int* compact_unitcell_edges ( )

Compute unitcell edges.

Returns: an array of unitcell edges of length NCUCEDGE*2, this is an index in vert[], which is calculated by calc_unitcell_vertices. The index consists of NCUCEDGE pairs of vertex indices. The index does not change, so it needs to be retrieved only once.

gmx_bool correct_box	(	FILE *	fplog,
		int64_t	step,
		tensor	box
	)

Corrects the box if necessary.

Checks for un-allowed box angles and corrects the box.

Parameters

[in]	fplog	File for debug output
[in]	step	The MD step number
[in]	box	The simulation cell

Returns: TRUE when the box was corrected.

void do_pbc_first_mtop	(	FILE *	fplog,
		PbcType	pbcType,
		bool	correctVelocitiesForBoxDeformation,
		const matrix	boxDeformation,
		const matrix	box,
		const gmx_mtop_t *	mtop,
		gmx::ArrayRef< gmx::RVec >	x,
		gmx::ArrayRef< gmx::RVec >	v
	)

Make all molecules whole by shifting positions.

Parameters

[in]	fplog	Log file
[in]	pbcType	The PBC type
[in]	correctVelocitiesForBoxDeformation	Whether to correct the velocities for continous box deformation
[in]	boxDeformation	The box deformation velocity
[in]	box	The simulation box
[in]	mtop	System topology definition
[in,out]	x	The coordinates of the atoms
[in,out]	v	The velocities of the atoms, needed with correctVelocitiesForBoxDeformation

void do_pbc_mtop	(	PbcType	pbcType,
		const matrix	box,
		const gmx_mtop_t *	mtop,
		rvec	x[]
	)

Make molecules consisting of multiple charge groups whole by shifting positions.

Parameters

[in]	pbcType	The PBC type
[in]	box	The simulation box
[in]	mtop	System topology definition
[in,out]	x	The coordinates of the atoms

void dump_pbc	(	FILE *	fp,
		t_pbc *	pbc
	)

Dump the contents of the pbc structure to the file.

Parameters

[in]	fp	The file pointer to write to
[in]	pbc	The periodic boundary condition information structure

PbcType guessPbcType ( const matrix box )

Guess PBC type.

Guesses the type of periodic boundary conditions using the box

Parameters

[in] box The box matrix

Returns: The pbc type identifier

static void low_do_pbc_mtop	(	FILE *	fplog,
		PbcType	pbcType,
		const bool	correctVelocitiesForBoxDeformation,
		const matrix	boxDeformation,
		const matrix	box,
		const gmx_mtop_t *	mtop,
		gmx::ArrayRef< gmx::RVec >	x,
		gmx::ArrayRef< gmx::RVec >	v,
		gmx_bool	bFirst
	)

static

Make molecules whole by shifting positions.

Parameters

[in]	fplog	Log file
[in]	pbcType	The PBC type
[in]	correctVelocitiesForBoxDeformation	Whether to correct the velocities for continuous box deformation
[in]	boxDeformation	The box deformation velocity
[in]	box	The simulation box
[in]	mtop	System topology definition
[in,out]	x	The coordinates of the atoms
[in]	v	The velocities of the atoms
[in]	bFirst	Specifier for first-time PBC removal

void matrix_convert	(	matrix	box,
		const rvec	vec,
		const rvec	angleInDegrees
	)

Creates box matrix from edge lengths and angles.

Parameters

[in,out]	box	The box matrix
[in]	vec	The edge lengths
[in]	angleInDegrees	The angles

real max_cutoff2	(	PbcType	pbcType,
		const matrix	box
	)

Compute the maximum cutoff for the box.

Returns the square of the maximum cut-off allowed for the box, taking into account that the grid neighborsearch code and pbc_dx only check combinations of single box-vector shifts.

Parameters

[in]	pbcType	The pbc identifier
[in]	box	The box matrix

Returns: the maximum cut-off.

int numPbcDimensions ( PbcType pbcType )

Returns the number of dimensions that use pbc.

Parameters

[in] pbcType The periodic boundary condition type

Returns: the number of dimensions that use pbc, starting at X

void pbc_dx	(	const t_pbc *	pbc,
		const rvec	x1,
		const rvec	x2,
		rvec	dx
	)

Compute distance with PBC.

Calculate the correct distance vector from x2 to x1 and put it in dx. set_pbc must be called before ever calling this routine.

Note that for triclinic boxes that do not obey the GROMACS unit-cell restrictions, pbc_dx and pbc_dx_aiuc will not correct for PBC.

Parameters

[in,out]	pbc	The pbc information structure
[in]	x1	Coordinates for particle 1
[in]	x2	Coordinates for particle 2
[out]	dx	Distance vector

int pbc_dx_aiuc	(	const t_pbc *	pbc,
		const rvec	x1,
		const rvec	x2,
		rvec	dx
	)

Compute distance vector for simple PBC types.

Calculate the correct distance vector from x2 to x1 and put it in dx, This function can only be used when all atoms are in the rectangular or triclinic unit-cell. set_pbc_dd or set_pbc must be called before ever calling this routine.

Parameters

[in,out]	pbc	The pbc information structure
[in]	x1	Coordinates for particle 1
[in]	x2	Coordinates for particle 2
[out]	dx	Distance vector

Returns: the ishift required to shift x1 at closest distance to x2; i.e. if 0<=ishift<c_numShiftVectors then x1 - x2 + shift_vec[ishift] = dx (see calc_shifts below on how to obtain shift_vec)

void pbc_dx_d	(	const t_pbc *	pbc,
		const dvec	x1,
		const dvec	x2,
		dvec	dx
	)

Compute distance vector in double precision.

Compute distance with PBC.

void put_atoms_in_box	(	PbcType	pbcType,
		const matrix	box,
		gmx::ArrayRef< gmx::RVec >	x
	)

Put atoms inside the simulations box.

These routines puts ONE or ALL atoms in the box, not caring about charge groups! Also works for triclinic cells.

Parameters

[in]	pbcType	The pbc type
[in]	box	The simulation box
[in,out]	x	The coordinates of the atoms

void put_atoms_in_box_omp	(	PbcType	pbcType,
		const matrix	box,
		bool	haveBoxDeformation,
		const matrix	boxDeformation,
		gmx::ArrayRef< gmx::RVec >	x,
		gmx::ArrayRef< gmx::RVec >	v,
		gmx_unused int	nth
	)

Parallellizes put_atoms_in_box()

This wrapper function around put_atoms_in_box() with the ugly manual workload splitting is needed to avoid silently introducing multithreading in tools.

Parameters

[in]	pbcType	The pbc type
[in]	box	The simulation box
[in]	haveBoxDeformation	Whether the box is being continously deformed
[in]	boxDeformation	The deformation speed of the box components in units of nm/ps
[in,out]	x	The coordinates of the atoms
[in,out]	v	The velocities of the atoms
[in]	nth	number of threads to be used in the given module

void put_atoms_in_compact_unitcell	(	PbcType	pbcType,
		int	ecenter,
		const matrix	box,
		gmx::ArrayRef< gmx::RVec >	x
	)

Put atoms inside the unitcell.

This puts ALL atoms at the closest distance for the center of the box as calculated by calc_box_center. When pbcType=PbcType::Unset, the type of pbc is guessed from the box matrix.

Parameters

[in]	pbcType	The pbc type
[in]	ecenter	The pbc center type
[in]	box	The simulation box
[in,out]	x	The coordinates of the atoms

void put_atoms_in_triclinic_unitcell	(	int	ecenter,
		const matrix	box,
		gmx::ArrayRef< gmx::RVec >	x
	)

Put atoms inside triclinic box.

This puts ALL atoms in the triclinic unit cell, centered around the box center as calculated by calc_box_center.

Parameters

[in]	ecenter	The pbc center type
[in]	box	The simulation box
[in,out]	x	The coordinates of the atoms

void set_pbc	(	t_pbc *	pbc,
		PbcType	pbcType,
		const matrix	box
	)

Initiate the periodic boundary condition algorithms.

pbc_dx will not use pbc and return the normal difference vector when one or more of the diagonal elements of box are zero. When pbcType=PbcType::Unset, the type of pbc is guessed from the box matrix.

Parameters

[in,out]	pbc	The pbc information structure
[in]	pbcType	The PBC identifier
[in]	box	The box tensor

t_pbc* set_pbc_dd	(	t_pbc *	pbc,
		PbcType	pbcType,
		const ivec	domdecCells,
		bool	bSingleDir,
		const matrix	box
	)

Initiate the periodic boundary condition algorithms.

As set_pbc, but additionally sets that correct distances can be obtained using (combinations of) single box-vector shifts. Should be used with pbc_dx_aiuc. If domdecCells!=NULL pbc is not used for directions with dd->nc[i]==1 with bSingleDir==TRUE or with dd->nc[i]<=2 with bSingleDir==FALSE. Note that when no PBC is required only pbc->pbcType is set, the rest of the struct will be invalid.

Parameters

[in,out]	pbc	The pbc information structure
[in]	pbcType	The PBC identifier
[in]	domdecCells	3D integer vector describing the number of DD cells or nullptr if not using DD.
[in]	bSingleDir	TRUE if DD communicates only in one direction along dimensions
[in]	box	The box tensor

Returns: the pbc structure when pbc operations are required, NULL otherwise.

void setBoxDeformationRate	(	const matrix	boxDeformation,
		const matrix	box,
		matrix	boxDeformationRate
	)

Sets the box deformation rate.

Parameters

[in]	boxDeformation	The deformation speed of the box components in units of nm/ps
[in]	box	The box
[out]	boxDeformationRate	The deformation rate of the box in units of 1/ps

Variable Documentation

const gmx::EnumerationArray<PbcType, std::string> c_pbcTypeNames

Initial value:

= {
    { "xyz", "no", "xy", "screw", "unset" }
}

Names for all values in PBC types enumeration.

constexpr real sc_boxSkewnessMarginForWarning = sc_skewnessMargin + 0.004

static

Margin factor for warning/error message.

The term 0.004 is just sufficient to not warn for a box that is deformed with a shear rate of 0.001 ps^-1 over a pairlist lifetime of 0.2 ps.

Note that there can be errors in periodic images for atom pairs close to half the box size. In nearly all cases the cut-off distance is more then 0.5% shorter than half the minimum periodic vector length, so there is no issue.

Description

Enumerations

Functions

Variables

Function Documentation

Variable Documentation