made SHAKE work again with particle decompostion

Making a graph not starting at atom zero was not allowed,
as this did not properly with shifting atoms (and was never
called that way for shifting), but SHAKE with PD requires this.
Now also shifting not starting at 0 works correctly.

Change-Id: Icfbf51db36f2dfd7c76ea7e4e2a4bb4b0b895e47

diff --git a/include/mshift.h b/include/mshift.h
index 0e600165d773c0c75943b37336ef7972df95e6b2..a692e31f2530c6b3191e59c2dc61257784a3acd2 100644 (file)
--- a/include/mshift.h
+++ b/include/mshift.h
@@ -51,6 +51,8 @@ t_graph *mk_graph(FILE *fplog,
                         gmx_bool bShakeOnly,gmx_bool bSettle);
 /* Build a graph from an idef description. The graph can be used
  * to generate mol-shift indices.
+ * at_start and at_end should coincide will molecule boundaries,
+ * for the whole system this is simply 0 and natoms.
  * If bShakeOnly, only the connections in the shake list are used.
  * If bSettle && bShakeOnly the settles are used too.
  */

diff --git a/include/types/graph.h b/include/types/graph.h
index 96c94ba30cb912518a7e32ef47ffa60c58cc1459..306c93d1df6f292ea09d3ad5a2d6e713f21645f3 100644 (file)
--- a/include/types/graph.h
+++ b/include/types/graph.h
@@ -49,9 +49,10 @@ extern "C" {
 typedef enum { egcolWhite, egcolGrey, egcolBlack, egcolNR } egCol;
 
 typedef struct {
+  int      at0;     /* The first atom the graph was constructed for */
+  int      at1;     /* The last atom the graph was constructed for */
   int      nnodes;     /* The number of nodes, nnodes=at_end-at_start  */
   int      nbound;     /* The number of nodes with edges               */
-  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            */

diff --git a/src/gmxlib/mshift.c b/src/gmxlib/mshift.c
index abf9b7acd33222f3a6c53fe1d43ac81bb51389cf..d7542c1dc4951838d6ee41808c3e836686f97593 100644 (file)
--- a/src/gmxlib/mshift.c
+++ b/src/gmxlib/mshift.c
@@ -324,10 +324,13 @@ 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;
+  /* The naming is somewhat confusing, but we need g->at0 and g->at1
+   * for shifthing coordinates to a new array (not in place) when
+   * some atoms are not connected by the graph, which runs from
+   * g->at_start (>= g->at0) to g->at_end (<= g->at1).
+   */
+  g->at0 = at_start;
+  g->at1 = at_end;
 
   snew(nbond,at_end);
   nbtot = calc_start_end(fplog,g,ilist,at_start,at_end,nbond);
@@ -393,7 +396,7 @@ void mk_graph_ilist(FILE *fplog,
 
   sfree(nbond);
 
-  snew(g->ishift,g->natoms);
+  snew(g->ishift,g->at1);
 
   if (gmx_debug_at)
     p_graph(debug,"graph",g);
@@ -612,7 +615,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->natoms); i++) {
+  for(i=g->at0; (i<g->at1); i++) {
       g->ishift[i][XX]=g->ishift[i][YY]=g->ishift[i][ZZ]=0;
   }
     
@@ -711,7 +714,7 @@ void shift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
   g1 = g->at_end;
   is = g->ishift;
   
-  for(j=0; j<g0; j++) {
+  for(j=g->at0; j<g0; j++) {
     copy_rvec(x[j],x_s[j]);
   }
 
@@ -754,7 +757,7 @@ void shift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
      }
   }       
 
-  for(j=g1; j<g->natoms; j++) {
+  for(j=g1; j<g->at1; j++) {
     copy_rvec(x[j],x_s[j]);
   }
 }
@@ -811,7 +814,7 @@ void unshift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
   g1 = g->at_end;
   is = g->ishift;
 
-  for(j=0; j<g0; j++) {
+  for(j=g->at0; j<g0; j++) {
     copy_rvec(x_s[j],x[j]);
   }
 
@@ -837,7 +840,7 @@ void unshift_x(t_graph *g,matrix box,rvec x[],rvec x_s[])
       }
   }
 
-  for(j=g1; j<g->natoms; j++) {
+  for(j=g1; j<g->at1; j++) {
     copy_rvec(x_s[j],x[j]);
   }
 }