Merge release-5-0 into master
[alexxy/gromacs.git] / src / gromacs / mdlib / shellfc.c
index 8ddbc94859f6c3b1150e41d49a5398d89a422c78..ffe4469567ba13ef7bc539a0567b45be8fc3e539 100644 (file)
@@ -60,6 +60,7 @@
 
 #include "gromacs/math/vec.h"
 #include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/smalloc.h"
 
@@ -233,7 +234,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)
 {
@@ -287,11 +287,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;
 
@@ -535,6 +530,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;
         }
     }
@@ -995,26 +996,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);