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