Merge origin/release-4-6 into master
[alexxy/gromacs.git] / src / gromacs / mdlib / ns.c
index 811a92f26ed737af1780c712458d7196c174a990..208773dc99a827fc13bb492cb7f90aed6e9caf2d 100644 (file)
@@ -17,7 +17,7 @@
  * modify it under the terms of the GNU General Public License
  * as published by the Free Software Foundation; either version 2
  * of the License, or (at your option) any later version.
- * 
+ *
  * If you want to redistribute modifications, please consider that
  * scientific software is very special. Version control is crucial -
  * bugs must be traceable. We will be happy to consider code for
@@ -94,11 +94,11 @@ static void reallocate_nblist(t_nblist *nl)
 {
     if (gmx_debug_at)
     {
-        fprintf(debug,"reallocating neigborlist il_code=%d, maxnri=%d\n",
-                nl->il_code,nl->maxnri); 
+        fprintf(debug,"reallocating neigborlist (ielec=%d, ivdw=%d, igeometry=%d, free_energy=%d), maxnri=%d\n",
+                nl->ielec,nl->ivdw,nl->igeometry,nl->free_energy,nl->maxnri);
     }
     srenew(nl->iinr,   nl->maxnri);
-    if (nl->enlist == enlistCG_CG)
+    if (nl->igeometry == GMX_NBLIST_GEOMETRY_CG_CG)
     {
         srenew(nl->iinr_end,nl->maxnri);
     }
@@ -107,60 +107,17 @@ static void reallocate_nblist(t_nblist *nl)
     srenew(nl->jindex, nl->maxnri+1);
 }
 
