Fixes potential bug in neighborsearching.
[alexxy/gromacs.git] / src / gromacs / mdlib / ns.c
index 4527ebec7571ba52524f61ba0ff9ff49693a9822..430e26d304780b5b3f828efdad076163127cff1d 100644 (file)
@@ -1,37 +1,38 @@
-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
- *
- *
- *                This source code is part of
- *
- *                 G   R   O   M   A   C   S
- *
- *          GROningen MAchine for Chemical Simulations
+/*
+ * This file is part of the GROMACS molecular simulation package.
  *
- *                        VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
+ * Copyright (c) 2001-2004, The GROMACS development team.
+ * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
  * 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
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
  *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  *
- * For more info, check our website at http://www.gromacs.org
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
  *
- * And Hey:
- * GROwing Monsters And Cloning Shrimps
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
  */
 #ifdef HAVE_CONFIG_H
 #include <config.h>
 #include <math.h>
 #include <string.h>
 #include "sysstuff.h"
-#include "smalloc.h"
+#include "gromacs/utility/smalloc.h"
 #include "macros.h"
-#include "maths.h"
+#include "gromacs/math/utilities.h"
 #include "vec.h"
+#include "types/commrec.h"
 #include "network.h"
 #include "nsgrid.h"
 #include "force.h"
@@ -103,7 +105,7 @@ round_up_to_simd_width(int length, int simd_width)
  *
  ************************************************/
 
