Fixed essential dynamics / flooding group PBC serial

In former versions, the PBC representation of essential dynamics /
flooding group atoms could be incorrect in serial runs if the ED group
contained more than a single molecule. In multi-molecule cases, the required
steps to choose the correct PBC image in communicate_group_positions()
therefore need to be performed also in serial runs. Since the PBC representation
can only change in neigborsearching steps, we only need to check the
shifts then. In parallel, NS is signalled by the bUpdateShifts
variable, which is set in dd_make_local_ed_indices(). The latter
function is however not called in serial runs; but still we can pass
the bNS boolean to do_flood() to signal the NS status. For essential
dynamics, unfortunately, since do_edsam() is called from constrain(), there
is no information about the NS status at that point. Until someone
comes up with a better idea, we therefore do the PBC check in every step
in serial essential dynamics - the performance impact will be negligible
anyway.

Change-Id: I86336a5e34131bdeac7e28f35b1ccb633450e54e

diff --git a/include/edsam.h b/include/edsam.h
index d33c52ad2446d2fb07e222055f75f74d9cb059c1..564d61128b09d99daa6df36b9b0b0eb33bb2618b 100644 (file)
--- a/include/edsam.h
+++ b/include/edsam.h
@@ -59,7 +59,7 @@ void dd_make_local_ed_indices(gmx_domdec_t *dd, gmx_edsam_t ed);
  * Should be called at every domain decomposition. */
  
 void do_flood(FILE *log, t_commrec *cr, rvec x[],rvec force[], gmx_edsam_t ed,
-        matrix box, gmx_large_int_t step);
+        matrix box, gmx_large_int_t step, gmx_bool bNS);
 /* Flooding - called from do_force() */
 
 #ifdef __cplusplus

