Project

General

Profile

mshift.c

fixed src/gmxlib/mshift.c - Berk Hess, 12/11/2008 03:26 PM

 
1
/*
2
 * $Id: mshift.c,v 1.62.2.3 2008/12/11 14:24:43 hess Exp $
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.2.0
11
 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12
 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13
 * Copyright (c) 2001-2004, The GROMACS development team,
14
 * check out http://www.gromacs.org for more information.
15

16
 * This program is free software; you can redistribute it and/or
17
 * modify it under the terms of the GNU General Public License
18
 * as published by the Free Software Foundation; either version 2
19
 * of the License, or (at your option) any later version.
20
 * 
21
 * If you want to redistribute modifications, please consider that
22
 * scientific software is very special. Version control is crucial -
23
 * bugs must be traceable. We will be happy to consider code for
24
 * inclusion in the official distribution, but derived work must not
25
 * be called official GROMACS. Details are found in the README & COPYING
26
 * files - if they are missing, get the official version at www.gromacs.org.
27
 * 
28
 * To help us fund GROMACS development, we humbly ask that you cite
29
 * the papers on the package - you can find them in the top README file.
30
 * 
31
 * For more info, check our website at http://www.gromacs.org
32
 * 
33
 * And Hey:
34
 * GROningen Mixture of Alchemy and Childrens' Stories
35
 */
36
#ifdef HAVE_CONFIG_H
37
#include <config.h>
38
#endif
39

    
40
#include <string.h>
41
#include "smalloc.h"
42
#include "gmx_fatal.h"
43
#include "macros.h"
44
#include "vec.h"
45
#include "futil.h"
46
#include "copyrite.h"
47
#include "mshift.h"
48
#include "main.h"
49
#include "pbc.h"
50

    
51
/************************************************************
52
 *
53
 *      S H I F T   U T I L I T I E S
54
 *
55
 ************************************************************/
56
 
57

    
58
/************************************************************
59
 *
60
 *      G R A P H   G E N E R A T I O N    C O D E
61
 *
62
 ************************************************************/
63

    
64
static void add_gbond(t_graph *g,atom_id a0,atom_id a1)
65
{
66
  int     i;
67
  atom_id inda0,inda1;
68
  bool    bFound;
69

    
70
  inda0 = a0 - g->start;
71
  inda1 = a1 - g->start;
72
  bFound = FALSE;
73
  /* Search for a direct edge between a0 and a1.
74
   * All egdes are bidirectional, so we only need to search one way.
75
   */
76
  for(i=0; (i<g->nedge[inda0] && !bFound); i++) {
77
    bFound = (g->edge[inda0][i] == a1);
78
  }
79

    
80
  if (!bFound) {
81
    g->edge[inda0][g->nedge[inda0]++] = a1;
82
    g->edge[inda1][g->nedge[inda1]++] = a0;
83
  }
84
}
85

    
86
static void mk_igraph(t_graph *g,int ftype,t_ilist *il,
87
                      int at_start,int at_end,
88
                      int *part)
