Made shells work with the Verlet scheme

The issue was that atoms iso charge group should be put in the box.
Additionally water_pol only worked within a charge group, now fixed.
Note that with shells working across charge groups, which is always
the case with the Verlet scheme, no shell prediction is done.
We should implement prediction, since it improves performance.

Fixes #1429.

Change-Id: I2ebfc2d91fc161167f8f2573b61e1f519cf11fd8

diff --git a/include/shellfc.h b/include/shellfc.h
index 8c11aedb155f2bcb8e7f7c420146ee1ac2b43ed7..87e291d2fa8a5bb63c93b25673215dce087c17f8 100644 (file)
--- a/include/shellfc.h
+++ b/include/shellfc.h
@@ -48,7 +48,6 @@ extern "C" {
  */
 GMX_LIBMD_EXPORT
 gmx_shellfc_t init_shell_flexcon(FILE *fplog,
-                                 gmx_bool bCutoffSchemeIsVerlet,
                                  gmx_mtop_t *mtop, int nflexcon,
                                  rvec *x);
 

diff --git a/src/gmxlib/bondfree.c b/src/gmxlib/bondfree.c
index 66a91d4ac0d26a68b58ad6193cb210964a474e86..2c25351ec20380e4eb27ad9e2c535f451ba0c5f4 100644 (file)
--- a/src/gmxlib/bondfree.c
+++ b/src/gmxlib/bondfree.c
@@ -736,7 +736,8 @@ real water_pol(int nbonds,
      * a shell connected to a dummy with spring constant that differ in the
      * three spatial dimensions in the molecular frame.
      */
-    int  i, m, aO, aH1, aH2, aD, aS, type, type0;
+    int  i, m, aO, aH1, aH2, aD, aS, type, type0, ki;
+    ivec dt;
     rvec dOH1, dOH2, dHH, dOD, dDS, nW, kk, dx, kdx, proj;
 #ifdef DEBUG
     rvec df;
@@ -779,11 +780,11 @@ real water_pol(int nbonds,
             aS   = forceatoms[i+5];
 
             /* Compute vectors describing the water frame */
-            rvec_sub(x[aH1], x[aO], dOH1);
-            rvec_sub(x[aH2], x[aO], dOH2);
-            rvec_sub(x[aH2], x[aH1], dHH);
-            rvec_sub(x[aD], x[aO], dOD);
-            rvec_sub(x[aS], x[aD], dDS);
+            pbc_rvec_sub(pbc, x[aH1], x[aO], dOH1);
+            pbc_rvec_sub(pbc, x[aH2], x[aO], dOH2);
+            pbc_rvec_sub(pbc, x[aH2], x[aH1], dHH);
+            pbc_rvec_sub(pbc, x[aD], x[aO], dOD);
+            ki = pbc_rvec_sub(pbc, x[aS], x[aD], dDS);
             cprod(dOH1, dOH2, nW);
 
             /* Compute inverse length of normal vector
@@ -838,6 +839,13 @@ real water_pol(int nbonds,
             kdx[YY] = kk[YY]*dx[YY];
             kdx[ZZ] = kk[ZZ]*dx[ZZ];
             vtot   += iprod(dx, kdx);
+
+            if (g)
+            {
+                ivec_sub(SHIFT_IVEC(g, aS), SHIFT_IVEC(g, aD), dt);
+                ki = IVEC2IS(dt);
+            }
+
             for (m = 0; (m < DIM); m++)
             {
                 /* This is a tensor operation but written out for speed */
@@ -848,8 +856,10 @@ real water_pol(int nbonds,
 #ifdef DEBUG
                 df[m] = fij;
 #endif
-                f[aS][m] += fij;
-                f[aD][m] -= fij;
+                f[aS][m]           += fij;
+                f[aD][m]           -= fij;
+                fshift[ki][m]      += fij;
+                fshift[CENTRAL][m] -= fij;
             }
 #ifdef DEBUG
             if (debug)

diff --git a/src/kernel/md.c b/src/kernel/md.c
index da68da10f2bc219b4bc67902be786f81529dc0a4..8909fb8ed759fa31d43cc474ecf6b3346e55f406 100644 (file)
--- a/src/kernel/md.c
+++ b/src/kernel/md.c
@@ -342,7 +342,7 @@ double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     debug_gmx();
 
     /* Check for polarizable models and flexible constraints */
-    shellfc = init_shell_flexcon(fplog, fr->cutoff_scheme == ecutsVERLET,
+    shellfc = init_shell_flexcon(fplog,
                                  top_global, n_flexible_constraints(constr),
                                  (ir->bContinuation ||
                                   (DOMAINDECOMP(cr) && !MASTER(cr))) ?

diff --git a/src/mdlib/shellfc.c b/src/mdlib/shellfc.c
index 68e6c3be5e6e3287f66912f16b53a482e59a0ae7..0acd97530215a7272ad0a650b4fdbf98fef04497 100644 (file)
--- a/src/mdlib/shellfc.c
+++ b/src/mdlib/shellfc.c
@@ -231,7 +231,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)
 {
@@ -285,11 +284,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;
 
@@ -533,6 +527,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;
         }
     }
@@ -1001,34 +1001,31 @@ int relax_shell_flexcon(FILE *fplog, t_commrec *cr, gmx_bool bVerbose,
         force[i] = shfc->f[i];
     }
 
-    /* With particle decomposition this code only works
-     * when all particles involved with each shell are in the same cg.
-     */
-
     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.
          */
-        if (PARTDECOMP(cr))
+        if (inputrec->cutoff_scheme == ecutsVERLET)
         {
-            pd_cg_range(cr, &cg0, &cg1);
+            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);
         }
-        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);