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