89
{
90
  t_iatom *ia;
91
  int     i,j,np;
92
  int     end;
93

    
94
  end=il->nr;
95
  ia=il->iatoms;
96

    
97
  i = 0;
98
  while (i < end) {
99
    np = interaction_function[ftype].nratoms;
100
    
101
    if (ia[1] >= at_start && ia[1] < at_end) {
102
      if (ia[np] >= at_end)
103
        gmx_fatal(FARGS,
104
                  "Molecule in topology has atom numbers below and "
105
                  "above natoms (%d).\n"
106
                  "You are probably trying to use a trajectory which does "
107
                  "not match the first %d atoms of the run input file.\n"
108
                  "You can make a matching run input file with tpbconv.",
109
                  at_end,at_end);
110
      if (ftype == F_SETTLE) {
111
        /* Bond all the atoms in the settle */
112
        add_gbond(g,ia[1],ia[1]+1);
113
        add_gbond(g,ia[1],ia[1]+2);
114
      } else if (part == NULL) {
115
        /* Simply add this bond */
116
        for(j=1; j<np; j++) {
117
          add_gbond(g,ia[j],ia[j+1]);
118
        }
119
      } else {
120
        /* Add this bond when it connects two unlinked parts of the graph */
121
        for(j=1; j<np; j++) {
122
          if (part[ia[j]] != part[ia[j+1]]) {
123
            add_gbond(g,ia[j],ia[j+1]);
124
          }
125
        }
126
      }
127
    }
128
    ia+=np+1;
129
    i+=np+1;
130
  }
131
}
132

    
133
static void g_error(int line,char *file)
134
{
135
  gmx_fatal(FARGS,"Tring to print non existant graph (file %s, line %d)",
136
              file,line);
137
}
138
#define GCHECK(g) if (g == NULL) g_error(__LINE__,__FILE__)
139

    
140
void p_graph(FILE *log,char *title,t_graph *g)
141
{
142
  int  i,j;
143
  char *cc[egcolNR] = { "W", "G", "B" };
144
  
145
  GCHECK(g);
146
  fprintf(log,"graph:  %s\n",title);
147
  fprintf(log,"nnodes: %d\n",g->nnodes);
148
  fprintf(log,"nbound: %d\n",g->nbound);
149
  fprintf(log,"start:  %d\n",g->start);
150
  fprintf(log,"end:    %d\n",g->end);
151
  fprintf(log," atom shiftx shifty shiftz C nedg    e1    e2 etc.\n");
152
  for(i=0; (i<g->nnodes); i++)
153
    if (g->nedge[i] > 0) {
154
      fprintf(log,"%5d%7d%7d%7d %1s%5d",g->start+i+1,
155
              g->ishift[i][XX],g->ishift[i][YY],
156
              g->ishift[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[iaa+1] += 1;
183
        nbond[iaa+2] += 1;
184
        g->start      = min(g->start,iaa);
185
        g->end        = max(g->end,iaa+2);
186
      }
187
    } else {
188
      for(k=1; (k<=nratoms); k++) {
189
        iaa=ia[k];
190
        if (iaa >= at_start && iaa < at_end) {
191
          g->start=min(g->start,iaa);
192
          g->end  =max(g->end,  iaa);
193
          /*if (interaction_function[tp].flags & IF_CHEMBOND)*/
194
          nbond[iaa]++;
195
        }
196
      }
197
    }
198
  }
199
}
200

    
201
static int calc_start_end(FILE *fplog,t_graph *g,t_ilist il[],
202
                          int at_start,int at_end,
203
                          int nbond[])
