src/gromacs/pbcutil/mshift.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 #include "gmxpre.h"
  39
  40 #include "mshift.h"
  41
  42 #include <cstring>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/pbcutil/pbc.h"
  48 #include "gromacs/topology/idef.h"
  49 #include "gromacs/topology/ifunc.h"
  50 #include "gromacs/topology/topology.h"
  51 #include "gromacs/utility/fatalerror.h"
  52 #include "gromacs/utility/strconvert.h"
  53 #include "gromacs/utility/stringutil.h"
  54
  55 /************************************************************
  56  *
  57  *      S H I F T   U T I L I T I E S
  58  *
  59  ************************************************************/
  60
  61
  62 /************************************************************
  63  *
  64  *      G R A P H   G E N E R A T I O N    C O D E
  65  *
  66  ************************************************************/
  67
  68 using gmx::ArrayRef;
  69 using gmx::IVec;
  70
  71 // Class for generating the edges of the graph
  72 class EdgesGenerator
  73 {
  74 public:
  75     EdgesGenerator(const int endAtom) : edges_(endAtom) {}
  76
  77     // Adds an edge, bi-directional
  78     void addEdge(int a0, int a1);
  79
  80     // Returns the edges
  81     const std::vector<std::vector<int>>& edges() const { return edges_; }
  82
  83 private:
  84     // The edges stored as list (for each atom starting at \p startAtom_) of lists of atoms
  85     std::vector<std::vector<int>> edges_;
  86 };
  87
  88 void EdgesGenerator::addEdge(const int a0, const int a1)
  89 {
  90     const auto& edges = edges_[a0];
  91     if (std::find(edges.begin(), edges.end(), a1) == edges.end())
  92     {
  93         edges_[a0].push_back(a1);
  94         edges_[a1].push_back(a0);
  95     }
  96 }
  97
  98 /* Make the actual graph for an ilist, returns whether an edge was added.
  99  *
 100  * When a non-null part array is supplied with part indices for each atom,
 101  * edges are only added when atoms have a different part index.
 102  */
 103 template<typename T>
 104 static bool mk_igraph(EdgesGenerator* edgesG, int ftype, const T& il, int at_end, ArrayRef<const int> part)
 105 {
 106     int  i, j, np;
 107     int  end;
 108     bool addedEdge = false;
 109
 110     end = il.size();
 111
 112     i = 0;
 113     while (i < end)
 114     {
 115         np = interaction_function[ftype].nratoms;
 116
 117         if (np > 1 && il.iatoms[i + 1] < at_end)
 118         {
 119             if (il.iatoms[i + np] >= at_end)
 120             {
 121                 gmx_fatal(FARGS,
 122                           "Molecule in topology has atom numbers below and "
 123                           "above natoms (%d).\n"
 124                           "You are probably trying to use a trajectory which does "
 125                           "not match the first %d atoms of the run input file.\n"
 126                           "You can make a matching run input file with gmx convert-tpr.",
 127                           at_end,
 128                           at_end);
 129             }
 130             if (ftype == F_SETTLE)
 131             {
 132                 /* Bond all the atoms in the settle */
 133                 edgesG->addEdge(il.iatoms[i + 1], il.iatoms[i + 2]);
 134                 edgesG->addEdge(il.iatoms[i + 1], il.iatoms[i + 3]);
 135                 addedEdge = true;
 136             }
 137             else if (part.empty())
 138             {
 139                 /* Simply add this bond */
 140                 for (j = 1; j < np; j++)
 141                 {
 142                     edgesG->addEdge(il.iatoms[i + j], il.iatoms[i + j + 1]);
 143                 }
 144                 addedEdge = true;
 145             }
 146             else
 147             {
 148                 /* Add this bond when it connects two unlinked parts of the graph */
 149                 for (j = 1; j < np; j++)
 150                 {
 151                     if (part[il.iatoms[i + j]] != part[il.iatoms[i + j + 1]])
 152                     {
 153                         edgesG->addEdge(il.iatoms[i + j], il.iatoms[i + j + 1]);
 154                         addedEdge = true;
 155                     }
 156                 }
 157             }
 158         }
 159
 160         i += np + 1;
 161     }
 162
 163     return addedEdge;
 164 }
 165
 166 [[noreturn]] static void g_error(int line, const char* file)
 167 {
 168     gmx_fatal(FARGS, "Trying to print nonexistent graph (file %s, line %d)", file, line);
 169 }
 170 #define GCHECK(g)    \
 171     if ((g) == NULL) \
 172     g_error(__LINE__, __FILE__)
 173
 174 void p_graph(FILE* log, const char* title, const t_graph* g)
 175 {
 176     int         i;
 177     const char* cc[egcolNR] = { "W", "G", "B" };
 178
 179     GCHECK(g);
 180     fprintf(log, "graph:  %s\n", title);
 181     fprintf(log, "nnodes: %d\n", g->numNodes());
 182     fprintf(log, "nbound: %d\n", g->numConnectedAtoms);
 183     fprintf(log, "start:  %d\n", g->edgeAtomBegin);
 184     fprintf(log, "end:    %d\n", g->edgeAtomEnd);
 185     fprintf(log, " atom shiftx shifty shiftz C nedg    e1    e2 etc.\n");
 186     for (i = 0; i < int(g->edges.size()); i++)
 187     {
 188         if (!g->edges[i].empty())
 189         {
 190             fprintf(log,
 191                     "%5d%7d%7d%7d %1s%5zu",
 192                     g->edgeAtomBegin + i + 1,
 193                     g->ishift[g->edgeAtomBegin + i][XX],
 194                     g->ishift[g->edgeAtomBegin + i][YY],
 195                     g->ishift[g->edgeAtomBegin + i][ZZ],
 196                     (!g->edgeColor.empty()) ? cc[g->edgeColor[i]] : " ",
 197                     g->edges[i].size());
 198             for (const int edge : g->edges[i])
 199             {
 200                 fprintf(log, " %5d", edge + 1);
 201             }
 202             fprintf(log, "\n");
 203         }
 204     }
 205     fflush(log);
 206 }
 207
 208 /* Converts the vector of vector of edges to ListOfLists
 209  * and removes leading and trailing uncoupled atoms
 210  */
 211 static gmx::ListOfLists<int> convertGraph(FILE*                 fplog,
 212                                           const EdgesGenerator& edgesG,
 213                                           int*                  firstConnectedAtom,
 214                                           int*                  numConnectedAtoms)
 215 {
 216     gmx::ListOfLists<int> edgesLists;
 217
 218     *firstConnectedAtom        = edgesG.edges().size();
 219     *numConnectedAtoms         = 0;
 220     int numEmptyEntriesSkipped = 0;
 221     int max_nedge              = 0;
 222     for (const auto& edges : edgesG.edges())
 223     {
 224         if (edges.empty())
 225         {
 226             numEmptyEntriesSkipped++;
 227         }
 228         else
 229         {
 230             if (edgesLists.empty())
 231             {
 232                 /* We ignore empty entries before the first connected entry */
 233                 *firstConnectedAtom = numEmptyEntriesSkipped;
 234             }
 235             else
 236             {
 237                 /* Push any empty entries we skipped */
 238                 for (int i = 0; i < numEmptyEntriesSkipped; i++)
 239                 {
 240                     edgesLists.pushBack({});
 241                 }
 242             }
 243             numEmptyEntriesSkipped = 0;
 244
 245             edgesLists.pushBack(edges);
 246
 247             (*numConnectedAtoms)++;
 248
 249             max_nedge = std::max(max_nedge, int(gmx::ssize(edges)));
 250         }
 251     }
 252
 253     if (fplog)
 254     {
 255         fprintf(fplog, "Max number of graph edges per atom is %d\n", max_nedge);
 256         fprintf(fplog, "Total number of graph edges is %d\n", edgesLists.numElements());
 257     }
 258
 259     return edgesLists;
 260 }
 261
 262 static gmx_bool determine_graph_parts(const EdgesGenerator& edgesG, ArrayRef<int> partNr)
 263 {
 264     /* Initialize the part array with all entries different */
 265     const int endAtom = edgesG.edges().size();
 266     for (int at_i = 0; at_i < endAtom; at_i++)
 267     {
 268         partNr[at_i] = at_i;
 269     }
 270
 271     /* Loop over the graph until the part array is fixed */
 272     bool haveMultipleParts = false;
 273     int  numAtomsChanged   = 0;
 274     do
 275     {
 276         haveMultipleParts = false;
 277         numAtomsChanged   = 0;
 278         for (gmx::index at_i = 0; at_i < gmx::ssize(edgesG.edges()); at_i++)
 279         {
 280             for (const int at_i2 : edgesG.edges()[at_i])
 281             {
 282                 /* Set part for both nodes to the minimum */
 283                 if (partNr[at_i2] > partNr[at_i])
 284                 {
 285                     partNr[at_i2] = partNr[at_i];
 286                     numAtomsChanged++;
 287                 }
 288                 else if (partNr[at_i2] < partNr[at_i])
 289                 {
 290                     partNr[at_i] = partNr[at_i2];
 291                     numAtomsChanged++;
 292                 }
 293             }
 294             if (partNr[at_i] != partNr[0])
 295             {
 296                 haveMultipleParts = true;
 297             }
 298         }
 299         if (debug)
 300         {
 301             fprintf(debug,
 302                     "graph partNr[] numAtomsChanged=%d, bMultiPart=%s\n",
 303                     numAtomsChanged,
 304                     gmx::boolToString(haveMultipleParts));
 305         }
 306     } while (numAtomsChanged > 0);
 307
 308     return haveMultipleParts;
 309 }
 310
 311 template<typename T>
 312 static t_graph mk_graph_ilist(FILE* fplog, const T* ilist, int at_end, gmx_bool bShakeOnly, gmx_bool bSettle)
 313 {
 314     EdgesGenerator edgesG(at_end);
 315
 316     t_graph::BondedParts parts = t_graph::BondedParts::Single;
 317
 318     if (at_end > 0)
 319     {
 320         if (!bShakeOnly)
 321         {
 322             /* First add all the real bonds: they should determine the molecular
 323              * graph.
 324              */
 325             for (int i = 0; (i < F_NRE); i++)
 326             {
 327                 if (interaction_function[i].flags & IF_CHEMBOND)
 328                 {
 329                     mk_igraph(&edgesG, i, ilist[i], at_end, {});
 330                 }
 331             }
 332
 333             /* Determine of which separated parts the IF_CHEMBOND graph consists.
 334              * Store the part numbers in the partNr array.
 335              */
 336             std::vector<int> partNr(at_end);
 337             const bool       bMultiPart = determine_graph_parts(edgesG, partNr);
 338
 339             if (bMultiPart)
 340             {
 341                 /* Then add all the other interactions in fixed lists,
 342                  * but only when they connect parts of the graph
 343                  * that are not connected through IF_CHEMBOND interactions.
 344                  */
 345                 bool addedEdge = false;
 346                 for (int i = 0; (i < F_NRE); i++)
 347                 {
 348                     if (!(interaction_function[i].flags & IF_CHEMBOND))
 349                     {
 350                         bool addedEdgeForType = mk_igraph(&edgesG, i, ilist[i], at_end, partNr);
 351                         addedEdge             = (addedEdge || addedEdgeForType);
 352                     }
 353                 }
 354
 355                 if (addedEdge)
 356                 {
 357                     parts = t_graph::BondedParts::MultipleConnected;
 358                 }
 359                 else
 360                 {
 361                     parts = t_graph::BondedParts::MultipleDisconnected;
 362                 }
 363             }
 364         }
 365         else
 366         {
 367             /* This is a special thing used in splitter.c to generate shake-blocks */
 368             mk_igraph(&edgesG, F_CONSTR, ilist[F_CONSTR], at_end, {});
 369             if (bSettle)
 370             {
 371                 mk_igraph(&edgesG, F_SETTLE, ilist[F_SETTLE], at_end, {});
 372             }
 373         }
 374     }
 375
 376     t_graph graph;
 377     /* The naming is somewhat confusing, but we need g->shiftAtomEnd
 378      * for shifthing coordinates to a new array (not in place) when
 379      * some atoms are not connected by the graph, which runs from
 380      * g->edgeAtomBegin (>= 0) to g->edgeAtomEnd (<= g->shiftAtomEnd).
 381      */
 382     graph.shiftAtomEnd  = at_end;
 383     graph.edgeAtomBegin = 0;
 384     graph.edgeAtomEnd   = at_end;
 385     graph.parts         = parts;
 386     if (at_end > 0)
 387     {
 388         /* Convert the vector of vector of edges to ListOfLists */
 389         graph.edges = convertGraph(fplog, edgesG, &graph.edgeAtomBegin, &graph.numConnectedAtoms);
 390
 391         graph.edgeAtomEnd = graph.edgeAtomBegin + graph.edges.ssize();
 392     }
 393
 394     graph.edgeColor.resize(graph.edges.size());
 395     graph.ishift.resize(graph.shiftAtomEnd);
 396
 397     if (gmx_debug_at)
 398     {
 399         p_graph(debug, "graph", &graph);
 400     }
 401
 402     return graph;
 403 }
 404
 405 t_graph mk_graph_moltype(const gmx_moltype_t& moltype)
 406 {
 407     return mk_graph_ilist(nullptr, moltype.ilist.data(), moltype.atoms.nr, FALSE, FALSE);
 408 }
 409
 410 t_graph mk_graph(const InteractionDefinitions& idef, const int numAtoms)
 411 {
 412     return mk_graph_ilist(nullptr, idef.il.data(), numAtoms, false, false);
 413 }
 414
 415 t_graph* mk_graph(FILE* fplog, const InteractionDefinitions& idef, int at_end, gmx_bool bShakeOnly, gmx_bool bSettle)
 416 {
 417     t_graph* g = new (t_graph);
 418
 419     *g = mk_graph_ilist(fplog, idef.il.data(), at_end, bShakeOnly, bSettle);
 420
 421     return g;
 422 }
 423
 424 t_graph* mk_graph(FILE* fplog, const t_idef* idef, int at_end, gmx_bool bShakeOnly, gmx_bool bSettle)
 425 {
 426     t_graph* g = new (t_graph);
 427
 428     *g = mk_graph_ilist(fplog, idef->il, at_end, bShakeOnly, bSettle);
 429
 430     return g;
 431 }
 432
 433 void done_graph(t_graph* g)
 434 {
 435     delete g;
 436 }
 437
 438 /************************************************************
 439  *
 440  *      S H I F T   C A L C U L A T I O N   C O D E
 441  *
 442  ************************************************************/
 443
 444 static void mk_1shift_tric(int          npbcdim,
 445                            const matrix box,
 446                            const rvec   hbox,
 447                            const rvec   xi,
 448                            const rvec   xj,
 449                            const int*   mi,
 450                            int*         mj)
 451 {
 452     /* Calculate periodicity for triclinic box... */
 453     int  m, d;
 454     rvec dx;
 455
 456     rvec_sub(xi, xj, dx);
 457
 458     mj[ZZ] = 0;
 459     for (m = npbcdim - 1; (m >= 0); m--)
 460     {
 461         /* If dx < hbox, then xj will be reduced by box, so that
 462          * xi - xj will be bigger
 463          */
 464         if (dx[m] < -hbox[m])
 465         {
 466             mj[m] = mi[m] - 1;
 467             for (d = m - 1; d >= 0; d--)
 468             {
 469                 dx[d] += box[m][d];
 470             }
 471         }
 472         else if (dx[m] >= hbox[m])
 473         {
 474             mj[m] = mi[m] + 1;
 475             for (d = m - 1; d >= 0; d--)
 476             {
 477                 dx[d] -= box[m][d];
 478             }
 479         }
 480         else
 481         {
 482             mj[m] = mi[m];
 483         }
 484     }
 485 }
 486
 487 static void mk_1shift(int npbcdim, const rvec hbox, const rvec xi, const rvec xj, const int* mi, int* mj)
 488 {
 489     /* Calculate periodicity for rectangular box... */
 490     int  m;
 491     rvec dx;
 492
 493     rvec_sub(xi, xj, dx);
 494
 495     mj[ZZ] = 0;
 496     for (m = 0; (m < npbcdim); m++)
 497     {
 498         /* If dx < hbox, then xj will be reduced by box, so that
 499          * xi - xj will be bigger
 500          */
 501         if (dx[m] < -hbox[m])
 502         {
 503             mj[m] = mi[m] - 1;
 504         }
 505         else if (dx[m] >= hbox[m])
 506         {
 507             mj[m] = mi[m] + 1;
 508         }
 509         else
 510         {
 511             mj[m] = mi[m];
 512         }
 513     }
 514 }
 515
 516 static void mk_1shift_screw(const matrix box, const rvec hbox, const rvec xi, const rvec xj, const int* mi, int* mj)
 517 {
 518     /* Calculate periodicity for rectangular box... */
 519     int  signi, m;
 520     rvec dx;
 521
 522     if ((mi[XX] > 0 && mi[XX] % 2 == 1) || (mi[XX] < 0 && -mi[XX] % 2 == 1))
 523     {
 524         signi = -1;
 525     }
 526     else
 527     {
 528         signi = 1;
 529     }
 530
 531     rvec_sub(xi, xj, dx);
 532
 533     if (dx[XX] < -hbox[XX])
 534     {
 535         mj[XX] = mi[XX] - 1;
 536     }
 537     else if (dx[XX] >= hbox[XX])
 538     {
 539         mj[XX] = mi[XX] + 1;
 540     }
 541     else
 542     {
 543         mj[XX] = mi[XX];
 544     }
 545     if (mj[XX] != mi[XX])
 546     {
 547         /* Rotate */
 548         dx[YY] = xi[YY] - (box[YY][YY] + box[ZZ][YY] - xj[YY]);
 549         dx[ZZ] = xi[ZZ] - (box[ZZ][ZZ] - xj[ZZ]);
 550     }
 551     for (m = 1; (m < DIM); m++)
 552     {
 553         /* The signs are taken such that we can first shift x and rotate
 554          * and then shift y and z.
 555          */
 556         if (dx[m] < -hbox[m])
 557         {
 558             mj[m] = mi[m] - signi;
 559         }
 560         else if (dx[m] >= hbox[m])
 561         {
 562             mj[m] = mi[m] + signi;
 563         }
 564         else
 565         {
 566             mj[m] = mi[m];
 567         }
 568     }
 569 }
 570
 571 static int mk_grey(ArrayRef<egCol> edgeColor,
 572                    t_graph*        g,
 573                    int*            AtomI,
 574                    int             npbcdim,
 575                    const matrix    box,
 576                    const rvec      x[],
 577                    int*            nerror)
 578 {
 579     int      m, ng, ai, g0;
 580     rvec     dx, hbox;
 581     gmx_bool bTriclinic;
 582     ivec     is_aj;
 583     t_pbc    pbc;
 584
 585     for (m = 0; (m < DIM); m++)
 586     {
 587         hbox[m] = box[m][m] * 0.5;
 588     }
 589     bTriclinic = TRICLINIC(box);
 590
 591     g0 = g->edgeAtomBegin;
 592     ng = 0;
 593     ai = g0 + *AtomI;
 594
 595     /* Loop over all the bonds */
 596     for (const int aj : g->edges[ai - g0])
 597     {
 598         /* If there is a white one, make it grey and set pbc */
 599         if (g->useScrewPbc)
 600         {
 601             mk_1shift_screw(box, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
 602         }
 603         else if (bTriclinic)
 604         {
 605             mk_1shift_tric(npbcdim, box, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
 606         }
 607         else
 608         {
 609             mk_1shift(npbcdim, hbox, x[ai], x[aj], g->ishift[ai], is_aj);
 610         }
 611
 612         if (edgeColor[aj - g0] == egcolWhite)
 613         {
 614             if (aj - g0 < *AtomI)
 615             {
 616                 *AtomI = aj - g0;
 617             }
 618             edgeColor[aj - g0] = egcolGrey;
 619
 620             copy_ivec(is_aj, g->ishift[aj]);
 621
 622             ng++;
 623         }
 624         else if ((is_aj[XX] != g->ishift[aj][XX]) || (is_aj[YY] != g->ishift[aj][YY])
 625                  || (is_aj[ZZ] != g->ishift[aj][ZZ]))
 626         {
 627             if (gmx_debug_at)
 628             {
 629                 set_pbc(&pbc, PbcType::Unset, box);
 630                 pbc_dx(&pbc, x[ai], x[aj], dx);
 631                 fprintf(debug,
 632                         "mk_grey: shifts for atom %d due to atom %d\n"
 633                         "are (%d,%d,%d), should be (%d,%d,%d)\n"
 634                         "dx = (%g,%g,%g)\n",
 635                         aj + 1,
 636                         ai + 1,
 637                         is_aj[XX],
 638                         is_aj[YY],
 639                         is_aj[ZZ],
 640                         g->ishift[aj][XX],
 641                         g->ishift[aj][YY],
 642                         g->ishift[aj][ZZ],
 643                         dx[XX],
 644                         dx[YY],
 645                         dx[ZZ]);
 646             }
 647             (*nerror)++;
 648         }
 649     }
 650     return ng;
 651 }
 652
 653 /* Return the first node/atom with colour Col starting at fC.
 654  * return -1 if none found.
 655  */
 656 static gmx::index first_colour(const int fC, const egCol Col, const t_graph* g, ArrayRef<const egCol> edgeColor)
 657 {
 658     for (gmx::index i = fC; i < gmx::ssize(g->edges); i++)
 659     {
 660         if (!g->edges[i].empty() && edgeColor[i] == Col)
 661         {
 662             return i;
 663         }
 664     }
 665
 666     return -1;
 667 }
 668
 669 /* Returns the maximum length of the graph edges for coordinates x */
 670 static real maxEdgeLength(const t_graph& g, PbcType pbcType, const matrix box, const rvec x[])
 671 {
 672     t_pbc pbc;
 673
 674     set_pbc(&pbc, pbcType, box);
 675
 676     real maxEdgeLength2 = 0;
 677
 678     for (int node = 0; node < int(g.edges.size()); node++)
 679     {
 680         for (const int nodeJ : g.edges[node])
 681         {
 682             rvec dx;
 683             pbc_dx(&pbc, x[node], x[nodeJ], dx);
 684             maxEdgeLength2 = std::max(maxEdgeLength2, norm2(dx));
 685         }
 686     }
 687
 688     return std::sqrt(maxEdgeLength2);
 689 }
 690
 691 void mk_mshift(FILE* log, t_graph* g, PbcType pbcType, const matrix box, const rvec x[])
 692 {
 693     static int nerror_tot = 0;
 694     int        npbcdim;
 695     int        ng, i;
 696     int        nW, nG, nB; /* Number of Grey, Black, White      */
 697     int        fW, fG;     /* First of each category    */
 698     int        nerror = 0;
 699
 700     g->useScrewPbc = (pbcType == PbcType::Screw);
 701
 702     if (pbcType == PbcType::XY)
 703     {
 704         npbcdim = 2;
 705     }
 706     else
 707     {
 708         npbcdim = 3;
 709     }
 710
 711     GCHECK(g);
 712     /* This puts everything in the central box, that is does not move it
 713      * at all. If we return without doing this for a system without bonds
 714      * (i.e. only settles) all water molecules are moved to the opposite octant
 715      */
 716     for (i = 0; i < g->shiftAtomEnd; i++)
 717     {
 718         g->ishift[i][XX] = g->ishift[i][YY] = g->ishift[i][ZZ] = 0;
 719     }
 720
 721     if (!g->numConnectedAtoms)
 722     {
 723         return;
 724     }
 725
 726     std::fill(g->edgeColor.begin(), g->edgeColor.end(), egcolWhite);
 727
 728     nW = g->numConnectedAtoms;
 729     nG = 0;
 730     nB = 0;
 731
 732     fW = 0;
 733
 734     /* We even have a loop invariant:
 735      * nW+nG+nB == g->nbound
 736      */
 737 #ifdef DEBUG2
 738     fprintf(log, "Starting W loop\n");
 739 #endif
 740     while (nW > 0)
 741     {
 742         /* Find the first white, this will allways be a larger
 743          * number than before, because no nodes are made white
 744          * in the loop
 745          */
 746         if ((fW = first_colour(fW, egcolWhite, g, g->edgeColor)) == -1)
 747         {
 748             gmx_fatal(FARGS, "No WHITE nodes found while nW=%d\n", nW);
 749         }
 750
 751         /* Make the first white node grey */
 752         g->edgeColor[fW] = egcolGrey;
 753         nG++;
 754         nW--;
 755
 756         /* Initial value for the first grey */
 757         fG = fW;
 758 #ifdef DEBUG2
 759         fprintf(log, "Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n", nW, nG, nB, nW + nG + nB);
 760 #endif
 761         while (nG > 0)
 762         {
 763             if ((fG = first_colour(fG, egcolGrey, g, g->edgeColor)) == -1)
 764             {
 765                 gmx_fatal(FARGS, "No GREY nodes found while nG=%d\n", nG);
 766             }
 767
 768             /* Make the first grey node black */
 769             g->edgeColor[fG] = egcolBlack;
 770             nB++;
 771             nG--;
 772
 773             /* Make all the neighbours of this black node grey
 774              * and set their periodicity
 775              */
 776             ng = mk_grey(g->edgeColor, g, &fG, npbcdim, box, x, &nerror);
 777             /* ng is the number of white nodes made grey */
 778             nG += ng;
 779             nW -= ng;
 780         }
 781     }
 782     if (nerror > 0)
 783     {
 784         /* We use a threshold of 0.25*boxSize for generating a fatal error
 785          * to allow removing PBC for systems with periodic molecules.
 786          *
 787          * TODO: Consider a better solution for systems with periodic
 788          *       molecules. Ideally analysis tools should only ask to make
 789          *       non-periodic molecules whole in a system with periodic mols.
 790          */
 791         constexpr real c_relativeDistanceThreshold = 0.25;
 792
 793         int npbcdim = numPbcDimensions(pbcType);
 794         GMX_RELEASE_ASSERT(npbcdim > 0, "Expect PBC with graph");
 795         real minBoxSize = norm(box[XX]);
 796         for (int d = 1; d < npbcdim; d++)
 797         {
 798             minBoxSize = std::min(minBoxSize, norm(box[d]));
 799         }
 800         real maxDistance = maxEdgeLength(*g, pbcType, box, x);
 801         if (maxDistance >= c_relativeDistanceThreshold * minBoxSize)
 802         {
 803             std::string mesg = gmx::formatString(
 804                     "There are inconsistent shifts over periodic boundaries in a molecule type "
 805                     "consisting of %d atoms. The longest distance involved in such interactions is "
 806                     "%.3f nm which is %s half the box length.",
 807                     g->shiftAtomEnd,
 808                     maxDistance,
 809                     maxDistance >= 0.5 * minBoxSize ? "above" : "close to");
 810
 811             switch (g->parts)
 812             {
 813                 case t_graph::BondedParts::MultipleConnected:
 814                     /* Ideally we should check if the long distances are
 815                      * actually between the parts, but that would require
 816                      * a lot of extra code.
 817                      */
 818                     mesg += " This molecule type consists of muliple parts, e.g. monomers, that "
 819                             "are connected by interactions that are not chemical bonds, e.g. "
 820                             "restraints. Such systems can not be treated. The only solution is "
 821                             "increasing the box size.";
 822                     break;
 823                 default:
 824                     mesg += " Either you have excessively large distances between atoms in bonded "
 825                             "interactions or your system is exploding.";
 826             }
 827             gmx_fatal(FARGS, "%s", mesg.c_str());
 828         }
 829
 830         /* The most likely reason for arriving here is a periodic molecule. */
 831
 832         nerror_tot++;
 833         if (nerror_tot <= 100)
 834         {
 835             fprintf(stderr, "There were %d inconsistent shifts. Check your topology\n", nerror);
 836             if (log)
 837             {
 838                 fprintf(log, "There were %d inconsistent shifts. Check your topology\n", nerror);
 839             }
 840         }
 841         if (nerror_tot == 100)
 842         {
 843             fprintf(stderr, "Will stop reporting inconsistent shifts\n");
 844             if (log)
 845             {
 846                 fprintf(log, "Will stop reporting inconsistent shifts\n");
 847             }
 848         }
 849     }
 850 }
 851
 852 /************************************************************
 853  *
 854  *      A C T U A L   S H I F T   C O D E
 855  *
 856  ************************************************************/
 857
 858 void shift_x(const t_graph* g, const matrix box, const rvec x[], rvec x_s[])
 859 {
 860     int j, tx, ty, tz;
 861
 862     GCHECK(g);
 863     const int            g0 = g->edgeAtomBegin;
 864     const int            g1 = g->edgeAtomEnd;
 865     ArrayRef<const IVec> is = g->ishift;
 866
 867     for (j = 0; j < g0; j++)
 868     {
 869         copy_rvec(x[j], x_s[j]);
 870     }
 871
 872     if (g->useScrewPbc)
 873     {
 874         for (j = g0; (j < g1); j++)
 875         {
 876             tx = is[j][XX];
 877             ty = is[j][YY];
 878             tz = is[j][ZZ];
 879
 880             if ((tx > 0 && tx % 2 == 1) || (tx < 0 && -tx % 2 == 1))
 881             {
 882                 x_s[j][XX] = x[j][XX] + tx * box[XX][XX];
 883                 x_s[j][YY] = box[YY][YY] + box[ZZ][YY] - x[j][YY];
 884                 x_s[j][ZZ] = box[ZZ][ZZ] - x[j][ZZ];
 885             }
 886             else
 887             {
 888                 x_s[j][XX] = x[j][XX];
 889             }
 890             x_s[j][YY] = x[j][YY] + ty * box[YY][YY] + tz * box[ZZ][YY];
 891             x_s[j][ZZ] = x[j][ZZ] + tz * box[ZZ][ZZ];
 892         }
 893     }
 894     else if (TRICLINIC(box))
 895     {
 896         for (j = g0; (j < g1); j++)
 897         {
 898             tx = is[j][XX];
 899             ty = is[j][YY];
 900             tz = is[j][ZZ];
 901
 902             x_s[j][XX] = x[j][XX] + tx * box[XX][XX] + ty * box[YY][XX] + tz * box[ZZ][XX];
 903             x_s[j][YY] = x[j][YY] + ty * box[YY][YY] + tz * box[ZZ][YY];
 904             x_s[j][ZZ] = x[j][ZZ] + tz * box[ZZ][ZZ];
 905         }
 906     }
 907     else
 908     {
 909         for (j = g0; (j < g1); j++)
 910         {
 911             tx = is[j][XX];
 912             ty = is[j][YY];
 913             tz = is[j][ZZ];
 914
 915             x_s[j][XX] = x[j][XX] + tx * box[XX][XX];
 916             x_s[j][YY] = x[j][YY] + ty * box[YY][YY];
 917             x_s[j][ZZ] = x[j][ZZ] + tz * box[ZZ][ZZ];
 918         }
 919     }
 920
 921     for (j = g1; j < g->shiftAtomEnd; j++)
 922     {
 923         copy_rvec(x[j], x_s[j]);
 924     }
 925 }
 926
 927 void shift_self(const t_graph& g, const matrix box, rvec x[])
 928 {
 929     int j, tx, ty, tz;
 930
 931     GMX_RELEASE_ASSERT(!g.useScrewPbc, "screw pbc not implemented for shift_self");
 932
 933     const int            g0 = g.edgeAtomBegin;
 934     const int            g1 = g.edgeAtomEnd;
 935     ArrayRef<const IVec> is = g.ishift;
 936
 937 #ifdef DEBUG
 938     fprintf(stderr, "Shifting atoms %d to %d\n", g0, g0 + gn);
 939 #endif
 940     if (TRICLINIC(box))
 941     {
 942         for (j = g0; (j < g1); j++)
 943         {
 944             tx = is[j][XX];
 945             ty = is[j][YY];
 946             tz = is[j][ZZ];
 947
 948             x[j][XX] = x[j][XX] + tx * box[XX][XX] + ty * box[YY][XX] + tz * box[ZZ][XX];
 949             x[j][YY] = x[j][YY] + ty * box[YY][YY] + tz * box[ZZ][YY];
 950             x[j][ZZ] = x[j][ZZ] + tz * box[ZZ][ZZ];
 951         }
 952     }
 953     else
 954     {
 955         for (j = g0; (j < g1); j++)
 956         {
 957             tx = is[j][XX];
 958             ty = is[j][YY];
 959             tz = is[j][ZZ];
 960
 961             x[j][XX] = x[j][XX] + tx * box[XX][XX];
 962             x[j][YY] = x[j][YY] + ty * box[YY][YY];
 963             x[j][ZZ] = x[j][ZZ] + tz * box[ZZ][ZZ];
 964         }
 965     }
 966 }
 967
 968 void shift_self(const t_graph* g, const matrix box, rvec x[])
 969 {
 970     shift_self(*g, box, x);
 971 }
 972
 973 void unshift_x(const t_graph* g, const matrix box, rvec x[], const rvec x_s[])
 974 {
 975     int j, tx, ty, tz;
 976
 977     if (g->useScrewPbc)
 978     {
 979         gmx_incons("screw pbc not implemented (yet) for unshift_x");
 980     }
 981
 982     const int            g0 = g->edgeAtomBegin;
 983     const int            g1 = g->edgeAtomEnd;
 984     ArrayRef<const IVec> is = g->ishift;
 985
 986     for (j = 0; j < g0; j++)
 987     {
 988         copy_rvec(x_s[j], x[j]);
 989     }
 990
 991     if (TRICLINIC(box))
 992     {
 993         for (j = g0; (j < g1); j++)
 994         {
 995             tx = is[j][XX];
 996             ty = is[j][YY];
 997             tz = is[j][ZZ];
 998
 999             x[j][XX] = x_s[j][XX] - tx * box[XX][XX] - ty * box[YY][XX] - tz * box[ZZ][XX];
1000             x[j][YY] = x_s[j][YY] - ty * box[YY][YY] - tz * box[ZZ][YY];
1001             x[j][ZZ] = x_s[j][ZZ] - tz * box[ZZ][ZZ];
1002         }
1003     }
1004     else
1005     {
1006         for (j = g0; (j < g1); j++)
1007         {
1008             tx = is[j][XX];
1009             ty = is[j][YY];
1010             tz = is[j][ZZ];
1011
1012             x[j][XX] = x_s[j][XX] - tx * box[XX][XX];
1013             x[j][YY] = x_s[j][YY] - ty * box[YY][YY];
1014             x[j][ZZ] = x_s[j][ZZ] - tz * box[ZZ][ZZ];
1015         }
1016     }
1017
1018     for (j = g1; j < g->shiftAtomEnd; j++)
1019     {
1020         copy_rvec(x_s[j], x[j]);
1021     }
1022 }
1023
1024 void unshift_self(const t_graph* g, const matrix box, rvec x[])
1025 {
1026     int j, tx, ty, tz;
1027
1028     if (g->useScrewPbc)
1029     {
1030         gmx_incons("screw pbc not implemented for unshift_self");
1031     }
1032
1033     const int            g0 = g->edgeAtomBegin;
1034     const int            g1 = g->edgeAtomEnd;
1035     ArrayRef<const IVec> is = g->ishift;
1036
1037     if (TRICLINIC(box))
1038     {
1039         for (j = g0; (j < g1); j++)
1040         {
1041             tx = is[j][XX];
1042             ty = is[j][YY];
1043             tz = is[j][ZZ];
1044
1045             x[j][XX] = x[j][XX] - tx * box[XX][XX] - ty * box[YY][XX] - tz * box[ZZ][XX];
1046             x[j][YY] = x[j][YY] - ty * box[YY][YY] - tz * box[ZZ][YY];
1047             x[j][ZZ] = x[j][ZZ] - tz * box[ZZ][ZZ];
1048         }
1049     }
1050     else
1051     {
1052         for (j = g0; (j < g1); j++)
1053         {
1054             tx = is[j][XX];
1055             ty = is[j][YY];
1056             tz = is[j][ZZ];
1057
1058             x[j][XX] = x[j][XX] - tx * box[XX][XX];
1059             x[j][YY] = x[j][YY] - ty * box[YY][YY];
1060             x[j][ZZ] = x[j][ZZ] - tz * box[ZZ][ZZ];
1061         }
1062     }
1063 }
1064 #undef GCHECK