-static void reallocate_nblist(t_nblist *nl)
+void reallocate_nblist(t_nblist *nl)
 {
     if (gmx_debug_at)
     {
@@ -125,7 +127,8 @@ static void init_nblist(FILE *log, t_nblist *nl_sr, t_nblist *nl_lr,
                         int maxsr, int maxlr,
                         int ivdw, int ivdwmod,
                         int ielec, int ielecmod,
-                        int igeometry, int type)
+                        int igeometry, int type,
+                        gmx_bool bElecAndVdwSwitchDiffers)
 {
     t_nblist *nl;
     int       homenr;
@@ -158,7 +161,7 @@ static void init_nblist(FILE *log, t_nblist *nl_sr, t_nblist *nl_lr,
         }
 
         /* This will also set the simd_padding_width field */
-        gmx_nonbonded_set_kernel_pointers( (i == 0) ? log : NULL, nl);
+        gmx_nonbonded_set_kernel_pointers( (i == 0) ? log : NULL, nl, bElecAndVdwSwitchDiffers);
 
         /* maxnri is influenced by the number of shifts (maximum is 8)
          * and the number of energy groups.
@@ -168,13 +171,14 @@ static void init_nblist(FILE *log, t_nblist *nl_sr, t_nblist *nl_lr,
          */
         nl->maxnri      = homenr*4;
         nl->maxnrj      = 0;
-        nl->maxlen      = 0;
         nl->nri         = -1;
         nl->nrj         = 0;
         nl->iinr        = NULL;
         nl->gid         = NULL;
         nl->shift       = NULL;
         nl->jindex      = NULL;
+        nl->jjnr        = NULL;
+        nl->excl_fep    = NULL;
         reallocate_nblist(nl);
         nl->jindex[0] = 0;
 
@@ -194,10 +198,11 @@ void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr)
      * cache trashing.
      */
     int        maxsr, maxsr_wat, maxlr, maxlr_wat;
-    int        ielec, ielecf, ivdw, ielecmod, ielecmodf, ivdwmod, type;
+    int        ielec, ivdw, ielecmod, ivdwmod, type;
     int        solvent;
     int        igeometry_def, igeometry_w, igeometry_ww;
     int        i;
+    gmx_bool   bElecAndVdwSwitchDiffers;
     t_nblists *nbl;
 
     /* maxsr     = homenr-fr->nWatMol*3; */
@@ -224,11 +229,12 @@ void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr)
     }
 
     /* 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;
-    type     = GMX_NBLIST_INTERACTION_STANDARD;
+    ielec                    = fr->nbkernel_elec_interaction;
+    ivdw                     = fr->nbkernel_vdw_interaction;
+    ielecmod                 = fr->nbkernel_elec_modifier;
+    ivdwmod                  = fr->nbkernel_vdw_modifier;
+    type                     = GMX_NBLIST_INTERACTION_STANDARD;
+    bElecAndVdwSwitchDiffers = ( (fr->rcoulomb_switch != fr->rvdw_switch) || (fr->rcoulomb != fr->rvdw));
 
     fr->ns.bCGlist = (getenv("GMX_NBLISTCG") != 0);
     if (!fr->ns.bCGlist)
@@ -264,19 +270,19 @@ void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr)
             type = GMX_NBLIST_INTERACTION_ADRESS;
         }
         init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ], &nbl->nlist_lr[eNL_VDWQQ],
-                    maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, igeometry_def, type);
+                    maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, igeometry_def, type, bElecAndVdwSwitchDiffers);
         init_nblist(log, &nbl->nlist_sr[eNL_VDW], &nbl->nlist_lr[eNL_VDW],
-                    maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, igeometry_def, type);
+                    maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, igeometry_def, type, bElecAndVdwSwitchDiffers);
         init_nblist(log, &nbl->nlist_sr[eNL_QQ], &nbl->nlist_lr[eNL_QQ],
-                    maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_def, type);
+                    maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_def, type, bElecAndVdwSwitchDiffers);
         init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_WATER], &nbl->nlist_lr[eNL_VDWQQ_WATER],
-                    maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_w, type);
+                    maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_w, type, bElecAndVdwSwitchDiffers);
         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, igeometry_w, type);
+                    maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_w, type, bElecAndVdwSwitchDiffers);
         init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_WATERWATER], &nbl->nlist_lr[eNL_VDWQQ_WATERWATER],
-                    maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_ww, type);
+                    maxsr_wat, maxlr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_ww, type, bElecAndVdwSwitchDiffers);
         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, igeometry_ww, type);
+                    maxsr_wat, maxlr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_ww, type, bElecAndVdwSwitchDiffers);
 
         /* Did we get the solvent loops so we can use optimized water kernels? */
         if (nbl->nlist_sr[eNL_VDWQQ_WATER].kernelptr_vf == NULL
@@ -288,35 +294,27 @@ void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr)
             )
         {
             fr->solvent_opt = esolNO;
-            fprintf(log, "Note: The available nonbonded kernels do not support water optimization - disabling.\n");
+            if (log != NULL)
+            {
+                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 */
-            {
-                ielecf    = GMX_NBKERNEL_ELEC_EWALD;
-                ielecmodf = eintmodNONE;
-            }
-            else
-            {
-                ielecf    = ielec;
-                ielecmodf = ielecmod;
-            }
-
             init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_FREE], &nbl->nlist_lr[eNL_VDWQQ_FREE],
-                        maxsr, maxlr, ivdw, ivdwmod, ielecf, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
+                        maxsr, maxlr, ivdw, ivdwmod, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
             init_nblist(log, &nbl->nlist_sr[eNL_VDW_FREE], &nbl->nlist_lr[eNL_VDW_FREE],
-                        maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
+                        maxsr, maxlr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
             init_nblist(log, &nbl->nlist_sr[eNL_QQ_FREE], &nbl->nlist_lr[eNL_QQ_FREE],
-                        maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielecf, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
+                        maxsr, maxlr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY, bElecAndVdwSwitchDiffers);
         }
     }
     /* QMMM MM list */
     if (fr->bQMMM && fr->qr->QMMMscheme != eQMMMschemeoniom)
     {
         init_nblist(log, &fr->QMMMlist, NULL,
-                    maxsr, maxlr, 0, 0, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_STANDARD);
+                    maxsr, maxlr, 0, 0, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_STANDARD, bElecAndVdwSwitchDiffers);
     }
 
     if (log != NULL)