204
{
205
  int   i,nnb,nbtot;
206
  
207
  g->start=at_end;
208
  g->end=0;
209

    
210
  /* First add all the real bonds: they should determine the molecular 
211
   * graph.
212
   */
213
  for(i=0; (i<F_NRE); i++)
214
    if (interaction_function[i].flags & IF_CHEMBOND)
215
      calc_1se(g,i,&il[i],nbond,at_start,at_end);
216
  /* Then add all the other interactions in fixed lists, but first
217
   * check to see what's there already.
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
    }
223
  }
224
  
225
  nnb   = 0;
226
  nbtot = 0;
227
  for(i=g->start; (i<=g->end); i++) {
228
    nbtot += nbond[i];
229
    nnb    = max(nnb,nbond[i]);
230
  }
231
  if (fplog) {
232
    fprintf(fplog,"Max number of connections per atom is %d\n",nnb);
233
    fprintf(fplog,"Total number of connections is %d\n",nbtot);
234
  }
235
  return nbtot;
236
}
237

    
238

    
239

    
240
static void compact_graph(FILE *fplog,t_graph *g)
241
{
242
  int i,j,n,max_nedge;
243
  atom_id *e;
244

    
245
  max_nedge = 0;
246
  n = 0;
247
  for(i=0; i<g->nnodes; i++) {
248
    srenew(g->edge[i],g->nedge[i]);
249
    max_nedge = max(max_nedge,g->nedge[i]);
250
    n+=g->nedge[i];
251
  }
252
  if (fplog) {
253
    fprintf(fplog,"Max number of graph edges per atom is %d\n",
254
            max_nedge);
255
    fprintf(fplog,"Total number of graph edges is %d\n",n);
256
  }
257
}
258

    
259
static bool determine_graph_parts(t_graph *g,int *part)
260
{
261
  int  i,e;
262
  int  nchanged;
263
  atom_id at_i,*at_i2;
264
  bool bMultiPart;
265

    
266
  /* Initialize the part array with all entries different */
267
  for(at_i=g->start; at_i<g->end; at_i++) {
268
    part[at_i] = at_i;
269
  }
270

    
271
  /* Loop over the graph until the part array is fixed */
272
  do {
273
    bMultiPart = FALSE;
274
    nchanged = 0;
275
    for(i=0; (i<g->nnodes); i++) {
276
      at_i  = g->start + i;
277
      at_i2 = g->edge[i];
278
      for(e=0; e<g->nedge[i]; e++) {
279
        /* Set part for both nodes to the minimum */
280
        if (part[at_i2[e]] > part[at_i]) {
281
          part[at_i2[e]] = part[at_i];
282
          nchanged++;
283
        } else if (part[at_i2[e]] < part[at_i]) {
284
          part[at_i] = part[at_i2[e]];
285
          nchanged++;
286
        }
287
      }
288
      if (part[at_i] != part[g->start]) {
289
        bMultiPart = TRUE;
290
      }
291
    }
292
    if (debug) {
293
      fprintf(debug,"graph part[] nchanged=%d, bMultiPart=%d\n",
294
              nchanged,bMultiPart);
295
    }
296
  } while (nchanged > 0);
297

    
298
  return bMultiPart;
299
}
300

    
301
void mk_graph_ilist(FILE *fplog,
302
                    t_ilist *ilist,int at_start,int at_end,
303
                    bool bShakeOnly,bool bSettle,
304
                    t_graph *g)
305
{
306
  int     *nbond;
307
  int     i,nbtot;
308
  bool    bMultiPart;
309

    
310
  snew(nbond,at_end);
311
  nbtot = calc_start_end(fplog,g,ilist,at_start,at_end,nbond);
312
  
313
  if (g->start >= g->end) {
314
    g->nnodes = 0;
315
    g->nbound = 0;
316
  }
317
  else {
318
    g->nnodes = g->end - g->start + 1;
319
    snew(g->ishift,g->nnodes);
320
    snew(g->nedge,g->nnodes);
321
    snew(g->edge,g->nnodes);
322
    for(i=0; (i<g->nnodes); i++)
323
      snew(g->edge[i],nbond[g->start+i]);
324

    
325
    if (!bShakeOnly) {
326
      /* First add all the real bonds: they should determine the molecular 
327
       * graph.
328
       */
329
      for(i=0; (i<F_NRE); i++)
330
        if (interaction_function[i].flags & IF_CHEMBOND)
331
          mk_igraph(g,i,&(ilist[i]),at_start,at_end,NULL);
332

    
333
      /* Determine of which separated parts the IF_CHEMBOND graph consists.
334
       * Store the parts in the nbond array.
335
       */
336
      bMultiPart = determine_graph_parts(g,nbond);
337

    
338
      if (bMultiPart) {
339
        /* Then add all the other interactions in fixed lists,
340
         * but only when they connect parts of the graph
341
         * that are not connected through IF_CHEMBOND interactions.
342
         */         
343
        for(i=0; (i<F_NRE); i++) {
344
          if (!(interaction_function[i].flags & IF_CHEMBOND)) {
345
            mk_igraph(g,i,&(ilist[i]),at_start,at_end,nbond);
346
          }
347
        }
348
      }
349
      
350
      /* Removed all the unused space from the edge array */
351
      compact_graph(fplog,g);
352
    }
353
    else {
354
      /* This is a special thing used in splitter.c to generate shake-blocks */
355
      mk_igraph(g,F_CONSTR,&(ilist[F_CONSTR]),at_start,at_end,NULL);
356
      if (bSettle)
357
        mk_igraph(g,F_SETTLE,&(ilist[F_SETTLE]),at_start,at_end,NULL);
358
    }
359
    g->nbound = 0;
360
    for(i=0; (i<g->nnodes); i++)
361
      if (g->nedge[i] > 0)
362
        g->nbound++;
363
  }
364

    
365
  g->negc = 0;
366
  g->egc = NULL;
367

    
368
  sfree(nbond);
369

    
370
  if (gmx_debug_at)
371
    p_graph(debug,"graph",g);
372
}
373

    
374
t_graph *mk_graph(FILE *fplog,
375
                  t_idef *idef,int at_start,int at_end,
376
                  bool bShakeOnly,bool bSettle)
