Merge "reduced PLINCS communication with H-bond constraints" into release-4-6

diff --git a/src/mdlib/clincs.c b/src/mdlib/clincs.c
index adb1a7c99306183526e6561f56bca6036122f569..d88241d6734bd24bb915fb13f209c568d432cff5 100644 (file)
--- a/src/mdlib/clincs.c
+++ b/src/mdlib/clincs.c
@@ -90,6 +90,7 @@ typedef struct gmx_lincsdata {
     int  *triangle;   /* the list of triangle constraints */
     int  *tri_bits;   /* the bits tell if the matrix element should be used */
     int  ncc_triangle;/* the number of constraint connections in triangles */
+    gmx_bool bCommIter; /* communicate before each LINCS interation */
     real *blmf;       /* matrix of mass factors for constraint connections */
     real *blmf1;      /* as blmf, but with all masses 1 */
     real *bllen;      /* the reference bond length */
@@ -589,7 +590,7 @@ static void do_lincs(rvec *x,rvec *xp,matrix box,t_pbc *pbc,
        
     for(iter=0; iter<lincsd->nIter; iter++)
     {
-        if ((DOMAINDECOMP(cr) && cr->dd->constraints) ||
+        if ((lincsd->bCommIter && DOMAINDECOMP(cr) && cr->dd->constraints) ||
             PARTDECOMP(cr))
         {
 #pragma omp barrier
@@ -850,6 +851,37 @@ static int count_triangle_constraints(t_ilist *ilist,t_blocka *at2con)
     return ncon_triangle;
 }
 
+static gmx_bool more_than_two_sequential_constraints(const t_ilist *ilist,
+                                                     const t_blocka *at2con)
+{
+    t_iatom  *ia1,*ia2,*iap;
+    int      ncon1,ncon_tot,c;
+    int      a1,a2;
+    gmx_bool bMoreThanTwoSequentialConstraints;
+
+    ncon1    = ilist[F_CONSTR].nr/3;
+    ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3;
+
+    ia1 = ilist[F_CONSTR].iatoms;
+    ia2 = ilist[F_CONSTRNC].iatoms;
+
+    bMoreThanTwoSequentialConstraints = FALSE;
+    for(c=0; c<ncon_tot && !bMoreThanTwoSequentialConstraints; c++)
+    {
+        iap = constr_iatomptr(ncon1,ia1,ia2,c);
+        a1 = iap[1];
+        a2 = iap[2];
+        /* Check if this constraint has constraints connected at both atoms */
+        if (at2con->index[a1+1] - at2con->index[a1] > 1 &&
+            at2con->index[a2+1] - at2con->index[a2] > 1)
+        {
+            bMoreThanTwoSequentialConstraints = TRUE;
+        }
+    }
+
+    return bMoreThanTwoSequentialConstraints;
+}
+
 static int int_comp(const void *a,const void *b)
 {
     return (*(int *)a) - (*(int *)b);
@@ -876,7 +908,11 @@ gmx_lincsdata_t init_lincs(FILE *fplog,gmx_mtop_t *mtop,
         gmx_mtop_ftype_count(mtop,F_CONSTRNC);
     li->ncg_flex = nflexcon_global;
     
+    li->nIter  = nIter;
+    li->nOrder = nProjOrder;
+
     li->ncg_triangle = 0;
+    li->bCommIter = FALSE;
     for(mb=0; mb<mtop->nmolblock; mb++)
     {
         molt = &mtop->moltype[mtop->molblock[mb].type];
@@ -884,10 +920,23 @@ gmx_lincsdata_t init_lincs(FILE *fplog,gmx_mtop_t *mtop,
             mtop->molblock[mb].nmol*
             count_triangle_constraints(molt->ilist,
                                        &at2con[mtop->molblock[mb].type]);
+        if (bPLINCS && li->bCommIter == FALSE)
+        {
+            /* Check if we need to communicate not only before LINCS,
+             * but also before each iteration.
+             * The check for only two sequential constraints is only
+             * useful for the common case of H-bond only constraints.
+             * With more effort we could also make it useful for small
+             * molecules with nr. sequential constraints <= nOrder-1.
+             */
+            li->bCommIter = (li->nOrder < 1 || more_than_two_sequential_constraints(molt->ilist,&at2con[mtop->molblock[mb].type]));
+        }
+    }
+    if (debug && bPLINCS)
+    {
+        fprintf(debug,"PLINCS communication before each iteration: %d\n",
+                li->bCommIter);
     }
-    
-    li->nIter  = nIter;
-    li->nOrder = nProjOrder;
 
     /* LINCS can run on any number of threads.
      * Currently the number is fixed for the whole simulation,