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