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