src/gromacs/pbcutil/pbc.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2012,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #ifndef GMX_PBCUTIL_PBC_H
  38 #define GMX_PBCUTIL_PBC_H
  39
  40 #include <stdio.h>
  41
  42 #include "gromacs/math/vectypes.h"
  43 #include "gromacs/utility/arrayref.h"
  44 #include "gromacs/utility/basedefinitions.h"
  45 #include "gromacs/utility/real.h"
  46
  47 struct gmx_domdec_t;
  48 struct gmx_mtop_t;
  49
  50 enum
  51 {
  52     epbcXYZ,
  53     epbcNONE,
  54     epbcXY,
  55     epbcSCREW,
  56     epbcNR
  57 };
  58
  59 //! Strings corresponding to epbc enum values.
  60 extern const char* epbc_names[epbcNR + 1];
  61
  62 /* Maximum number of combinations of single triclinic box vectors
  63  * required to shift atoms that are within a brick of the size of
  64  * the diagonal of the box to within the maximum cut-off distance.
  65  */
  66 #define MAX_NTRICVEC 12
  67
  68 /*! \brief Structure containing info on periodic boundary conditions */
  69 typedef struct t_pbc
  70 {
  71     //! The PBC type
  72     int ePBC;
  73     //! Number of dimensions in which PBC is exerted
  74     int ndim_ePBC;
  75     /*! \brief Determines how to compute distance vectors.
  76      *
  77      *  Indicator of how to compute distance vectors, depending
  78      *  on PBC type (depends on ePBC and dimensions with(out) DD)
  79      *  and the box angles.
  80      */
  81     int ePBCDX;
  82     /*! \brief Used for selecting which dimensions to use in PBC.
  83      *
  84      *  In case of 1-D PBC this indicates which dimension is used,
  85      *  in case of 2-D PBC this indicates the opposite
  86      */
  87     int dim;
  88     //! The simulation box
  89     matrix box;
  90     //! The lengths of the diagonal of the full box
  91     rvec fbox_diag;
  92     //! Halve of the above
  93     rvec hbox_diag;
  94     //! Negative of the above
  95     rvec mhbox_diag;
  96     //! Maximum allowed cutoff squared for the box and PBC used
  97     real max_cutoff2;
  98     /*! \brief Number of triclinic shift vectors.
  99      *
 100      *  Number of triclinic shift vectors depends on the skewedness
 101      *  of the box, that is mostly on the angles. For triclinic boxes
 102      *  we first take the closest image along each Cartesian dimension
 103      *  independently. When the resulting distance^2 is larger than
 104      *  max_cutoff2, up to ntric_vec triclinic shift vectors need to
 105      *  be tried. Because of the restrictions imposed on the unit-cell
 106      *  by GROMACS, ntric_vec <= MAX_NTRICVEC = 12.
 107      */
 108     int ntric_vec;
 109     //! The triclinic shift vectors in grid cells. Internal use only.
 110     ivec tric_shift[MAX_NTRICVEC];
 111     //!  The triclinic shift vectors in length units
 112     rvec tric_vec[MAX_NTRICVEC];
 113 } t_pbc;
 114
 115 #define TRICLINIC(box) ((box)[YY][XX] != 0 || (box)[ZZ][XX] != 0 || (box)[ZZ][YY] != 0)
 116
 117 #define NTRICIMG 14
 118 #define NCUCVERT 24
 119 #define NCUCEDGE 36
 120
 121 enum
 122 {
 123     ecenterTRIC, /* 0.5*(a+b+c)                  */
 124     ecenterRECT, /* (0.5*a[x],0.5*b[y],0.5*c[z]) */
 125     ecenterZERO, /* (0,0,0)                      */
 126     ecenterDEF = ecenterTRIC
 127 };
 128
 129 struct t_graph;
 130
 131 /*! \brief Returns the number of dimensions that use pbc
 132  *
 133  * \param[in] ePBC The periodic boundary condition type
 134  * \return the number of dimensions that use pbc, starting at X
 135  */
 136 int ePBC2npbcdim(int ePBC);
 137
 138 /*! \brief Dump the contents of the pbc structure to the file
 139  *
 140  * \param[in] fp  The file pointer to write to
 141  * \param[in] pbc The periodic boundary condition information structure
 142  */
 143 void dump_pbc(FILE* fp, t_pbc* pbc);
 144
 145 /*! \brief Check the box for consistency
 146  *
 147  * \param[in] ePBC The pbc identifier
 148  * \param[in] box  The box matrix
 149  * \return NULL if the box is supported by Gromacs.
 150  * Otherwise returns a string with the problem.
 151  * When ePBC=-1, the type of pbc is guessed from the box matrix.
 152  */
 153 const char* check_box(int ePBC, const matrix box);
 154
 155 /*! \brief Creates box matrix from edge lengths and angles.
 156  *
 157  * \param[in,out] box         The box matrix
 158  * \param[in] vec            The edge lengths
 159  * \param[in] angleInDegrees The angles
 160  */
 161 void matrix_convert(matrix box, const rvec vec, const rvec angleInDegrees);
 162
 163 /*! \brief Compute the maximum cutoff for the box
 164
 165  * Returns the square of the maximum cut-off allowed for the box,
 166  * taking into account that the grid neighborsearch code and pbc_dx
 167  * only check combinations of single box-vector shifts.
 168  * \param[in] ePBC The pbc identifier
 169  * \param[in] box  The box matrix
 170  * \return the maximum cut-off.
 171  */
 172 real max_cutoff2(int ePBC, const matrix box);
 173
 174 /*! \brief Guess PBC typr
 175  *
 176  * Guesses the type of periodic boundary conditions using the box
 177  * \param[in] box  The box matrix
 178  * \return The pbc identifier
 179  */
 180 int guess_ePBC(const matrix box);
 181
 182 /*! \brief Corrects the box if necessary
 183  *
 184  * Checks for un-allowed box angles and corrects the box
 185  * and the integer shift vectors in the graph (if graph!=NULL) if necessary.
 186  * \param[in] fplog File for debug output
 187  * \param[in] step  The MD step number
 188  * \param[in] box   The simulation cell
 189  * \param[in] graph Information about molecular connectivity
 190  * \return TRUE when the box was corrected.
 191  */
 192 gmx_bool correct_box(FILE* fplog, int step, tensor box, struct t_graph* graph);
 193
 194 /*! \brief Initiate the periodic boundary condition algorithms.
 195  *
 196  * pbc_dx will not use pbc and return the normal difference vector
 197  * when one or more of the diagonal elements of box are zero.
 198  * When ePBC=-1, the type of pbc is guessed from the box matrix.
 199  * \param[in,out] pbc The pbc information structure
 200  * \param[in] ePBC The PBC identifier
 201  * \param[in] box  The box tensor
 202  */
 203 void set_pbc(t_pbc* pbc, int ePBC, const matrix box);
 204
 205 /*! \brief Initiate the periodic boundary condition algorithms.
 206  *
 207  * As set_pbc, but additionally sets that correct distances can
 208  * be obtained using (combinations of) single box-vector shifts.
 209  * Should be used with pbc_dx_aiuc.
 210  * If domdecCells!=NULL pbc is not used for directions
 211  * with dd->nc[i]==1 with bSingleDir==TRUE or
 212  * with dd->nc[i]<=2 with bSingleDir==FALSE.
 213  * Note that when no PBC is required only pbc->ePBC is set,
 214  * the rest of the struct will be invalid.
 215  * \param[in,out] pbc The pbc information structure
 216  * \param[in] ePBC        The PBC identifier
 217  * \param[in] domdecCells 3D integer vector describing the number of DD cells
 218  *                        or nullptr if not using DD.
 219  * \param[in] bSingleDir  TRUE if DD communicates only in one direction along dimensions
 220  * \param[in] box         The box tensor
 221  * \return the pbc structure when pbc operations are required, NULL otherwise.
 222  */
 223 t_pbc* set_pbc_dd(t_pbc* pbc, int ePBC, const ivec domdecCells, gmx_bool bSingleDir, const matrix box);
 224
 225 /*! \brief Compute distance with PBC
 226  *
 227  * Calculate the correct distance vector from x2 to x1 and put it in dx.
 228  * set_pbc must be called before ever calling this routine.
 229  *
 230  * Note that for triclinic boxes that do not obey the GROMACS unit-cell
 231  * restrictions, pbc_dx and pbc_dx_aiuc will not correct for PBC.
 232  * \param[in,out] pbc The pbc information structure
 233  * \param[in]    x1  Coordinates for particle 1
 234  * \param[in]    x2  Coordinates for particle 2
 235  * \param[out]   dx  Distance vector
 236  */
 237 void pbc_dx(const t_pbc* pbc, const rvec x1, const rvec x2, rvec dx);
 238
 239 /*! \brief Compute distance vector for simple PBC types
 240  *
 241  * Calculate the correct distance vector from x2 to x1 and put it in dx,
 242  * This function can only be used when all atoms are in the rectangular
 243  * or triclinic unit-cell.
 244  * set_pbc_dd or set_pbc must be called before ever calling this routine.
 245  * \param[in,out] pbc The pbc information structure
 246  * \param[in]    x1  Coordinates for particle 1
 247  * \param[in]    x2  Coordinates for particle 2
 248  * \param[out]   dx  Distance vector
 249  * \return the ishift required to shift x1 at closest distance to x2;
 250  * i.e. if 0<=ishift<SHIFTS then x1 - x2 + shift_vec[ishift] = dx
 251  * (see calc_shifts below on how to obtain shift_vec)
 252  */
 253 int pbc_dx_aiuc(const t_pbc* pbc, const rvec x1, const rvec x2, rvec dx);
 254
 255 /*! \brief Compute distance with PBC
 256  *
 257  * As pbc_dx, but for double precision vectors.
 258  * set_pbc must be called before ever calling this routine.
 259  * \param[in,out] pbc The pbc information structure
 260  * \param[in]    x1  Coordinates for particle 1
 261  * \param[in]    x2  Coordinates for particle 2
 262  * \param[out]   dx  Distance vector
 263  */
 264 void pbc_dx_d(const t_pbc* pbc, const dvec x1, const dvec x2, dvec dx);
 265
 266 /*! \brief Computes shift vectors
 267  *
 268  * This routine calculates ths shift vectors necessary to use the
 269  * neighbor searching routine.
 270  * \param[in]  box       The simulation box
 271  * \param[out] shift_vec The shifting vectors
 272  */
 273 void calc_shifts(const matrix box, rvec shift_vec[]);
 274
 275 /*! \brief Calculates the center of the box.
 276  *
 277  * See the description for the enum ecenter above.
 278  * \param[in]  ecenter    Description of center type
 279  * \param[in]  box        The simulation box
 280  * \param[out] box_center The center of the box
 281  */
 282 void calc_box_center(int ecenter, const matrix box, rvec box_center);
 283
 284 /*! \brief Calculates the NTRICIMG box images
 285  *
 286  * \param[in]  box The simulation box
 287  * \param[out] img The triclinic box images
 288  */
 289 void calc_triclinic_images(const matrix box, rvec img[]);
 290
 291 /*! \brief Calculates the NCUCVERT vertices of a compact unitcell
 292  *
 293  * \param[in]  ecenter The center type
 294  * \param[in]  box     The simulation box
 295  * \param[out] vert    The vertices
 296  */
 297 void calc_compact_unitcell_vertices(int ecenter, const matrix box, rvec vert[]);
 298
 299 /*! \brief Compute unitcell edges
 300  *
 301  * \return an array of unitcell edges of length NCUCEDGE*2,
 302  * this is an index in vert[], which is calculated by calc_unitcell_vertices.
 303  * The index consists of NCUCEDGE pairs of vertex indices.
 304  * The index does not change, so it needs to be retrieved only once.
 305  */
 306 int* compact_unitcell_edges();
 307
 308 /*! \brief Put atoms inside the simulations box
 309  *
 310  * These routines puts ONE or ALL atoms in the box, not caring
 311  * about charge groups!
 312  * Also works for triclinic cells.
 313  * \param[in]     ePBC   The pbc type
 314  * \param[in]     box    The simulation box
 315  * \param[in,out] x      The coordinates of the atoms
 316  */
 317 void put_atoms_in_box(int ePBC, const matrix box, gmx::ArrayRef<gmx::RVec> x);
 318
 319 /*! \brief Parallellizes put_atoms_in_box()
 320  *
 321  * This wrapper function around put_atoms_in_box() with the ugly manual
 322  * workload splitting is needed to avoid silently introducing multithreading
 323  * in tools.
 324  * \param[in]     ePBC       The pbc type
 325  * \param[in]     box        The simulation box
 326  * \param[in,out] x          The coordinates of the atoms
 327  * \param[in]     nth        number of threads to be used in the given module
 328  */
 329 void put_atoms_in_box_omp(int ePBC, const matrix box, gmx::ArrayRef<gmx::RVec> x, gmx_unused int nth);
 330
 331 /*! \brief Put atoms inside triclinic box
 332  *
 333  * This puts ALL atoms in the triclinic unit cell, centered around the
 334  * box center as calculated by calc_box_center.
 335  * \param[in]    ecenter The pbc center type
 336  * \param[in]    box     The simulation box
 337  * \param[in,out] x       The coordinates of the atoms
 338  */
 339 void put_atoms_in_triclinic_unitcell(int ecenter, const matrix box, gmx::ArrayRef<gmx::RVec> x);
 340
 341 /*! \brief Put atoms inside the unitcell
 342  *
 343  * This puts ALL atoms at the closest distance for the center of the box
 344  * as calculated by calc_box_center.
 345  * When ePBC=-1, the type of pbc is guessed from the box matrix.
 346  * \param[in]    ePBC    The pbc type
 347  * \param[in]    ecenter The pbc center type
 348  * \param[in]    box     The simulation box
 349  * \param[in,out] x       The coordinates of the atoms
 350  */
 351 void put_atoms_in_compact_unitcell(int ePBC, int ecenter, const matrix box, gmx::ArrayRef<gmx::RVec> x);
 352
 353 /*! \brief Make all molecules whole by shifting positions
 354  *
 355  * \param[in]     fplog     Log file
 356  * \param[in]     ePBC      The PBC type
 357  * \param[in]     box       The simulation box
 358  * \param[in]     mtop      System topology definition
 359  * \param[in,out] x         The coordinates of the atoms
 360  */
 361 void do_pbc_first_mtop(FILE* fplog, int ePBC, const matrix box, const gmx_mtop_t* mtop, rvec x[]);
 362
 363 /*! \brief Make molecules consisting of multiple charge groups whole by shifting positions
 364  *
 365  * \param[in]     ePBC      The PBC type
 366  * \param[in]     box       The simulation box
 367  * \param[in]     mtop      System topology definition
 368  * \param[in,out] x         The coordinates of the atoms
 369  */
 370 void do_pbc_mtop(int ePBC, const matrix box, const gmx_mtop_t* mtop, rvec x[]);
 371
 372 #endif