Merge release-5-0 into master
authorRoland Schulz <roland@utk.edu>
Fri, 5 Dec 2014 23:31:17 +0000 (18:31 -0500)
committerRoland Schulz <roland@utk.edu>
Sun, 7 Dec 2014 00:08:19 +0000 (19:08 -0500)
Conflicts (trivial):
src/gromacs/gmxlib/gmx_thread_affinity.c
src/gromacs/mdlib/clincs.cpp

Manual changes:
Removed config.h because not necessary anymore:
src/gromacs/mdlib/nbnxn_search_simd_4xn.h
src/gromacs/mdlib/nbnxn_search.c
Removed config.h supression:
docs/doxygen/suppressions.txt
uncrustify:
src/gromacs/mdlib/clincs.cpp

Change-Id: I00064120f12f609fabecac9b3c823d33c015e59c

12 files changed:
docs/doxygen/suppressions.txt
share/top/amber99sb-ildn.ff/aminoacids.rtp
src/gromacs/domdec/domdec.cpp
src/gromacs/gmxlib/gmx_thread_affinity.c
src/gromacs/gmxpreprocess/readir.c
src/gromacs/gmxpreprocess/readpull.c
src/gromacs/mdlib/clincs.cpp
src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_common.h
src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_inner.h
src/gromacs/mdlib/nbnxn_pairlist.h
src/gromacs/mdlib/nbnxn_search.c
src/gromacs/mdlib/nbnxn_search_simd_4xn.h

index d1337b46dab06f6afe85173a414a23a45efa0ee5..71d9f29d23e6ac0dedbf67ee83ac07c9ed0f61ae 100644 (file)
@@ -24,7 +24,6 @@ src/gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_ref.c: warning: includes "config.h"
 src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn.c: warning: includes "config.h" unnecessarily
 src/gromacs/mdlib/nbnxn_kernels/simd_2xnn/nbnxn_kernel_simd_2xnn_common.h: warning: should include "config.h"
 src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn.c: warning: includes "config.h" unnecessarily
-src/gromacs/mdlib/nbnxn_kernels/simd_4xn/nbnxn_kernel_simd_4xn_inner.h: warning: should include "config.h"
 
 # These are specific to Folding@Home, and easiest to suppress here
 *: warning: includes non-local file as "corewrap.h"
index efb90779472040a2615229629968322e39fedb7e..70370089903a4617761bc70b3bca386228b5b2e2 100644 (file)
     -C     N
  [ dihedrals ]
      N    CA    CB   CG2        torsion_ILE_N_CA_CB_CG2_mult1
-     B    CA    CB   CG2        torsion_ILE_N_CA_CB_CG2_mult2
+     N    CA    CB   CG2        torsion_ILE_N_CA_CB_CG2_mult2
  [ impropers ]
-    -C    CA     N     H    
-    CA   OC1     C   OC2    
-                        
+    -C    CA     N     H
+    CA   OC1     C   OC2
+
 [ CVAL ]
  [ atoms ]
      N    N           -0.38210     1
      C    +N
  [ dihedrals ]
      N    CA    CB   CG2        torsion_ILE_N_CA_CB_CG2_mult1
-     B    CA    CB   CG2        torsion_ILE_N_CA_CB_CG2_mult2
+     N    CA    CB   CG2        torsion_ILE_N_CA_CB_CG2_mult2
  [ impropers ]
     CA    +N     C     O
-                        
+
 [ NVAL ]
  [ atoms ]
      N    N3           0.05770     1
