src/gromacs/gmxlib/mshift.c

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