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