src/gromacs/pbcutil/pbc.cpp

   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) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /*! \internal \file
  39  * \brief
  40  * Implements routines in pbc.h.
  41  *
  42  * Utility functions for handling periodic boundary conditions.
  43  * Mainly used in analysis tools.
  44  */
  45 #include "gmxpre.h"
  46
  47 #include "pbc.h"
  48
  49 #include <cmath>
  50
  51 #include <algorithm>
  52
  53 #include "gromacs/math/functions.h"
  54 #include "gromacs/math/units.h"
  55 #include "gromacs/math/utilities.h"
  56 #include "gromacs/math/vec.h"
  57 #include "gromacs/math/vecdump.h"
  58 #include "gromacs/pbcutil/ishift.h"
  59 #include "gromacs/pbcutil/mshift.h"
  60 #include "gromacs/topology/topology.h"
  61 #include "gromacs/utility/exceptions.h"
  62 #include "gromacs/utility/fatalerror.h"
  63 #include "gromacs/utility/gmxassert.h"
  64 #include "gromacs/utility/smalloc.h"
  65
  66 const gmx::EnumerationArray<PbcType, std::string> c_pbcTypeNames = { { "xyz", "no", "xy", "screw",
  67                                                                        "unset" } };
  68
  69 /* Skip 0 so we have more chance of detecting if we forgot to call set_pbc. */
  70 enum
  71 {
  72     epbcdxRECTANGULAR = 1,
  73     epbcdxTRICLINIC,
  74     epbcdx2D_RECT,
  75     epbcdx2D_TRIC,
  76     epbcdx1D_RECT,
  77     epbcdx1D_TRIC,
  78     epbcdxSCREW_RECT,
  79     epbcdxSCREW_TRIC,
  80     epbcdxNOPBC,
  81     epbcdxUNSUPPORTED
  82 };
  83
  84 //! Margin factor for error message
  85 #define BOX_MARGIN 1.0010
  86 //! Margin correction if the box is too skewed
  87 #define BOX_MARGIN_CORRECT 1.0005
  88
  89 int numPbcDimensions(PbcType pbcType)
  90 {
  91     int npbcdim = 0;
  92
  93     switch (pbcType)
  94     {
  95         case PbcType::Unset:
  96             gmx_fatal(FARGS, "Number of PBC dimensions was requested before the PBC type set.");
  97         case PbcType::Xyz: npbcdim = 3; break;
  98         case PbcType::XY: npbcdim = 2; break;
  99         case PbcType::Screw: npbcdim = 3; break;
 100         case PbcType::No: npbcdim = 0; break;
 101         default:
 102             gmx_fatal(FARGS, "Unknown pbcType=%s in numPbcDimensions", c_pbcTypeNames[pbcType].c_str());
 103     }
 104
 105     return npbcdim;
 106 }
 107
 108 void dump_pbc(FILE* fp, t_pbc* pbc)
 109 {
 110     rvec sum_box;
 111
 112     fprintf(fp, "pbcTypeDX = %d\n", pbc->pbcTypeDX);
 113     pr_rvecs(fp, 0, "box", pbc->box, DIM);
 114     pr_rvecs(fp, 0, "fbox_diag", &pbc->fbox_diag, 1);
 115     pr_rvecs(fp, 0, "hbox_diag", &pbc->hbox_diag, 1);
 116     pr_rvecs(fp, 0, "mhbox_diag", &pbc->mhbox_diag, 1);
 117     rvec_add(pbc->hbox_diag, pbc->mhbox_diag, sum_box);
 118     pr_rvecs(fp, 0, "sum of the above two", &sum_box, 1);
 119     fprintf(fp, "max_cutoff2 = %g\n", pbc->max_cutoff2);
 120     fprintf(fp, "ntric_vec = %d\n", pbc->ntric_vec);
 121     if (pbc->ntric_vec > 0)
 122     {
 123         pr_ivecs(fp, 0, "tric_shift", pbc->tric_shift, pbc->ntric_vec, FALSE);
 124         pr_rvecs(fp, 0, "tric_vec", pbc->tric_vec, pbc->ntric_vec);
 125     }
 126 }
 127
 128 const char* check_box(PbcType pbcType, const matrix box)
 129 {
 130     const char* ptr;
 131
 132     if (pbcType == PbcType::Unset)
 133     {
 134         pbcType = guessPbcType(box);
 135     }
 136
 137     if (pbcType == PbcType::No)
 138     {
 139         return nullptr;
 140     }
 141
 142     if ((box[XX][YY] != 0) || (box[XX][ZZ] != 0) || (box[YY][ZZ] != 0))
 143     {
 144         ptr = "Only triclinic boxes with the first vector parallel to the x-axis and the second "
 145               "vector in the xy-plane are supported.";
 146     }
 147     else if (pbcType == PbcType::Screw && (box[YY][XX] != 0 || box[ZZ][XX] != 0))
 148     {
 149         ptr = "The unit cell can not have off-diagonal x-components with screw pbc";
 150     }
 151     else if (std::fabs(box[YY][XX]) > BOX_MARGIN * 0.5 * box[XX][XX]
 152              || (pbcType != PbcType::XY
 153                  && (std::fabs(box[ZZ][XX]) > BOX_MARGIN * 0.5 * box[XX][XX]
 154                      || std::fabs(box[ZZ][YY]) > BOX_MARGIN * 0.5 * box[YY][YY])))
 155     {
 156         ptr = "Triclinic box is too skewed.";
 157     }
 158     else
 159     {
 160         ptr = nullptr;
 161     }
 162
 163     return ptr;
 164 }
 165
 166 void matrix_convert(matrix box, const rvec vec, const rvec angleInDegrees)
 167 {
 168     rvec angle;
 169     svmul(DEG2RAD, angleInDegrees, angle);
 170     box[XX][XX] = vec[XX];
 171     box[YY][XX] = vec[YY] * cos(angle[ZZ]);
 172     box[YY][YY] = vec[YY] * sin(angle[ZZ]);
 173     box[ZZ][XX] = vec[ZZ] * cos(angle[YY]);
 174     box[ZZ][YY] = vec[ZZ] * (cos(angle[XX]) - cos(angle[YY]) * cos(angle[ZZ])) / sin(angle[ZZ]);
 175     box[ZZ][ZZ] =
 176             std::sqrt(gmx::square(vec[ZZ]) - box[ZZ][XX] * box[ZZ][XX] - box[ZZ][YY] * box[ZZ][YY]);
 177 }
 178
 179 real max_cutoff2(PbcType pbcType, const matrix box)
 180 {
 181     real       min_hv2, min_ss;
 182     const real oneFourth = 0.25;
 183
 184     /* Physical limitation of the cut-off
 185      * by half the length of the shortest box vector.
 186      */
 187     min_hv2 = oneFourth * std::min(norm2(box[XX]), norm2(box[YY]));
 188     if (pbcType != PbcType::XY)
 189     {
 190         min_hv2 = std::min(min_hv2, oneFourth * norm2(box[ZZ]));
 191     }
 192
 193     /* Limitation to the smallest diagonal element due to optimizations:
 194      * checking only linear combinations of single box-vectors (2 in x)
 195      * in the grid search and pbc_dx is a lot faster
 196      * than checking all possible combinations.
 197      */
 198     if (pbcType == PbcType::XY)
 199     {
 200         min_ss = std::min(box[XX][XX], box[YY][YY]);
 201     }
 202     else
 203     {
 204         min_ss = std::min(box[XX][XX], std::min(box[YY][YY] - std::fabs(box[ZZ][YY]), box[ZZ][ZZ]));
 205     }
 206
 207     return std::min(min_hv2, min_ss * min_ss);
 208 }
 209
 210 //! Set to true if warning has been printed
 211 static gmx_bool bWarnedGuess = FALSE;
 212
 213 PbcType guessPbcType(const matrix box)
 214 {
 215     PbcType pbcType;
 216
 217     if (box[XX][XX] > 0 && box[YY][YY] > 0 && box[ZZ][ZZ] > 0)
 218     {
 219         pbcType = PbcType::Xyz;
 220     }
 221     else if (box[XX][XX] > 0 && box[YY][YY] > 0 && box[ZZ][ZZ] == 0)
 222     {
 223         pbcType = PbcType::XY;
 224     }
 225     else if (box[XX][XX] == 0 && box[YY][YY] == 0 && box[ZZ][ZZ] == 0)
 226     {
 227         pbcType = PbcType::No;
 228     }
 229     else
 230     {
 231         if (!bWarnedGuess)
 232         {
 233             fprintf(stderr,
 234                     "WARNING: Unsupported box diagonal %f %f %f, "
 235                     "will not use periodic boundary conditions\n\n",
 236                     box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
 237             bWarnedGuess = TRUE;
 238         }
 239         pbcType = PbcType::No;
 240     }
 241
 242     if (debug)
 243     {
 244         fprintf(debug, "Guessed pbc = %s from the box matrix\n", c_pbcTypeNames[pbcType].c_str());
 245     }
 246
 247     return pbcType;
 248 }
 249
 250 //! Check if the box still obeys the restrictions, if not, correct it
 251 static int correct_box_elem(FILE* fplog, int step, tensor box, int v, int d)
 252 {
 253     int shift, maxshift = 10;
 254
 255     shift = 0;
 256
 257     /* correct elem d of vector v with vector d */
 258     while (box[v][d] > BOX_MARGIN_CORRECT * 0.5 * box[d][d])
 259     {
 260         if (fplog)
 261         {
 262             fprintf(fplog, "Step %d: correcting invalid box:\n", step);
 263             pr_rvecs(fplog, 0, "old box", box, DIM);
 264         }
 265         rvec_dec(box[v], box[d]);
 266         shift--;
 267         if (fplog)
 268         {
 269             pr_rvecs(fplog, 0, "new box", box, DIM);
 270         }
 271         if (shift <= -maxshift)
 272         {
 273             gmx_fatal(FARGS, "Box was shifted at least %d times. Please see log-file.", maxshift);
 274         }
 275     }
 276     while (box[v][d] < -BOX_MARGIN_CORRECT * 0.5 * box[d][d])
 277     {
 278         if (fplog)
 279         {
 280             fprintf(fplog, "Step %d: correcting invalid box:\n", step);
 281             pr_rvecs(fplog, 0, "old box", box, DIM);
 282         }
 283         rvec_inc(box[v], box[d]);
 284         shift++;
 285         if (fplog)
 286         {
 287             pr_rvecs(fplog, 0, "new box", box, DIM);
 288         }
 289         if (shift >= maxshift)
 290         {
 291             gmx_fatal(FARGS, "Box was shifted at least %d times. Please see log-file.", maxshift);
 292         }
 293     }
 294
 295     return shift;
 296 }
 297
 298 gmx_bool correct_box(FILE* fplog, int step, tensor box)
 299 {
 300     int      zy, zx, yx;
 301     gmx_bool bCorrected;
 302
 303     zy = correct_box_elem(fplog, step, box, ZZ, YY);
 304     zx = correct_box_elem(fplog, step, box, ZZ, XX);
 305     yx = correct_box_elem(fplog, step, box, YY, XX);
 306
 307     bCorrected = ((zy != 0) || (zx != 0) || (yx != 0));
 308
 309     return bCorrected;
 310 }
 311
 312 //! Do the real arithmetic for filling the pbc struct
 313 static void low_set_pbc(t_pbc* pbc, PbcType pbcType, const ivec dd_pbc, const matrix box)
 314 {
 315     int         order[3] = { 0, -1, 1 };
 316     ivec        bPBC;
 317     const char* ptr;
 318
 319     pbc->pbcType   = pbcType;
 320     pbc->ndim_ePBC = numPbcDimensions(pbcType);
 321
 322     if (pbc->pbcType == PbcType::No)
 323     {
 324         pbc->pbcTypeDX = epbcdxNOPBC;
 325
 326         return;
 327     }
 328
 329     copy_mat(box, pbc->box);
 330     pbc->max_cutoff2 = 0;
 331     pbc->dim         = -1;
 332     pbc->ntric_vec   = 0;
 333
 334     for (int i = 0; (i < DIM); i++)
 335     {
 336         pbc->fbox_diag[i]  = box[i][i];
 337         pbc->hbox_diag[i]  = pbc->fbox_diag[i] * 0.5;
 338         pbc->mhbox_diag[i] = -pbc->hbox_diag[i];
 339     }
 340
 341     ptr = check_box(pbcType, box);
 342     if (ptr)
 343     {
 344         fprintf(stderr, "Warning: %s\n", ptr);
 345         pr_rvecs(stderr, 0, "         Box", box, DIM);
 346         fprintf(stderr, "         Can not fix pbc.\n\n");
 347         pbc->pbcTypeDX = epbcdxUNSUPPORTED;
 348     }
 349     else
 350     {
 351         if (pbcType == PbcType::Screw && nullptr != dd_pbc)
 352         {
 353             /* This combinated should never appear here */
 354             gmx_incons("low_set_pbc called with screw pbc and dd_nc != NULL");
 355         }
 356
 357         int npbcdim = 0;
 358         for (int i = 0; i < DIM; i++)
 359         {
 360             if ((dd_pbc && dd_pbc[i] == 0) || (pbcType == PbcType::XY && i == ZZ))
 361             {
 362                 bPBC[i] = 0;
 363             }
 364             else
 365             {
 366                 bPBC[i] = 1;
 367                 npbcdim++;
 368             }
 369         }
 370         switch (npbcdim)
 371         {
 372             case 1:
 373                 /* 1D pbc is not an mdp option and it is therefore only used
 374                  * with single shifts.
 375                  */
 376                 pbc->pbcTypeDX = epbcdx1D_RECT;
 377                 for (int i = 0; i < DIM; i++)
 378                 {
 379                     if (bPBC[i])
 380                     {
 381                         pbc->dim = i;
 382                     }
 383                 }
 384                 GMX_ASSERT(pbc->dim < DIM, "Dimension for PBC incorrect");
 385                 for (int i = 0; i < pbc->dim; i++)
 386                 {
 387                     if (pbc->box[pbc->dim][i] != 0)
 388                     {
 389                         pbc->pbcTypeDX = epbcdx1D_TRIC;
 390                     }
 391                 }
 392                 break;
 393             case 2:
 394                 pbc->pbcTypeDX = epbcdx2D_RECT;
 395                 for (int i = 0; i < DIM; i++)
 396                 {
 397                     if (!bPBC[i])
 398                     {
 399                         pbc->dim = i;
 400                     }
 401                 }
 402                 for (int i = 0; i < DIM; i++)
 403                 {
 404                     if (bPBC[i])
 405                     {
 406                         for (int j = 0; j < i; j++)
 407                         {
 408                             if (pbc->box[i][j] != 0)
 409                             {
 410                                 pbc->pbcTypeDX = epbcdx2D_TRIC;
 411                             }
 412                         }
 413                     }
 414                 }
 415                 break;
 416             case 3:
 417                 if (pbcType != PbcType::Screw)
 418                 {
 419                     if (TRICLINIC(box))
 420                     {
 421                         pbc->pbcTypeDX = epbcdxTRICLINIC;
 422                     }
 423                     else
 424                     {
 425                         pbc->pbcTypeDX = epbcdxRECTANGULAR;
 426                     }
 427                 }
 428                 else
 429                 {
 430                     pbc->pbcTypeDX = (box[ZZ][YY] == 0 ? epbcdxSCREW_RECT : epbcdxSCREW_TRIC);
 431                     if (pbc->pbcTypeDX == epbcdxSCREW_TRIC)
 432                     {
 433                         fprintf(stderr,
 434                                 "Screw pbc is not yet implemented for triclinic boxes.\n"
 435                                 "Can not fix pbc.\n");
 436                         pbc->pbcTypeDX = epbcdxUNSUPPORTED;
 437                     }
 438                 }
 439                 break;
 440             default: gmx_fatal(FARGS, "Incorrect number of pbc dimensions with DD: %d", npbcdim);
 441         }
 442         pbc->max_cutoff2 = max_cutoff2(pbcType, box);
 443
 444         if (pbc->pbcTypeDX == epbcdxTRICLINIC || pbc->pbcTypeDX == epbcdx2D_TRIC
 445             || pbc->pbcTypeDX == epbcdxSCREW_TRIC)
 446         {
 447             if (debug)
 448             {
 449                 pr_rvecs(debug, 0, "Box", box, DIM);
 450                 fprintf(debug, "max cutoff %.3f\n", sqrt(pbc->max_cutoff2));
 451             }
 452             /* We will only need single shifts here */
 453             for (int kk = 0; kk < 3; kk++)
 454             {
 455                 int k = order[kk];
 456                 if (!bPBC[ZZ] && k != 0)
 457                 {
 458                     continue;
 459                 }
 460                 for (int jj = 0; jj < 3; jj++)
 461                 {
 462                     int j = order[jj];
 463                     if (!bPBC[YY] && j != 0)
 464                     {
 465                         continue;
 466                     }
 467                     for (int ii = 0; ii < 3; ii++)
 468                     {
 469                         int i = order[ii];
 470                         if (!bPBC[XX] && i != 0)
 471                         {
 472                             continue;
 473                         }
 474                         /* A shift is only useful when it is trilinic */
 475                         if (j != 0 || k != 0)
 476                         {
 477                             rvec trial;
 478                             rvec pos;
 479                             real d2old = 0;
 480                             real d2new = 0;
 481
 482                             for (int d = 0; d < DIM; d++)
 483                             {
 484                                 trial[d] = i * box[XX][d] + j * box[YY][d] + k * box[ZZ][d];
 485                                 /* Choose the vector within the brick around 0,0,0 that
 486                                  * will become the shortest due to shift try.
 487                                  */
 488                                 if (d == pbc->dim)
 489                                 {
 490                                     trial[d] = 0;
 491                                     pos[d]   = 0;
 492                                 }
 493                                 else
 494                                 {
 495                                     if (trial[d] < 0)
 496                                     {
 497                                         pos[d] = std::min(pbc->hbox_diag[d], -trial[d]);
 498                                     }
 499                                     else
 500                                     {
 501                                         pos[d] = std::max(-pbc->hbox_diag[d], -trial[d]);
 502                                     }
 503                                 }
 504                                 d2old += gmx::square(pos[d]);
 505                                 d2new += gmx::square(pos[d] + trial[d]);
 506                             }
 507                             if (BOX_MARGIN * d2new < d2old)
 508                             {
 509                                 /* Check if shifts with one box vector less do better */
 510                                 gmx_bool bUse = TRUE;
 511                                 for (int dd = 0; dd < DIM; dd++)
 512                                 {
 513                                     int shift = (dd == 0 ? i : (dd == 1 ? j : k));
 514                                     if (shift)
 515                                     {
 516                                         real d2new_c = 0;
 517                                         for (int d = 0; d < DIM; d++)
 518                                         {
 519                                             d2new_c += gmx::square(pos[d] + trial[d] - shift * box[dd][d]);
 520                                         }
 521                                         if (d2new_c <= BOX_MARGIN * d2new)
 522                                         {
 523                                             bUse = FALSE;
 524                                         }
 525                                     }
 526                                 }
 527                                 if (bUse)
 528                                 {
 529                                     /* Accept this shift vector. */
 530                                     if (pbc->ntric_vec >= MAX_NTRICVEC)
 531                                     {
 532                                         fprintf(stderr,
 533                                                 "\nWARNING: Found more than %d triclinic "
 534                                                 "correction vectors, ignoring some.\n"
 535                                                 "  There is probably something wrong with your "
 536                                                 "box.\n",
 537                                                 MAX_NTRICVEC);
 538                                         pr_rvecs(stderr, 0, "         Box", box, DIM);
 539                                     }
 540                                     else
 541                                     {
 542                                         copy_rvec(trial, pbc->tric_vec[pbc->ntric_vec]);
 543                                         pbc->tric_shift[pbc->ntric_vec][XX] = i;
 544                                         pbc->tric_shift[pbc->ntric_vec][YY] = j;
 545                                         pbc->tric_shift[pbc->ntric_vec][ZZ] = k;
 546                                         pbc->ntric_vec++;
 547
 548                                         if (debug)
 549                                         {
 550                                             fprintf(debug,
 551                                                     "  tricvec %2d = %2d %2d %2d  %5.2f %5.2f  "
 552                                                     "%5.2f %5.2f %5.2f  %5.2f %5.2f %5.2f\n",
 553                                                     pbc->ntric_vec, i, j, k, sqrt(d2old),
 554                                                     sqrt(d2new), trial[XX], trial[YY], trial[ZZ],
 555                                                     pos[XX], pos[YY], pos[ZZ]);
 556                                         }
 557                                     }
 558                                 }
 559                             }
 560                         }
 561                     }
 562                 }
 563             }
 564         }
 565     }
 566 }
 567
 568 void set_pbc(t_pbc* pbc, PbcType pbcType, const matrix box)
 569 {
 570     if (pbcType == PbcType::Unset)
 571     {
 572         pbcType = guessPbcType(box);
 573     }
 574
 575     low_set_pbc(pbc, pbcType, nullptr, box);
 576 }
 577
 578 t_pbc* set_pbc_dd(t_pbc* pbc, PbcType pbcType, const ivec domdecCells, gmx_bool bSingleDir, const matrix box)
 579 {
 580     if (pbcType == PbcType::No)
 581     {
 582         pbc->pbcType = pbcType;
 583
 584         return nullptr;
 585     }
 586
 587     if (nullptr == domdecCells)
 588     {
 589         low_set_pbc(pbc, pbcType, nullptr, box);
 590     }
 591     else
 592     {
 593         if (pbcType == PbcType::Screw && domdecCells[XX] > 1)
 594         {
 595             /* The rotation has been taken care of during coordinate communication */
 596             pbcType = PbcType::Xyz;
 597         }
 598
 599         ivec usePBC;
 600         int  npbcdim = 0;
 601         for (int i = 0; i < DIM; i++)
 602         {
 603             usePBC[i] = 0;
 604             if (domdecCells[i] <= (bSingleDir ? 1 : 2) && !(pbcType == PbcType::XY && i == ZZ))
 605             {
 606                 usePBC[i] = 1;
 607                 npbcdim++;
 608             }
 609         }
 610
 611         if (npbcdim > 0)
 612         {
 613             low_set_pbc(pbc, pbcType, usePBC, box);
 614         }
 615         else
 616         {
 617             pbc->pbcType = PbcType::No;
 618         }
 619     }
 620
 621     return (pbc->pbcType != PbcType::No ? pbc : nullptr);
 622 }
 623
 624 void pbc_dx(const t_pbc* pbc, const rvec x1, const rvec x2, rvec dx)
 625 {
 626     int      i, j;
 627     rvec     dx_start, trial;
 628     real     d2min, d2trial;
 629     gmx_bool bRot;
 630
 631     rvec_sub(x1, x2, dx);
 632
 633     switch (pbc->pbcTypeDX)
 634     {
 635         case epbcdxRECTANGULAR:
 636             for (i = 0; i < DIM; i++)
 637             {
 638                 while (dx[i] > pbc->hbox_diag[i])
 639                 {
 640                     dx[i] -= pbc->fbox_diag[i];
 641                 }
 642                 while (dx[i] <= pbc->mhbox_diag[i])
 643                 {
 644                     dx[i] += pbc->fbox_diag[i];
 645                 }
 646             }
 647             break;
 648         case epbcdxTRICLINIC:
 649             for (i = DIM - 1; i >= 0; i--)
 650             {
 651                 while (dx[i] > pbc->hbox_diag[i])
 652                 {
 653                     for (j = i; j >= 0; j--)
 654                     {
 655                         dx[j] -= pbc->box[i][j];
 656                     }
 657                 }
 658                 while (dx[i] <= pbc->mhbox_diag[i])
 659                 {
 660                     for (j = i; j >= 0; j--)
 661                     {
 662                         dx[j] += pbc->box[i][j];
 663                     }
 664                 }
 665             }
 666             /* dx is the distance in a rectangular box */
 667             d2min = norm2(dx);
 668             if (d2min > pbc->max_cutoff2)
 669             {
 670                 copy_rvec(dx, dx_start);
 671                 d2min = norm2(dx);
 672                 /* Now try all possible shifts, when the distance is within max_cutoff
 673                  * it must be the shortest possible distance.
 674                  */
 675                 i = 0;
 676                 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
 677                 {
 678                     rvec_add(dx_start, pbc->tric_vec[i], trial);
 679                     d2trial = norm2(trial);
 680                     if (d2trial < d2min)
 681                     {
 682                         copy_rvec(trial, dx);
 683                         d2min = d2trial;
 684                     }
 685                     i++;
 686                 }
 687             }
 688             break;
 689         case epbcdx2D_RECT:
 690             for (i = 0; i < DIM; i++)
 691             {
 692                 if (i != pbc->dim)
 693                 {
 694                     while (dx[i] > pbc->hbox_diag[i])
 695                     {
 696                         dx[i] -= pbc->fbox_diag[i];
 697                     }
 698                     while (dx[i] <= pbc->mhbox_diag[i])
 699                     {
 700                         dx[i] += pbc->fbox_diag[i];
 701                     }
 702                 }
 703             }
 704             break;
 705         case epbcdx2D_TRIC:
 706             d2min = 0;
 707             for (i = DIM - 1; i >= 0; i--)
 708             {
 709                 if (i != pbc->dim)
 710                 {
 711                     while (dx[i] > pbc->hbox_diag[i])
 712                     {
 713                         for (j = i; j >= 0; j--)
 714                         {
 715                             dx[j] -= pbc->box[i][j];
 716                         }
 717                     }
 718                     while (dx[i] <= pbc->mhbox_diag[i])
 719                     {
 720                         for (j = i; j >= 0; j--)
 721                         {
 722                             dx[j] += pbc->box[i][j];
 723                         }
 724                     }
 725                     d2min += dx[i] * dx[i];
 726                 }
 727             }
 728             if (d2min > pbc->max_cutoff2)
 729             {
 730                 copy_rvec(dx, dx_start);
 731                 d2min = norm2(dx);
 732                 /* Now try all possible shifts, when the distance is within max_cutoff
 733                  * it must be the shortest possible distance.
 734                  */
 735                 i = 0;
 736                 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
 737                 {
 738                     rvec_add(dx_start, pbc->tric_vec[i], trial);
 739                     d2trial = 0;
 740                     for (j = 0; j < DIM; j++)
 741                     {
 742                         if (j != pbc->dim)
 743                         {
 744                             d2trial += trial[j] * trial[j];
 745                         }
 746                     }
 747                     if (d2trial < d2min)
 748                     {
 749                         copy_rvec(trial, dx);
 750                         d2min = d2trial;
 751                     }
 752                     i++;
 753                 }
 754             }
 755             break;
 756         case epbcdxSCREW_RECT:
 757             /* The shift definition requires x first */
 758             bRot = FALSE;
 759             while (dx[XX] > pbc->hbox_diag[XX])
 760             {
 761                 dx[XX] -= pbc->fbox_diag[XX];
 762                 bRot = !bRot;
 763             }
 764             while (dx[XX] <= pbc->mhbox_diag[XX])
 765             {
 766                 dx[XX] += pbc->fbox_diag[YY];
 767                 bRot = !bRot;
 768             }
 769             if (bRot)
 770             {
 771                 /* Rotate around the x-axis in the middle of the box */
 772                 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
 773                 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
 774             }
 775             /* Normal pbc for y and z */
 776             for (i = YY; i <= ZZ; i++)
 777             {
 778                 while (dx[i] > pbc->hbox_diag[i])
 779                 {
 780                     dx[i] -= pbc->fbox_diag[i];
 781                 }
 782                 while (dx[i] <= pbc->mhbox_diag[i])
 783                 {
 784                     dx[i] += pbc->fbox_diag[i];
 785                 }
 786             }
 787             break;
 788         case epbcdxNOPBC:
 789         case epbcdxUNSUPPORTED: break;
 790         default: gmx_fatal(FARGS, "Internal error in pbc_dx, set_pbc has not been called");
 791     }
 792 }
 793
 794 int pbc_dx_aiuc(const t_pbc* pbc, const rvec x1, const rvec x2, rvec dx)
 795 {
 796     int  i, j, is;
 797     rvec dx_start, trial;
 798     real d2min, d2trial;
 799     ivec ishift, ishift_start;
 800
 801     rvec_sub(x1, x2, dx);
 802     clear_ivec(ishift);
 803
 804     switch (pbc->pbcTypeDX)
 805     {
 806         case epbcdxRECTANGULAR:
 807             for (i = 0; i < DIM; i++)
 808             {
 809                 if (dx[i] > pbc->hbox_diag[i])
 810                 {
 811                     dx[i] -= pbc->fbox_diag[i];
 812                     ishift[i]--;
 813                 }
 814                 else if (dx[i] <= pbc->mhbox_diag[i])
 815                 {
 816                     dx[i] += pbc->fbox_diag[i];
 817                     ishift[i]++;
 818                 }
 819             }
 820             break;
 821         case epbcdxTRICLINIC:
 822             /* For triclinic boxes the performance difference between
 823              * if/else and two while loops is negligible.
 824              * However, the while version can cause extreme delays
 825              * before a simulation crashes due to large forces which
 826              * can cause unlimited displacements.
 827              * Also allowing multiple shifts would index fshift beyond bounds.
 828              */
 829             for (i = DIM - 1; i >= 1; i--)
 830             {
 831                 if (dx[i] > pbc->hbox_diag[i])
 832                 {
 833                     for (j = i; j >= 0; j--)
 834                     {
 835                         dx[j] -= pbc->box[i][j];
 836                     }
 837                     ishift[i]--;
 838                 }
 839                 else if (dx[i] <= pbc->mhbox_diag[i])
 840                 {
 841                     for (j = i; j >= 0; j--)
 842                     {
 843                         dx[j] += pbc->box[i][j];
 844                     }
 845                     ishift[i]++;
 846                 }
 847             }
 848             /* Allow 2 shifts in x */
 849             if (dx[XX] > pbc->hbox_diag[XX])
 850             {
 851                 dx[XX] -= pbc->fbox_diag[XX];
 852                 ishift[XX]--;
 853                 if (dx[XX] > pbc->hbox_diag[XX])
 854                 {
 855                     dx[XX] -= pbc->fbox_diag[XX];
 856                     ishift[XX]--;
 857                 }
 858             }
 859             else if (dx[XX] <= pbc->mhbox_diag[XX])
 860             {
 861                 dx[XX] += pbc->fbox_diag[XX];
 862                 ishift[XX]++;
 863                 if (dx[XX] <= pbc->mhbox_diag[XX])
 864                 {
 865                     dx[XX] += pbc->fbox_diag[XX];
 866                     ishift[XX]++;
 867                 }
 868             }
 869             /* dx is the distance in a rectangular box */
 870             d2min = norm2(dx);
 871             if (d2min > pbc->max_cutoff2)
 872             {
 873                 copy_rvec(dx, dx_start);
 874                 copy_ivec(ishift, ishift_start);
 875                 d2min = norm2(dx);
 876                 /* Now try all possible shifts, when the distance is within max_cutoff
 877                  * it must be the shortest possible distance.
 878                  */
 879                 i = 0;
 880                 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
 881                 {
 882                     rvec_add(dx_start, pbc->tric_vec[i], trial);
 883                     d2trial = norm2(trial);
 884                     if (d2trial < d2min)
 885                     {
 886                         copy_rvec(trial, dx);
 887                         ivec_add(ishift_start, pbc->tric_shift[i], ishift);
 888                         d2min = d2trial;
 889                     }
 890                     i++;
 891                 }
 892             }
 893             break;
 894         case epbcdx2D_RECT:
 895             for (i = 0; i < DIM; i++)
 896             {
 897                 if (i != pbc->dim)
 898                 {
 899                     if (dx[i] > pbc->hbox_diag[i])
 900                     {
 901                         dx[i] -= pbc->fbox_diag[i];
 902                         ishift[i]--;
 903                     }
 904                     else if (dx[i] <= pbc->mhbox_diag[i])
 905                     {
 906                         dx[i] += pbc->fbox_diag[i];
 907                         ishift[i]++;
 908                     }
 909                 }
 910             }
 911             break;
 912         case epbcdx2D_TRIC:
 913             d2min = 0;
 914             for (i = DIM - 1; i >= 1; i--)
 915             {
 916                 if (i != pbc->dim)
 917                 {
 918                     if (dx[i] > pbc->hbox_diag[i])
 919                     {
 920                         for (j = i; j >= 0; j--)
 921                         {
 922                             dx[j] -= pbc->box[i][j];
 923                         }
 924                         ishift[i]--;
 925                     }
 926                     else if (dx[i] <= pbc->mhbox_diag[i])
 927                     {
 928                         for (j = i; j >= 0; j--)
 929                         {
 930                             dx[j] += pbc->box[i][j];
 931                         }
 932                         ishift[i]++;
 933                     }
 934                     d2min += dx[i] * dx[i];
 935                 }
 936             }
 937             if (pbc->dim != XX)
 938             {
 939                 /* Allow 2 shifts in x */
 940                 if (dx[XX] > pbc->hbox_diag[XX])
 941                 {
 942                     dx[XX] -= pbc->fbox_diag[XX];
 943                     ishift[XX]--;
 944                     if (dx[XX] > pbc->hbox_diag[XX])
 945                     {
 946                         dx[XX] -= pbc->fbox_diag[XX];
 947                         ishift[XX]--;
 948                     }
 949                 }
 950                 else if (dx[XX] <= pbc->mhbox_diag[XX])
 951                 {
 952                     dx[XX] += pbc->fbox_diag[XX];
 953                     ishift[XX]++;
 954                     if (dx[XX] <= pbc->mhbox_diag[XX])
 955                     {
 956                         dx[XX] += pbc->fbox_diag[XX];
 957                         ishift[XX]++;
 958                     }
 959                 }
 960                 d2min += dx[XX] * dx[XX];
 961             }
 962             if (d2min > pbc->max_cutoff2)
 963             {
 964                 copy_rvec(dx, dx_start);
 965                 copy_ivec(ishift, ishift_start);
 966                 /* Now try all possible shifts, when the distance is within max_cutoff
 967                  * it must be the shortest possible distance.
 968                  */
 969                 i = 0;
 970                 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
 971                 {
 972                     rvec_add(dx_start, pbc->tric_vec[i], trial);
 973                     d2trial = 0;
 974                     for (j = 0; j < DIM; j++)
 975                     {
 976                         if (j != pbc->dim)
 977                         {
 978                             d2trial += trial[j] * trial[j];
 979                         }
 980                     }
 981                     if (d2trial < d2min)
 982                     {
 983                         copy_rvec(trial, dx);
 984                         ivec_add(ishift_start, pbc->tric_shift[i], ishift);
 985                         d2min = d2trial;
 986                     }
 987                     i++;
 988                 }
 989             }
 990             break;
 991         case epbcdx1D_RECT:
 992             i = pbc->dim;
 993             if (dx[i] > pbc->hbox_diag[i])
 994             {
 995                 dx[i] -= pbc->fbox_diag[i];
 996                 ishift[i]--;
 997             }
 998             else if (dx[i] <= pbc->mhbox_diag[i])
 999             {
1000                 dx[i] += pbc->fbox_diag[i];
1001                 ishift[i]++;
1002             }
1003             break;
1004         case epbcdx1D_TRIC:
1005             i = pbc->dim;
1006             if (dx[i] > pbc->hbox_diag[i])
1007             {
1008                 rvec_dec(dx, pbc->box[i]);
1009                 ishift[i]--;
1010             }
1011             else if (dx[i] <= pbc->mhbox_diag[i])
1012             {
1013                 rvec_inc(dx, pbc->box[i]);
1014                 ishift[i]++;
1015             }
1016             break;
1017         case epbcdxSCREW_RECT:
1018             /* The shift definition requires x first */
1019             if (dx[XX] > pbc->hbox_diag[XX])
1020             {
1021                 dx[XX] -= pbc->fbox_diag[XX];
1022                 ishift[XX]--;
1023             }
1024             else if (dx[XX] <= pbc->mhbox_diag[XX])
1025             {
1026                 dx[XX] += pbc->fbox_diag[XX];
1027                 ishift[XX]++;
1028             }
1029             if (ishift[XX] == 1 || ishift[XX] == -1)
1030             {
1031                 /* Rotate around the x-axis in the middle of the box */
1032                 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
1033                 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
1034             }
1035             /* Normal pbc for y and z */
1036             for (i = YY; i <= ZZ; i++)
1037             {
1038                 if (dx[i] > pbc->hbox_diag[i])
1039                 {
1040                     dx[i] -= pbc->fbox_diag[i];
1041                     ishift[i]--;
1042                 }
1043                 else if (dx[i] <= pbc->mhbox_diag[i])
1044                 {
1045                     dx[i] += pbc->fbox_diag[i];
1046                     ishift[i]++;
1047                 }
1048             }
1049             break;
1050         case epbcdxNOPBC:
1051         case epbcdxUNSUPPORTED: break;
1052         default:
1053             gmx_fatal(FARGS,
1054                       "Internal error in pbc_dx_aiuc, set_pbc_dd or set_pbc has not been called");
1055     }
1056
1057     is = IVEC2IS(ishift);
1058     if (debug)
1059     {
1060         range_check_mesg(is, 0, SHIFTS, "PBC shift vector index range check.");
1061     }
1062
1063     return is;
1064 }
1065
1066 //! Compute distance vector in double precision
1067 void pbc_dx_d(const t_pbc* pbc, const dvec x1, const dvec x2, dvec dx)
1068 {
1069     int      i, j;
1070     dvec     dx_start, trial;
1071     double   d2min, d2trial;
1072     gmx_bool bRot;
1073
1074     dvec_sub(x1, x2, dx);
1075
1076     switch (pbc->pbcTypeDX)
1077     {
1078         case epbcdxRECTANGULAR:
1079         case epbcdx2D_RECT:
1080             for (i = 0; i < DIM; i++)
1081             {
1082                 if (i != pbc->dim)
1083                 {
1084                     while (dx[i] > pbc->hbox_diag[i])
1085                     {
1086                         dx[i] -= pbc->fbox_diag[i];
1087                     }
1088                     while (dx[i] <= pbc->mhbox_diag[i])
1089                     {
1090                         dx[i] += pbc->fbox_diag[i];
1091                     }
1092                 }
1093             }
1094             break;
1095         case epbcdxTRICLINIC:
1096         case epbcdx2D_TRIC:
1097             d2min = 0;
1098             for (i = DIM - 1; i >= 0; i--)
1099             {
1100                 if (i != pbc->dim)
1101                 {
1102                     while (dx[i] > pbc->hbox_diag[i])
1103                     {
1104                         for (j = i; j >= 0; j--)
1105                         {
1106                             dx[j] -= pbc->box[i][j];
1107                         }
1108                     }
1109                     while (dx[i] <= pbc->mhbox_diag[i])
1110                     {
1111                         for (j = i; j >= 0; j--)
1112                         {
1113                             dx[j] += pbc->box[i][j];
1114                         }
1115                     }
1116                     d2min += dx[i] * dx[i];
1117                 }
1118             }
1119             if (d2min > pbc->max_cutoff2)
1120             {
1121                 copy_dvec(dx, dx_start);
1122                 /* Now try all possible shifts, when the distance is within max_cutoff
1123                  * it must be the shortest possible distance.
1124                  */
1125                 i = 0;
1126                 while ((d2min > pbc->max_cutoff2) && (i < pbc->ntric_vec))
1127                 {
1128                     for (j = 0; j < DIM; j++)
1129                     {
1130                         trial[j] = dx_start[j] + pbc->tric_vec[i][j];
1131                     }
1132                     d2trial = 0;
1133                     for (j = 0; j < DIM; j++)
1134                     {
1135                         if (j != pbc->dim)
1136                         {
1137                             d2trial += trial[j] * trial[j];
1138                         }
1139                     }
1140                     if (d2trial < d2min)
1141                     {
1142                         copy_dvec(trial, dx);
1143                         d2min = d2trial;
1144                     }
1145                     i++;
1146                 }
1147             }
1148             break;
1149         case epbcdxSCREW_RECT:
1150             /* The shift definition requires x first */
1151             bRot = FALSE;
1152             while (dx[XX] > pbc->hbox_diag[XX])
1153             {
1154                 dx[XX] -= pbc->fbox_diag[XX];
1155                 bRot = !bRot;
1156             }
1157             while (dx[XX] <= pbc->mhbox_diag[XX])
1158             {
1159                 dx[XX] += pbc->fbox_diag[YY];
1160                 bRot = !bRot;
1161             }
1162             if (bRot)
1163             {
1164                 /* Rotate around the x-axis in the middle of the box */
1165                 dx[YY] = pbc->box[YY][YY] - x1[YY] - x2[YY];
1166                 dx[ZZ] = pbc->box[ZZ][ZZ] - x1[ZZ] - x2[ZZ];
1167             }
1168             /* Normal pbc for y and z */
1169             for (i = YY; i <= ZZ; i++)
1170             {
1171                 while (dx[i] > pbc->hbox_diag[i])
1172                 {
1173                     dx[i] -= pbc->fbox_diag[i];
1174                 }
1175                 while (dx[i] <= pbc->mhbox_diag[i])
1176                 {
1177                     dx[i] += pbc->fbox_diag[i];
1178                 }
1179             }
1180             break;
1181         case epbcdxNOPBC:
1182         case epbcdxUNSUPPORTED: break;
1183         default: gmx_fatal(FARGS, "Internal error in pbc_dx, set_pbc has not been called");
1184     }
1185 }
1186
1187 void calc_shifts(const matrix box, rvec shift_vec[])
1188 {
1189     for (int n = 0, m = -D_BOX_Z; m <= D_BOX_Z; m++)
1190     {
1191         for (int l = -D_BOX_Y; l <= D_BOX_Y; l++)
1192         {
1193             for (int k = -D_BOX_X; k <= D_BOX_X; k++, n++)
1194             {
1195                 for (int d = 0; d < DIM; d++)
1196                 {
1197                     shift_vec[n][d] = k * box[XX][d] + l * box[YY][d] + m * box[ZZ][d];
1198                 }
1199             }
1200         }
1201     }
1202 }
1203
1204 void calc_box_center(int ecenter, const matrix box, rvec box_center)
1205 {
1206     int d, m;
1207
1208     clear_rvec(box_center);
1209     switch (ecenter)
1210     {
1211         case ecenterTRIC:
1212             for (m = 0; (m < DIM); m++)
1213             {
1214                 for (d = 0; d < DIM; d++)
1215                 {
1216                     box_center[d] += 0.5 * box[m][d];
1217                 }
1218             }
1219             break;
1220         case ecenterRECT:
1221             for (d = 0; d < DIM; d++)
1222             {
1223                 box_center[d] = 0.5 * box[d][d];
1224             }
1225             break;
1226         case ecenterZERO: break;
1227         default: gmx_fatal(FARGS, "Unsupported value %d for ecenter", ecenter);
1228     }
1229 }
1230
1231 void calc_triclinic_images(const matrix box, rvec img[])
1232 {
1233     int i;
1234
1235     /* Calculate 3 adjacent images in the xy-plane */
1236     copy_rvec(box[0], img[0]);
1237     copy_rvec(box[1], img[1]);
1238     if (img[1][XX] < 0)
1239     {
1240         svmul(-1, img[1], img[1]);
1241     }
1242     rvec_sub(img[1], img[0], img[2]);
1243
1244     /* Get the next 3 in the xy-plane as mirror images */
1245     for (i = 0; i < 3; i++)
1246     {
1247         svmul(-1, img[i], img[3 + i]);
1248     }
1249
1250     /* Calculate the first 4 out of xy-plane images */
1251     copy_rvec(box[2], img[6]);
1252     if (img[6][XX] < 0)
1253     {
1254         svmul(-1, img[6], img[6]);
1255     }
1256     for (i = 0; i < 3; i++)
1257     {
1258         rvec_add(img[6], img[i + 1], img[7 + i]);
1259     }
1260
1261     /* Mirror the last 4 from the previous in opposite rotation */
1262     for (i = 0; i < 4; i++)
1263     {
1264         svmul(-1, img[6 + (2 + i) % 4], img[10 + i]);
1265     }
1266 }
1267
1268 void calc_compact_unitcell_vertices(int ecenter, const matrix box, rvec vert[])
1269 {
1270     rvec       img[NTRICIMG], box_center;
1271     int        n, i, j, tmp[4], d;
1272     const real oneFourth = 0.25;
1273
1274     calc_triclinic_images(box, img);
1275
1276     n = 0;
1277     for (i = 2; i <= 5; i += 3)
1278     {
1279         tmp[0] = i - 1;
1280         if (i == 2)
1281         {
1282             tmp[1] = 8;
1283         }
1284         else
1285         {
1286             tmp[1] = 6;
1287         }
1288         tmp[2] = (i + 1) % 6;
1289         tmp[3] = tmp[1] + 4;
1290         for (j = 0; j < 4; j++)
1291         {
1292             for (d = 0; d < DIM; d++)
1293             {
1294                 vert[n][d] = img[i][d] + img[tmp[j]][d] + img[tmp[(j + 1) % 4]][d];
1295             }
1296             n++;
1297         }
1298     }
1299     for (i = 7; i <= 13; i += 6)
1300     {
1301         tmp[0] = (i - 7) / 2;
1302         tmp[1] = tmp[0] + 1;
1303         if (i == 7)
1304         {
1305             tmp[2] = 8;
1306         }
1307         else
1308         {
1309             tmp[2] = 10;
1310         }
1311         tmp[3] = i - 1;
1312         for (j = 0; j < 4; j++)
1313         {
1314             for (d = 0; d < DIM; d++)
1315             {
1316                 vert[n][d] = img[i][d] + img[tmp[j]][d] + img[tmp[(j + 1) % 4]][d];
1317             }
1318             n++;
1319         }
1320     }
1321     for (i = 9; i <= 11; i += 2)
1322     {
1323         if (i == 9)
1324         {
1325             tmp[0] = 3;
1326         }
1327         else
1328         {
1329             tmp[0] = 0;
1330         }
1331         tmp[1] = tmp[0] + 1;
1332         if (i == 9)
1333         {
1334             tmp[2] = 6;
1335         }
1336         else
1337         {
1338             tmp[2] = 12;
1339         }
1340         tmp[3] = i - 1;
1341         for (j = 0; j < 4; j++)
1342         {
1343             for (d = 0; d < DIM; d++)
1344             {
1345                 vert[n][d] = img[i][d] + img[tmp[j]][d] + img[tmp[(j + 1) % 4]][d];
1346             }
1347             n++;
1348         }
1349     }
1350
1351     calc_box_center(ecenter, box, box_center);
1352     for (i = 0; i < NCUCVERT; i++)
1353     {
1354         for (d = 0; d < DIM; d++)
1355         {
1356             vert[i][d] = vert[i][d] * oneFourth + box_center[d];
1357         }
1358     }
1359 }
1360
1361 int* compact_unitcell_edges()
1362 {
1363     /* this is an index in vert[] (see calc_box_vertices) */
1364     /*static int edge[NCUCEDGE*2];*/
1365     int*             edge;
1366     static const int hexcon[24] = { 0, 9,  1, 19, 2, 15, 3,  21, 4,  17, 5,  11,
1367                                     6, 23, 7, 13, 8, 20, 10, 18, 12, 16, 14, 22 };
1368     int              e, i, j;
1369
1370     snew(edge, NCUCEDGE * 2);
1371
1372     e = 0;
1373     for (i = 0; i < 6; i++)
1374     {
1375         for (j = 0; j < 4; j++)
1376         {
1377             edge[e++] = 4 * i + j;
1378             edge[e++] = 4 * i + (j + 1) % 4;
1379         }
1380     }
1381     for (i = 0; i < 12 * 2; i++)
1382     {
1383         edge[e++] = hexcon[i];
1384     }
1385
1386     return edge;
1387 }
1388
1389 void put_atoms_in_box(PbcType pbcType, const matrix box, gmx::ArrayRef<gmx::RVec> x)
1390 {
1391     int npbcdim, m, d;
1392
1393     if (pbcType == PbcType::Screw)
1394     {
1395         gmx_fatal(FARGS, "Sorry, %s pbc is not yet supported", c_pbcTypeNames[pbcType].c_str());
1396     }
1397
1398     if (pbcType == PbcType::XY)
1399     {
1400         npbcdim = 2;
1401     }
1402     else
1403     {
1404         npbcdim = 3;
1405     }
1406
1407     if (TRICLINIC(box))
1408     {
1409         for (gmx::index i = 0; (i < x.ssize()); ++i)
1410         {
1411             for (m = npbcdim - 1; m >= 0; m--)
1412             {
1413                 while (x[i][m] < 0)
1414                 {
1415                     for (d = 0; d <= m; d++)
1416                     {
1417                         x[i][d] += box[m][d];
1418                     }
1419                 }
1420                 while (x[i][m] >= box[m][m])
1421                 {
1422                     for (d = 0; d <= m; d++)
1423                     {
1424                         x[i][d] -= box[m][d];
1425                     }
1426                 }
1427             }
1428         }
1429     }
1430     else
1431     {
1432         for (gmx::index i = 0; (i < x.ssize()); ++i)
1433         {
1434             for (d = 0; d < npbcdim; d++)
1435             {
1436                 while (x[i][d] < 0)
1437                 {
1438                     x[i][d] += box[d][d];
1439                 }
1440                 while (x[i][d] >= box[d][d])
1441                 {
1442                     x[i][d] -= box[d][d];
1443                 }
1444             }
1445         }
1446     }
1447 }
1448
1449 void put_atoms_in_box_omp(PbcType pbcType, const matrix box, gmx::ArrayRef<gmx::RVec> x, gmx_unused int nth)
1450 {
1451 #pragma omp parallel for num_threads(nth) schedule(static)
1452     for (int t = 0; t < nth; t++)
1453     {
1454         try
1455         {
1456             size_t natoms = x.size();
1457             size_t offset = (natoms * t) / nth;
1458             size_t len    = (natoms * (t + 1)) / nth - offset;
1459             put_atoms_in_box(pbcType, box, x.subArray(offset, len));
1460         }
1461         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1462     }
1463 }
1464
1465 void put_atoms_in_triclinic_unitcell(int ecenter, const matrix box, gmx::ArrayRef<gmx::RVec> x)
1466 {
1467     rvec box_center, shift_center;
1468     real shm01, shm02, shm12, shift;
1469     int  m, d;
1470
1471     calc_box_center(ecenter, box, box_center);
1472
1473     /* The product of matrix shm with a coordinate gives the shift vector
1474        which is required determine the periodic cell position */
1475     shm01 = box[1][0] / box[1][1];
1476     shm02 = (box[1][1] * box[2][0] - box[2][1] * box[1][0]) / (box[1][1] * box[2][2]);
1477     shm12 = box[2][1] / box[2][2];
1478
1479     clear_rvec(shift_center);
1480     for (d = 0; d < DIM; d++)
1481     {
1482         rvec_inc(shift_center, box[d]);
1483     }
1484     svmul(0.5, shift_center, shift_center);
1485     rvec_sub(box_center, shift_center, shift_center);
1486
1487     shift_center[0] = shm01 * shift_center[1] + shm02 * shift_center[2];
1488     shift_center[1] = shm12 * shift_center[2];
1489     shift_center[2] = 0;
1490
1491     for (gmx::index i = 0; (i < x.ssize()); ++i)
1492     {
1493         for (m = DIM - 1; m >= 0; m--)
1494         {
1495             shift = shift_center[m];
1496             if (m == 0)
1497             {
1498                 shift += shm01 * x[i][1] + shm02 * x[i][2];
1499             }
1500             else if (m == 1)
1501             {
1502                 shift += shm12 * x[i][2];
1503             }
1504             while (x[i][m] - shift < 0)
1505             {
1506                 for (d = 0; d <= m; d++)
1507                 {
1508                     x[i][d] += box[m][d];
1509                 }
1510             }
1511             while (x[i][m] - shift >= box[m][m])
1512             {
1513                 for (d = 0; d <= m; d++)
1514                 {
1515                     x[i][d] -= box[m][d];
1516                 }
1517             }
1518         }
1519     }
1520 }
1521
1522 void put_atoms_in_compact_unitcell(PbcType pbcType, int ecenter, const matrix box, gmx::ArrayRef<gmx::RVec> x)
1523 {
1524     t_pbc pbc;
1525     rvec  box_center, dx;
1526
1527     set_pbc(&pbc, pbcType, box);
1528
1529     if (pbc.pbcTypeDX == epbcdxUNSUPPORTED)
1530     {
1531         gmx_fatal(FARGS, "Can not put atoms in compact unitcell with unsupported PBC");
1532     }
1533
1534     calc_box_center(ecenter, box, box_center);
1535     for (gmx::index i = 0; (i < x.ssize()); ++i)
1536     {
1537         pbc_dx(&pbc, x[i], box_center, dx);
1538         rvec_add(box_center, dx, x[i]);
1539     }
1540 }
1541
1542 /*! \brief Make molecules whole by shifting positions
1543  *
1544  * \param[in]     fplog     Log file
1545  * \param[in]     pbcType   The PBC type
1546  * \param[in]     box       The simulation box
1547  * \param[in]     mtop      System topology definition
1548  * \param[in,out] x         The coordinates of the atoms
1549  * \param[in]     bFirst    Specifier for first-time PBC removal
1550  */
1551 static void
1552 low_do_pbc_mtop(FILE* fplog, PbcType pbcType, const matrix box, const gmx_mtop_t* mtop, rvec x[], gmx_bool bFirst)
1553 {
1554     int as, mol;
1555
1556     if (bFirst && fplog)
1557     {
1558         fprintf(fplog, "Removing pbc first time\n");
1559     }
1560
1561     as = 0;
1562     for (const gmx_molblock_t& molb : mtop->molblock)
1563     {
1564         const gmx_moltype_t& moltype = mtop->moltype[molb.type];
1565         if (moltype.atoms.nr == 1 || (!bFirst && moltype.atoms.nr == 1))
1566         {
1567             /* Just one atom or charge group in the molecule, no PBC required */
1568             as += molb.nmol * moltype.atoms.nr;
1569         }
1570         else
1571         {
1572             t_graph graph = mk_graph_moltype(moltype);
1573
1574             for (mol = 0; mol < molb.nmol; mol++)
1575             {
1576                 mk_mshift(fplog, &graph, pbcType, box, x + as);
1577
1578                 shift_self(graph, box, x + as);
1579                 /* The molecule is whole now.
1580                  * We don't need the second mk_mshift call as in do_pbc_first,
1581                  * since we no longer need this graph.
1582                  */
1583
1584                 as += moltype.atoms.nr;
1585             }
1586         }
1587     }
1588 }
1589
1590 void do_pbc_first_mtop(FILE* fplog, PbcType pbcType, const matrix box, const gmx_mtop_t* mtop, rvec x[])
1591 {
1592     low_do_pbc_mtop(fplog, pbcType, box, mtop, x, TRUE);
1593 }
1594
1595 void do_pbc_mtop(PbcType pbcType, const matrix box, const gmx_mtop_t* mtop, rvec x[])
1596 {
1597     low_do_pbc_mtop(nullptr, pbcType, box, mtop, x, FALSE);
1598 }