index ba974d348f92b7afff6aa1703fc7e63de343fcb0..54fb8f7c35120334dcf036902d8ee1d5cd5872c4 100644 (file)
@@ -8737,21 +8737,22 @@ static void set_zones_size(gmx_domdec_t *dd,
             {
                 corner[ZZ] = zones->size[z].x1[ZZ];
             }
-            if (dd->ndim == 1 && box[ZZ][YY] != 0)
-            {
-                /* With 1D domain decomposition the cg's are not in
-                 * the triclinic box, but triclinic x-y and rectangular y-z.
-                 * Shift y back, so it will later end up at 0.
-                 */
-                corner[YY] -= corner[ZZ]*box[ZZ][YY]/box[ZZ][ZZ];
-            }
             /* Apply the triclinic couplings */
             assert(ddbox->npbcdim <= DIM);
             for (i = YY; i < ddbox->npbcdim; i++)
             {
                 for (j = XX; j < i; j++)
                 {
-                    corner[j] += corner[i]*box[i][j]/box[i][i];
+                    /* With 1D domain decomposition the cg's are not in
+                     * a triclinic box, but triclinic x-y and rectangular y/x-z.
+                     * So we should ignore the coupling for the non
+                     * domain-decomposed dimension of the pair x and y.
+                     */
+                    if (!(dd->ndim == 1 && ((dd->dim[0] == XX && j == YY) ||
+                                            (dd->dim[0] == YY && j == XX))))
+                    {
+                        corner[j] += corner[i]*box[i][j]/box[i][i];
+                    }
                 }
             }
             if (c == 0)
index 75e277999be1fcef991d5d14a30d2456a8449214..e4cc7386e1324c5e825e251ec8f25460f38cb1f8 100644 (file)
@@ -384,6 +384,9 @@ gmx_check_thread_affinity_set(FILE            *fplog,
     int       i, ret, cpu_count, cpu_set;
     gmx_bool  bAllSet;
 #endif
+#ifdef GMX_LIB_MPI
+    gmx_bool  bAllSet_All;
+#endif
 
     assert(hw_opt);
     if (!bAfterOpenmpInit)
@@ -458,6 +461,11 @@ gmx_check_thread_affinity_set(FILE            *fplog,
         bAllSet = bAllSet && (CPU_ISSET(i, &mask_current) != 0);
     }
 
+#ifdef GMX_LIB_MPI
+    MPI_Allreduce(&bAllSet, &bAllSet_All, 1, MPI_INT, MPI_LAND, MPI_COMM_WORLD);
+    bAllSet = bAllSet_All;
+#endif
+
     if (!bAllSet)
     {
         if (hw_opt->thread_affinity == threadaffAUTO)
index 8f55b6c6609358dcc0d19c872e9da8e9fdb006d3..9a9cc3ed3c23c06b56ea454806cf1dad7218aab2 100644 (file)
@@ -455,7 +455,7 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
         warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
     }
 
-    if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
+    if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rlist && ir->nstcalclr > 1)
     {
         warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
     }
@@ -1225,6 +1225,28 @@ void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
         sprintf(err_buf, "wall-ewald-zfac should be >= 2");
         CHECK(ir->wall_ewald_zfac < 2);
     }
