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, 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
  49 enum {
  50     epbcXYZ, epbcNONE, epbcXY, epbcSCREW, epbcNR
  51 };
  52
  53 //! Strings corresponding to epbc enum values.
  54 extern const char *epbc_names[epbcNR+1];
  55
  56 /* Maximum number of combinations of single triclinic box vectors
  57  * required to shift atoms that are within a brick of the size of
  58  * the diagonal of the box to within the maximum cut-off distance.
  59  */
  60 #define MAX_NTRICVEC 12
  61
  62 /*! \brief Structure containing info on periodic boundary conditions */
  63 typedef struct t_pbc {
  64     //! The PBC type
  65     int        ePBC;
  66     //! Number of dimensions in which PBC is exerted
  67     int        ndim_ePBC;
  68     /*! \brief Determines how to compute distance vectors.
  69      *
  70      *  Indicator of how to compute distance vectors, depending
  71      *  on PBC type (depends on ePBC and dimensions with(out) DD)
  72      *  and the box angles.
  73      */
  74     int        ePBCDX;
  75     /*! \brief Used for selecting which dimensions to use in PBC.
  76      *
  77      *  In case of 1-D PBC this indicates which dimension is used,
  78      *  in case of 2-D PBC this indicates the opposite
  79      */
  80     int        dim;
  81     //! The simulation box
  82     matrix     box;
  83     //! The lengths of the diagonal of the full box
  84     rvec       fbox_diag;
  85     //! Halve of the above
  86     rvec       hbox_diag;
  87     //! Negative of the above
  88     rvec       mhbox_diag;
  89     //! Maximum allowed cutoff squared for the box and PBC used
  90     real       max_cutoff2;
  91     /*! \brief Number of triclinic shift vectors.
  92      *
  93      *  Number of triclinic shift vectors depends on the skewedness
  94      *  of the box, that is mostly on the angles. For triclinic boxes
  95      *  we first take the closest image along each Cartesian dimension
  96      *  independently. When the resulting distance^2 is larger than
  97      *  max_cutoff2, up to ntric_vec triclinic shift vectors need to
  98      *  be tried. Because of the restrictions imposed on the unit-cell
  99      *  by GROMACS, ntric_vec <= MAX_NTRICVEC = 12.
 100      */
 101     int        ntric_vec;
 102     //! The triclinic shift vectors in grid cells. Internal use only.
 103     ivec       tric_shift[MAX_NTRICVEC];
 104     //!  The triclinic shift vectors in length units
 105     rvec       tric_vec[MAX_NTRICVEC];
 106 } t_pbc;
 107
 108 #define TRICLINIC(box) ((box)[YY][XX] != 0 || (box)[ZZ][XX] != 0 || (box)[ZZ][YY] != 0)
 109
 110 #define NTRICIMG 14
 111 #define NCUCVERT 24
 112 #define NCUCEDGE 36
 113
 114 enum {
 115     ecenterTRIC, /* 0.5*(a+b+c)                  */
 116     ecenterRECT, /* (0.5*a[x],0.5*b[y],0.5*c[z]) */
 117     ecenterZERO, /* (0,0,0)                      */
 118     ecenterDEF = ecenterTRIC
 119 };
 120
 121 struct t_graph;
 122
 123 /*! \brief Returns the number of dimensions that use pbc
 124  *
 125  * \param[in] ePBC The periodic boundary condition type
 126  * \return the number of dimensions that use pbc, starting at X
 127  */
 128 int ePBC2npbcdim(int ePBC);
 129
 130 /*! \brief Dump the contents of the pbc structure to the file
 131  *
 132  * \param[in] fp  The file pointer to write to
 133  * \param[in] pbc The periodic boundary condition information structure
 134  */
 135 void dump_pbc(FILE *fp, t_pbc *pbc);
 136
 137 /*! \brief Check the box for consistency
 138  *
 139  * \param[in] ePBC The pbc identifier
 140  * \param[in] box  The box matrix
 141  * \return NULL if the box is supported by Gromacs.
 142  * Otherwise returns a string with the problem.
 143  * When ePBC=-1, the type of pbc is guessed from the box matrix.
 144  */
 145 const char *check_box(int ePBC, const matrix box);
 146
 147 /*! \brief Creates box matrix from edge lengths and angles.
 148  *
 149  * \param[in,out] box         The box matrix
 150  * \param[in] vec            The edge lengths
 151  * \param[in] angleInDegrees The angles
 152  */
 153 void matrix_convert(matrix box, const rvec vec, const rvec angleInDegrees);
 154
 155 /*! \brief Compute the maximum cutoff for the box
 156
 157  * Returns the square of the maximum cut-off allowed for the box,
 158  * taking into account that the grid neighborsearch code and pbc_dx
 159  * only check combinations of single box-vector shifts.
 160  * \param[in] ePBC The pbc identifier
 161  * \param[in] box  The box matrix
 162  * \return the maximum cut-off.
 163  */
 164 real max_cutoff2(int ePBC, const matrix box);
 165
 166 /*! \brief Guess PBC typr
 167  *
 168  * Guesses the type of periodic boundary conditions using the box
 169  * \param[in] box  The box matrix
 170  * \return The pbc identifier
 171  */
 172 int guess_ePBC(const matrix box);
 173
 174 /*! \brief Corrects the box if necessary
 175  *
 176  * Checks for un-allowed box angles and corrects the box
 177  * and the integer shift vectors in the graph (if graph!=NULL) if necessary.
 178  * \param[in] fplog File for debug output
 179  * \param[in] step  The MD step number
 180  * \param[in] box   The simulation cell
 181  * \param[in] graph Information about molecular connectivity
 182  * \return TRUE when the box was corrected.
 183  */
 184 gmx_bool correct_box(FILE *fplog, int step, tensor box, struct t_graph *graph);
 185
 186 /*! \brief Initiate the periodic boundary condition algorithms.
 187  *
 188  * pbc_dx will not use pbc and return the normal difference vector
 189  * when one or more of the diagonal elements of box are zero.
 190  * When ePBC=-1, the type of pbc is guessed from the box matrix.
 191  * \param[in,out] pbc The pbc information structure
 192  * \param[in] ePBC The PBC identifier
 193  * \param[in] box  The box tensor
 194  */
 195 void set_pbc(t_pbc *pbc, int ePBC, const matrix box);
 196
 197 /*! \brief Initiate the periodic boundary condition algorithms.
 198  *
 199  * As set_pbc, but additionally sets that correct distances can
 200  * be obtained using (combinations of) single box-vector shifts.
 201  * Should be used with pbc_dx_aiuc.
 202  * If domdecCells!=NULL pbc is not used for directions
 203  * with dd->nc[i]==1 with bSingleDir==TRUE or
 204  * with dd->nc[i]<=2 with bSingleDir==FALSE.
 205  * Note that when no PBC is required only pbc->ePBC is set,
 206  * the rest of the struct will be invalid.
 207  * \param[in,out] pbc The pbc information structure
 208  * \param[in] ePBC        The PBC identifier
 209  * \param[in] domdecCells 3D integer vector describing the number of DD cells
 210  *                        or nullptr if not using DD.
 211  * \param[in] bSingleDir  TRUE if DD communicates only in one direction along dimensions
 212  * \param[in] box         The box tensor
 213  * \return the pbc structure when pbc operations are required, NULL otherwise.
 214  */
 215 t_pbc *set_pbc_dd(t_pbc *pbc, int ePBC,
 216                   const ivec domdecCells, gmx_bool bSingleDir,
 217                   const matrix box);
 218
 219 /*! \brief Compute distance with PBC
 220  *
 221  * Calculate the correct distance vector from x2 to x1 and put it in dx.
 222  * set_pbc must be called before ever calling this routine.
 223  *
 224  * Note that for triclinic boxes that do not obey the GROMACS unit-cell
 225  * restrictions, pbc_dx and pbc_dx_aiuc will not correct for PBC.
 226  * \param[in,out] pbc The pbc information structure
 227  * \param[in]    x1  Coordinates for particle 1
 228  * \param[in]    x2  Coordinates for particle 2
 229  * \param[out]   dx  Distance vector
 230  */
 231 void pbc_dx(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx);
 232
 233 /*! \brief Compute distance vector for simple PBC types
 234  *
 235  * Calculate the correct distance vector from x2 to x1 and put it in dx,
 236  * This function can only be used when all atoms are in the rectangular
 237  * or triclinic unit-cell.
 238  * set_pbc_dd or set_pbc must be called before ever calling this routine.
 239  * \param[in,out] pbc The pbc information structure
 240  * \param[in]    x1  Coordinates for particle 1
 241  * \param[in]    x2  Coordinates for particle 2
 242  * \param[out]   dx  Distance vector
 243  * \return the ishift required to shift x1 at closest distance to x2;
 244  * i.e. if 0<=ishift<SHIFTS then x1 - x2 + shift_vec[ishift] = dx
 245  * (see calc_shifts below on how to obtain shift_vec)
 246  */
 247 int pbc_dx_aiuc(const t_pbc *pbc, const rvec x1, const rvec x2, rvec dx);
 248
 249 /*! \brief Compute distance with PBC
 250  *
 251  * As pbc_dx, but for double precision vectors.
 252  * set_pbc must be called before ever calling this routine.
 253  * \param[in,out] pbc The pbc information structure
 254  * \param[in]    x1  Coordinates for particle 1
 255  * \param[in]    x2  Coordinates for particle 2
 256  * \param[out]   dx  Distance vector
 257  */
 258 void pbc_dx_d(const t_pbc *pbc, const dvec x1, const dvec x2, dvec dx);
 259
 260 /*! \brief Computes shift vectors
 261  *
 262  * This routine calculates ths shift vectors necessary to use the
 263  * neighbor searching routine.
 264  * \param[in]  box       The simulation box
 265  * \param[out] shift_vec The shifting vectors
 266  */
 267 void calc_shifts(const matrix box, rvec shift_vec[]);
 268
 269 /*! \brief Calculates the center of the box.
 270  *
 271  * See the description for the enum ecenter above.
 272  * \param[in]  ecenter    Description of center type
 273  * \param[in]  box        The simulation box
 274  * \param[out] box_center The center of the box
 275  */
 276 void calc_box_center(int ecenter, const matrix box, rvec box_center);
 277
 278 /*! \brief Calculates the NTRICIMG box images
 279  *
 280  * \param[in]  box The simulation box
 281  * \param[out] img The triclinic box images
 282  */
 283 void calc_triclinic_images(const matrix box, rvec img[]);
 284
 285 /*! \brief Calculates the NCUCVERT vertices of a compact unitcell
 286  *
 287  * \param[in]  ecenter The center type
 288  * \param[in]  box     The simulation box
 289  * \param[out] vert    The vertices
 290  */
 291 void calc_compact_unitcell_vertices(int ecenter, const matrix box,
 292                                     rvec vert[]);
 293
 294 /*! \brief Compute unitcell edges
 295  *
 296  * \return an array of unitcell edges of length NCUCEDGE*2,
 297  * this is an index in vert[], which is calculated by calc_unitcell_vertices.
 298  * The index consists of NCUCEDGE pairs of vertex indices.
 299  * The index does not change, so it needs to be retrieved only once.
 300  */
 301 int *compact_unitcell_edges(void);
 302
 303 /*! \brief Put atoms inside the simulations box
 304  *
 305  * These routines puts ONE or ALL atoms in the box, not caring
 306  * about charge groups!
 307  * Also works for triclinic cells.
 308  * \param[in]    ePBC   The pbc type
 309  * \param[in]    box    The simulation box
 310  * \param[in,out] x      The coordinates of the atoms
 311  */
 312 void put_atoms_in_box(int ePBC, const matrix box, gmx::ArrayRef<gmx::RVec> x);
 313
 314 /*! \brief Put atoms inside triclinic box
 315  *
 316  * This puts ALL atoms in the triclinic unit cell, centered around the
 317  * box center as calculated by calc_box_center.
 318  * \param[in]    ecenter The pbc center type
 319  * \param[in]    box     The simulation box
 320  * \param[in,out] x       The coordinates of the atoms
 321  */
 322 void put_atoms_in_triclinic_unitcell(int ecenter, const matrix box,
 323                                      gmx::ArrayRef<gmx::RVec> x);
 324
 325 /*! \brief Put atoms inside the unitcell
 326  *
 327  * This puts ALL atoms at the closest distance for the center of the box
 328  * as calculated by calc_box_center.
 329  * When ePBC=-1, the type of pbc is guessed from the box matrix.
 330  * \param[in]    ePBC    The pbc type
 331  * \param[in]    ecenter The pbc center type
 332  * \param[in]    box     The simulation box
 333  * \param[in,out] x       The coordinates of the atoms
 334  */
 335 void put_atoms_in_compact_unitcell(int ePBC, int ecenter,
 336                                    const matrix box,
 337                                    gmx::ArrayRef<gmx::RVec> x);
 338
 339 #endif