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