377
{
378
  t_graph *g;
379

    
380
  snew(g,1);
381

    
382
  mk_graph_ilist(fplog,idef->il,at_start,at_end,bShakeOnly,bSettle,g);
383

    
384
  return g;
385
}
386

    
387
void done_graph(t_graph *g)
388
{
389
  int i;
390
  
391
  GCHECK(g);
392
  if (g->nnodes > 0) {
393
    sfree(g->ishift);
394
    sfree(g->nedge);
395
    for(i=0; (i<g->nnodes); i++)
396
      sfree(g->edge[i]);
397
    sfree(g->edge);
398
    sfree(g->egc);
399
  }
400
}
401

    
402
/************************************************************
403
 *
404
 *      S H I F T   C A L C U L A T I O N   C O D E
405
 *
406
 ************************************************************/
407

    
408
static void mk_1shift_tric(int npbcdim,matrix box,rvec hbox,
409
                           rvec xi,rvec xj,int *mi,int *mj)
410
{
411
  /* Calculate periodicity for triclinic box... */
412
  int  m,d;
413
  rvec dx;
414
  
415
  rvec_sub(xi,xj,dx);
416

    
417
  mj[ZZ] = 0;
418
  for(m=npbcdim-1; (m>=0); m--) {
419
    /* If dx < hbox, then xj will be reduced by box, so that
420
     * xi - xj will be bigger
421
     */
422
    if (dx[m] < -hbox[m]) {
423
      mj[m]=mi[m]-1;
424
      for(d=m-1; d>=0; d--)
425
        dx[d]+=box[m][d];
426
    } else if (dx[m] >= hbox[m]) {
427
      mj[m]=mi[m]+1;
428
      for(d=m-1; d>=0; d--)
429
        dx[d]-=box[m][d];
430
    } else
431
      mj[m]=mi[m];
432
  }
433
}
434

    
435
static void mk_1shift(int npbcdim,rvec hbox,rvec xi,rvec xj,int *mi,int *mj)
436
{
437
  /* Calculate periodicity for rectangular box... */
438
  int  m;
439
  rvec dx;
440
  
441
  rvec_sub(xi,xj,dx);
442

    
443
  mj[ZZ] = 0;
444
  for(m=0; (m<npbcdim); m++) {
445
    /* If dx < hbox, then xj will be reduced by box, so that
446
     * xi - xj will be bigger
447
     */
448
    if (dx[m] < -hbox[m])
449
      mj[m]=mi[m]-1;
450
    else if (dx[m] >= hbox[m])
451
      mj[m]=mi[m]+1;
452
    else
453
      mj[m]=mi[m];
454
  }
455
}
456

    
457
static void mk_1shift_screw(matrix box,rvec hbox,
458
                            rvec xi,rvec xj,int *mi,int *mj)
