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