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