diff --git a/src/mdlib/edsam.c b/src/mdlib/edsam.c
index 7ddb1ca444b1f375c24222cfd5e1f60e95c02003..c629ec60d39a173c251522cc2628f14d2eaaa548 100644 (file)
--- a/src/mdlib/edsam.c
+++ b/src/mdlib/edsam.c
@@ -857,7 +857,8 @@ static void do_single_flood(
         t_edpar *edi,
         gmx_large_int_t step,
         matrix box,
-        t_commrec *cr)
+        t_commrec *cr,
+        gmx_bool bNS)       /* Are we in a neighbor searching step? */
 {
     int i;
     matrix  rotmat;         /* rotation matrix */
@@ -868,15 +869,16 @@ static void do_single_flood(
 
     buf=edi->buf->do_edsam;
 
+
     /* Broadcast the positions of the AVERAGE structure such that they are known on
      * every processor. Each node contributes its local positions x and stores them in
      * the collective ED array buf->xcoll */
-    communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, buf->bUpdateShifts, x,
+    communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, bNS, x,
                     edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old, box);
 
     /* Only assembly REFERENCE positions if their indices differ from the average ones */
     if (!edi->bRefEqAv)
-        communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, buf->bUpdateShifts, x,
+        communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, bNS, x,
                 edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
 
     /* If bUpdateShifts was TRUE, the shifts have just been updated in get_positions.
@@ -935,7 +937,8 @@ extern void do_flood(
         rvec            force[], /* forcefield forces, to these the flooding forces are added */
         gmx_edsam_t     ed,      /* ed data structure contains all ED and flooding datasets */
         matrix          box,     /* the box */
-        gmx_large_int_t step)    /* The relative time step since ir->init_step is already subtracted */
+        gmx_large_int_t step,    /* The relative time step since ir->init_step is already subtracted */
+        gmx_bool        bNS)     /* Are we in a neighbor searching step? */
 {
     t_edpar *edi;
 
@@ -948,7 +951,7 @@ extern void do_flood(
     {
         /* Call flooding for one matrix */
         if (edi->flood.vecs.neig)
-            do_single_flood(ed->edo,x,force,edi,step,box,cr);
+            do_single_flood(ed->edo,x,force,edi,step,box,cr,bNS);
         edi = edi->next_edi;
     }
 }
@@ -1534,8 +1537,7 @@ static int read_edi(FILE* in, gmx_edsam_t ed,t_edpar *edi,int nr_mdatoms, int ed
     /* allocate space for reference positions and read them */
     snew(edi->sref.anrs,edi->sref.nr);
     snew(edi->sref.x   ,edi->sref.nr);
-    if (PAR(cr))
-        snew(edi->sref.x_old,edi->sref.nr);
+    snew(edi->sref.x_old,edi->sref.nr);
     edi->sref.sqrtm    =NULL;
     read_edx(in,edi->sref.nr,edi->sref.anrs,edi->sref.x);
 
@@ -1543,8 +1545,7 @@ static int read_edi(FILE* in, gmx_edsam_t ed,t_edpar *edi,int nr_mdatoms, int ed
     edi->sav.nr=read_checked_edint(in,"NAV");
     snew(edi->sav.anrs,edi->sav.nr);
     snew(edi->sav.x   ,edi->sav.nr);
-    if (PAR(cr))
-        snew(edi->sav.x_old,edi->sav.nr);
+    snew(edi->sav.x_old,edi->sav.nr);
     read_edx(in,edi->sav.nr,edi->sav.anrs,edi->sav.x);
 
     /* Check if the same atom indices are used for reference and average positions */
@@ -2214,10 +2215,9 @@ void init_edsam(gmx_mtop_t  *mtop,   /* global topology                    */
             {
                 copy_rvec(x_pbc[edi->sref.anrs[i]], xfit[i]);
 
-                /* Save the sref positions such that in the next time step the molecule can
-                 * be made whole again (in the parallel case) */
-                if (PAR(cr))
-                    copy_rvec(xfit[i], edi->sref.x_old[i]);
+                /* Save the sref positions such that in the next time step we can make the ED group whole
+                 * in case any of the atoms do not have the correct PBC representation */
+                copy_rvec(xfit[i], edi->sref.x_old[i]);
             }
 
             /* Extract the positions of the atoms subject to ED sampling */
@@ -2225,10 +2225,9 @@ void init_edsam(gmx_mtop_t  *mtop,   /* global topology                    */
             {
                 copy_rvec(x_pbc[edi->sav.anrs[i]], xstart[i]);
 
-                /* Save the sav positions such that in the next time step the molecule can
-                 * be made whole again (in the parallel case) */
-                if (PAR(cr))
-                    copy_rvec(xstart[i], edi->sav.x_old[i]);
+                /* Save the sav positions such that in the next time step we can make the ED group whole
+                 * in case any of the atoms do not have the correct PBC representation */
+                copy_rvec(xstart[i], edi->sav.x_old[i]);
             }
 
             /* Make the fit to the REFERENCE structure, get translation and rotation */
@@ -2245,6 +2244,14 @@ void init_edsam(gmx_mtop_t  *mtop,   /* global topology                    */
             /* calculate initial projections */
             project(xstart, edi);
 
+            /* For the target and origin structure both a reference (fit) and an
+             * average structure can be provided in make_edi. If both structures
+             * are the same, make_edi only stores one of them in the .edi file.
+             * If they differ, first the fit and then the average structure is stored
+             * in star (or sor), thus the number of entries in star/sor is
+             * (n_fit + n_av) with n_fit the size of the fitting group and n_av
+             * the size of the average group. */
+
             /* process target structure, if required */
             if (edi->star.nr > 0)
             {
@@ -2498,7 +2505,7 @@ void do_edsam(t_inputrec  *ir,
              * the collective buf->xcoll array. Note that for edinr > 1
              * xs could already have been modified by an earlier ED */
 
-            communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, buf->bUpdateShifts, xs,
+            communicate_group_positions(cr, buf->xcoll, buf->shifts_xcoll, buf->extra_shifts_xcoll, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
                     edi->sav.nr, edi->sav.nr_loc, edi->sav.anrs_loc, edi->sav.c_ind, edi->sav.x_old,  box);
 
 #ifdef DEBUG_ED
@@ -2506,11 +2513,12 @@ void do_edsam(t_inputrec  *ir,
 #endif
             /* Only assembly reference positions if their indices differ from the average ones */
             if (!edi->bRefEqAv)
-                communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, buf->bUpdateShifts, xs,
+                communicate_group_positions(cr, buf->xc_ref, buf->shifts_xc_ref, buf->extra_shifts_xc_ref, PAR(cr) ? buf->bUpdateShifts : TRUE, xs,
                         edi->sref.nr, edi->sref.nr_loc, edi->sref.anrs_loc, edi->sref.c_ind, edi->sref.x_old, box);
 
-            /* If bUpdateShifts was TRUE then the shifts have just been updated in get_positions.
-             * We do not need to uptdate the shifts until the next NS step */
+            /* If bUpdateShifts was TRUE then the shifts have just been updated in communicate_group_positions.
+             * We do not need to update the shifts until the next NS step. Note that dd_make_local_ed_indices
+             * set bUpdateShifts=TRUE in the parallel case. */
             buf->bUpdateShifts = FALSE;
 
             /* Now all nodes have all of the ED positions in edi->sav->xcoll,

diff --git a/src/mdlib/groupcoord.c b/src/mdlib/groupcoord.c
index 6ec672badb0e2c175ed28e9336fb1b9ca746d263..2f46d5c60fbcd609ea19daba028dddf0e21b8b51 100644 (file)
--- a/src/mdlib/groupcoord.c
+++ b/src/mdlib/groupcoord.c
@@ -195,7 +195,7 @@ extern void communicate_group_positions(
         rvec       *xcoll,        /* OUT: Collective array of positions */
         ivec       *shifts,       /* IN+OUT: Collective array of shifts for xcoll */
         ivec       *extra_shifts, /* BUF: Extra shifts since last time step */
-        const gmx_bool bNS,           /* IN:  NS step, the shifts have changed */
+        const gmx_bool bNS,       /* IN:  NS step, the shifts have changed */
         rvec       *x_loc,        /* IN:  Local positions on this node */ 
         const int  nr,            /* IN:  Total number of atoms in the group */
         const int  nr_loc,        /* IN:  Local number of atoms in the group */
@@ -221,35 +221,34 @@ extern void communicate_group_positions(
     {
         /* Add the arrays from all nodes together */
         gmx_sum(nr*3, xcoll[0], cr);
+    }
+    /* To make the group whole, start with a whole group and each
+     * step move the assembled positions at closest distance to the positions
+     * from the last step. First shift the positions with the saved shift
+     * vectors (these are 0 when this routine is called for the first time!) */
+    shift_positions_group(box, xcoll, shifts, nr);
+
+    /* Now check if some shifts changed since the last step.
+     * This only needs to be done when the shifts are expected to have changed,
+     * i.e. after neighboursearching */
+    if (bNS)
+    {
+        get_shifts_group(3, box, xcoll, nr, xcoll_old, extra_shifts);
 
-        /* To make the group whole, start with a whole group and each
-         * step move the assembled positions at closest distance to the positions 
-         * from the last step. First shift the positions with the saved shift 
-         * vectors (these are 0 when this routine is called for the first time!) */
-        shift_positions_group(box, xcoll, shifts, nr);
-        
-        /* Now check if some shifts changed since the last step.
-         * This only needs to be done when the shifts are expected to have changed,
-         * i.e. after neighboursearching */
-        if (bNS) 
+        /* Shift with the additional shifts such that we get a whole group now */
+        shift_positions_group(box, xcoll, extra_shifts, nr);
+
+        /* Add the shift vectors together for the next time step */
+        for (i=0; i<nr; i++)
         {
-            get_shifts_group(3, box, xcoll, nr, xcoll_old, extra_shifts);
-            
-            /* Shift with the additional shifts such that we get a whole group now */
-            shift_positions_group(box, xcoll, extra_shifts, nr);
-            
-            /* Add the shift vectors together for the next time step */
-            for (i=0; i<nr; i++)
-            {
-                shifts[i][XX] += extra_shifts[i][XX];
-                shifts[i][YY] += extra_shifts[i][YY];
-                shifts[i][ZZ] += extra_shifts[i][ZZ];
-            }
-            
-            /* Store current correctly-shifted positions for comparison in the next NS time step */
-            for (i=0; i<nr; i++)
-                copy_rvec(xcoll[i],xcoll_old[i]);   
+            shifts[i][XX] += extra_shifts[i][XX];
+            shifts[i][YY] += extra_shifts[i][YY];
+            shifts[i][ZZ] += extra_shifts[i][ZZ];
         }
+
+        /* Store current correctly-shifted positions for comparison in the next NS time step */
+        for (i=0; i<nr; i++)
+            copy_rvec(xcoll[i],xcoll_old[i]);
     }
     
     GMX_MPE_LOG(ev_get_group_x_finish);

diff --git a/src/mdlib/sim_util.c b/src/mdlib/sim_util.c
index 2d7fd82a742210b9b24fa90bb9382f07e89f6ed4..047dc9d4c2b1cef66333fd09f3c2909991ca803c 100644 (file)
--- a/src/mdlib/sim_util.c
+++ b/src/mdlib/sim_util.c
@@ -734,7 +734,7 @@ void do_force(FILE *fplog,t_commrec *cr,
     
     if (ed)
     {
-        do_flood(fplog,cr,x,f,ed,box,step);
+        do_flood(fplog,cr,x,f,ed,box,step,bNS);
     }
        
     if (DOMAINDECOMP(cr))