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