-/* ivdw/icoul are used to determine the type of interaction, so we
- * can set an innerloop index here. The obvious choice for this would have
- * been the vdwtype/coultype values in the forcerecord, but unfortunately 
- * those types are braindead - for instance both Buckingham and normal 
- * Lennard-Jones use the same value (evdwCUT), and a separate gmx_boolean variable
- * to determine which interaction is used. There is further no special value
- * for 'no interaction'. For backward compatibility with old TPR files we won't
- * change this in the 3.x series, so when calling this routine you should use:
- *
- * icoul=0 no coulomb interaction
- * icoul=1 cutoff standard coulomb
- * icoul=2 reaction-field coulomb
- * icoul=3 tabulated coulomb
- *
- * ivdw=0 no vdw interaction
- * ivdw=1 standard L-J interaction
- * ivdw=2 Buckingham
- * ivdw=3 tabulated vdw.
- *
- * Kind of ugly, but it works.
- */
-static void init_nblist(t_nblist *nl_sr,t_nblist *nl_lr,
+
+static void init_nblist(FILE *log, t_nblist *nl_sr,t_nblist *nl_lr,
                         int maxsr,int maxlr,
-                        int ivdw, int icoul, 
-                        gmx_bool bfree, int enlist)
+                        int ivdw, int ivdwmod,
+                        int ielec, int ielecmod,
+                        gmx_bool bfree, int igeometry)
 {
     t_nblist *nl;
     int      homenr;
     int      i,nn;
     
-    int inloop[20] =
-    { 
-        eNR_NBKERNEL_NONE,
-        eNR_NBKERNEL010,
-        eNR_NBKERNEL020,
-        eNR_NBKERNEL030,
-        eNR_NBKERNEL100,
-        eNR_NBKERNEL110,
-        eNR_NBKERNEL120,
-        eNR_NBKERNEL130,
-        eNR_NBKERNEL200,
-        eNR_NBKERNEL210,
-        eNR_NBKERNEL220,
-        eNR_NBKERNEL230,
-        eNR_NBKERNEL300,
-        eNR_NBKERNEL310,
-        eNR_NBKERNEL320,
-        eNR_NBKERNEL330,
-        eNR_NBKERNEL400,
-        eNR_NBKERNEL410,
-        eNR_NBKERNEL_NONE,
-        eNR_NBKERNEL430
-    };
-  
     for(i=0; (i<2); i++)
     {
         nl     = (i == 0) ? nl_sr : nl_lr;
@@ -170,47 +127,24 @@ static void init_nblist(t_nblist *nl_sr,t_nblist *nl_lr,
         {
             continue;
         }
-        
+
+
         /* Set coul/vdw in neighborlist, and for the normal loops we determine
          * an index of which one to call.
          */
-        nl->ivdw  = ivdw;
-        nl->icoul = icoul;
+        nl->ivdw        = ivdw;
+        nl->ivdwmod     = ivdwmod;
+        nl->ielec       = ielec;
+        nl->ielecmod    = ielecmod;
         nl->free_energy = bfree;
-    
+        nl->igeometry   = igeometry;
+
         if (bfree)
         {
-            nl->enlist  = enlistATOM_ATOM;
-            nl->il_code = eNR_NBKERNEL_FREE_ENERGY;
+            nl->igeometry  = GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE;
         }
-        else
-        {
-            nl->enlist = enlist;
-
-            nn = inloop[4*icoul + ivdw];
-            
-            /* solvent loops follow directly after the corresponding
-            * ordinary loops, in the order:
-            *
-            * SPC, SPC-SPC, TIP4p, TIP4p-TIP4p
-            *   
-            */
-            switch (enlist) {
-            case enlistATOM_ATOM:
-            case enlistCG_CG:
-                break;
-            case enlistSPC_ATOM:     nn += 1; break;
-            case enlistSPC_SPC:      nn += 2; break;
-            case enlistTIP4P_ATOM:   nn += 3; break;
-            case enlistTIP4P_TIP4P:  nn += 4; break;
-            }
-            
-            nl->il_code = nn;
-        }
-
-        if (debug)
-            fprintf(debug,"Initiating neighbourlist type %d for %s interactions,\nwith %d SR, %d LR atoms.\n",
-                    nl->il_code,ENLISTTYPE(enlist),maxsr,maxlr);
+        
+        gmx_nonbonded_set_kernel_pointers( (i==0) ? log : NULL,nl);
         
         /* maxnri is influenced by the number of shifts (maximum is 8)
          * and the number of energy groups.
@@ -229,6 +163,13 @@ static void init_nblist(t_nblist *nl_sr,t_nblist *nl_lr,
         nl->jindex      = NULL;
         reallocate_nblist(nl);
         nl->jindex[0] = 0;
+
+        if(debug)
+        {
+            fprintf(debug,"Initiating neighbourlist (ielec=%d, ivdw=%d, free=%d) for %s interactions,\nwith %d SR, %d LR atoms.\n",
+                    nl->ielec,nl->ivdw,nl->free_energy,gmx_nblist_geometry_names[nl->igeometry],maxsr,maxlr);
+        }
+
 #ifdef GMX_THREAD_SHM_FDECOMP
         nl->counter = 0;
         snew(nl->mtx,1);
@@ -245,9 +186,9 @@ void init_neighbor_list(FILE *log,t_forcerec *fr,int homenr)
     * cache trashing.
     */
    int maxsr,maxsr_wat,maxlr,maxlr_wat;
-   int icoul,icoulf,ivdw;
+   int ielec,ielecf,ivdw,ielecmod,ielecmodf,ivdwmod;
    int solvent;
-   int enlist_def,enlist_w,enlist_ww;
+   int igeometry_def,igeometry_w,igeometry_ww;
    int i;
    t_nblists *nbl;
 
@@ -274,46 +215,20 @@ void init_neighbor_list(FILE *log,t_forcerec *fr,int homenr)
      maxlr = maxlr_wat = 0;
    }  
 
-   /* Determine the values for icoul/ivdw. */
-   /* Start with GB */
-   if(fr->bGB)
-   {
-       icoul=enbcoulGB;
-   }
-   else if (fr->bcoultab)
-   {
-       icoul = enbcoulTAB;
-   }
-   else if (EEL_RF(fr->eeltype))
-   {
-       icoul = enbcoulRF;
-   }
-   else 
-   {
-       icoul = enbcoulOOR;
-   }
-   
-   if (fr->bvdwtab)
-   {
-       ivdw = enbvdwTAB;
-   }
-   else if (fr->bBHAM)
-   {
-       ivdw = enbvdwBHAM;
-   }
-   else 
-   {
-       ivdw = enbvdwLJ;
-   }
+   /* Determine the values for ielec/ivdw. */
+   ielec = fr->nbkernel_elec_interaction;
+   ivdw  = fr->nbkernel_vdw_interaction;
+   ielecmod = fr->nbkernel_elec_modifier;
+   ivdwmod  = fr->nbkernel_vdw_modifier;
 
    fr->ns.bCGlist = (getenv("GMX_NBLISTCG") != 0);
    if (!fr->ns.bCGlist)
    {
-       enlist_def = enlistATOM_ATOM;
+       igeometry_def = GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE;
    }
    else
    {
-       enlist_def = enlistCG_CG;
+       igeometry_def = GMX_NBLIST_GEOMETRY_CG_CG;
        if (log != NULL)
        {
            fprintf(log,"\nUsing charge-group - charge-group neighbor lists and kernels\n\n");
@@ -321,55 +236,75 @@ void init_neighbor_list(FILE *log,t_forcerec *fr,int homenr)
    }
    
    if (fr->solvent_opt == esolTIP4P) {
-       enlist_w  = enlistTIP4P_ATOM;
-       enlist_ww = enlistTIP4P_TIP4P;
+       igeometry_w  = GMX_NBLIST_GEOMETRY_WATER4_PARTICLE;
+       igeometry_ww = GMX_NBLIST_GEOMETRY_WATER4_WATER4;
    } else {
-       enlist_w  = enlistSPC_ATOM;
-       enlist_ww = enlistSPC_SPC;
+       igeometry_w  = GMX_NBLIST_GEOMETRY_WATER3_PARTICLE;
+       igeometry_ww = GMX_NBLIST_GEOMETRY_WATER3_WATER3;
    }
 
    for(i=0; i<fr->nnblists; i++) 
    {
        nbl = &(fr->nblists[i]);
-       init_nblist(&nbl->nlist_sr[eNL_VDWQQ],&nbl->nlist_lr[eNL_VDWQQ],
-                   maxsr,maxlr,ivdw,icoul,FALSE,enlist_def);
-       init_nblist(&nbl->nlist_sr[eNL_VDW],&nbl->nlist_lr[eNL_VDW],
-                   maxsr,maxlr,ivdw,0,FALSE,enlist_def);
-       init_nblist(&nbl->nlist_sr[eNL_QQ],&nbl->nlist_lr[eNL_QQ],
-                   maxsr,maxlr,0,icoul,FALSE,enlist_def);
-       init_nblist(&nbl->nlist_sr[eNL_VDWQQ_WATER],&nbl->nlist_lr[eNL_VDWQQ_WATER],
-                   maxsr_wat,maxlr_wat,ivdw,icoul, FALSE,enlist_w);
-       init_nblist(&nbl->nlist_sr[eNL_QQ_WATER],&nbl->nlist_lr[eNL_QQ_WATER],
-                   maxsr_wat,maxlr_wat,0,icoul, FALSE,enlist_w);
-       init_nblist(&nbl->nlist_sr[eNL_VDWQQ_WATERWATER],&nbl->nlist_lr[eNL_VDWQQ_WATERWATER],
-                   maxsr_wat,maxlr_wat,ivdw,icoul, FALSE,enlist_ww);
-       init_nblist(&nbl->nlist_sr[eNL_QQ_WATERWATER],&nbl->nlist_lr[eNL_QQ_WATERWATER],
-                   maxsr_wat,maxlr_wat,0,icoul, FALSE,enlist_ww);
+       init_nblist(log,&nbl->nlist_sr[eNL_VDWQQ],&nbl->nlist_lr[eNL_VDWQQ],
+                   maxsr,maxlr,ivdw,ivdwmod,ielec,ielecmod,FALSE,igeometry_def);
+       init_nblist(log,&nbl->nlist_sr[eNL_VDW],&nbl->nlist_lr[eNL_VDW],
+                   maxsr,maxlr,ivdw,ivdwmod,GMX_NBKERNEL_ELEC_NONE,eintmodNONE,FALSE,igeometry_def);
+       init_nblist(log,&nbl->nlist_sr[eNL_QQ],&nbl->nlist_lr[eNL_QQ],
+                   maxsr,maxlr,GMX_NBKERNEL_VDW_NONE,eintmodNONE,ielec,ielecmod,FALSE,igeometry_def);
+       init_nblist(log,&nbl->nlist_sr[eNL_VDWQQ_WATER],&nbl->nlist_lr[eNL_VDWQQ_WATER],
+                   maxsr_wat,maxlr_wat,ivdw,ivdwmod,ielec,ielecmod, FALSE,igeometry_w);
+       init_nblist(log,&nbl->nlist_sr[eNL_QQ_WATER],&nbl->nlist_lr[eNL_QQ_WATER],
+                   maxsr_wat,maxlr_wat,GMX_NBKERNEL_VDW_NONE,eintmodNONE,ielec,ielecmod, FALSE,igeometry_w);
+       init_nblist(log,&nbl->nlist_sr[eNL_VDWQQ_WATERWATER],&nbl->nlist_lr[eNL_VDWQQ_WATERWATER],
+                   maxsr_wat,maxlr_wat,ivdw,ivdwmod,ielec,ielecmod, FALSE,igeometry_ww);
+       init_nblist(log,&nbl->nlist_sr[eNL_QQ_WATERWATER],&nbl->nlist_lr[eNL_QQ_WATERWATER],
+                   maxsr_wat,maxlr_wat,GMX_NBKERNEL_VDW_NONE,eintmodNONE,ielec,ielecmod, FALSE,igeometry_ww);
+
+       /* Did we get the solvent loops so we can use optimized water kernels? */
+       if(nbl->nlist_sr[eNL_VDWQQ_WATER].kernelptr_vf==NULL
+          || nbl->nlist_sr[eNL_QQ_WATER].kernelptr_vf==NULL
+#ifndef DISABLE_WATERWATER_NLIST
+          || nbl->nlist_sr[eNL_VDWQQ_WATERWATER].kernelptr_vf==NULL
+          || nbl->nlist_sr[eNL_QQ_WATERWATER].kernelptr_vf==NULL
+#endif
+          )
+       {
+           fr->solvent_opt = esolNO;
+           fprintf(log,"Note: The available nonbonded kernels do not support water optimization - disabling.\n");
+       }
        
        if (fr->efep != efepNO) 
        {
            if ((fr->bEwald) && (fr->sc_alphacoul > 0)) /* need to handle long range differently if using softcore */
            {
-               icoulf = enbcoulFEWALD;
+               ielecf = GMX_NBKERNEL_ELEC_EWALD;
+               ielecmodf = eintmodNONE;
            }
            else
            {
-               icoulf = icoul;
+               ielecf = ielec;
+               ielecmodf = ielecmod;
            }
 
-           init_nblist(&nbl->nlist_sr[eNL_VDWQQ_FREE],&nbl->nlist_lr[eNL_VDWQQ_FREE],
-                       maxsr,maxlr,ivdw,icoulf,TRUE,enlistATOM_ATOM);
-           init_nblist(&nbl->nlist_sr[eNL_VDW_FREE],&nbl->nlist_lr[eNL_VDW_FREE],
-                       maxsr,maxlr,ivdw,0,TRUE,enlistATOM_ATOM);
-           init_nblist(&nbl->nlist_sr[eNL_QQ_FREE],&nbl->nlist_lr[eNL_QQ_FREE],
-                       maxsr,maxlr,0,icoulf,TRUE,enlistATOM_ATOM);
+           init_nblist(log,&nbl->nlist_sr[eNL_VDWQQ_FREE],&nbl->nlist_lr[eNL_VDWQQ_FREE],
+                       maxsr,maxlr,ivdw,ivdwmod,ielecf,ielecmod,TRUE,GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE);
+           init_nblist(log,&nbl->nlist_sr[eNL_VDW_FREE],&nbl->nlist_lr[eNL_VDW_FREE],
+                       maxsr,maxlr,ivdw,ivdwmod,GMX_NBKERNEL_ELEC_NONE,eintmodNONE,TRUE,GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE);
+           init_nblist(log,&nbl->nlist_sr[eNL_QQ_FREE],&nbl->nlist_lr[eNL_QQ_FREE],
+                       maxsr,maxlr,GMX_NBKERNEL_VDW_NONE,eintmodNONE,ielecf,ielecmod,TRUE,GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE);
        }  
    }
    /* QMMM MM list */
    if (fr->bQMMM && fr->qr->QMMMscheme != eQMMMschemeoniom)
    {
-       init_nblist(&fr->QMMMlist,NULL,
-                   maxsr,maxlr,0,icoul,FALSE,enlistATOM_ATOM);
+       init_nblist(log,&fr->QMMMlist,NULL,
+                   maxsr,maxlr,0,0,ielec,ielecmod,FALSE,GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE);
+   }
+
+   if(log!=NULL)
+   {
+       fprintf(log,"\n");
    }
 
    fr->ns.nblist_initialized=TRUE;
@@ -386,27 +321,28 @@ static void reset_nblist(t_nblist *nl)
      }
 }
 
-static void reset_neighbor_list(t_forcerec *fr,gmx_bool bLR,int nls,int eNL)
+static void reset_neighbor_lists(t_forcerec *fr,gmx_bool bResetSR, gmx_bool bResetLR)
 {
     int n,i;
   
-    if (bLR) 
+    if (fr->bQMMM)
     {
-        reset_nblist(&(fr->nblists[nls].nlist_lr[eNL]));
+        /* only reset the short-range nblist */
+        reset_nblist(&(fr->QMMMlist));
     }
-    else 
+
+    for(n=0; n<fr->nnblists; n++)
     {
-        for(n=0; n<fr->nnblists; n++)
+        for(i=0; i<eNL_NR; i++)
         {
-            for(i=0; i<eNL_NR; i++)
+            if(bResetSR)
             {
-                reset_nblist(&(fr->nblists[n].nlist_sr[i]));
+                reset_nblist( &(fr->nblists[n].nlist_sr[i]) );
+            }
+            if(bResetLR)
+            {
+                reset_nblist( &(fr->nblists[n].nlist_lr[i]) );
             }
-        }
-        if (fr->bQMMM)
-        { 
-            /* only reset the short-range nblist */
-            reset_nblist(&(fr->QMMMlist));
         }
     }
 }
@@ -494,36 +430,26 @@ static inline void close_nblist(t_nblist *nlist)
     }
 }
 