459
{
460
  /* Calculate periodicity for rectangular box... */
461
  int  signi,m;
462
  rvec dx;
463

    
464
  if ((mi[XX] > 0 &&  mi[XX] % 2 == 1) ||
465
      (mi[XX] < 0 && -mi[XX] % 2 == 1)) {
466
    signi = -1;
467
  } else {
468
    signi =  1;
469
  }
470

    
471
  rvec_sub(xi,xj,dx);
472

    
473
  if (dx[XX] < -hbox[XX])
474
    mj[XX] = mi[XX] - 1;
475
  else if (dx[XX] >= hbox[XX])
476
    mj[XX] = mi[XX] + 1;
477
  else
478
    mj[XX] = mi[XX];
479
  if (mj[XX] != mi[XX]) {
480
    /* Rotate */
481
    dx[YY] = xi[YY] - (box[YY][YY] + box[ZZ][YY] - xj[YY]);
482
    dx[ZZ] = xi[ZZ] - (box[ZZ][ZZ]               - xj[ZZ]);
483
  }
484
  for(m=1; (m<DIM); m++) {
485
    /* The signs are taken such that we can first shift x and rotate
486
     * and then shift y and z.
487
     */
488
    if (dx[m] < -hbox[m])
489
      mj[m] = mi[m] - signi;
490
    else if (dx[m] >= hbox[m])
491
      mj[m] = mi[m] + signi;
492
    else
493
      mj[m] = mi[m];
494
  }
495
}
496

    
497
static int mk_grey(FILE *log,int nnodes,egCol egc[],t_graph *g,int *AtomI,
498
                   int npbcdim,matrix box,rvec x[],int *nerror)
499
{
500
  int      m,j,ng,ai,aj,g0;
501
  rvec     dx,hbox;
502
  bool     bTriclinic;
503
  ivec     is_aj;
504
  t_pbc    pbc;
505
   
506
  for(m=0; (m<DIM); m++)
507
    hbox[m]=box[m][m]*0.5;
508
  bTriclinic = TRICLINIC(box);
509
  
510
  ng=0;
511
  ai=*AtomI;
512
  
513
  g0=g->start;
514
  /* Loop over all the bonds */
515
  for(j=0; (j<g->nedge[ai]); j++) {
516
    aj=g->edge[ai][j]-g0;
517
    /* If there is a white one, make it grey and set pbc */
518
    if (g->bScrewPBC)
519
      mk_1shift_screw(box,hbox,x[g0+ai],x[g0+aj],g->ishift[ai],is_aj);
520
    else if (bTriclinic)
521
      mk_1shift_tric(npbcdim,box,hbox,x[g0+ai],x[g0+aj],g->ishift[ai],is_aj);
522
    else
523
      mk_1shift(npbcdim,hbox,x[g0+ai],x[g0+aj],g->ishift[ai],is_aj);
524
    
525
    if (egc[aj] == egcolWhite) {
526
      if (aj < *AtomI)
527
        *AtomI = aj;
528
      egc[aj] = egcolGrey;
529
      
530
      copy_ivec(is_aj,g->ishift[aj]);
531

    
532
      ng++;
533
    }
534
    else if ((is_aj[XX] != g->ishift[aj][XX]) ||
535
             (is_aj[YY] != g->ishift[aj][YY]) ||
536
             (is_aj[ZZ] != g->ishift[aj][ZZ])) {
537
      if (gmx_debug_at) {
538
        set_pbc(&pbc,-1,box);
539
        pbc_dx(&pbc,x[g0+ai],x[g0+aj],dx);
540
        fprintf(debug,"mk_grey: shifts for atom %d due to atom %d\n"
541
                "are (%d,%d,%d), should be (%d,%d,%d)\n"
542
                "dx = (%g,%g,%g)\n",
543
                aj+g0+1,ai+g0+1,is_aj[XX],is_aj[YY],is_aj[ZZ],
544
                g->ishift[aj][XX],g->ishift[aj][YY],g->ishift[aj][ZZ],
545
                dx[XX],dx[YY],dx[ZZ]);
546
      }
547
      (*nerror)++;
548
    }
549
  }
550
  return ng;
551
}
552

    
553
static int first_colour(int fC,egCol Col,t_graph *g,egCol egc[])
554
/* Return the first node with colour Col starting at fC.
555
 * return -1 if none found.
556
 */