+    if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) &&
+        EEL_FULL(ir->coulombtype))
+    {
+        sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
+                eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
+        warning(wi, warn_buf);
+    }
+    if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
+    {
+        if (ir->cutoff_scheme == ecutsVERLET)
+        {
+            sprintf(warn_buf, "Since molecules/charge groups are broken using the Verlet scheme, you can not use a dipole correction to the %s electrostatics.",
+                    eel_names[ir->coulombtype]);
+            warning(wi, warn_buf);
+        }
+        else
+        {
+            sprintf(warn_buf, "Dipole corrections to %s electrostatics only work if all charge groups that can cross PBC boundaries are dipoles. If this is not the case set epsilon_surface to 0",
+                    eel_names[ir->coulombtype]);
+            warning_note(wi, warn_buf);
+        }
+    }
 
     if (ir_vdw_switched(ir))
     {
@@ -2399,6 +2421,22 @@ void get_ir(const char *mdparin, const char *mdparout,
         {
             do_simtemp_params(ir);
         }
+
+        /* Because sc-coul (=FALSE by default) only acts on the lambda state
+         * setup and not on the old way of specifying the free-energy setup,
+         * we should check for using soft-core when not needed, since that
+         * can complicate the sampling significantly.
+         * Note that we only check for the automated coupling setup.
+         * If the (advanced) user does FEP through manual topology changes,
+         * this check will not be triggered.
+         */
+        if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 &&
+            ir->fepvals->sc_alpha != 0 &&
+            ((opts->couple_lam0 == ecouplamVDW  && opts->couple_lam0 == ecouplamVDWQ) ||
+             (opts->couple_lam1 == ecouplamVDWQ && opts->couple_lam1 == ecouplamVDW)))
+        {
+            warning(wi, "You are using soft-core interactions while the Van der Waals interactions are not decoupled (note that the sc-coul option is only active when using lambda states). Although this will not lead to errors, you will need much more sampling than without soft-core interactions. Consider using sc-alpha=0.");
+        }
     }
     else
     {
index 2e78022c5b1d32b148f31348616f69ab83af7e1a..e181244bf516fb04d18096bced06e5bf4dfdb890 100644 (file)
@@ -204,7 +204,8 @@ char **read_pullparams(int *ninp_p, t_inpfile **inp_p,
         nscan = sscanf(groups, "%d %d %d", &pcrd->group[0], &pcrd->group[1], &idum);
         if (nscan != 2)
         {
-            fprintf(stderr, "ERROR: %s should have %d components\n", buf, 2);
+            fprintf(stderr, "ERROR: %s should contain %d pull group indices\n",
+                    buf, 2);
             nerror++;
         }
         sprintf(buf, "pull-coord%d-origin", i);
@@ -248,7 +249,7 @@ void make_pull_groups(t_pull *pull,
 
         if (strcmp(pgnames[g], "") == 0)
         {
-            gmx_fatal(FARGS, "Group pull_group%d required by grompp was undefined.", g);
+            gmx_fatal(FARGS, "Pull option pull_group%d required by grompp has not been set.", g);
         }
 
         ig        = search_string(pgnames[g], grps->nr, gnames);
index 57ff2d7b8fb86f5924a748bdba5f6f0a6d92725b..aa4e633f77e0729350f3b990760068ca7f8bafa8 100644 (file)
@@ -67,6 +67,7 @@ typedef struct {
     int   *ind_r;      /* constraint index for updating atom data */
     int    ind_nalloc; /* allocation size of ind and ind_r */
     tensor vir_r_m_dr; /* temporary variable for virial calculation */
+    real   dhdlambda;  /* temporary variable for lambda derivative */
 } lincs_thread_t;
 
 typedef struct gmx_lincsdata {
@@ -360,7 +361,7 @@ static void lincs_update_atoms(struct gmx_lincsdata *li, int th,
 static void do_lincsp(rvec *x, rvec *f, rvec *fp, t_pbc *pbc,
                       struct gmx_lincsdata *lincsd, int th,
                       real *invmass,
-                      int econq, real *dvdlambda,
+                      int econq, gmx_bool bCalcDHDL,
                       gmx_bool bCalcVir, tensor rmdf)
 {
     int      b0, b1, b, i, j, k, n;
@@ -463,14 +464,17 @@ static void do_lincsp(rvec *x, rvec *f, rvec *fp, t_pbc *pbc,
     lincs_update_atoms(lincsd, th, 1.0, sol, r,
                        (econq != econqForce) ? invmass : NULL, fp);
 
-    if (dvdlambda != NULL)
+    if (bCalcDHDL)
     {
-#pragma omp barrier
+        real dhdlambda;
+
+        dhdlambda = 0;
         for (b = b0; b < b1; b++)
         {
-            *dvdlambda -= sol[b]*lincsd->ddist[b];
+            dhdlambda -= sol[b]*lincsd->ddist[b];
         }
-        /* 10 ncons flops */
+
+        lincsd->th[th].dhdlambda = dhdlambda;
     }
 
     if (bCalcVir)
@@ -499,7 +503,7 @@ static void do_lincs(rvec *x, rvec *xp, matrix box, t_pbc *pbc,
                      struct gmx_lincsdata *lincsd, int th,
                      real *invmass,
                      t_commrec *cr,
-                     gmx_bool bCalcLambda,
+                     gmx_bool bCalcDHDL,
                      real wangle, int *warn,
                      real invdt, rvec *v,
                      gmx_bool bCalcVir, tensor vir_r_m_dr)
@@ -685,9 +689,9 @@ static void do_lincs(rvec *x, rvec *xp, matrix box, t_pbc *pbc,
         /* 16 ncons flops */
     }
 
-    if (nlocat != NULL && bCalcLambda)
+    if (nlocat != NULL && (bCalcDHDL || bCalcVir))
     {
-        /* In lincs_update_atoms thread might cross-read mlambda */
+        /* In lincs_update_atoms threads might cross-read mlambda */
 #pragma omp barrier
 
         /* Only account for local atoms */
@@ -697,6 +701,21 @@ static void do_lincs(rvec *x, rvec *xp, matrix box, t_pbc *pbc,
         }
     }
 
+    if (bCalcDHDL)
+    {
+        real dhdl;
+
+        dhdl = 0;
+        for (b = b0; b < b1; b++)
+        {
+            /* Note that this this is dhdl*dt^2, the dt^2 factor is corrected
+             * later after the contributions are reduced over the threads.
+             */
+            dhdl -= lincsd->mlambda[b]*lincsd->ddist[b];
+        }
+        lincsd->th[th].dhdlambda = dhdl;
+    }
+
     if (bCalcVir)
     {
         /* Constraint virial */
@@ -1462,6 +1481,7 @@ gmx_bool constrain_lincs(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
                          t_nrnb *nrnb,
                          int maxwarn, int *warncount)
 {
+    gmx_bool  bCalcDHDL;
     char      buf[STRLEN], buf2[22], buf3[STRLEN];
     int       i, warn, p_imax;
     real      ncons_loc, p_ssd, p_max = 0;
@@ -1470,6 +1490,12 @@ gmx_bool constrain_lincs(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
 
     bOK = TRUE;
 
+    /* This boolean should be set by a flag passed to this routine.
+     * We can also easily check if any constraint length is changed,
+     * if not dH/dlambda=0 and we can also set the boolean to FALSE.
+     */
+    bCalcDHDL = (ir->efep != efepNO && dvdlambda != NULL);
+
     if (lincsd->nc == 0 && cr->dd == NULL)
     {
         if (bLog || bEner)
@@ -1491,6 +1517,9 @@ gmx_bool constrain_lincs(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
 
     if (econq == econqCoord)
     {
+        /* We can't use bCalcDHDL here, since NULL can be passed for dvdlambda
+         * also with efep!=fepNO.
+         */
         if (ir->efep != efepNO)
         {
             if (md->nMassPerturbed && lincsd->matlam != md->lambda)
@@ -1553,25 +1582,12 @@ gmx_bool constrain_lincs(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
 
             do_lincs(x, xprime, box, pbc, lincsd, th,
                      md->invmass, cr,
-                     bCalcVir || (ir->efep != efepNO),
+                     bCalcDHDL,
                      ir->LincsWarnAngle, &warn,
                      invdt, v, bCalcVir,
                      th == 0 ? vir_r_m_dr : lincsd->th[th].vir_r_m_dr);
         }
 
-        if (ir->efep != efepNO)
-        {
-            real dt_2, dvdl = 0;
-
-            /* TODO This should probably use invdt, so that sd integrator scaling works properly */
-            dt_2 = 1.0/(ir->delta_t*ir->delta_t);
-            for (i = 0; (i < lincsd->nc); i++)
-            {
-                dvdl -= lincsd->mlambda[i]*dt_2*lincsd->ddist[i];
-            }
-            *dvdlambda += dvdl;
-        }
-
         if (bLog && fplog && lincsd->nc > 0)
         {
             fprintf(fplog, "   Rel. Constraint Deviation:  RMS         MAX     between atoms\n");
@@ -1664,11 +1680,31 @@ gmx_bool constrain_lincs(FILE *fplog, gmx_bool bLog, gmx_bool bEner,
             int th = gmx_omp_get_thread_num();
 
             do_lincsp(x, xprime, min_proj, pbc, lincsd, th,
-                      md->invmass, econq, ir->efep != efepNO ? dvdlambda : NULL,
+                      md->invmass, econq, bCalcDHDL,
                       bCalcVir, th == 0 ? vir_r_m_dr : lincsd->th[th].vir_r_m_dr);
         }
     }
 
+    if (bCalcDHDL)
+    {
+        /* Reduce the dH/dlambda contributions over the threads */
+        real dhdlambda;
+        int  th;
+
+        dhdlambda = 0;
+        for (th = 0; th < lincsd->nth; th++)
+        {
+            dhdlambda += lincsd->th[th].dhdlambda;
+        }
+        if (econqCoord)
+        {
+            /* dhdlambda contains dH/dlambda*dt^2, correct for this */
+            /* TODO This should probably use invdt, so that sd integrator scaling works properly */
+            dhdlambda /= ir->delta_t*ir->delta_t;
+        }
+        *dvdlambda += dhdlambda;
+    }
+
     if (bCalcVir && lincsd->nth > 1)
     {
         for (i = 1; i < lincsd->nth; i++)
index 146c95855be950c76df0f85aa910b34dd6994de0..70d6edff67ce522c8c678f1f7a8cab28bd5163e7 100644 (file)
 #include "gromacs/mdlib/nbnxn_kernels/nbnxn_kernel_simd_utils.h"
 
 static gmx_inline void gmx_simdcall
-gmx_load_simd_4xn_interactions(int gmx_unused             excl,
+gmx_load_simd_4xn_interactions(int                        excl,
                                gmx_exclfilter gmx_unused  filter_S0,
                                gmx_exclfilter gmx_unused  filter_S1,
                                gmx_exclfilter gmx_unused  filter_S2,
                                gmx_exclfilter gmx_unused  filter_S3,
-                               const char gmx_unused     *interaction_mask_indices,
                                real gmx_unused           *simd_interaction_array,
                                gmx_simd_bool_t           *interact_S0,
                                gmx_simd_bool_t           *interact_S1,
@@ -79,13 +78,14 @@ gmx_load_simd_4xn_interactions(int gmx_unused             excl,
     *interact_S1  = gmx_checkbitmask_pb(mask_pr_S, filter_S1);
     *interact_S2  = gmx_checkbitmask_pb(mask_pr_S, filter_S2);
     *interact_S3  = gmx_checkbitmask_pb(mask_pr_S, filter_S3);
-#endif
-#ifdef GMX_SIMD_IBM_QPX
+#elif defined GMX_SIMD_IBM_QPX
     const int size = GMX_SIMD_REAL_WIDTH * sizeof(real);
-    *interact_S0  = gmx_load_interaction_mask_pb(size*interaction_mask_indices[0], simd_interaction_array);
-    *interact_S1  = gmx_load_interaction_mask_pb(size*interaction_mask_indices[1], simd_interaction_array);
-    *interact_S2  = gmx_load_interaction_mask_pb(size*interaction_mask_indices[2], simd_interaction_array);
-    *interact_S3  = gmx_load_interaction_mask_pb(size*interaction_mask_indices[3], simd_interaction_array);
+    *interact_S0  = gmx_load_interaction_mask_pb(size*((excl >> (0 * UNROLLJ)) & 0xF), simd_interaction_array);
+    *interact_S1  = gmx_load_interaction_mask_pb(size*((excl >> (1 * UNROLLJ)) & 0xF), simd_interaction_array);
+    *interact_S2  = gmx_load_interaction_mask_pb(size*((excl >> (2 * UNROLLJ)) & 0xF), simd_interaction_array);
+    *interact_S3  = gmx_load_interaction_mask_pb(size*((excl >> (3 * UNROLLJ)) & 0xF), simd_interaction_array);
+#else
+#error "Need implementation of gmx_load_simd_4xn_interactions"
 #endif
 }
 
index 64467a121c95bca9f9a22629816ac9b27e03435a..b465954dc4246cd1b66277a554d8c56e01510905 100644 (file)
     gmx_load_simd_4xn_interactions(l_cj[cjind].excl,
                                    filter_S0, filter_S1,
                                    filter_S2, filter_S3,
-#ifdef GMX_SIMD_IBM_QPX
-                                   l_cj[cjind].interaction_mask_indices,
                                    nbat->simd_interaction_array,
-#else
-                                   /* The struct fields do not exist
-                                      except on BlueGene/Q */
-                                   NULL,
-                                   NULL,
-#endif
                                    &interact_S0, &interact_S1,
                                    &interact_S2, &interact_S3);
 #endif /* CHECK_EXCLS */
index 4558266bc39a2c899a9ad0c763ee9edea58510b8..09d429a27d79ce456d69d0933ff5a038f3ec941d 100644 (file)
@@ -83,8 +83,6 @@ typedef void nbnxn_free_t (void *ptr);
 typedef struct {
     int          cj;    /* The j-cluster                    */
     unsigned int excl;  /* The exclusion (interaction) bits */
-    /* Indices into the arrays of SIMD interaction masks. */
-    char         interaction_mask_indices[4];
 } nbnxn_cj_t;
 
 /* In nbnxn_ci_t the integer shift contains the shift in the lower 7 bits.
index 856ae32d0f0d8c00a97f736ab14d5d66b64026d1..f7c27be2f28fd8f310acd44252ceb5ba9527a79b 100644 (file)
@@ -37,8 +37,6 @@
 
 #include "nbnxn_search.h"
 
-#include "config.h"
-
 #include <assert.h>
 #include <math.h>
 #include <string.h>
@@ -3258,12 +3256,6 @@ static void set_ci_top_excls(const nbnxn_search_t nbs,
                         inner_e = ge - (se << na_cj_2log);
 
                         nbl->cj[found].excl &= ~(1U<<((inner_i<<na_cj_2log) + inner_e));
-/* The next code line is usually not needed. We do not want to version
- * away the above line, because there is logic that relies on being
- * able to detect easily whether any exclusions exist. */
-#if (defined GMX_SIMD_IBM_QPX)
-                        nbl->cj[found].interaction_mask_indices[inner_i] &= ~(1U << inner_e);
-#endif
                     }
                 }
             }
index 33a307b58456ab8080b37e6c0309915d0f148217..4e1c15a0ee65c8b19e4bc59d12345db7971d06f4 100644 (file)
@@ -33,8 +33,6 @@
  * the research papers on the package. Check out http://www.gromacs.org.
  */
 
-#include "config.h"
-
 #if GMX_SIMD_REAL_WIDTH >= NBNXN_CPU_CLUSTER_I_SIZE
 #define STRIDE_S  (GMX_SIMD_REAL_WIDTH)
 #else
@@ -265,12 +263,6 @@ make_cluster_list_simd_4xn(const nbnxn_grid_t *gridj,
             /* Store cj and the interaction mask */
             nbl->cj[nbl->ncj].cj   = CI_TO_CJ_SIMD_4XN(gridj->cell0) + cj;
             nbl->cj[nbl->ncj].excl = get_imask_simd_4xn(remove_sub_diag, ci, cj);
-#ifdef GMX_SIMD_IBM_QPX
-            nbl->cj[nbl->ncj].interaction_mask_indices[0] = (nbl->cj[nbl->ncj].excl & 0x000F) >> (0 * 4);
-            nbl->cj[nbl->ncj].interaction_mask_indices[1] = (nbl->cj[nbl->ncj].excl & 0x00F0) >> (1 * 4);
-            nbl->cj[nbl->ncj].interaction_mask_indices[2] = (nbl->cj[nbl->ncj].excl & 0x0F00) >> (2 * 4);
-            nbl->cj[nbl->ncj].interaction_mask_indices[3] = (nbl->cj[nbl->ncj].excl & 0xF000) >> (3 * 4);
-#endif
             nbl->ncj++;
         }
         /* Increase the closing index in i super-cell list */