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