src/gmxlib/mshift.c

   1 /*
   2  * $Id$
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * Gyas ROwers Mature At Cryogenic Speed
  31  */
  32
  33 #include <string.h>
  34 #include "assert.h"
  35 #include "smalloc.h"
  36 #include "fatal.h"
  37 #include "macros.h"
  38 #include "vec.h"
  39 #include "futil.h"
  40 #include "copyrite.h"
  41 #include "mshift.h"
  42 #include "main.h"
  43 #include "pbc.h"
  44
  45 /************************************************************
  46  *
  47  *      S H I F T   U T I L I T I E S
  48  *
  49  ************************************************************/
  50
  51
  52 /************************************************************
  53  *
  54  *      G R A P H   G E N E R A T I O N    C O D E
  55  *
  56  ************************************************************/
  57
  58 static void add_gbond(t_graph *g,t_iatom ia[],int np)
  59 {
  60   int     j,k,l;
  61   atom_id aa,inda;
  62
  63   for(j=0; (j<np); j++) {
  64     aa=ia[j];
  65     inda=aa-g->start;
  66     for(k=0; (k<np); k++)
  67       if (j != k) {
  68         for(l=0; (l<g->nedge[inda]); l++)
  69           if (g->edge[inda][l] == ia[k])
  70             break;
  71         if (l == g->nedge[inda]) {
  72           if (g->nedge[inda] == g->maxedge)
  73             fatal_error(0,"More than %d graph edges per atom (atom %d)\n",
  74                         g->maxedge,aa+1);
  75           g->edge[inda][g->nedge[inda]++]=ia[k];
  76         }
  77       }
  78   }
  79 }
  80
  81 static void mk_igraph(t_graph *g,t_functype ftype[],t_ilist *il,
  82                       int natoms,bool bAll)
  83 {
  84   t_iatom *ia,waterh[3],*iap;
  85   t_iatom tp;
  86   int     i,j,np,nbonded;
  87   int     end;
  88
  89   end=il->nr;
  90   ia=il->iatoms;
  91
  92   i = 0;
  93   while (i < end) {
  94     tp=ftype[ia[0]];
  95     np=interaction_function[tp].nratoms;
  96
  97     if ((ia[1] < natoms) && (interaction_function[tp].flags & IF_GRAPH)) {
  98       if (ia[np] >= natoms)
  99         fatal_error(0,"Molecule in topology has atom numbers below and "
 100                     "above natoms (%d).\n"
 101                     "You are probably trying to use a trajectory which does "
 102                     "not match the first %d atoms of the run input file.\n"
 103                     "You can make a matching run input file with tpbconv.",
 104                     natoms,natoms);
 105       if (tp == F_SETTLE) {
 106         /* Bond all the atoms in the settle */
 107         nbonded = 3;
 108         waterh[0] = ia[1];
 109         waterh[1] = ia[1]+1;
 110         waterh[2] = ia[1]+2;
 111         iap = waterh;
 112       } else {
 113         if (interaction_function[tp].flags & IF_DUMMY)
 114           /* Bond a dummy only to the first constructing atom */
 115           nbonded = 2;
 116         else
 117           nbonded = np;
 118         iap = &(ia[1]);
 119       }
 120       if (bAll)
 121         add_gbond(g,iap,nbonded);
 122       else {
 123         /* Check whether all atoms are bonded now! */
 124         for(j=0; j<nbonded; j++)
 125           if (g->nedge[iap[j]-g->start] == 0)
 126             break;
 127         if (j < nbonded)
 128           add_gbond(g,iap,nbonded);
 129       }
 130     }
 131     ia+=np+1;
 132     i+=np+1;
 133   }
 134 }
 135
 136 void p_graph(FILE *log,char *title,t_graph *g)
 137 {
 138   int i,j;
 139
 140   fprintf(log,"graph:  %s\n",title);
 141   fprintf(log,"maxedge:%d\n",g->maxedge);
 142   fprintf(log,"nnodes: %d\n",g->nnodes);
 143   fprintf(log,"nbound: %d\n",g->nbound);
 144   fprintf(log,"start:  %d\n",g->start);
 145   fprintf(log,"end:    %d\n",g->end);
 146   fprintf(log," atom shiftx shifty shiftz nedg   e1   e2 etc.\n");
 147   for(i=0; (i<g->nnodes); i++)
 148     if (g->nedge[i] > 0) {
 149         fprintf(log,"%5d%7d%7d%7d%5d",g->start+i+1,g->ishift[i][XX],g->ishift[i][YY],
 150                 g->ishift[i][ZZ],g->nedge[i]);
 151       for(j=0; (j<g->nedge[i]); j++)
 152         fprintf(log,"%5u",g->edge[i][j]+1);
 153       fprintf(log,"\n");
 154     }
 155   fflush(log);
 156 }
 157
 158 static void calc_1se(t_graph *g,t_ilist *il,t_functype ftype[],
 159                      short nbond[],int natoms)
 160 {
 161   int     k,nratoms,nbonded,tp,end,j;
 162   t_iatom *ia,iaa;
 163
 164   end=il->nr;
 165
 166   ia=il->iatoms;
 167   for(j=0; (j<end); j+=nratoms+1,ia+=nratoms+1) {
 168     tp      = ftype[ia[0]];
 169     nratoms = interaction_function[tp].nratoms;
 170
 171     if (tp == F_SETTLE) {
 172       iaa          = ia[1];
 173       if (iaa<natoms) {
 174         nbond[iaa]   = 2;
 175         nbond[iaa+1] = 2;
 176         nbond[iaa+2] = 2;
 177         g->start     = min(g->start,iaa);
 178         g->end       = max(g->end,iaa+2);
 179       }
 180     } else {
 181       if (interaction_function[tp].flags & IF_DUMMY)
 182         /* Bond a dummy only to the first constructing atom */
 183         nbonded = 2;
 184       else
 185         nbonded = nratoms;
 186       for(k=0; (k<nbonded); k++) {
 187         iaa=ia[k+1];
 188         if (iaa<natoms) {
 189           g->start=min(g->start,iaa);
 190           g->end  =max(g->end,  iaa);
 191           if (interaction_function[tp].flags & IF_GRAPH)
 192             nbond[iaa]++;
 193         }
 194       }
 195     }
 196   }
 197 }
 198
 199 static void calc_start_end(t_graph *g,t_idef *idef,int natoms)
 200 {
 201   short *nbond;
 202   int   i,nnb;
 203
 204   g->start=natoms;
 205   g->end=0;
 206
 207   snew(nbond,natoms);
 208   for(i=0; (i<F_NRE); i++) {
 209     if (interaction_function[i].flags & IF_GRAPH)
 210       calc_1se(g,&idef->il[i],idef->functype,nbond,natoms);
 211   }
 212
 213   nnb=0;
 214   for(i=g->start; (i<=g->end); i++)
 215     nnb=max(nnb,nbond[i]);
 216   if (stdlog)
 217     fprintf(stdlog,"Max number of graph edges per atom is %d\n",nnb);
 218
 219   sfree(nbond);
 220
 221   g->maxedge=nnb+6;
 222 }
 223
 224 t_graph *mk_graph(t_idef *idef,int natoms,bool bShakeOnly,bool bSettle)
 225 {
 226   t_graph *g;
 227   int     i;
 228
 229   snew(g,1);
 230
 231   calc_start_end(g,idef,natoms);
 232
 233   if (g->start >= g->end) {
 234     g->nnodes=0;
 235   }
 236   else {
 237     g->nnodes=g->end-g->start+1;
 238     snew(g->ishift,g->nnodes);
 239     snew(g->nedge,g->nnodes);
 240
 241     /* To prevent malloc problems with many small arrays, using realloc,
 242      * we allocate some more memory, and divide it ourselves.
 243      * We calculate pointers... (Yuck Yuck)
 244      */
 245     if (debug)
 246       fprintf(debug,"MSHIFT: nnodes=%d, maxedge=%d\n",g->nnodes,g->maxedge);
 247     snew(g->edge,g->nnodes);
 248     snew(g->edge[0],g->maxedge*g->nnodes);
 249
 250     for(i=1; (i<g->nnodes); i++)
 251       g->edge[i]=g->edge[i-1]+g->maxedge;
 252
 253     if (!bShakeOnly) {
 254       /* First add all the real bonds: they should determine the molecular
 255        * graph.
 256        */
 257       for(i=0; (i<F_NRE); i++)
 258         if (interaction_function[i].flags & IF_GRAPH)
 259           mk_igraph(g,idef->functype,&(idef->il[i]),natoms,TRUE);
 260       /* Then add all the other interactions in fixed lists, but first
 261        * check to see what's there already.
 262        */
 263       for(i=0; (i<F_NRE); i++) {
 264         if (!(interaction_function[i].flags & IF_GRAPH)) {
 265           mk_igraph(g,idef->functype,&(idef->il[i]),natoms,FALSE);
 266         }
 267       }
 268     }
 269     else {
 270       /* This is a special thing used in grompp to generate shake-blocks */
 271       mk_igraph(g,idef->functype,&(idef->il[F_SHAKE]),natoms,TRUE);
 272       if (bSettle)
 273         mk_igraph(g,idef->functype,&(idef->il[F_SETTLE]),natoms,TRUE);
 274     }
 275     g->nbound=0;
 276     for(i=0; (i<g->nnodes); i++)
 277       if (g->nedge[i] > 0)
 278         g->nbound++;
 279   }
 280 #ifdef DEBUG
 281   p_graph(stdlog,"graph",g);
 282 #endif
 283
 284   g->negc = 0;
 285   g->egc = NULL;
 286
 287   return g;
 288 }
 289
 290 void done_graph(t_graph *g)
 291 {
 292   if (g->nnodes > 0) {
 293     sfree(g->ishift);
 294     sfree(g->nedge);
 295     /* This is malloced in a NASTY way, see above */
 296     sfree(g->edge[0]);
 297     sfree(g->edge);
 298     sfree(g->egc);
 299   }
 300 }
 301
 302 /************************************************************
 303  *
 304  *      S H I F T   C A L C U L A T I O N   C O D E
 305  *
 306  ************************************************************/
 307
 308 static void mk_1shift_tric(matrix box,rvec hbox,rvec xi,rvec xj,int *mi,int *mj)
 309 {
 310   /* Calculate periodicity for triclinic box... */
 311   int  m,d;
 312   rvec dx;
 313
 314   rvec_sub(xi,xj,dx);
 315
 316   for(m=DIM-1; (m>=0); m--) {
 317     /* If dx < hbox, then xj will be reduced by box, so that
 318      * xi - xj will be bigger
 319      */
 320     if (dx[m] < -hbox[m]) {
 321       mj[m]=mi[m]-1;
 322       for(d=m-1; d>=0; d--)
 323         dx[d]+=box[m][d];
 324     } else if (dx[m] >= hbox[m]) {
 325       mj[m]=mi[m]+1;
 326       for(d=m-1; d>=0; d--)
 327         dx[d]-=box[m][d];
 328     } else
 329       mj[m]=mi[m];
 330   }
 331 }
 332
 333 static void mk_1shift(rvec hbox,rvec xi,rvec xj,int *mi,int *mj)
 334 {
 335   /* Calculate periodicity for rectangular box... */
 336   int  m;
 337   rvec dx;
 338
 339   rvec_sub(xi,xj,dx);
 340
 341   for(m=0; (m<DIM); m++) {
 342     /* If dx < hbox, then xj will be reduced by box, so that
 343      * xi - xj will be bigger
 344      */
 345     if (dx[m] < -hbox[m])
 346       mj[m]=mi[m]-1;
 347     else if (dx[m] >= hbox[m])
 348       mj[m]=mi[m]+1;
 349     else
 350       mj[m]=mi[m];
 351   }
 352 }
 353
 354 static int mk_grey(FILE *log,int nnodes,egCol egc[],t_graph *g,int *AtomI,
 355                    matrix box,rvec x[],int *nerror)
 356 {
 357   int      m,j,ng,ai,aj,g0;
 358   rvec     hbox;
 359   bool     bTriclinic;
 360   ivec     is_aj;
 361
 362   for(m=0; (m<DIM); m++)
 363     hbox[m]=box[m][m]*0.5;
 364   bTriclinic = TRICLINIC(box);
 365
 366   ng=0;
 367   ai=*AtomI;
 368
 369   g0=g->start;
 370   /* Loop over all the bonds */
 371   for(j=0; (j<g->nedge[ai]); j++) {
 372     aj=g->edge[ai][j]-g0;
 373     /* If there is a white one, make it gray and set pbc */
 374     if (bTriclinic)
 375       mk_1shift_tric(box,hbox,x[g0+ai],x[g0+aj],g->ishift[ai],is_aj);
 376     else
 377       mk_1shift(hbox,x[g0+ai],x[g0+aj],g->ishift[ai],is_aj);
 378
 379     if (egc[aj] == egcolWhite) {
 380       if (aj < *AtomI)
 381         *AtomI = aj;
 382       egc[aj] = egcolGrey;
 383
 384       copy_ivec(is_aj,g->ishift[aj]);
 385
 386       ng++;
 387     }
 388     else if ((is_aj[XX] != g->ishift[aj][XX]) ||
 389              (is_aj[YY] != g->ishift[aj][YY]) ||
 390              (is_aj[ZZ] != g->ishift[aj][ZZ]))
 391              {
 392       (*nerror)++;
 393     }
 394   }
 395   return ng;
 396 }
 397
 398 static int first_colour(int fC,egCol Col,t_graph *g,egCol egc[])
 399 /* Return the first node with colour Col starting at fC.
 400  * return -1 if none found.
 401  */
 402 {
 403   int i;
 404
 405   for(i=fC; (i<g->nnodes); i++)
 406     if ((g->nedge[i] > 0) && (egc[i]==Col))
 407       return i;
 408
 409   return -1;
 410 }
 411
 412 void mk_mshift(FILE *log,t_graph *g,matrix box,rvec x[])
 413 {
 414   int    ng,nnodes,i;
 415   int    nW,nG,nB;              /* Number of Grey, Black, White */
 416   int    fW,fG;                 /* First of each category       */
 417   int    nerror=0;
 418
 419   /* This puts everything in the central box, that is does not move it
 420    * at all. If we return without doing this for a system without bonds
 421    * (i.e. only settles) all water molecules are moved to the opposite octant
 422    */
 423   for(i=0; (i<g->nnodes); i++) {
 424       g->ishift[i][XX]=g->ishift[i][YY]=g->ishift[i][ZZ]=0;
 425   }
 426
 427   if (!g->nbound)
 428     return;
 429
 430   nnodes=g->nnodes;
 431   if (nnodes > g->negc) {
 432     g->negc = nnodes;
 433     srenew(g->egc,g->negc);
 434   }
 435   memset(g->egc,0,(size_t)(nnodes*sizeof(g->egc[0])));
 436
 437   nW=g->nbound;
 438   nG=0;
 439   nB=0;
 440
 441   fW=0;
 442
 443   /* We even have a loop invariant:
 444    * nW+nG+nB == g->nbound
 445    */
 446 #ifdef DEBUG2
 447   fprintf(log,"Starting W loop\n");
 448 #endif
 449   while (nW > 0) {
 450     /* Find the first white, this will allways be a larger
 451      * number than before, because no nodes are made white
 452      * in the loop
 453      */
 454     if ((fW=first_colour(fW,egcolWhite,g,g->egc)) == -1)
 455       fatal_error(0,"No WHITE nodes found while nW=%d\n",nW);
 456
 457     /* Make the first white node grey */
 458     g->egc[fW]=egcolGrey;
 459     nG++;
 460     nW--;
 461
 462     /* Initial value for the first grey */
 463     fG=fW;
 464 #ifdef DEBUG2
 465     fprintf(log,"Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n",
 466             nW,nG,nB,nW+nG+nB);
 467 #endif
 468     while (nG > 0) {
 469       if ((fG=first_colour(fG,egcolGrey,g,g->egc)) == -1)
 470         fatal_error(0,"No GREY nodes found while nG=%d\n",nG);
 471
 472       /* Make the first grey node black */
 473       g->egc[fG]=egcolBlack;
 474       nB++;
 475       nG--;
 476
 477       /* Make all the neighbours of this black node grey
 478        * and set their periodicity
 479        */
 480       ng=mk_grey(log,nnodes,g->egc,g,&fG,box,x,&nerror);
 481       /* ng is the number of white nodes made grey */
 482       nG+=ng;
 483       nW-=ng;
 484     }
 485   }
 486   if (nerror > 0)
 487     fprintf(log,"There were %d inconsistent shifts. Check your topology\n",
 488             nerror);
 489 }
 490
 491 /************************************************************
 492  *
 493  *      A C T U A L   S H I F T   C O D E
 494  *
 495  ************************************************************/
 496
 497 void shift_atom(t_graph *g,matrix box,rvec x[],rvec x_s[],atom_id ai)
 498 {
 499   int tx,ty,tz;
 500   tx=(g->ishift[ai-g->start])[XX];
 501   ty=(g->ishift[ai-g->start])[YY];
 502   tz=(g->ishift[ai-g->start])[ZZ];
 503
 504   x_s[ai][XX]=x[ai][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
 505   x_s[ai][YY]=x[ai][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
 506   x_s[ai][ZZ]=x[ai][ZZ]+tz*box[ZZ][ZZ];
 507 }
 508
 509 void unshift_atom(t_graph *g,matrix box,rvec x[],rvec x_s[],atom_id ai)
 510 {
 511   int tx,ty,tz;
 512   tx=(g->ishift[ai-g->start])[XX];
 513   ty=(g->ishift[ai-g->start])[YY];
 514   tz=(g->ishift[ai-g->start])[ZZ];
 515
 516   x_s[ai][XX]=x[ai][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
 517   x_s[ai][YY]=x[ai][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
 518   x_s[ai][ZZ]=x[ai][ZZ]-tz*box[ZZ][ZZ];
 519 }
 520
 521 void shift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
 522 {
 523   ivec *is;
 524   int      g0,gn;
 525   int      i,j,tx,ty,tz;
 526
 527   g0=g->start;
 528   gn=g->nnodes;
 529   is=g->ishift;
 530
 531   if(TRICLINIC(box)) {
 532      for(i=0,j=g0; (i<gn); i++,j++) {
 533          tx=is[i][XX];
 534          ty=is[i][YY];
 535          tz=is[i][ZZ];
 536
 537          x_s[j][XX]=x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
 538          x_s[j][YY]=x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
 539          x_s[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
 540      }
 541   } else {
 542      for(i=0,j=g0; (i<gn); i++,j++) {
 543          tx=is[i][XX];
 544          ty=is[i][YY];
 545          tz=is[i][ZZ];
 546
 547          x_s[j][XX]=x[j][XX]+tx*box[XX][XX];
 548          x_s[j][YY]=x[j][YY]+ty*box[YY][YY];
 549          x_s[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
 550      }
 551   }
 552
 553 }
 554
 555 void shift_self(t_graph *g,matrix box,rvec x[])
 556 {
 557   ivec *is;
 558   int      g0,gn;
 559   int      i,j,tx,ty,tz;
 560
 561   g0=g->start;
 562   gn=g->nnodes;
 563   is=g->ishift;
 564
 565 #ifdef DEBUG
 566   fprintf(stderr,"Shifting atoms %d to %d\n",g0,g0+gn);
 567 #endif
 568   if(TRICLINIC(box)) {
 569       for(i=0,j=g0; (i<gn); i++,j++) {
 570           tx=is[i][XX];
 571           ty=is[i][YY];
 572           tz=is[i][ZZ];
 573
 574           x[j][XX]=x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
 575           x[j][YY]=x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
 576           x[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
 577       }
 578   } else {
 579       for(i=0,j=g0; (i<gn); i++,j++) {
 580           tx=is[i][XX];
 581           ty=is[i][YY];
 582           tz=is[i][ZZ];
 583
 584           x[j][XX]=x[j][XX]+tx*box[XX][XX];
 585           x[j][YY]=x[j][YY]+ty*box[YY][YY];
 586           x[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
 587       }
 588   }
 589
 590 }
 591
 592 void unshift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
 593 {
 594   ivec *is;
 595   int      g0,gn;
 596   int      i,j,tx,ty,tz;
 597
 598   g0=g->start;
 599   gn=g->nnodes;
 600   is=g->ishift;
 601   if(TRICLINIC(box)) {
 602       for(i=0,j=g0; (i<gn); i++,j++) {
 603           tx=is[i][XX];
 604           ty=is[i][YY];
 605           tz=is[i][ZZ];
 606
 607           x[j][XX]=x_s[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
 608           x[j][YY]=x_s[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
 609           x[j][ZZ]=x_s[j][ZZ]-tz*box[ZZ][ZZ];
 610       }
 611   } else {
 612       for(i=0,j=g0; (i<gn); i++,j++) {
 613           tx=is[i][XX];
 614           ty=is[i][YY];
 615           tz=is[i][ZZ];
 616
 617           x[j][XX]=x_s[j][XX]-tx*box[XX][XX];
 618           x[j][YY]=x_s[j][YY]-ty*box[YY][YY];
 619           x[j][ZZ]=x_s[j][ZZ]-tz*box[ZZ][ZZ];
 620       }
 621   }
 622 }
 623
 624 void unshift_self(t_graph *g,matrix box,rvec x[])
 625 {
 626   ivec *is;
 627   int g0,gn;
 628   int i,j,tx,ty,tz;
 629
 630   g0=g->start;
 631   gn=g->nnodes;
 632   is=g->ishift;
 633   if(TRICLINIC(box)) {
 634       for(i=0,j=g0; (i<gn); i++,j++) {
 635           tx=is[i][XX];
 636           ty=is[i][YY];
 637           tz=is[i][ZZ];
 638
 639           x[j][XX]=x[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
 640           x[j][YY]=x[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
 641           x[j][ZZ]=x[j][ZZ]-tz*box[ZZ][ZZ];
 642       }
 643   } else {
 644       for(i=0,j=g0; (i<gn); i++,j++) {
 645           tx=is[i][XX];
 646           ty=is[i][YY];
 647           tz=is[i][ZZ];
 648
 649           x[j][XX]=x[j][XX]-tx*box[XX][XX];
 650           x[j][YY]=x[j][YY]-ty*box[YY][YY];
 651           x[j][ZZ]=x[j][ZZ]-tz*box[ZZ][ZZ];
 652       }
 653   }
 654 }
 655
 656 #ifdef DEBUGMSHIFT
 657 void main(int argc,char *argv[])
 658 {
 659   FILE         *out;
 660   t_args       targ;
 661   t_topology   top;
 662   t_statheader sh;
 663   rvec         *x;
 664   ivec         *mshift;
 665   matrix       box;
 666
 667   t_graph      *g;
 668   int          i,idum,pid;
 669   real         rdum;
 670
 671   CopyRight(stderr,argv[0]);
 672   parse_common_args(&argc,argv,&targ,PCA_NEED_INOUT,NULL);
 673   if (argc > 1)
 674     pid=atoi(argv[1]);
 675   else
 676     pid=0;
 677
 678   read_status_header(targ.infile,&sh);
 679   snew(x,sh.natoms);
 680   snew(mshift,sh.natoms);
 681
 682   fprintf(stderr,"Reading Status %s\n",targ.infile);
 683   read_status(targ.infile,&idum,&rdum,&rdum,NULL,
 684               box,NULL,NULL,&idum,x,NULL,NULL,&idum,NULL,&top);
 685
 686   fprintf(stderr,"Making Graph Structure...\n");
 687   g=mk_graph(&(top.idef),top.atoms.nr,FALSE,FALSE);
 688
 689   out=ffopen(targ.outfile,"w");
 690
 691   fprintf(stderr,"Making Shift...\n");
 692   mk_mshift(out,g,box,x,mshift);
 693
 694   p_graph(out,"In Den Haag daar woont een graaf...",g);
 695   fclose(out);
 696 }
 697 #endif
 698