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