src/gromacs/pbcutil/pbc.cpp

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