557
{
558
  int i;
559
  
560
  for(i=fC; (i<g->nnodes); i++)
561
    if ((g->nedge[i] > 0) && (egc[i]==Col))
562
      return i;
563
  
564
  return -1;
565
}
566

    
567
void mk_mshift(FILE *log,t_graph *g,int ePBC,matrix box,rvec x[])
568
{
569
  static int nerror_tot = 0;
570
  int    npbcdim;
571
  int    ng,nnodes,i;
572
  int    nW,nG,nB;                /* Number of Grey, Black, White        */
573
  int    fW,fG;                        /* First of each category        */
574
  int    nerror=0;
575

    
576
  g->bScrewPBC = (ePBC == epbcSCREW);
577

    
578
  if (ePBC == epbcXY)
579
    npbcdim = 2;
580
  else
581
    npbcdim = 3;
582

    
583
  GCHECK(g);
584
  /* This puts everything in the central box, that is does not move it 
585
   * at all. If we return without doing this for a system without bonds
586
   * (i.e. only settles) all water molecules are moved to the opposite octant
587
   */
588
  for(i=0; (i<g->nnodes); i++) {
589
      g->ishift[i][XX]=g->ishift[i][YY]=g->ishift[i][ZZ]=0;
590
  }
591
    
592
  if (!g->nbound)
593
    return;
594

    
595
  nnodes=g->nnodes;
596
  if (nnodes > g->negc) {
597
    g->negc = nnodes;
598
    srenew(g->egc,g->negc);
599
  }
600
  memset(g->egc,0,(size_t)(nnodes*sizeof(g->egc[0])));
601

    
602
  nW=g->nbound;
603
  nG=0;
604
  nB=0;
605

    
606
  fW=0;
607

    
608
  /* We even have a loop invariant:
609
   * nW+nG+nB == g->nbound
610
   */
611
#ifdef DEBUG2
612
  fprintf(log,"Starting W loop\n");
613
#endif
614
  while (nW > 0) {
615
    /* Find the first white, this will allways be a larger
616
     * number than before, because no nodes are made white
617
     * in the loop
618
     */
619
    if ((fW=first_colour(fW,egcolWhite,g,g->egc)) == -1) 
620
      gmx_fatal(FARGS,"No WHITE nodes found while nW=%d\n",nW);
621
    
622
    /* Make the first white node grey */
623
    g->egc[fW]=egcolGrey;
624
    nG++;
625
    nW--;
626

    
627
    /* Initial value for the first grey */
628
    fG=fW;
629
#ifdef DEBUG2
630
    fprintf(log,"Starting G loop (nW=%d, nG=%d, nB=%d, total %d)\n",
631
            nW,nG,nB,nW+nG+nB);
632
#endif
633
    while (nG > 0) {
634
      if ((fG=first_colour(fG,egcolGrey,g,g->egc)) == -1)
635
        gmx_fatal(FARGS,"No GREY nodes found while nG=%d\n",nG);
636
      
637
      /* Make the first grey node black */
638
      g->egc[fG]=egcolBlack;
639
      nB++;
640
      nG--;
641

    
642
      /* Make all the neighbours of this black node grey
643
       * and set their periodicity 
644
       */
645
      ng=mk_grey(log,nnodes,g->egc,g,&fG,npbcdim,box,x,&nerror);
646
      /* ng is the number of white nodes made grey */
647
      nG+=ng;
648
      nW-=ng;
649
    }
650
  }
651
  if (nerror > 0) {
652
    nerror_tot++;
653
    if (nerror_tot <= 100) {
654
      fprintf(stderr,"There were %d inconsistent shifts. Check your topology\n",
655
              nerror);
656
      if (log) {
657
        fprintf(log,"There were %d inconsistent shifts. Check your topology\n",
658
                nerror);
659
      }
660
    }
661
    if (nerror_tot == 100) {
662
      fprintf(stderr,"Will stop reporting inconsistent shifts\n");
663
      if (log) {
664
        fprintf(log,"Will stop reporting inconsistent shifts\n");
665
      }
666
    }
667
  }
668
}
669

    
670
/************************************************************
671
 *
672
 *      A C T U A L   S H I F T   C O D E
673
 *
674
 ************************************************************/
