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