-static inline void close_neighbor_list(t_forcerec *fr,gmx_bool bLR,int nls,int eNL, 
-                                       gmx_bool bMakeQMMMnblist)
+static inline void close_neighbor_lists(t_forcerec *fr,gmx_bool bMakeQMMMnblist)
 {
     int n,i;
     
-    if (bMakeQMMMnblist) {
-        if (!bLR)
-        {
+    if (bMakeQMMMnblist)
+    {
             close_nblist(&(fr->QMMMlist));
-        }
     }
-    else 
+
+    for(n=0; n<fr->nnblists; n++)
     {
-        if (bLR)
+        for(i=0; (i<eNL_NR); i++)
         {
-            close_nblist(&(fr->nblists[nls].nlist_lr[eNL]));
-        }
-        else
-        { 
-            for(n=0; n<fr->nnblists; n++)
-            {
-                for(i=0; (i<eNL_NR); i++)
-                {
-                    close_nblist(&(fr->nblists[n].nlist_sr[i]));
-                }
-            }
+            close_nblist(&(fr->nblists[n].nlist_sr[i]));
+            close_nblist(&(fr->nblists[n].nlist_lr[i]));
         }
     }
 }
 
+
 static inline void add_j_to_nblist(t_nblist *nlist,atom_id j_atom,gmx_bool bLR)
 {
     int nrj=nlist->nrj;
@@ -532,8 +458,8 @@ static inline void add_j_to_nblist(t_nblist *nlist,atom_id j_atom,gmx_bool bLR)
     {
         nlist->maxnrj = over_alloc_small(nlist->nrj + 1);
         if (gmx_debug_at)
-            fprintf(debug,"Increasing %s nblist %s j size to %d\n",
-                    bLR ? "LR" : "SR",nrnb_str(nlist->il_code),nlist->maxnrj);
+            fprintf(debug,"Increasing %s nblist (ielec=%d,ivdw=%d,free=%d,igeometry=%d) j size to %d\n",
+                    bLR ? "LR" : "SR",nlist->ielec,nlist->ivdw,nlist->free_energy,nlist->igeometry,nlist->maxnrj);
         
         srenew(nlist->jjnr,nlist->maxnrj);
     }
@@ -554,8 +480,8 @@ static inline void add_j_to_nblist_cg(t_nblist *nlist,
     {
         nlist->maxnrj = over_alloc_small(nlist->nrj + 1);
         if (gmx_debug_at)
-            fprintf(debug,"Increasing %s nblist %s j size to %d\n",
-                    bLR ? "LR" : "SR",nrnb_str(nlist->il_code),nlist->maxnrj);
+            fprintf(debug,"Increasing %s nblist (ielec=%d,ivdw=%d,free=%d,igeometry=%d) j size to %d\n",
+                    bLR ? "LR" : "SR",nlist->ielec,nlist->ivdw,nlist->free_energy,nlist->igeometry,nlist->maxnrj);
         
         srenew(nlist->jjnr    ,nlist->maxnrj);
         srenew(nlist->jjnr_end,nlist->maxnrj);
@@ -599,9 +525,10 @@ put_in_list_t(gmx_bool              bHaveVdW[],
               t_excl            bExcl[],
               int               shift,
               t_forcerec *      fr,
-              gmx_bool              bLR,
-              gmx_bool              bDoVdW,
-              gmx_bool              bDoCoul);
+              gmx_bool          bLR,
+              gmx_bool          bDoVdW,
+              gmx_bool          bDoCoul,
+              int               solvent_opt);
 
 static void 
 put_in_list_at(gmx_bool              bHaveVdW[],
@@ -615,9 +542,10 @@ put_in_list_at(gmx_bool              bHaveVdW[],
                t_excl            bExcl[],
                int               shift,
                t_forcerec *      fr,
-               gmx_bool              bLR,
-               gmx_bool              bDoVdW,
-               gmx_bool              bDoCoul)
+               gmx_bool          bLR,
+               gmx_bool          bDoVdW,
+               gmx_bool          bDoCoul,
+               int               solvent_opt)
 {
     /* The a[] index has been removed,
      * to put it back in i_atom should be a[i0] and jj should be a[jj].
@@ -658,7 +586,8 @@ put_in_list_at(gmx_bool              bHaveVdW[],
     
     /* Get the i charge group info */
     igid   = GET_CGINFO_GID(cginfo[icg]);
-    iwater = GET_CGINFO_SOLOPT(cginfo[icg]);
+
+    iwater = (solvent_opt!=esolNO) ? GET_CGINFO_SOLOPT(cginfo[icg]) : esolNO;
     
     bFreeEnergy = FALSE;
     if (md->nPerturbed) 
@@ -1152,9 +1081,10 @@ put_in_list_qmmm(gmx_bool              bHaveVdW[],
                  t_excl            bExcl[],
                  int               shift,
                  t_forcerec *      fr,
-                 gmx_bool              bLR,
-                 gmx_bool              bDoVdW,
-                 gmx_bool              bDoCoul)
+                 gmx_bool          bLR,
+                 gmx_bool          bDoVdW,
+                 gmx_bool          bDoCoul,
+                 int               solvent_opt)
 {
     t_nblist *   coul;
     int          i,j,jcg,igid,gid;
@@ -1214,9 +1144,10 @@ put_in_list_cg(gmx_bool              bHaveVdW[],
                t_excl            bExcl[],
                int               shift,
                t_forcerec *      fr,
-               gmx_bool              bLR,
-               gmx_bool              bDoVdW,
-               gmx_bool              bDoCoul)
+               gmx_bool          bLR,
+               gmx_bool          bDoVdW,
+               gmx_bool          bDoCoul,
+               int               solvent_opt)
 {
     int          cginfo;
     int          igid,gid,nbl_ind;
@@ -1444,7 +1375,7 @@ static void add_simple(t_ns_buf *nsbuf,int nrj,atom_id cg_j,
     if (nsbuf->nj + nrj > MAX_CG)
     {
         put_in_list(bHaveVdW,ngid,md,icg,jgid,nsbuf->ncg,nsbuf->jcg,
-                    cgs->index,bexcl,shift,fr,FALSE,TRUE,TRUE);
+                    cgs->index,bexcl,shift,fr,FALSE,TRUE,TRUE,fr->solvent_opt);
         /* Reset buffer contents */
         nsbuf->ncg = nsbuf->nj = 0;
     }
@@ -1518,7 +1449,7 @@ static void ns_inner_rect(rvec x[],int icg,int *i_egp_flags,
                 }
             }
         }
-    } 
+    }
     else
     {
         for(j=0; (j<njcg); j++)
@@ -1616,7 +1547,7 @@ static int ns_simple_core(t_forcerec *fr,
                 if (nsbuf->ncg > 0)
                 {
                     put_in_list(bHaveVdW,ngid,md,icg,nn,nsbuf->ncg,nsbuf->jcg,
-                                cgs->index,bexcl,k,fr,FALSE,TRUE,TRUE);
+                                cgs->index,bexcl,k,fr,FALSE,TRUE,TRUE,fr->solvent_opt);
                     nsbuf->ncg=nsbuf->nj=0;
                 }
             }