675

    
676
void shift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
677
{
678
  ivec *is;
679
  int      g0,gn;
680
  int      i,j,tx,ty,tz;
681

    
682
  GCHECK(g);
683
  g0=g->start;
684
  gn=g->nnodes;
685
  is=g->ishift;
686
  
687
  if (g->bScrewPBC) {
688
    for(i=0,j=g0; (i<gn); i++,j++) { 
689
      tx=is[i][XX];
690
      ty=is[i][YY];
691
      tz=is[i][ZZ];
692
      
693
      if ((tx > 0 && tx % 2 == 1) ||
694
          (tx < 0 && -tx %2 == 1)) {
695
        x_s[j][XX] = x[j][XX] + tx*box[XX][XX];
696
        x_s[j][YY] = box[YY][YY] + box[ZZ][YY] - x[j][YY];
697
        x_s[j][ZZ] = box[ZZ][ZZ]               - x[j][ZZ];
698
      } else {
699
        x_s[j][XX] = x[j][XX];
700
      }
701
      x_s[j][YY] = x[j][YY] + ty*box[YY][YY] + tz*box[ZZ][YY];
702
      x_s[j][ZZ] = x[j][ZZ] + tz*box[ZZ][ZZ];
703
    }
704
  } else if (TRICLINIC(box)) {
705
     for(i=0,j=g0; (i<gn); i++,j++) { 
706
         tx=is[i][XX];
707
         ty=is[i][YY];
708
         tz=is[i][ZZ];
709
         
710
         x_s[j][XX]=x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
711
         x_s[j][YY]=x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
712
         x_s[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
713
     }
714
  } else {
715
     for(i=0,j=g0; (i<gn); i++,j++) { 
716
         tx=is[i][XX];
717
         ty=is[i][YY];
718
         tz=is[i][ZZ];
719
         
720
         x_s[j][XX]=x[j][XX]+tx*box[XX][XX];
721
         x_s[j][YY]=x[j][YY]+ty*box[YY][YY];
722
         x_s[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
723
     }
724
  }       
725
     
726
}
727

    
728
void shift_self(t_graph *g,matrix box,rvec x[])
729
{
730
  ivec *is;
731
  int      g0,gn;
732
  int      i,j,tx,ty,tz;
733

    
734
  if (g->bScrewPBC)
735
    gmx_incons("screw pbc not implemented for shift_self");
736

    
737
  g0=g->start;
738
  gn=g->nnodes;
739
  is=g->ishift;
740

    
741
#ifdef DEBUG
742
  fprintf(stderr,"Shifting atoms %d to %d\n",g0,g0+gn);
743
#endif
744
  if(TRICLINIC(box)) {
745
      for(i=0,j=g0; (i<gn); i++,j++) { 
746
          tx=is[i][XX];
747
          ty=is[i][YY];
748
          tz=is[i][ZZ];
749
          
750
          x[j][XX]=x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
751
          x[j][YY]=x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
752
          x[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
753
      }
754
  } else {
755
      for(i=0,j=g0; (i<gn); i++,j++) { 
756
          tx=is[i][XX];
757
          ty=is[i][YY];
758
          tz=is[i][ZZ];
759
          
760
          x[j][XX]=x[j][XX]+tx*box[XX][XX];
761
          x[j][YY]=x[j][YY]+ty*box[YY][YY];
762
          x[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
763
      }
764
  }       
765
  
766
}
767

    
768
void unshift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
769
{
770
  ivec *is;
771
  int      g0,gn;
772
  int      i,j,tx,ty,tz;
773

    
774
  if (g->bScrewPBC)
775
    gmx_incons("screw pbc not implemented for unshift_x");
776

    
777
  g0=g->start;
778
  gn=g->nnodes;
779
  is=g->ishift;
780
  if(TRICLINIC(box)) {
781
      for(i=0,j=g0; (i<gn); i++,j++) {
782
          tx=is[i][XX];
783
          ty=is[i][YY];
784
          tz=is[i][ZZ];
785
          
786
          x[j][XX]=x_s[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
787
          x[j][YY]=x_s[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
788
          x[j][ZZ]=x_s[j][ZZ]-tz*box[ZZ][ZZ];
789
      }
790
  } else {
791
      for(i=0,j=g0; (i<gn); i++,j++) {
792
          tx=is[i][XX];
793
          ty=is[i][YY];
794
          tz=is[i][ZZ];
795
          
796
          x[j][XX]=x_s[j][XX]-tx*box[XX][XX];
797
          x[j][YY]=x_s[j][YY]-ty*box[YY][YY];
798
          x[j][ZZ]=x_s[j][ZZ]-tz*box[ZZ][ZZ];
799
      }
800
  }
801
}
802

    
803
void unshift_self(t_graph *g,matrix box,rvec x[])
804
{
805
  ivec *is;
806
  int g0,gn;
807
  int i,j,tx,ty,tz;
808

    
809
  if (g->bScrewPBC)
810
    gmx_incons("screw pbc not implemented for unshift_self");
811

    
812
  g0=g->start;
813
  gn=g->nnodes;
814
  is=g->ishift;
815
  if(TRICLINIC(box)) {
816
      for(i=0,j=g0; (i<gn); i++,j++) {
817
          tx=is[i][XX];
818
          ty=is[i][YY];
819
          tz=is[i][ZZ];
820
          
821
          x[j][XX]=x[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
822
          x[j][YY]=x[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
823
          x[j][ZZ]=x[j][ZZ]-tz*box[ZZ][ZZ];
824
      }
825
  } else {
826
      for(i=0,j=g0; (i<gn); i++,j++) {
827
          tx=is[i][XX];
828
          ty=is[i][YY];
829
          tz=is[i][ZZ];
830
          
831
          x[j][XX]=x[j][XX]-tx*box[XX][XX];
832
          x[j][YY]=x[j][YY]-ty*box[YY][YY];
833
          x[j][ZZ]=x[j][ZZ]-tz*box[ZZ][ZZ];
834
      }
835
  }
836
}
837
#undef GCHECK
838

    
839
#ifdef DEBUGMSHIFT
840
void main(int argc,char *argv[])
841
{
842
  FILE         *out;
843
  t_args       targ;
844
  t_topology   top;
845
  t_statheader sh;
846
  rvec         *x;
847
  ivec         *mshift;
848
  matrix       box;
849

    
850
  t_graph      *g;
851
  int          i,idum,pid;
852
  real         rdum;
853

    
854
  CopyRight(stderr,argv[0]);
855
  parse_common_args(&argc,argv,&targ,PCA_NEED_INOUT,NULL);
856
  if (argc > 1)
857
    pid=atoi(argv[1]);
858
  else
859
    pid=0;
860
  
861
  read_status_header(targ.infile,&sh);
862
  snew(x,sh.natoms);
863
  snew(mshift,sh.natoms);
864

    
865
  fprintf(stderr,"Reading Status %s\n",targ.infile);
866
  read_status(targ.infile,&idum,&rdum,&rdum,NULL,
867
              box,NULL,NULL,&idum,x,NULL,NULL,&idum,NULL,&top);
868

    
869
  fprintf(stderr,"Making Graph Structure...\n");
870
  g=mk_graph(&(top.idef),top.atoms.nr,FALSE,FALSE);
871

    
872
  out=gmx_fio_fopen(targ.outfile,"w");
873

    
874
  fprintf(stderr,"Making Shift...\n");
875
  mk_mshift(out,g,box,x,mshift);
876

    
877
  p_graph(out,"In Den Haag daar woont een graaf...",g);
878
  gmx_fio_fclose(out);
879
}
880
#endif
881