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