@@ -331,7 +329,6 @@ static void reset_nblist(t_nblist *nl)
 {
     nl->nri       = -1;
     nl->nrj       = 0;
-    nl->maxlen    = 0;
     if (nl->jindex)
     {
         nl->jindex[0] = 0;
@@ -367,7 +364,7 @@ static void reset_neighbor_lists(t_forcerec *fr, gmx_bool bResetSR, gmx_bool bRe
 
 
 
-static inline void new_i_nblist(t_nblist *nlist, atom_id i_atom, int shift, int gid)
+static gmx_inline void new_i_nblist(t_nblist *nlist, atom_id i_atom, int shift, int gid)
 {
     int    i, k, nri, nshift;
 
@@ -401,9 +398,20 @@ static inline void new_i_nblist(t_nblist *nlist, atom_id i_atom, int shift, int
         nlist->gid[nri]      = gid;
         nlist->shift[nri]    = shift;
     }
+    else
+    {
+        /* Adding to previous list. First remove possible previous padding */
+        if (nlist->simd_padding_width > 1)
+        {
+            while (nlist->nrj > 0 && nlist->jjnr[nlist->nrj-1] < 0)
+            {
+                nlist->nrj--;
+            }
+        }
+    }
 }
 
-static inline void close_i_nblist(t_nblist *nlist)
+static gmx_inline void close_i_nblist(t_nblist *nlist)
 {
     int nri = nlist->nri;
     int len;
@@ -422,18 +430,18 @@ static inline void close_i_nblist(t_nblist *nlist)
         nlist->jindex[nri+1] = nlist->nrj;
 
         len = nlist->nrj -  nlist->jindex[nri];
-
-        /* nlist length for water i molecules is treated statically
-         * in the innerloops
+        /* If there are no j-particles we have to reduce the
+         * number of i-particles again, to prevent errors in the
+         * kernel functions.
          */
-        if (len > nlist->maxlen)
+        if ((len == 0) && (nlist->nri > 0))
         {
-            nlist->maxlen = len;
+            nlist->nri--;
         }
     }
 }
 
-static inline void close_nblist(t_nblist *nlist)
+static gmx_inline void close_nblist(t_nblist *nlist)
 {
     /* Only close this nblist when it has been initialized.
      * Avoid the creation of i-lists with no j-particles.
@@ -455,7 +463,7 @@ static inline void close_nblist(t_nblist *nlist)
     }
 }
 
-static inline void close_neighbor_lists(t_forcerec *fr, gmx_bool bMakeQMMMnblist)
+static gmx_inline void close_neighbor_lists(t_forcerec *fr, gmx_bool bMakeQMMMnblist)
 {
     int n, i;
 
@@ -475,7 +483,7 @@ static inline void close_neighbor_lists(t_forcerec *fr, gmx_bool bMakeQMMMnblist
 }
 
 
-static inline void add_j_to_nblist(t_nblist *nlist, atom_id j_atom, gmx_bool bLR)
+static gmx_inline void add_j_to_nblist(t_nblist *nlist, atom_id j_atom, gmx_bool bLR)
 {
     int nrj = nlist->nrj;
 
@@ -496,10 +504,10 @@ static inline void add_j_to_nblist(t_nblist *nlist, atom_id j_atom, gmx_bool bLR
     nlist->nrj++;
 }
 
-static inline void add_j_to_nblist_cg(t_nblist *nlist,
-                                      atom_id j_start, int j_end,
-                                      t_excl *bexcl, gmx_bool i_is_j,
-                                      gmx_bool bLR)
+static gmx_inline void add_j_to_nblist_cg(t_nblist *nlist,
+                                          atom_id j_start, int j_end,
+                                          t_excl *bexcl, gmx_bool i_is_j,
+                                          gmx_bool bLR)
 {
     int nrj = nlist->nrj;
     int j;
@@ -1873,8 +1881,8 @@ static int ns_simple_core(t_forcerec *fr,
  *
  ************************************************/
 
-static inline void get_dx(int Nx, real gridx, real rc2, int xgi, real x,
-                          int *dx0, int *dx1, real *dcx2)
+static gmx_inline void get_dx(int Nx, real gridx, real rc2, int xgi, real x,
+                              int *dx0, int *dx1, real *dcx2)
 {
     real dcx, tmp;
     int  xgi0, xgi1, i;
@@ -1926,9 +1934,9 @@ static inline void get_dx(int Nx, real gridx, real rc2, int xgi, real x,
     }
 }
 
-static inline void get_dx_dd(int Nx, real gridx, real rc2, int xgi, real x,
-                             int ncpddc, int shift_min, int shift_max,
-                             int *g0, int *g1, real *dcx2)
+static gmx_inline void get_dx_dd(int Nx, real gridx, real rc2, int xgi, real x,
+                                 int ncpddc, int shift_min, int shift_max,
+                                 int *g0, int *g1, real *dcx2)
 {
     real dcx, tmp;
     int  g_min, g_max, shift_home;
@@ -2034,26 +2042,36 @@ static void get_cutoff2(t_forcerec *fr, gmx_bool bDoLongRange,
 
     if (bDoLongRange && fr->bTwinRange)
     {
-        /* The VdW and elec. LR cut-off's could be different,
+        /* With plain cut-off or RF we need to make the list exactly
+         * up to the cut-off and the cut-off's can be different,
          * so we can not simply set them to rlistlong.
+         * To keep this code compatible with (exotic) old cases,
+         * we also create lists up to rvdw/rcoulomb for PME and Ewald.
+         * The interaction check should correspond to:
+         * !ir_vdw/coulomb_might_be_zero_at_cutoff from inputrec.c.
          */
-        if (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(fr->vdwtype) &&
-            fr->rvdw > fr->rlist)
+        if (((fr->vdwtype == evdwCUT || fr->vdwtype == evdwPME) &&
+             fr->vdw_modifier == eintmodNONE) ||
+            fr->rvdw <= fr->rlist)
         {
-            *rvdw2  = sqr(fr->rlistlong);
+            *rvdw2  = sqr(fr->rvdw);
         }
         else
         {
-            *rvdw2  = sqr(fr->rvdw);
+            *rvdw2  = sqr(fr->rlistlong);
         }
-        if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(fr->eeltype) &&
-            fr->rcoulomb > fr->rlist)
+        if (((fr->eeltype == eelCUT ||
+              (EEL_RF(fr->eeltype) && fr->eeltype != eelRF_ZERO) ||
+              fr->eeltype == eelPME ||
+              fr->eeltype == eelEWALD) &&
+             fr->coulomb_modifier == eintmodNONE) ||
+            fr->rcoulomb <= fr->rlist)
         {
-            *rcoul2 = sqr(fr->rlistlong);
+            *rcoul2 = sqr(fr->rcoulomb);
         }
         else
         {
-            *rcoul2 = sqr(fr->rcoulomb);
+            *rcoul2 = sqr(fr->rlistlong);
         }
     }
     else
@@ -2655,7 +2673,6 @@ void init_ns(FILE *fplog, const t_commrec *cr,
     ns->bexcl     = NULL;
     if (!DOMAINDECOMP(cr))
     {
-        /* This could be reduced with particle decomposition */
         ns_realloc_natoms(ns, mtop->natoms);
     }
 
@@ -2759,14 +2776,8 @@ int search_neighbours(FILE *log, t_forcerec *fr,
         }
         debug_gmx();
 
-        /* Don't know why this all is... (DvdS 3/99) */
-#ifndef SEGV
         start = 0;
         end   = cgs->nr;
-#else
-        start = fr->cg0;
-        end   = (cgs->nr+1)/2;
-#endif
 
         if (DOMAINDECOMP(cr))
         {
@@ -2781,12 +2792,6 @@ int search_neighbours(FILE *log, t_forcerec *fr,
             grid->icg0 = fr->cg0;
             grid->icg1 = fr->hcg;
             debug_gmx();
-
-            if (PARTDECOMP(cr))
-            {
-                mv_grid(cr, grid);
-            }
-            debug_gmx();
         }
 
         calc_elemnr(grid, start, end, cgs->nr);