@@ -1624,7 +1555,7 @@ static int ns_simple_core(t_forcerec *fr,
         /* setexcl(nri,i_atoms,excl,FALSE,bexcl); */
         setexcl(cgs->index[icg],cgs->index[icg+1],excl,FALSE,bexcl);
     }
-    close_neighbor_list(fr,FALSE,-1,-1,FALSE);
+    close_neighbor_lists(fr,FALSE);
     
     return nsearch;
 }
@@ -1785,56 +1716,13 @@ static inline void get_dx_dd(int Nx,real gridx,real rc2,int xgi,real x,
  *
  ****************************************************/
 
-static void do_longrange(t_commrec *cr,gmx_localtop_t *top,t_forcerec *fr,
-                         int ngid,t_mdatoms *md,int icg,
-                         int jgid,int nlr,
-                         atom_id lr[],t_excl bexcl[],int shift,
-                         rvec x[],rvec box_size,t_nrnb *nrnb,
-                         real *lambda,real *dvdlambda,
-                         gmx_grppairener_t *grppener,
-                         gmx_bool bDoVdW,gmx_bool bDoCoul,
-                         gmx_bool bEvaluateNow,put_in_list_t *put_in_list,
-                         gmx_bool bHaveVdW[],
-                         gmx_bool bDoForces,rvec *f)
-{
-    int n,i;
-    t_nblist *nl;
-    
-    for(n=0; n<fr->nnblists; n++)
-    {
-        for(i=0; (i<eNL_NR); i++)
-        {
-            nl = &fr->nblists[n].nlist_lr[i];
-            if ((nl->nri > nl->maxnri-32) || bEvaluateNow)
-            {
-                close_neighbor_list(fr,TRUE,n,i,FALSE);
-                /* Evaluate the energies and forces */
-                do_nonbonded(cr,fr,x,f,md,NULL,
-                             grppener->ener[fr->bBHAM ? egBHAMLR : egLJLR],
-                             grppener->ener[egCOULLR],
-                                                        grppener->ener[egGB],box_size,
-                             nrnb,lambda,dvdlambda,n,i,
-                             GMX_DONB_LR | GMX_DONB_FORCES);
-                
-                reset_neighbor_list(fr,TRUE,n,i);
-            }
-        }
-    }
-    
-    if (!bEvaluateNow)
-    {  
-        /* Put the long range particles in a list */
-        /* do_longrange is never called for QMMM  */
-        put_in_list(bHaveVdW,ngid,md,icg,jgid,nlr,lr,top->cgs.index,
-                    bexcl,shift,fr,TRUE,bDoVdW,bDoCoul);
-    }
-}
 
 static void get_cutoff2(t_forcerec *fr,gmx_bool bDoLongRange,
                         real *rvdw2,real *rcoul2,
                         real *rs2,real *rm2,real *rl2)
 {
     *rs2 = sqr(fr->rlist);
+
     if (bDoLongRange && fr->bTwinRange)
     {
         /* The VdW and elec. LR cut-off's could be different,
@@ -1883,26 +1771,16 @@ static void init_nsgrid_lists(t_forcerec *fr,int ngid,gmx_ns_t *ns)
     snew(ns->nlr_ljc,ngid);
     snew(ns->nlr_one,ngid);
     
-    if (rm2 > rs2)
-    {
-            /* Long range VdW and Coul buffers */
-        snew(ns->nl_lr_ljc,ngid);
-    }
-    if (rl2 > rm2)
-    {
-        /* Long range VdW or Coul only buffers */
-        snew(ns->nl_lr_one,ngid);
-    }
+    /* Always allocate both list types, since rcoulomb might now change with PME load balancing */
+    /* Long range VdW and Coul buffers */
+    snew(ns->nl_lr_ljc,ngid);
+    /* Long range VdW or Coul only buffers */
+    snew(ns->nl_lr_one,ngid);
+
     for(j=0; (j<ngid); j++) {
         snew(ns->nl_sr[j],MAX_CG);
-        if (rm2 > rs2)
-        {
-            snew(ns->nl_lr_ljc[j],MAX_CG);
-        }
-        if (rl2 > rm2)
-        {
-            snew(ns->nl_lr_one[j],MAX_CG);
-        }
+        snew(ns->nl_lr_ljc[j],MAX_CG);
+        snew(ns->nl_lr_one[j],MAX_CG);
     }
     if (debug)
     {
@@ -1922,8 +1800,7 @@ static int nsgrid_core(FILE *log,t_commrec *cr,t_forcerec *fr,
                        gmx_grppairener_t *grppener,
                        put_in_list_t *put_in_list,
                        gmx_bool bHaveVdW[],
-                       gmx_bool bDoLongRange,gmx_bool bDoForces,rvec *f,
-                       gmx_bool bMakeQMMMnblist)
+                       gmx_bool bDoLongRange,gmx_bool bMakeQMMMnblist)
 {
     gmx_ns_t *ns;
     atom_id **nl_lr_ljc,**nl_lr_one,**nl_sr;
@@ -2273,10 +2150,11 @@ static int nsgrid_core(FILE *log,t_commrec *cr,t_forcerec *fr,
                                                         {
                                                             if (nsr[jgid] >= MAX_CG)
                                                             {
+                                                                /* Add to short-range list */
                                                                 put_in_list(bHaveVdW,ngid,md,icg,jgid,
                                                                             nsr[jgid],nl_sr[jgid],
                                                                             cgs->index,/* cgsatoms, */ bexcl,
-                                                                            shift,fr,FALSE,TRUE,TRUE);
+                                                                            shift,fr,FALSE,TRUE,TRUE,fr->solvent_opt);
                                                                 nsr[jgid]=0;
                                                             }
                                                             nl_sr[jgid][nsr[jgid]++]=jjcg;
@@ -2285,33 +2163,22 @@ static int nsgrid_core(FILE *log,t_commrec *cr,t_forcerec *fr,
                                                         {
                                                             if (nlr_ljc[jgid] >= MAX_CG)
                                                             {
-                                                                do_longrange(cr,top,fr,ngid,md,icg,jgid,
-                                                                             nlr_ljc[jgid],
-                                                                             nl_lr_ljc[jgid],bexcl,shift,x,
-                                                                             box_size,nrnb,
-                                                                             lambda,dvdlambda,
-                                                                             grppener,
-                                                                             TRUE,TRUE,FALSE,
-                                                                             put_in_list,
-                                                                             bHaveVdW,
-                                                                             bDoForces,f);
+                                                                /* Add to LJ+coulomb long-range list */
+                                                                put_in_list(bHaveVdW,ngid,md,icg,jgid,
+                                                                            nlr_ljc[jgid],nl_lr_ljc[jgid],top->cgs.index,
+                                                                            bexcl,shift,fr,TRUE,TRUE,TRUE,fr->solvent_opt);
                                                                 nlr_ljc[jgid]=0;
                                                             }
                                                             nl_lr_ljc[jgid][nlr_ljc[jgid]++]=jjcg;
                                                         }
                                                         else
                                                         {
-                                                            if (nlr_one[jgid] >= MAX_CG) {
-                                                                do_longrange(cr,top,fr,ngid,md,icg,jgid,
-                                                                             nlr_one[jgid],
-                                                                             nl_lr_one[jgid],bexcl,shift,x,
-                                                                             box_size,nrnb,
-                                                                             lambda,dvdlambda,
-                                                                             grppener,
-                                                                             rvdw_lt_rcoul,rcoul_lt_rvdw,FALSE,
-                                                                             put_in_list,
-                                                                             bHaveVdW,
-                                                                             bDoForces,f);
+                                                            if (nlr_one[jgid] >= MAX_CG)
+                                                            {
+                                                                /* Add to long-range list with only coul, or only LJ */
+                                                                put_in_list(bHaveVdW,ngid,md,icg,jgid,
+                                                                            nlr_one[jgid],nl_lr_one[jgid],top->cgs.index,
+                                                                            bexcl,shift,fr,TRUE,rvdw_lt_rcoul,rcoul_lt_rvdw,fr->solvent_opt);
                                                                 nlr_one[jgid]=0;
                                                             }
                                                             nl_lr_one[jgid][nlr_one[jgid]++]=jjcg;
@@ -2333,24 +2200,21 @@ static int nsgrid_core(FILE *log,t_commrec *cr,t_forcerec *fr,
                         {
                             put_in_list(bHaveVdW,ngid,md,icg,nn,nsr[nn],nl_sr[nn],
                                         cgs->index, /* cgsatoms, */ bexcl,
-                                        shift,fr,FALSE,TRUE,TRUE);
+                                        shift,fr,FALSE,TRUE,TRUE,fr->solvent_opt);
                         }
                         
                         if (nlr_ljc[nn] > 0)
                         {
-                            do_longrange(cr,top,fr,ngid,md,icg,nn,nlr_ljc[nn],
-                                         nl_lr_ljc[nn],bexcl,shift,x,box_size,nrnb,
-                                         lambda,dvdlambda,grppener,TRUE,TRUE,FALSE,
-                                         put_in_list,bHaveVdW,bDoForces,f);
+                            put_in_list(bHaveVdW,ngid,md,icg,nn,nlr_ljc[nn],
+                                        nl_lr_ljc[nn],top->cgs.index,
+                                        bexcl,shift,fr,TRUE,TRUE,TRUE,fr->solvent_opt);
                         }
                         
                         if (nlr_one[nn] > 0)
                         {
-                            do_longrange(cr,top,fr,ngid,md,icg,nn,nlr_one[nn],
-                                         nl_lr_one[nn],bexcl,shift,x,box_size,nrnb,
-                                         lambda,dvdlambda,grppener,
-                                         rvdw_lt_rcoul,rcoul_lt_rvdw,FALSE,
-                                         put_in_list,bHaveVdW,bDoForces,f);
+                            put_in_list(bHaveVdW,ngid,md,icg,nn,nlr_one[nn],
+                                        nl_lr_one[nn],top->cgs.index,
+                                        bexcl,shift,fr,TRUE,rvdw_lt_rcoul,rcoul_lt_rvdw,fr->solvent_opt);
                         }
                     }
                 }
@@ -2359,28 +2223,14 @@ static int nsgrid_core(FILE *log,t_commrec *cr,t_forcerec *fr,
         /* setexcl(nri,i_atoms,&top->atoms.excl,FALSE,bexcl); */
         setexcl(cgs->index[icg],cgs->index[icg+1],&top->excls,FALSE,bexcl);
     }
-    /* Perform any left over force calculations */
-    for (nn=0; (nn<ngid); nn++)
-    {
-        if (rm2 > rs2)
-        {
-            do_longrange(cr,top,fr,0,md,icg,nn,nlr_ljc[nn],
-                         nl_lr_ljc[nn],bexcl,shift,x,box_size,nrnb,
-                         lambda,dvdlambda,grppener,
-                         TRUE,TRUE,TRUE,put_in_list,bHaveVdW,bDoForces,f);
-        }
-        if (rl2 > rm2) {
-            do_longrange(cr,top,fr,0,md,icg,nn,nlr_one[nn],
-                         nl_lr_one[nn],bexcl,shift,x,box_size,nrnb,
-                         lambda,dvdlambda,grppener,
-                         rvdw_lt_rcoul,rcoul_lt_rvdw,
-                         TRUE,put_in_list,bHaveVdW,bDoForces,f);
-        }
-    }
+    /* No need to perform any left-over force calculations anymore (as we used to do here)
+     * since we now save the proper long-range lists for later evaluation.
+     */
+
     debug_gmx();
-    
-    /* Close off short range neighbourlists */
-    close_neighbor_list(fr,FALSE,-1,-1,bMakeQMMMnblist);
+     
+    /* Close neighbourlists */
+    close_neighbor_lists(fr,bMakeQMMMnblist);
     
     return nns;
 }
@@ -2522,8 +2372,7 @@ int search_neighbours(FILE *log,t_forcerec *fr,
                       real *lambda,real *dvdlambda,
                       gmx_grppairener_t *grppener,
                       gmx_bool bFillGrid,
-                      gmx_bool bDoLongRange,
-                      gmx_bool bDoForces,rvec *f)
+                      gmx_bool bDoLongRangeNS)
 {
     t_block  *cgs=&(top->cgs);
     rvec     box_size,grid_x0,grid_x1;
@@ -2538,7 +2387,7 @@ int search_neighbours(FILE *log,t_forcerec *fr,
     t_grid   *grid;
     gmx_domdec_zones_t *dd_zones;
     put_in_list_t *put_in_list;
-       
+
     ns = &fr->ns;
 
     /* Set some local variables */
@@ -2571,7 +2420,7 @@ int search_neighbours(FILE *log,t_forcerec *fr,
     debug_gmx();
     
     /* Reset the neighbourlists */
-    reset_neighbor_list(fr,FALSE,-1,-1);
+    reset_neighbor_lists(fr,TRUE,TRUE);
     
     if (bGrid && bFillGrid)
     {
@@ -2657,8 +2506,7 @@ int search_neighbours(FILE *log,t_forcerec *fr,
                               grid,x,ns->bexcl,ns->bExcludeAlleg,
                               nrnb,md,lambda,dvdlambda,grppener,
                               put_in_list,ns->bHaveVdW,
-                              bDoLongRange,bDoForces,f,
-                              FALSE);
+                              bDoLongRangeNS,FALSE);
         
         /* neighbour searching withouth QMMM! QM atoms have zero charge in
          * the classical calculation. The charge-charge interaction
@@ -2675,8 +2523,7 @@ int search_neighbours(FILE *log,t_forcerec *fr,
                                    grid,x,ns->bexcl,ns->bExcludeAlleg,
                                    nrnb,md,lambda,dvdlambda,grppener,
                                    put_in_list_qmmm,ns->bHaveVdW,
-                                   bDoLongRange,bDoForces,f,
-                                   TRUE);
+                                   bDoLongRangeNS,TRUE);
         }
     }
     else 
@@ -2686,7 +2533,7 @@ int search_neighbours(FILE *log,t_forcerec *fr,
                                  ngid,ns->ns_buf,put_in_list,ns->bHaveVdW);
     }
     debug_gmx();
-    
+
 #ifdef DEBUG
     pr_nsblock(log);
 #endif