Merge release-4-6 into release-5-0

[alexxy/gromacs.git] / src / gromacs / mdlib / shellfc.c

diff --git a/src/gromacs/mdlib/shellfc.c b/src/gromacs/mdlib/shellfc.c
index 531b88b61673dc519f888287f881b43600db2a5b..7c1f8679f44741ccb1d2d16ad8c5842106623e89 100644 (file)
--- a/src/gromacs/mdlib/shellfc.c
+++ b/src/gromacs/mdlib/shellfc.c
@@ -230,7 +230,6 @@ static void predict_shells(FILE *fplog, rvec x[], rvec v[], real dt,
 }
 
 gmx_shellfc_t init_shell_flexcon(FILE *fplog,
-                                 gmx_bool bCutoffSchemeIsVerlet,
                                  gmx_mtop_t *mtop, int nflexcon,
                                  rvec *x)
 {
@@ -284,11 +283,6 @@ gmx_shellfc_t init_shell_flexcon(FILE *fplog,
         return NULL;
     }
 
-    if (bCutoffSchemeIsVerlet)
-    {
-        gmx_fatal(FARGS, "The shell code does not work with the Verlet cut-off scheme.\n");
-    }
-
     snew(shfc, 1);
     shfc->nflexcon = nflexcon;
 
@@ -532,6 +526,12 @@ gmx_shellfc_t init_shell_flexcon(FILE *fplog,
             {
                 fprintf(fplog, "\nNOTE: there all shells that are connected to particles outside thier own charge group, will not predict shells positions during the run\n\n");
             }
+            /* Prediction improves performance, so we should implement either:
+             * 1. communication for the atoms needed for prediction
+             * 2. prediction using the velocities of shells; currently the
+             *    shell velocities are zeroed, it's a bit tricky to keep
+             *    track of the shell displacements and thus the velocity.
+             */
             shfc->bPredict = FALSE;
         }
     }
@@ -992,26 +992,31 @@ int relax_shell_flexcon(FILE *fplog, t_commrec *cr, gmx_bool bVerbose,
         force[i] = shfc->f[i];
     }
 
-    /* When we had particle decomposition, this code only worked with
-     * PD when all particles involved with each shell were in the same
-     * charge group. Not sure if this is still relevant. */
     if (bDoNS && inputrec->ePBC != epbcNONE && !DOMAINDECOMP(cr))
     {
         /* This is the only time where the coordinates are used
          * before do_force is called, which normally puts all
          * charge groups in the box.
          */
-        cg0 = 0;
-        cg1 = top->cgs.nr;
-        put_charge_groups_in_box(fplog, cg0, cg1, fr->ePBC, state->box,
-                                 &(top->cgs), state->x, fr->cg_cm);
+        if (inputrec->cutoff_scheme == ecutsVERLET)
+        {
+            put_atoms_in_box_omp(fr->ePBC, state->box, md->homenr, state->x);
+        }
+        else
+        {
+            cg0 = 0;
+            cg1 = top->cgs.nr;
+            put_charge_groups_in_box(fplog, cg0, cg1, fr->ePBC, state->box,
+                                     &(top->cgs), state->x, fr->cg_cm);
+        }
+
         if (graph)
         {
             mk_mshift(fplog, graph, fr->ePBC, state->box, state->x);
         }
     }
 
-    /* After this all coordinate arrays will contain whole molecules */
+    /* After this all coordinate arrays will contain whole charge groups */
     if (graph)
     {
         shift_self(graph, state->box, state->x);