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