added support for graphs not starting at 0

It has always been implicitly assumed that graph->start==0.
This assumption would only fail for a topology starting with
single atom molecules followed by multi-atom molecules.
Furthermore, for mk_graph at_start could only be 0, a check
has been added for this, and at_end was not stored.

In t_graph start is now replaced by at_start and end by at_end,
where at_end is end+1, now the first atom not part of the graph.

Change-Id: I0d4c09c90cbba1a7f4a1d84cf1f90a7643c17cd8

diff --git a/include/types/graph.h b/include/types/graph.h
index b9b19355227b8722b7887319788ad53c3f8cf70d..5d80ee4a60420b13fd7403bbd25ccf6588b54dd0 100644 (file)
--- a/include/types/graph.h
+++ b/include/types/graph.h
@@ -43,10 +43,11 @@ extern "C" {
 typedef enum { egcolWhite, egcolGrey, egcolBlack, egcolNR } egCol;
 
 typedef struct {
-  int      nnodes;     /* The number of nodes                          */
+  int      nnodes;     /* The number of nodes, nnodes=at_end-at_start  */
   int      nbound;     /* The number of nodes with edges               */
-  int      start;      /* The first atom in this graph                 */
-  int      end;                /* The last atom in this graph                  */
+  int      natoms;      /* Total range for this graph: 0 to natoms      */
+  int      at_start;   /* The first connected atom in this graph       */
+  int      at_end;     /* The last+1 connected atom in this graph      */
   int      *nedge;     /* For each node the number of edges            */
   atom_id  **edge;     /* For each node, the actual edges (bidirect.)  */
   gmx_bool     bScrewPBC;   /* Screw boundary conditions                    */
@@ -56,7 +57,7 @@ typedef struct {
 } t_graph;
 
 
-#define SHIFT_IVEC(g,i) ((g)->ishift[(i)-(g)->start])
+#define SHIFT_IVEC(g,i) ((g)->ishift[i])
 
 #ifdef __cplusplus
 }

diff --git a/src/gmxlib/mshift.c b/src/gmxlib/mshift.c
index 0aae993a3275e6e38e26f3434223c93e6f1ff150..0e9e102539a6f3da260b9fb48b37ce2d7e4052e9 100644 (file)
--- a/src/gmxlib/mshift.c
+++ b/src/gmxlib/mshift.c
@@ -66,8 +66,8 @@ static void add_gbond(t_graph *g,atom_id a0,atom_id a1)
   atom_id inda0,inda1;
   gmx_bool    bFound;
 
-  inda0 = a0 - g->start;
-  inda1 = a1 - g->start;
+  inda0 = a0 - g->at_start;
+  inda1 = a1 - g->at_start;
   bFound = FALSE;
   /* Search for a direct edge between a0 and a1.
    * All egdes are bidirectional, so we only need to search one way.
@@ -145,14 +145,15 @@ void p_graph(FILE *log,const char *title,t_graph *g)
   fprintf(log,"graph:  %s\n",title);
   fprintf(log,"nnodes: %d\n",g->nnodes);
   fprintf(log,"nbound: %d\n",g->nbound);
-  fprintf(log,"start:  %d\n",g->start);
-  fprintf(log,"end:    %d\n",g->end);
+  fprintf(log,"start:  %d\n",g->at_start);
+  fprintf(log,"end:    %d\n",g->at_end);
   fprintf(log," atom shiftx shifty shiftz C nedg    e1    e2 etc.\n");
   for(i=0; (i<g->nnodes); i++)
     if (g->nedge[i] > 0) {
-      fprintf(log,"%5d%7d%7d%7d %1s%5d",g->start+i+1,
-             g->ishift[i][XX],g->ishift[i][YY],
-             g->ishift[i][ZZ],
+      fprintf(log,"%5d%7d%7d%7d %1s%5d",g->at_start+i+1,
+             g->ishift[g->at_start+i][XX],
+             g->ishift[g->at_start+i][YY],
+             g->ishift[g->at_start+i][ZZ],
              (g->negc > 0) ? cc[g->egc[i]] : " ",
              g->nedge[i]);
       for(j=0; (j<g->nedge[i]); j++)
@@ -180,15 +181,15 @@ static void calc_1se(t_graph *g,int ftype,t_ilist *il,
        nbond[iaa]   += 2;
        nbond[iaa+1] += 1;
        nbond[iaa+2] += 1;
-       g->start      = min(g->start,iaa);
-       g->end        = max(g->end,iaa+2);
+       g->at_start   = min(g->at_start,iaa);
+       g->at_end     = max(g->at_end,iaa+2+1);
       }
     } else {
       for(k=1; (k<=nratoms); k++) {
        iaa=ia[k];
        if (iaa >= at_start && iaa < at_end) {
-         g->start=min(g->start,iaa);
-         g->end  =max(g->end,  iaa);
+         g->at_start = min(g->at_start,iaa);
+         g->at_end   = max(g->at_end,  iaa+1);
          /* When making the graph we (might) link all atoms in an interaction
           * sequentially. Therefore the end atoms add 1 to the count,
           * the middle atoms 2.
@@ -210,8 +211,8 @@ static int calc_start_end(FILE *fplog,t_graph *g,t_ilist il[],
 {
   int   i,nnb,nbtot;
   
-  g->start=at_end;
-  g->end=0;
+  g->at_start = at_end;
+  g->at_end   = 0;
 
   /* First add all the real bonds: they should determine the molecular 
    * graph.
@@ -230,7 +231,7 @@ static int calc_start_end(FILE *fplog,t_graph *g,t_ilist il[],
   
   nnb   = 0;
   nbtot = 0;
-  for(i=g->start; (i<=g->end); i++) {
+  for(i=g->at_start; (i<g->at_end); i++) {
     nbtot += nbond[i];
     nnb    = max(nnb,nbond[i]);
   }
@@ -277,7 +278,7 @@ static gmx_bool determine_graph_parts(t_graph *g,int *part)
   gmx_bool bMultiPart;
 
   /* Initialize the part array with all entries different */
-  for(at_i=g->start; at_i<g->end; at_i++) {
+  for(at_i=g->at_start; at_i<g->at_end; at_i++) {
     part[at_i] = at_i;
   }
 
@@ -286,7 +287,7 @@ static gmx_bool determine_graph_parts(t_graph *g,int *part)
     bMultiPart = FALSE;
     nchanged = 0;
     for(i=0; (i<g->nnodes); i++) {
-      at_i  = g->start + i;
+      at_i  = g->at_start + i;
       at_i2 = g->edge[i];
       for(e=0; e<g->nedge[i]; e++) {
        /* Set part for both nodes to the minimum */
@@ -298,7 +299,7 @@ static gmx_bool determine_graph_parts(t_graph *g,int *part)
          nchanged++;
        }
       }
-      if (part[at_i] != part[g->start]) {
+      if (part[at_i] != part[g->at_start]) {
        bMultiPart = TRUE;
       }
     }
@@ -320,22 +321,28 @@ void mk_graph_ilist(FILE *fplog,
   int     i,nbtot;
   gmx_bool    bMultiPart;
 
+  if (at_start != 0) {
+    gmx_incons("In mk_graph_ilist at_start can not be != 0");
+  }
+  g->natoms = at_end;
+
   snew(nbond,at_end);
   nbtot = calc_start_end(fplog,g,ilist,at_start,at_end,nbond);
   
-  if (g->start >= g->end) {
-    g->nnodes = 0;
-    g->nbound = 0;
+  if (g->at_start >= g->at_end) {
+    g->at_start = at_start;
+    g->at_end   = at_end;
+    g->nnodes   = 0;
+    g->nbound   = 0;
   }
   else {
-    g->nnodes = g->end - g->start + 1;
-    snew(g->ishift,g->nnodes);
+    g->nnodes = g->at_end - g->at_start;
     snew(g->nedge,g->nnodes);
     snew(g->edge,g->nnodes);
     /* Allocate a single array and set pointers into it */
     snew(g->edge[0],nbtot);
     for(i=1; (i<g->nnodes); i++) {
-      g->edge[i] = g->edge[i-1] + nbond[g->start+i-1];
+      g->edge[i] = g->edge[i-1] + nbond[g->at_start+i-1];
     }
 
     if (!bShakeOnly) {
@@ -383,6 +390,8 @@ void mk_graph_ilist(FILE *fplog,
 
   sfree(nbond);
 
+  snew(g->ishift,g->natoms);
+
   if (gmx_debug_at)
     p_graph(debug,"graph",g);
 }
@@ -406,12 +415,12 @@ void done_graph(t_graph *g)
   
   GCHECK(g);
   if (g->nnodes > 0) {
-    sfree(g->ishift);
     sfree(g->nedge);
     sfree(g->edge[0]);
     sfree(g->edge);
     sfree(g->egc);
   }
+  sfree(g->ishift);
 }
 
 /************************************************************
@@ -522,25 +531,25 @@ static int mk_grey(FILE *log,int nnodes,egCol egc[],t_graph *g,int *AtomI,
     hbox[m]=box[m][m]*0.5;
   bTriclinic = TRICLINIC(box);
   
-  ng=0;
-  ai=*AtomI;
+  g0 = g->at_start;
+  ng = 0;
+  ai = g0 + *AtomI;
   
-  g0=g->start;
   /* Loop over all the bonds */
-  for(j=0; (j<g->nedge[ai]); j++) {
-    aj=g->edge[ai][j]-g0;
+  for(j=0; (j<g->nedge[ai-g0]); j++) {
+    aj = g->edge[ai-g0][j];
     /* If there is a white one, make it grey and set pbc */
     if (g->bScrewPBC)
-      mk_1shift_screw(box,hbox,x[g0+ai],x[g0+aj],g->ishift[ai],is_aj);
+      mk_1shift_screw(box,hbox,x[ai],x[aj],g->ishift[ai],is_aj);
     else if (bTriclinic)
-      mk_1shift_tric(npbcdim,box,hbox,x[g0+ai],x[g0+aj],g->ishift[ai],is_aj);
+      mk_1shift_tric(npbcdim,box,hbox,x[ai],x[aj],g->ishift[ai],is_aj);
     else
-      mk_1shift(npbcdim,hbox,x[g0+ai],x[g0+aj],g->ishift[ai],is_aj);
+      mk_1shift(npbcdim,hbox,x[ai],x[aj],g->ishift[ai],is_aj);
     
-    if (egc[aj] == egcolWhite) {
-      if (aj < *AtomI)
-       *AtomI = aj;
-      egc[aj] = egcolGrey;
+    if (egc[aj-g0] == egcolWhite) {
+      if (aj - g0 < *AtomI)
+       *AtomI = aj - g0;
+      egc[aj-g0] = egcolGrey;
       
       copy_ivec(is_aj,g->ishift[aj]);
 
@@ -551,11 +560,11 @@ static int mk_grey(FILE *log,int nnodes,egCol egc[],t_graph *g,int *AtomI,
             (is_aj[ZZ] != g->ishift[aj][ZZ])) {
       if (gmx_debug_at) {
        set_pbc(&pbc,-1,box);
-       pbc_dx(&pbc,x[g0+ai],x[g0+aj],dx);
+       pbc_dx(&pbc,x[ai],x[aj],dx);
        fprintf(debug,"mk_grey: shifts for atom %d due to atom %d\n"
                "are (%d,%d,%d), should be (%d,%d,%d)\n"
                "dx = (%g,%g,%g)\n",
-               aj+g0+1,ai+g0+1,is_aj[XX],is_aj[YY],is_aj[ZZ],
+               aj+1,ai+1,is_aj[XX],is_aj[YY],is_aj[ZZ],
                g->ishift[aj][XX],g->ishift[aj][YY],g->ishift[aj][ZZ],
                dx[XX],dx[YY],dx[ZZ]);
       }
@@ -600,7 +609,7 @@ void mk_mshift(FILE *log,t_graph *g,int ePBC,matrix box,rvec x[])
    * at all. If we return without doing this for a system without bonds
    * (i.e. only settles) all water molecules are moved to the opposite octant
    */
-  for(i=0; (i<g->nnodes); i++) {
+  for(i=0; (i<g->natoms); i++) {
       g->ishift[i][XX]=g->ishift[i][YY]=g->ishift[i][ZZ]=0;
   }
     
@@ -691,19 +700,23 @@ void mk_mshift(FILE *log,t_graph *g,int ePBC,matrix box,rvec x[])
 void shift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
 {
   ivec *is;
-  int      g0,gn;
-  int      i,j,tx,ty,tz;
+  int      g0,g1;
+  int      j,tx,ty,tz;
 
   GCHECK(g);
-  g0=g->start;
-  gn=g->nnodes;
-  is=g->ishift;
+  g0 = g->at_start;
+  g1 = g->at_end;
+  is = g->ishift;
   
+  for(j=0; j<g0; j++) {
+    copy_rvec(x[j],x_s[j]);
+  }
+
   if (g->bScrewPBC) {
-    for(i=0,j=g0; (i<gn); i++,j++) { 
-      tx=is[i][XX];
-      ty=is[i][YY];
-      tz=is[i][ZZ];
+    for(j=g0; (j<g1); j++) { 
+      tx=is[j][XX];
+      ty=is[j][YY];
+      tz=is[j][ZZ];
       
       if ((tx > 0 && tx % 2 == 1) ||
          (tx < 0 && -tx %2 == 1)) {
@@ -717,131 +730,143 @@ void shift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
       x_s[j][ZZ] = x[j][ZZ] + tz*box[ZZ][ZZ];
     }
   } else if (TRICLINIC(box)) {
-     for(i=0,j=g0; (i<gn); i++,j++) { 
-        tx=is[i][XX];
-        ty=is[i][YY];
-        tz=is[i][ZZ];
+     for(j=g0; (j<g1); j++) { 
+        tx=is[j][XX];
+        ty=is[j][YY];
+        tz=is[j][ZZ];
         
         x_s[j][XX]=x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
         x_s[j][YY]=x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
         x_s[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
      }
   } else {
-     for(i=0,j=g0; (i<gn); i++,j++) { 
-        tx=is[i][XX];
-        ty=is[i][YY];
-        tz=is[i][ZZ];
+     for(j=g0; (j<g1); j++) { 
+        tx=is[j][XX];
+        ty=is[j][YY];
+        tz=is[j][ZZ];
         
         x_s[j][XX]=x[j][XX]+tx*box[XX][XX];
         x_s[j][YY]=x[j][YY]+ty*box[YY][YY];
         x_s[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
      }
   }       
-     
+
+  for(j=g1; j<g->natoms; j++) {
+    copy_rvec(x[j],x_s[j]);
+  }
 }
 
 void shift_self(t_graph *g,matrix box,rvec x[])
 {
   ivec *is;
-  int      g0,gn;
-  int      i,j,tx,ty,tz;
+  int      g0,g1;
+  int      j,tx,ty,tz;
 
   if (g->bScrewPBC)
     gmx_incons("screw pbc not implemented for shift_self");
 
-  g0=g->start;
-  gn=g->nnodes;
-  is=g->ishift;
+  g0 = g->at_start;
+  g1 = g->at_end;
+  is = g->ishift;
 
 #ifdef DEBUG
   fprintf(stderr,"Shifting atoms %d to %d\n",g0,g0+gn);
 #endif
   if(TRICLINIC(box)) {
-      for(i=0,j=g0; (i<gn); i++,j++) { 
-         tx=is[i][XX];
-         ty=is[i][YY];
-         tz=is[i][ZZ];
+      for(j=g0; (j<g1); j++) { 
+         tx=is[j][XX];
+         ty=is[j][YY];
+         tz=is[j][ZZ];
          
          x[j][XX]=x[j][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
          x[j][YY]=x[j][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
          x[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
       }
   } else {
-      for(i=0,j=g0; (i<gn); i++,j++) { 
-         tx=is[i][XX];
-         ty=is[i][YY];
-         tz=is[i][ZZ];
+      for(j=g0; (j<g1); j++) { 
+         tx=is[j][XX];
+         ty=is[j][YY];
+         tz=is[j][ZZ];
          
          x[j][XX]=x[j][XX]+tx*box[XX][XX];
          x[j][YY]=x[j][YY]+ty*box[YY][YY];
          x[j][ZZ]=x[j][ZZ]+tz*box[ZZ][ZZ];
       }
   }       
-  
 }
 
 void unshift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
 {
   ivec *is;
-  int      g0,gn;
-  int      i,j,tx,ty,tz;
+  int      g0,g1;
+  int      j,tx,ty,tz;
 
   if (g->bScrewPBC)
-    gmx_incons("screw pbc not implemented for unshift_x");
+    gmx_incons("screw pbc not implemented (yet) for unshift_x");
+
+  g0 = g->at_start;
+  g1 = g->at_end;
+  is = g->ishift;
+
+  for(j=0; j<g0; j++) {
+    copy_rvec(x_s[j],x[j]);
+  }
 
-  g0=g->start;
-  gn=g->nnodes;
-  is=g->ishift;
   if(TRICLINIC(box)) {
-      for(i=0,j=g0; (i<gn); i++,j++) {
-         tx=is[i][XX];
-         ty=is[i][YY];
-         tz=is[i][ZZ];
+      for(j=g0; (j<g1); j++) {
+         tx=is[j][XX];
+         ty=is[j][YY];
+         tz=is[j][ZZ];
          
          x[j][XX]=x_s[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
          x[j][YY]=x_s[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
          x[j][ZZ]=x_s[j][ZZ]-tz*box[ZZ][ZZ];
       }
   } else {
-      for(i=0,j=g0; (i<gn); i++,j++) {
-         tx=is[i][XX];
-         ty=is[i][YY];
-         tz=is[i][ZZ];
+      for(j=g0; (j<g1); j++) {
+         tx=is[j][XX];
+         ty=is[j][YY];
+         tz=is[j][ZZ];
          
          x[j][XX]=x_s[j][XX]-tx*box[XX][XX];
          x[j][YY]=x_s[j][YY]-ty*box[YY][YY];
          x[j][ZZ]=x_s[j][ZZ]-tz*box[ZZ][ZZ];
       }
   }
+
+  for(j=g1; j<g->natoms; j++) {
+    copy_rvec(x_s[j],x[j]);
+  }
 }
 
 void unshift_self(t_graph *g,matrix box,rvec x[])
 {
   ivec *is;
-  int g0,gn;
-  int i,j,tx,ty,tz;
+  int g0,g1;
+  int j,tx,ty,tz;
 
   if (g->bScrewPBC)
     gmx_incons("screw pbc not implemented for unshift_self");
 
-  g0=g->start;
-  gn=g->nnodes;
-  is=g->ishift;
+  g0 = g->at_start;
+  g1 = g->at_end;
+  is = g->ishift;
+
   if(TRICLINIC(box)) {
-      for(i=0,j=g0; (i<gn); i++,j++) {
-         tx=is[i][XX];
-         ty=is[i][YY];
-         tz=is[i][ZZ];
+      for(j=g0; (j<g1); j++) {
+         tx=is[j][XX];
+         ty=is[j][YY];
+         tz=is[j][ZZ];
          
          x[j][XX]=x[j][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
          x[j][YY]=x[j][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
          x[j][ZZ]=x[j][ZZ]-tz*box[ZZ][ZZ];
       }
   } else {
-      for(i=0,j=g0; (i<gn); i++,j++) {
-         tx=is[i][XX];
-         ty=is[i][YY];
-         tz=is[i][ZZ];
+      for(j=g0; (j<g1); j++) {
+         tx=is[j][XX];
+         ty=is[j][YY];
+         tz=is[j][ZZ];
          
          x[j][XX]=x[j][XX]-tx*box[XX][XX];
          x[j][YY]=x[j][YY]-ty*box[YY][YY];

diff --git a/src/gmxlib/pbc.c b/src/gmxlib/pbc.c
index af63fc6d047219aaaae81050ad343c75d115b221..21655152c09a27e66c0d6e2c7738b262811b583f 100644 (file)
--- a/src/gmxlib/pbc.c
+++ b/src/gmxlib/pbc.c
@@ -243,7 +243,7 @@ gmx_bool correct_box(FILE *fplog,int step,tensor box,t_graph *graph)
 
     if (bCorrected && graph) {
         /* correct the graph */
-        for(i=0; i<graph->nnodes; i++) {
+        for(i=graph->at_start; i<graph->at_end; i++) {
             graph->ishift[i][YY] -= graph->ishift[i][ZZ]*zy;
             graph->ishift[i][XX] -= graph->ishift[i][ZZ]*zx;
             graph->ishift[i][XX] -= graph->ishift[i][YY]*yx;

diff --git a/src/gmxlib/rmpbc.c b/src/gmxlib/rmpbc.c
index fb2843f02125666e3dcaf5530bd617aede479a69..41aa1636d235a1b2c23d2e59f9a0a10f17d2e58e 100644 (file)
--- a/src/gmxlib/rmpbc.c
+++ b/src/gmxlib/rmpbc.c
@@ -168,7 +168,7 @@ void gmx_rmpbc(gmx_rmpbc_t gpbc,int natoms,matrix box,rvec x[])
     if (gr != NULL)
     {
         mk_mshift(stdout,gr,ePBC,box,x);
-        shift_x(gr,box,x,x);
+        shift_self(gr,box,x);
     }
 }
 
@@ -206,7 +206,7 @@ void gmx_rmpbc_trxfr(gmx_rmpbc_t gpbc,t_trxframe *fr)
         if (gr != NULL)
         {
             mk_mshift(stdout,gr,ePBC,fr->box,fr->x);
-            shift_x(gr,fr->box,fr->x,fr->x);
+            shift_self(gr,fr->box,fr->x);
         }
     }
 }

diff --git a/src/gmxlib/splitter.c b/src/gmxlib/splitter.c
index de954fbc059a7c2b0f36592af12477190c0367b0..284aa1c3b9f59ff3a0274e38215cef67b85b4a43 100644 (file)
--- a/src/gmxlib/splitter.c
+++ b/src/gmxlib/splitter.c
@@ -518,7 +518,7 @@ static int mk_grey(int nnodes,egCol egc[],t_graph *g,int *AtomI,
   ng=0;
   ai=*AtomI;
   
-  g0=g->start;
+  g0=g->at_start;
   /* Loop over all the bonds */
   for(j=0; (j<g->nedge[ai]); j++) {
     aj=g->edge[ai][j]-g0;
@@ -574,7 +574,7 @@ static int mk_sblocks(FILE *fp,t_graph *g,int maxsid,t_sid sid[])
   nnodes=g->nnodes;
   snew(egc,nnodes);
   
-  g0=g->start;
+  g0=g->at_start;
   nW=g->nbound;
   nG=0;
   nB=0;

diff --git a/src/mdlib/calcvir.c b/src/mdlib/calcvir.c
index 2cdc6dad72ffb6d0a9356f5b769aa735d7a997bf..6a15a76dcd090f61b4f34eaff1025a85bbdc8895 100644 (file)
--- a/src/mdlib/calcvir.c
+++ b/src/mdlib/calcvir.c
@@ -88,21 +88,20 @@ void calc_vir(FILE *log,int nxf,rvec x[],rvec f[],tensor vir,
 }
 
 
-static void lo_fcv(int i0,int i1,int g0,
+static void lo_fcv(int i0,int i1,
                   real x[],real f[],tensor vir,
                   int is[],real box[], gmx_bool bTriclinic)
 {
-  int      i,i3,gg,g3,tx,ty,tz;
+  int      i,i3,tx,ty,tz;
   real     xx,yy,zz;
   real     dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
 
   if(bTriclinic) {
-      for(i=i0,gg=g0; (i<i1); i++,gg++) {
+      for(i=i0; (i<i1); i++) {
          i3=DIM*i;
-         g3=DIM*gg;
-         tx=is[g3+XX];
-         ty=is[g3+YY];
-         tz=is[g3+ZZ];
+         tx=is[i3+XX];
+         ty=is[i3+YY];
+         tz=is[i3+ZZ];
          
          xx=x[i3+XX]-tx*box[XXXX]-ty*box[YYXX]-tz*box[ZZXX];
          dvxx+=xx*f[i3+XX];
@@ -120,12 +119,11 @@ static void lo_fcv(int i0,int i1,int g0,
          dvzz+=zz*f[i3+ZZ];
       }
   } else {
-      for(i=i0,gg=g0; (i<i1); i++,gg++) {
+      for(i=i0; (i<i1); i++) {
          i3=DIM*i;
-         g3=DIM*gg;
-         tx=is[g3+XX];
-         ty=is[g3+YY];
-         tz=is[g3+ZZ];
+         tx=is[i3+XX];
+         ty=is[i3+YY];
+         tz=is[i3+ZZ];
          
          xx=x[i3+XX]-tx*box[XXXX];
          dvxx+=xx*f[i3+XX];
@@ -215,12 +213,12 @@ void f_calc_vir(FILE *log,int i0,int i1,rvec x[],rvec f[],tensor vir,
     /* Calculate virial for bonded forces only when they belong to
      * this node.
      */
-    start = max(i0,g->start);
-    end   = min(i1,g->end+1);
+    start = max(i0,g->at_start);
+    end   = min(i1,g->at_end);
 #ifdef SAFE
     lo_fcv2(start,end,x,f,vir,g->ishift,box,TRICLINIC(box));
 #else
-    lo_fcv(start,end,0,x[0],f[0],vir,g->ishift[0],box[0],TRICLINIC(box));
+    lo_fcv(start,end,x[0],f[0],vir,g->ishift[0],box[0],TRICLINIC(box));
 #endif
     
     /* If not all atoms are bonded, calculate their virial contribution 

diff --git a/src/mdlib/update.c b/src/mdlib/update.c
index 96ed94565a556604823129d7a1447d42a645ceac..09c501cac86c1d6bfae6f82b521186685ee9c218 100644 (file)
--- a/src/mdlib/update.c
+++ b/src/mdlib/update.c
@@ -1427,8 +1427,6 @@ void update_constraints(FILE         *fplog,
             {
                 inc_nrnb(nrnb,eNR_SHIFTX,graph->nnodes);    
             }
-            copy_rvecn(upd->xp,state->x,start,graph->start);
-            copy_rvecn(upd->xp,state->x,graph->start+graph->nnodes,nrend);
         }
         else 
         {