-/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+/*
+ * This file is part of the GROMACS molecular simulation package.
*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * 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.
+ *
+ * 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.
+ *
+ * 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.
+ *
* 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.
- *
- * For more info, check our website at http://www.gromacs.org
- *
- * And Hey:
- * GROwing Monsters And Cloning Shrimps
+ * the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
+
+#include "gromacs/legacyheaders/ns.h"
#include <math.h>
+#include <stdlib.h>
#include <string.h>
-#include "sysstuff.h"
-#include "smalloc.h"
-#include "macros.h"
-#include "maths.h"
-#include "vec.h"
-#include "network.h"
-#include "nsgrid.h"
-#include "force.h"
-#include "nonbonded.h"
-#include "ns.h"
-#include "pbc.h"
-#include "names.h"
-#include "gmx_fatal.h"
-#include "nrnb.h"
-#include "txtdump.h"
-#include "mtop_util.h"
-
-#include "domdec.h"
-#include "adress.h"
+#include "gromacs/legacyheaders/domdec.h"
+#include "gromacs/legacyheaders/force.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/network.h"
+#include "gromacs/legacyheaders/nonbonded.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/nsgrid.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/smalloc.h"
-/*
+#include "adress.h"
+
+/*
* E X C L U S I O N H A N D L I N G
*/
#ifdef DEBUG
-static void SETEXCL_(t_excl e[],atom_id i,atom_id j)
-{ e[j] = e[j] | (1<<i); }
-static void RMEXCL_(t_excl e[],atom_id i,atom_id j)
-{ e[j]=e[j] & ~(1<<i); }
-static gmx_bool ISEXCL_(t_excl e[],atom_id i,atom_id j)
-{ return (gmx_bool)(e[j] & (1<<i)); }
-static gmx_bool NOTEXCL_(t_excl e[],atom_id i,atom_id j)
-{ return !(ISEXCL(e,i,j)); }
+static void SETEXCL_(t_excl e[], atom_id i, atom_id j)
+{
+ e[j] = e[j] | (1<<i);
+}
+static void RMEXCL_(t_excl e[], atom_id i, atom_id j)
+{
+ e[j] = e[j] & ~(1<<i);
+}
+static gmx_bool ISEXCL_(t_excl e[], atom_id i, atom_id j)
+{
+ return (gmx_bool)(e[j] & (1<<i));
+}
+static gmx_bool NOTEXCL_(t_excl e[], atom_id i, atom_id j)
+{
+ return !(ISEXCL(e, i, j));
+}
#else
-#define SETEXCL(e,i,j) (e)[((atom_id) (j))] |= (1<<((atom_id) (i)))
-#define RMEXCL(e,i,j) (e)[((atom_id) (j))] &= (~(1<<((atom_id) (i))))
-#define ISEXCL(e,i,j) (gmx_bool) ((e)[((atom_id) (j))] & (1<<((atom_id) (i))))
-#define NOTEXCL(e,i,j) !(ISEXCL(e,i,j))
+#define SETEXCL(e, i, j) (e)[((atom_id) (j))] |= (1<<((atom_id) (i)))
+#define RMEXCL(e, i, j) (e)[((atom_id) (j))] &= (~(1<<((atom_id) (i))))
+#define ISEXCL(e, i, j) (gmx_bool) ((e)[((atom_id) (j))] & (1<<((atom_id) (i))))
+#define NOTEXCL(e, i, j) !(ISEXCL(e, i, j))
#endif
+static int
+round_up_to_simd_width(int length, int simd_width)
+{
+ int offset, newlength;
+
+ offset = (simd_width > 0) ? length % simd_width : 0;
+
+ return (offset == 0) ? length : length-offset+simd_width;
+}
/************************************************
*
* U T I L I T I E S F O R N S
*
************************************************/
-static void reallocate_nblist(t_nblist *nl)
+void reallocate_nblist(t_nblist *nl)
{
if (gmx_debug_at)
{
- 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);
+ fprintf(debug, "reallocating neigborlist (ielec=%d, ivdw=%d, igeometry=%d, type=%d), maxnri=%d\n",
+ nl->ielec, nl->ivdw, nl->igeometry, nl->type, nl->maxnri);
}
srenew(nl->iinr, nl->maxnri);
if (nl->igeometry == GMX_NBLIST_GEOMETRY_CG_CG)
{
- srenew(nl->iinr_end,nl->maxnri);
+ srenew(nl->iinr_end, nl->maxnri);
}
srenew(nl->gid, nl->maxnri);
srenew(nl->shift, nl->maxnri);
}
-static void init_nblist(FILE *log, t_nblist *nl_sr,t_nblist *nl_lr,
- int maxsr,int maxlr,
+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,
- gmx_bool bfree, int igeometry)
+ int igeometry, int type,
+ gmx_bool bElecAndVdwSwitchDiffers)
{
t_nblist *nl;
- int homenr;
- int i,nn;
-
- for(i=0; (i<2); i++)
+ int homenr;
+ int i, nn;
+
+ for (i = 0; (i < 2); i++)
{
nl = (i == 0) ? nl_sr : nl_lr;
homenr = (i == 0) ? maxsr : maxlr;
nl->ivdwmod = ivdwmod;
nl->ielec = ielec;
nl->ielecmod = ielecmod;
- nl->free_energy = bfree;
+ nl->type = type;
nl->igeometry = igeometry;
- if (bfree)
+ if (nl->type == GMX_NBLIST_INTERACTION_FREE_ENERGY)
{
nl->igeometry = GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE;
}
-
- gmx_nonbonded_set_kernel_pointers( (i==0) ? log : NULL,nl);
-
+
+ /* This will also set the simd_padding_width field */
+ 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.
* If it is not enough, nl memory will be reallocated during the run.
*/
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;
- if(debug)
+ 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);
+ fprintf(debug, "Initiating neighbourlist (ielec=%d, ivdw=%d, type=%d) for %s interactions,\nwith %d SR, %d LR atoms.\n",
+ nl->ielec, nl->ivdw, nl->type, gmx_nblist_geometry_names[nl->igeometry], maxsr, maxlr);
}
}
}
-void init_neighbor_list(FILE *log,t_forcerec *fr,int homenr)
+void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr)
{
- /* Make maxlr tunable! (does not seem to be a big difference though)
- * This parameter determines the number of i particles in a long range
- * neighbourlist. Too few means many function calls, too many means
- * cache trashing.
- */
- int maxsr,maxsr_wat,maxlr,maxlr_wat;
- int ielec,ielecf,ivdw,ielecmod,ielecmodf,ivdwmod;
- int solvent;
- int igeometry_def,igeometry_w,igeometry_ww;
- int i;
- t_nblists *nbl;
-
- /* maxsr = homenr-fr->nWatMol*3; */
- maxsr = homenr;
-
- if (maxsr < 0)
- {
- gmx_fatal(FARGS,"%s, %d: Negative number of short range atoms.\n"
- "Call your Gromacs dealer for assistance.",__FILE__,__LINE__);
- }
- /* This is just for initial allocation, so we do not reallocate
- * all the nlist arrays many times in a row.
- * The numbers seem very accurate, but they are uncritical.
- */
- maxsr_wat = min(fr->nWatMol,(homenr+2)/3);
- if (fr->bTwinRange)
- {
- maxlr = 50;
- maxlr_wat = min(maxsr_wat,maxlr);
- }
- else
- {
- maxlr = maxlr_wat = 0;
- }
-
- /* 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)
- {
- igeometry_def = GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE;
- }
- else
- {
- igeometry_def = GMX_NBLIST_GEOMETRY_CG_CG;
- if (log != NULL)
- {
- fprintf(log,"\nUsing charge-group - charge-group neighbor lists and kernels\n\n");
- }
- }
-
- if (fr->solvent_opt == esolTIP4P) {
- igeometry_w = GMX_NBLIST_GEOMETRY_WATER4_PARTICLE;
- igeometry_ww = GMX_NBLIST_GEOMETRY_WATER4_WATER4;
- } else {
- 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(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 */
- {
- 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,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(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;
+ /* Make maxlr tunable! (does not seem to be a big difference though)
+ * This parameter determines the number of i particles in a long range
+ * neighbourlist. Too few means many function calls, too many means
+ * cache trashing.
+ */
+ int maxsr, maxsr_wat, maxlr, maxlr_wat;
+ 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; */
+ maxsr = homenr;
+
+ if (maxsr < 0)
+ {
+ gmx_fatal(FARGS, "%s, %d: Negative number of short range atoms.\n"
+ "Call your Gromacs dealer for assistance.", __FILE__, __LINE__);
+ }
+ /* This is just for initial allocation, so we do not reallocate
+ * all the nlist arrays many times in a row.
+ * The numbers seem very accurate, but they are uncritical.
+ */
+ maxsr_wat = min(fr->nWatMol, (homenr+2)/3);
+ if (fr->bTwinRange)
+ {
+ maxlr = 50;
+ maxlr_wat = min(maxsr_wat, maxlr);
+ }
+ else
+ {
+ maxlr = maxlr_wat = 0;
+ }
+
+ /* 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;
+ bElecAndVdwSwitchDiffers = ( (fr->rcoulomb_switch != fr->rvdw_switch) || (fr->rcoulomb != fr->rvdw));
+
+ fr->ns.bCGlist = (getenv("GMX_NBLISTCG") != 0);
+ if (!fr->ns.bCGlist)
+ {
+ igeometry_def = GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE;
+ }
+ else
+ {
+ igeometry_def = GMX_NBLIST_GEOMETRY_CG_CG;
+ if (log != NULL)
+ {
+ fprintf(log, "\nUsing charge-group - charge-group neighbor lists and kernels\n\n");
+ }
+ }
+
+ if (fr->solvent_opt == esolTIP4P)
+ {
+ igeometry_w = GMX_NBLIST_GEOMETRY_WATER4_PARTICLE;
+ igeometry_ww = GMX_NBLIST_GEOMETRY_WATER4_WATER4;
+ }
+ else
+ {
+ 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]);
+
+ if ((fr->adress_type != eAdressOff) && (i >= fr->nnblists/2))
+ {
+ 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, 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, 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, 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, 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, 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, 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, bElecAndVdwSwitchDiffers);
+
+ /* 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
+ || nbl->nlist_sr[eNL_VDWQQ_WATERWATER].kernelptr_vf == NULL
+ || nbl->nlist_sr[eNL_QQ_WATERWATER].kernelptr_vf == NULL)
+ {
+ fr->solvent_opt = esolNO;
+ if (log != NULL)
+ {
+ fprintf(log, "Note: The available nonbonded kernels do not support water optimization - disabling.\n");
+ }
+ }
+
+ if (fr->efep != efepNO)
+ {
+ init_nblist(log, &nbl->nlist_sr[eNL_VDWQQ_FREE], &nbl->nlist_lr[eNL_VDWQQ_FREE],
+ 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, bElecAndVdwSwitchDiffers);
+ init_nblist(log, &nbl->nlist_sr[eNL_QQ_FREE], &nbl->nlist_lr[eNL_QQ_FREE],
+ 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, bElecAndVdwSwitchDiffers);
+ }
+
+ if (log != NULL)
+ {
+ fprintf(log, "\n");
+ }
+
+ fr->ns.nblist_initialized = TRUE;
}
static void reset_nblist(t_nblist *nl)
{
- nl->nri = -1;
- nl->nrj = 0;
- nl->maxlen = 0;
- if (nl->jindex)
- {
- nl->jindex[0] = 0;
- }
+ nl->nri = -1;
+ nl->nrj = 0;
+ if (nl->jindex)
+ {
+ nl->jindex[0] = 0;
+ }
}
-static void reset_neighbor_lists(t_forcerec *fr,gmx_bool bResetSR, gmx_bool bResetLR)
+static void reset_neighbor_lists(t_forcerec *fr, gmx_bool bResetSR, gmx_bool bResetLR)
{
- int n,i;
-
+ int n, i;
+
if (fr->bQMMM)
{
/* only reset the short-range nblist */
reset_nblist(&(fr->QMMMlist));
}
- 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)
+ if (bResetSR)
{
reset_nblist( &(fr->nblists[n].nlist_sr[i]) );
}
- if(bResetLR)
+ if (bResetLR)
{
reset_nblist( &(fr->nblists[n].nlist_lr[i]) );
}
-static inline void new_i_nblist(t_nblist *nlist,
- gmx_bool bLR,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;
-
+ int i, k, nri, nshift;
+
nri = nlist->nri;
-
+
/* Check whether we have to increase the i counter */
if ((nri == -1) ||
- (nlist->iinr[nri] != i_atom) ||
- (nlist->shift[nri] != shift) ||
+ (nlist->iinr[nri] != i_atom) ||
+ (nlist->shift[nri] != shift) ||
(nlist->gid[nri] != gid))
{
- /* This is something else. Now see if any entries have
+ /* This is something else. Now see if any entries have
* been added in the list of the previous atom.
*/
if ((nri == -1) ||
- ((nlist->jindex[nri+1] > nlist->jindex[nri]) &&
+ ((nlist->jindex[nri+1] > nlist->jindex[nri]) &&
(nlist->gid[nri] != -1)))
{
/* If so increase the counter */
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;
-
+
if (nri >= 0)
{
- nlist->jindex[nri+1] = nlist->nrj;
-
- len=nlist->nrj - nlist->jindex[nri];
-
- /* nlist length for water i molecules is treated statically
- * in the innerloops
+ /* Add elements up to padding. Since we allocate memory in units
+ * of the simd_padding width, we do not have to check for possible
+ * list reallocation here.
*/
- if (len > nlist->maxlen)
+ while ((nlist->nrj % nlist->simd_padding_width) != 0)
{
- nlist->maxlen = len;
+ /* Use -4 here, so we can write forces for 4 atoms before real data */
+ nlist->jjnr[nlist->nrj++] = -4;
}
+ nlist->jindex[nri+1] = nlist->nrj;
+
+ len = nlist->nrj - nlist->jindex[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.
}
}
-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;
-
+ int n, i;
+
if (bMakeQMMMnblist)
{
- close_nblist(&(fr->QMMMlist));
+ close_nblist(&(fr->QMMMlist));
}
- 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++)
{
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)
+static gmx_inline void add_j_to_nblist(t_nblist *nlist, atom_id j_atom, gmx_bool bLR)
{
- int nrj=nlist->nrj;
-
+ int nrj = nlist->nrj;
+
if (nlist->nrj >= nlist->maxnrj)
{
- nlist->maxnrj = over_alloc_small(nlist->nrj + 1);
+ nlist->maxnrj = round_up_to_simd_width(over_alloc_small(nlist->nrj + 1), nlist->simd_padding_width);
+
if (gmx_debug_at)
- 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);
+ {
+ fprintf(debug, "Increasing %s nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
+ bLR ? "LR" : "SR", nlist->ielec, nlist->ivdw, nlist->type, nlist->igeometry, nlist->maxnrj);
+ }
+
+ srenew(nlist->jjnr, nlist->maxnrj);
}
nlist->jjnr[nrj] = j_atom;
- nlist->nrj ++;
+ 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 nrj = nlist->nrj;
int j;
if (nlist->nrj >= nlist->maxnrj)
{
nlist->maxnrj = over_alloc_small(nlist->nrj + 1);
if (gmx_debug_at)
- 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);
- srenew(nlist->excl ,nlist->maxnrj*MAX_CGCGSIZE);
+ {
+ fprintf(debug, "Increasing %s nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
+ bLR ? "LR" : "SR", nlist->ielec, nlist->ivdw, nlist->type, nlist->igeometry, nlist->maxnrj);
+ }
+
+ srenew(nlist->jjnr, nlist->maxnrj);
+ srenew(nlist->jjnr_end, nlist->maxnrj);
+ srenew(nlist->excl, nlist->maxnrj*MAX_CGCGSIZE);
}
nlist->jjnr[nrj] = j_start;
if (j_end - j_start > MAX_CGCGSIZE)
{
- gmx_fatal(FARGS,"The charge-group - charge-group neighborlist do not support charge groups larger than %d, found a charge group of size %d",MAX_CGCGSIZE,j_end-j_start);
+ gmx_fatal(FARGS, "The charge-group - charge-group neighborlist do not support charge groups larger than %d, found a charge group of size %d", MAX_CGCGSIZE, j_end-j_start);
}
/* Set the exclusions */
- for(j=j_start; j<j_end; j++)
+ for (j = j_start; j < j_end; j++)
{
nlist->excl[nrj*MAX_CGCGSIZE + j - j_start] = bexcl[j];
}
if (i_is_j)
{
/* Avoid double counting of intra-cg interactions */
- for(j=1; j<j_end-j_start; j++)
+ for (j = 1; j < j_end-j_start; j++)
{
nlist->excl[nrj*MAX_CGCGSIZE + j] |= (1<<j) - 1;
}
}
- nlist->nrj ++;
+ nlist->nrj++;
}
typedef void
-put_in_list_t(gmx_bool bHaveVdW[],
- int ngid,
- t_mdatoms * md,
- int icg,
- int jgid,
- int nj,
- atom_id jjcg[],
- atom_id index[],
- t_excl bExcl[],
- int shift,
- t_forcerec * fr,
- gmx_bool bLR,
- gmx_bool bDoVdW,
- gmx_bool bDoCoul,
- int solvent_opt);
-
-static void
+ put_in_list_t (gmx_bool bHaveVdW[],
+ int ngid,
+ t_mdatoms * md,
+ int icg,
+ int jgid,
+ int nj,
+ atom_id jjcg[],
+ atom_id index[],
+ t_excl bExcl[],
+ int shift,
+ t_forcerec * fr,
+ gmx_bool bLR,
+ gmx_bool bDoVdW,
+ gmx_bool bDoCoul,
+ int solvent_opt);
+
+static void
put_in_list_at(gmx_bool bHaveVdW[],
- int ngid,
- t_mdatoms * md,
- int icg,
- int jgid,
- int nj,
- atom_id jjcg[],
- atom_id index[],
- t_excl bExcl[],
- int shift,
- t_forcerec * fr,
- gmx_bool bLR,
- gmx_bool bDoVdW,
- gmx_bool bDoCoul,
- int solvent_opt)
+ int ngid,
+ t_mdatoms * md,
+ int icg,
+ int jgid,
+ int nj,
+ atom_id jjcg[],
+ atom_id index[],
+ t_excl bExcl[],
+ int shift,
+ t_forcerec * fr,
+ 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].
*/
- t_nblist * vdwc;
- t_nblist * vdw;
- t_nblist * coul;
- t_nblist * vdwc_free = NULL;
- t_nblist * vdw_free = NULL;
- t_nblist * coul_free = NULL;
- t_nblist * vdwc_ww = NULL;
- t_nblist * coul_ww = NULL;
-
- int i,j,jcg,igid,gid,nbl_ind,ind_ij;
- atom_id jj,jj0,jj1,i_atom;
- int i0,nicg,len;
-
- int *cginfo;
- int *type,*typeB;
- real *charge,*chargeB;
- real qi,qiB,qq,rlj;
- gmx_bool bFreeEnergy,bFree,bFreeJ,bNotEx,*bPert;
- gmx_bool bDoVdW_i,bDoCoul_i,bDoCoul_i_sol;
- int iwater,jwater;
- t_nblist *nlist;
-
+ t_nblist * vdwc;
+ t_nblist * vdw;
+ t_nblist * coul;
+ t_nblist * vdwc_free = NULL;
+ t_nblist * vdw_free = NULL;
+ t_nblist * coul_free = NULL;
+ t_nblist * vdwc_ww = NULL;
+ t_nblist * coul_ww = NULL;
+
+ int i, j, jcg, igid, gid, nbl_ind, ind_ij;
+ atom_id jj, jj0, jj1, i_atom;
+ int i0, nicg, len;
+
+ int *cginfo;
+ int *type, *typeB;
+ real *charge, *chargeB;
+ real qi, qiB, qq, rlj;
+ gmx_bool bFreeEnergy, bFree, bFreeJ, bNotEx, *bPert;
+ gmx_bool bDoVdW_i, bDoCoul_i, bDoCoul_i_sol;
+ int iwater, jwater;
+ t_nblist *nlist;
+
/* Copy some pointers */
cginfo = fr->cginfo;
charge = md->chargeA;
type = md->typeA;
typeB = md->typeB;
bPert = md->bPerturbed;
-
+
/* Get atom range */
i0 = index[icg];
nicg = index[icg+1]-i0;
-
+
/* Get the i charge group info */
igid = GET_CGINFO_GID(cginfo[icg]);
- iwater = (solvent_opt!=esolNO) ? GET_CGINFO_SOLOPT(cginfo[icg]) : esolNO;
-
+ iwater = (solvent_opt != esolNO) ? GET_CGINFO_SOLOPT(cginfo[icg]) : esolNO;
+
bFreeEnergy = FALSE;
- if (md->nPerturbed)
+ if (md->nPerturbed)
{
- /* Check if any of the particles involved are perturbed.
+ /* Check if any of the particles involved are perturbed.
* If not we can do the cheaper normal put_in_list
* and use more solvent optimization.
*/
- for(i=0; i<nicg; i++)
+ for (i = 0; i < nicg; i++)
{
bFreeEnergy |= bPert[i0+i];
}
/* Loop over the j charge groups */
- for(j=0; (j<nj && !bFreeEnergy); j++)
+ for (j = 0; (j < nj && !bFreeEnergy); j++)
{
jcg = jjcg[j];
jj0 = index[jcg];
jj1 = index[jcg+1];
- /* Finally loop over the atoms in the j-charge group */
- for(jj=jj0; jj<jj1; jj++)
+ /* Finally loop over the atoms in the j-charge group */
+ for (jj = jj0; jj < jj1; jj++)
{
bFreeEnergy |= bPert[jj];
}
}
}
-
+
/* Unpack pointers to neighbourlist structs */
if (fr->nnblists == 1)
{
}
else
{
- nbl_ind = fr->gid2nblists[GID(igid,jgid,ngid)];
+ nbl_ind = fr->gid2nblists[GID(igid, jgid, ngid)];
}
if (bLR)
{
{
nlist = fr->nblists[nbl_ind].nlist_sr;
}
-
+
if (iwater != esolNO)
{
- vdwc = &nlist[eNL_VDWQQ_WATER];
- vdw = &nlist[eNL_VDW];
- coul = &nlist[eNL_QQ_WATER];
-#ifndef DISABLE_WATERWATER_NLIST
+ vdwc = &nlist[eNL_VDWQQ_WATER];
+ vdw = &nlist[eNL_VDW];
+ coul = &nlist[eNL_QQ_WATER];
vdwc_ww = &nlist[eNL_VDWQQ_WATERWATER];
coul_ww = &nlist[eNL_QQ_WATERWATER];
-#endif
- }
- else
+ }
+ else
{
vdwc = &nlist[eNL_VDWQQ];
vdw = &nlist[eNL_VDW];
coul = &nlist[eNL_QQ];
}
-
- if (!bFreeEnergy)
+
+ if (!bFreeEnergy)
{
- if (iwater != esolNO)
+ if (iwater != esolNO)
{
- /* Loop over the atoms in the i charge group */
+ /* Loop over the atoms in the i charge group */
i_atom = i0;
- gid = GID(igid,jgid,ngid);
+ gid = GID(igid, jgid, ngid);
/* Create new i_atom for each energy group */
if (bDoCoul && bDoVdW)
{
- new_i_nblist(vdwc,bLR,i_atom,shift,gid);
-#ifndef DISABLE_WATERWATER_NLIST
- new_i_nblist(vdwc_ww,bLR,i_atom,shift,gid);
-#endif
+ new_i_nblist(vdwc, i_atom, shift, gid);
+ new_i_nblist(vdwc_ww, i_atom, shift, gid);
}
if (bDoVdW)
{
- new_i_nblist(vdw,bLR,i_atom,shift,gid);
+ new_i_nblist(vdw, i_atom, shift, gid);
}
- if (bDoCoul)
+ if (bDoCoul)
{
- new_i_nblist(coul,bLR,i_atom,shift,gid);
-#ifndef DISABLE_WATERWATER_NLIST
- new_i_nblist(coul_ww,bLR,i_atom,shift,gid);
-#endif
- }
- /* Loop over the j charge groups */
- for(j=0; (j<nj); j++)
+ new_i_nblist(coul, i_atom, shift, gid);
+ new_i_nblist(coul_ww, i_atom, shift, gid);
+ }
+ /* Loop over the j charge groups */
+ for (j = 0; (j < nj); j++)
{
- jcg=jjcg[j];
-
+ jcg = jjcg[j];
+
if (jcg == icg)
{
continue;
}
-
- jj0 = index[jcg];
+
+ jj0 = index[jcg];
jwater = GET_CGINFO_SOLOPT(cginfo[jcg]);
-
+
if (iwater == esolSPC && jwater == esolSPC)
{
/* Interaction between two SPC molecules */
if (!bDoCoul)
{
/* VdW only - only first atoms in each water interact */
- add_j_to_nblist(vdw,jj0,bLR);
+ add_j_to_nblist(vdw, jj0, bLR);
}
- else
+ else
{
-#ifdef DISABLE_WATERWATER_NLIST
- /* Add entries for the three atoms - only do VdW if we need to */
- if (!bDoVdW)
- {
- add_j_to_nblist(coul,jj0,bLR);
- }
- else
- {
- add_j_to_nblist(vdwc,jj0,bLR);
- }
- add_j_to_nblist(coul,jj0+1,bLR);
- add_j_to_nblist(coul,jj0+2,bLR);
-#else
/* One entry for the entire water-water interaction */
if (!bDoVdW)
{
- add_j_to_nblist(coul_ww,jj0,bLR);
+ add_j_to_nblist(coul_ww, jj0, bLR);
}
else
{
- add_j_to_nblist(vdwc_ww,jj0,bLR);
+ add_j_to_nblist(vdwc_ww, jj0, bLR);
}
-#endif
- }
- }
- else if (iwater == esolTIP4P && jwater == esolTIP4P)
+ }
+ }
+ else if (iwater == esolTIP4P && jwater == esolTIP4P)
{
/* Interaction between two TIP4p molecules */
if (!bDoCoul)
{
/* VdW only - only first atoms in each water interact */
- add_j_to_nblist(vdw,jj0,bLR);
+ add_j_to_nblist(vdw, jj0, bLR);
}
- else
+ else
{
-#ifdef DISABLE_WATERWATER_NLIST
- /* Add entries for the four atoms - only do VdW if we need to */
- if (bDoVdW)
- {
- add_j_to_nblist(vdw,jj0,bLR);
- }
- add_j_to_nblist(coul,jj0+1,bLR);
- add_j_to_nblist(coul,jj0+2,bLR);
- add_j_to_nblist(coul,jj0+3,bLR);
-#else
/* One entry for the entire water-water interaction */
if (!bDoVdW)
{
- add_j_to_nblist(coul_ww,jj0,bLR);
+ add_j_to_nblist(coul_ww, jj0, bLR);
}
else
{
- add_j_to_nblist(vdwc_ww,jj0,bLR);
+ add_j_to_nblist(vdwc_ww, jj0, bLR);
}
-#endif
- }
+ }
}
- else
+ else
{
/* j charge group is not water, but i is.
* Add entries to the water-other_atom lists; the geometry of the water
* molecule doesn't matter - that is taken care of in the nonbonded kernel,
* so we don't care if it is SPC or TIP4P...
*/
-
+
jj1 = index[jcg+1];
-
- if (!bDoVdW)
+
+ if (!bDoVdW)
{
- for(jj=jj0; (jj<jj1); jj++)
+ for (jj = jj0; (jj < jj1); jj++)
{
if (charge[jj] != 0)
{
- add_j_to_nblist(coul,jj,bLR);
+ add_j_to_nblist(coul, jj, bLR);
}
}
}
else if (!bDoCoul)
{
- for(jj=jj0; (jj<jj1); jj++)
+ for (jj = jj0; (jj < jj1); jj++)
{
if (bHaveVdW[type[jj]])
{
- add_j_to_nblist(vdw,jj,bLR);
+ add_j_to_nblist(vdw, jj, bLR);
}
}
}
- else
+ else
{
/* _charge_ _groups_ interact with both coulomb and LJ */
/* Check which atoms we should add to the lists! */
- for(jj=jj0; (jj<jj1); jj++)
+ for (jj = jj0; (jj < jj1); jj++)
{
- if (bHaveVdW[type[jj]])
+ if (bHaveVdW[type[jj]])
{
if (charge[jj] != 0)
{
- add_j_to_nblist(vdwc,jj,bLR);
+ add_j_to_nblist(vdwc, jj, bLR);
}
else
{
- add_j_to_nblist(vdw,jj,bLR);
+ add_j_to_nblist(vdw, jj, bLR);
}
}
else if (charge[jj] != 0)
{
- add_j_to_nblist(coul,jj,bLR);
+ add_j_to_nblist(coul, jj, bLR);
}
}
}
}
}
- close_i_nblist(vdw);
- close_i_nblist(coul);
- close_i_nblist(vdwc);
-#ifndef DISABLE_WATERWATER_NLIST
+ close_i_nblist(vdw);
+ close_i_nblist(coul);
+ close_i_nblist(vdwc);
close_i_nblist(coul_ww);
- close_i_nblist(vdwc_ww);
-#endif
- }
+ close_i_nblist(vdwc_ww);
+ }
else
- {
+ {
/* no solvent as i charge group */
- /* Loop over the atoms in the i charge group */
- for(i=0; i<nicg; i++)
+ /* Loop over the atoms in the i charge group */
+ for (i = 0; i < nicg; i++)
{
i_atom = i0+i;
- gid = GID(igid,jgid,ngid);
+ gid = GID(igid, jgid, ngid);
qi = charge[i_atom];
-
+
/* Create new i_atom for each energy group */
if (bDoVdW && bDoCoul)
{
- new_i_nblist(vdwc,bLR,i_atom,shift,gid);
+ new_i_nblist(vdwc, i_atom, shift, gid);
}
if (bDoVdW)
{
- new_i_nblist(vdw,bLR,i_atom,shift,gid);
+ new_i_nblist(vdw, i_atom, shift, gid);
}
if (bDoCoul)
{
- new_i_nblist(coul,bLR,i_atom,shift,gid);
+ new_i_nblist(coul, i_atom, shift, gid);
}
bDoVdW_i = (bDoVdW && bHaveVdW[type[i_atom]]);
- bDoCoul_i = (bDoCoul && qi!=0);
-
- if (bDoVdW_i || bDoCoul_i)
+ bDoCoul_i = (bDoCoul && qi != 0);
+
+ if (bDoVdW_i || bDoCoul_i)
{
/* Loop over the j charge groups */
- for(j=0; (j<nj); j++)
+ for (j = 0; (j < nj); j++)
{
- jcg=jjcg[j];
-
+ jcg = jjcg[j];
+
/* Check for large charge groups */
if (jcg == icg)
{
{
jj0 = index[jcg];
}
-
- jj1=index[jcg+1];
- /* Finally loop over the atoms in the j-charge group */
- for(jj=jj0; jj<jj1; jj++)
+
+ jj1 = index[jcg+1];
+ /* Finally loop over the atoms in the j-charge group */
+ for (jj = jj0; jj < jj1; jj++)
{
- bNotEx = NOTEXCL(bExcl,i,jj);
-
- if (bNotEx)
+ bNotEx = NOTEXCL(bExcl, i, jj);
+
+ if (bNotEx)
{
- if (!bDoVdW_i)
- {
+ if (!bDoVdW_i)
+ {
if (charge[jj] != 0)
{
- add_j_to_nblist(coul,jj,bLR);
+ add_j_to_nblist(coul, jj, bLR);
}
}
- else if (!bDoCoul_i)
+ else if (!bDoCoul_i)
{
if (bHaveVdW[type[jj]])
{
- add_j_to_nblist(vdw,jj,bLR);
+ add_j_to_nblist(vdw, jj, bLR);
}
}
- else
+ else
{
- if (bHaveVdW[type[jj]])
+ if (bHaveVdW[type[jj]])
{
if (charge[jj] != 0)
{
- add_j_to_nblist(vdwc,jj,bLR);
+ add_j_to_nblist(vdwc, jj, bLR);
}
else
{
- add_j_to_nblist(vdw,jj,bLR);
+ add_j_to_nblist(vdw, jj, bLR);
}
- }
+ }
else if (charge[jj] != 0)
{
- add_j_to_nblist(coul,jj,bLR);
+ add_j_to_nblist(coul, jj, bLR);
}
}
}
vdwc_free = &nlist[eNL_VDWQQ_FREE];
vdw_free = &nlist[eNL_VDW_FREE];
coul_free = &nlist[eNL_QQ_FREE];
- /* Loop over the atoms in the i charge group */
- for(i=0; i<nicg; i++)
+ /* Loop over the atoms in the i charge group */
+ for (i = 0; i < nicg; i++)
{
i_atom = i0+i;
- gid = GID(igid,jgid,ngid);
+ gid = GID(igid, jgid, ngid);
qi = charge[i_atom];
qiB = chargeB[i_atom];
-
+
/* Create new i_atom for each energy group */
- if (bDoVdW && bDoCoul)
- new_i_nblist(vdwc,bLR,i_atom,shift,gid);
- if (bDoVdW)
- new_i_nblist(vdw,bLR,i_atom,shift,gid);
- if (bDoCoul)
- new_i_nblist(coul,bLR,i_atom,shift,gid);
-
- new_i_nblist(vdw_free,bLR,i_atom,shift,gid);
- new_i_nblist(coul_free,bLR,i_atom,shift,gid);
- new_i_nblist(vdwc_free,bLR,i_atom,shift,gid);
-
+ if (bDoVdW && bDoCoul)
+ {
+ new_i_nblist(vdwc, i_atom, shift, gid);
+ }
+ if (bDoVdW)
+ {
+ new_i_nblist(vdw, i_atom, shift, gid);
+ }
+ if (bDoCoul)
+ {
+ new_i_nblist(coul, i_atom, shift, gid);
+ }
+
+ new_i_nblist(vdw_free, i_atom, shift, gid);
+ new_i_nblist(coul_free, i_atom, shift, gid);
+ new_i_nblist(vdwc_free, i_atom, shift, gid);
+
bDoVdW_i = (bDoVdW &&
(bHaveVdW[type[i_atom]] || bHaveVdW[typeB[i_atom]]));
- bDoCoul_i = (bDoCoul && (qi!=0 || qiB!=0));
+ bDoCoul_i = (bDoCoul && (qi != 0 || qiB != 0));
/* For TIP4P the first atom does not have a charge,
* but the last three do. So we should still put an atom
* without LJ but with charge in the water-atom neighborlist
{
bDoCoul_i_sol = bDoCoul;
}
-
- if (bDoVdW_i || bDoCoul_i_sol)
+
+ if (bDoVdW_i || bDoCoul_i_sol)
{
/* Loop over the j charge groups */
- for(j=0; (j<nj); j++)
+ for (j = 0; (j < nj); j++)
{
- jcg=jjcg[j];
-
+ jcg = jjcg[j];
+
/* Check for large charge groups */
if (jcg == icg)
{
{
jj0 = index[jcg];
}
-
- jj1=index[jcg+1];
- /* Finally loop over the atoms in the j-charge group */
+
+ jj1 = index[jcg+1];
+ /* Finally loop over the atoms in the j-charge group */
bFree = bPert[i_atom];
- for(jj=jj0; (jj<jj1); jj++)
+ for (jj = jj0; (jj < jj1); jj++)
{
bFreeJ = bFree || bPert[jj];
/* Complicated if, because the water H's should also
* see perturbed j-particles
*/
- if (iwater==esolNO || i==0 || bFreeJ)
+ if (iwater == esolNO || i == 0 || bFreeJ)
{
- bNotEx = NOTEXCL(bExcl,i,jj);
-
- if (bNotEx)
+ bNotEx = NOTEXCL(bExcl, i, jj);
+
+ if (bNotEx)
{
if (bFreeJ)
{
- if (!bDoVdW_i)
+ if (!bDoVdW_i)
{
- if (charge[jj]!=0 || chargeB[jj]!=0)
+ if (charge[jj] != 0 || chargeB[jj] != 0)
{
- add_j_to_nblist(coul_free,jj,bLR);
+ add_j_to_nblist(coul_free, jj, bLR);
}
}
- else if (!bDoCoul_i)
+ else if (!bDoCoul_i)
{
if (bHaveVdW[type[jj]] || bHaveVdW[typeB[jj]])
{
- add_j_to_nblist(vdw_free,jj,bLR);
+ add_j_to_nblist(vdw_free, jj, bLR);
}
}
- else
+ else
{
- if (bHaveVdW[type[jj]] || bHaveVdW[typeB[jj]])
+ if (bHaveVdW[type[jj]] || bHaveVdW[typeB[jj]])
{
- if (charge[jj]!=0 || chargeB[jj]!=0)
+ if (charge[jj] != 0 || chargeB[jj] != 0)
{
- add_j_to_nblist(vdwc_free,jj,bLR);
+ add_j_to_nblist(vdwc_free, jj, bLR);
}
else
{
- add_j_to_nblist(vdw_free,jj,bLR);
+ add_j_to_nblist(vdw_free, jj, bLR);
}
}
- else if (charge[jj]!=0 || chargeB[jj]!=0)
- add_j_to_nblist(coul_free,jj,bLR);
+ else if (charge[jj] != 0 || chargeB[jj] != 0)
+ {
+ add_j_to_nblist(coul_free, jj, bLR);
+ }
}
}
- else if (!bDoVdW_i)
- {
+ else if (!bDoVdW_i)
+ {
/* This is done whether or not bWater is set */
if (charge[jj] != 0)
{
- add_j_to_nblist(coul,jj,bLR);
+ add_j_to_nblist(coul, jj, bLR);
}
}
- else if (!bDoCoul_i_sol)
- {
+ else if (!bDoCoul_i_sol)
+ {
if (bHaveVdW[type[jj]])
{
- add_j_to_nblist(vdw,jj,bLR);
+ add_j_to_nblist(vdw, jj, bLR);
}
}
- else
+ else
{
- if (bHaveVdW[type[jj]])
+ if (bHaveVdW[type[jj]])
{
if (charge[jj] != 0)
{
- add_j_to_nblist(vdwc,jj,bLR);
+ add_j_to_nblist(vdwc, jj, bLR);
}
else
{
- add_j_to_nblist(vdw,jj,bLR);
+ add_j_to_nblist(vdw, jj, bLR);
}
- }
+ }
else if (charge[jj] != 0)
{
- add_j_to_nblist(coul,jj,bLR);
+ add_j_to_nblist(coul, jj, bLR);
}
}
}
}
}
-static void
-put_in_list_qmmm(gmx_bool bHaveVdW[],
- int ngid,
- t_mdatoms * md,
- int icg,
- int jgid,
- int nj,
- atom_id jjcg[],
- atom_id index[],
- t_excl bExcl[],
- int shift,
- t_forcerec * fr,
- gmx_bool bLR,
- gmx_bool bDoVdW,
- gmx_bool bDoCoul,
- int solvent_opt)
+static void
+put_in_list_adress(gmx_bool bHaveVdW[],
+ int ngid,
+ t_mdatoms * md,
+ int icg,
+ int jgid,
+ int nj,
+ atom_id jjcg[],
+ atom_id index[],
+ t_excl bExcl[],
+ int shift,
+ t_forcerec * fr,
+ 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].
+ */
+ t_nblist * vdwc;
+ t_nblist * vdw;
+ t_nblist * coul;
+ t_nblist * vdwc_adress = NULL;
+ t_nblist * vdw_adress = NULL;
+ t_nblist * coul_adress = NULL;
+ t_nblist * vdwc_ww = NULL;
+ t_nblist * coul_ww = NULL;
+
+ int i, j, jcg, igid, gid, nbl_ind, nbl_ind_adress;
+ atom_id jj, jj0, jj1, i_atom;
+ int i0, nicg, len;
+
+ int *cginfo;
+ int *type, *typeB;
+ real *charge, *chargeB;
+ real *wf;
+ real qi, qiB, qq, rlj;
+ gmx_bool bFreeEnergy, bFree, bFreeJ, bNotEx, *bPert;
+ gmx_bool bDoVdW_i, bDoCoul_i, bDoCoul_i_sol;
+ gmx_bool b_hybrid;
+ gmx_bool j_all_atom;
+ int iwater, jwater;
+ t_nblist *nlist, *nlist_adress;
+ gmx_bool bEnergyGroupCG;
+
+ /* Copy some pointers */
+ cginfo = fr->cginfo;
+ charge = md->chargeA;
+ chargeB = md->chargeB;
+ type = md->typeA;
+ typeB = md->typeB;
+ bPert = md->bPerturbed;
+ wf = md->wf;
+
+ /* Get atom range */
+ i0 = index[icg];
+ nicg = index[icg+1]-i0;
+
+ /* Get the i charge group info */
+ igid = GET_CGINFO_GID(cginfo[icg]);
+
+ iwater = (solvent_opt != esolNO) ? GET_CGINFO_SOLOPT(cginfo[icg]) : esolNO;
+
+ if (md->nPerturbed)
+ {
+ gmx_fatal(FARGS, "AdResS does not support free energy pertubation\n");
+ }
+
+ /* Unpack pointers to neighbourlist structs */
+ if (fr->nnblists == 2)
+ {
+ nbl_ind = 0;
+ nbl_ind_adress = 1;
+ }
+ else
+ {
+ nbl_ind = fr->gid2nblists[GID(igid, jgid, ngid)];
+ nbl_ind_adress = nbl_ind+fr->nnblists/2;
+ }
+ if (bLR)
+ {
+ nlist = fr->nblists[nbl_ind].nlist_lr;
+ nlist_adress = fr->nblists[nbl_ind_adress].nlist_lr;
+ }
+ else
+ {
+ nlist = fr->nblists[nbl_ind].nlist_sr;
+ nlist_adress = fr->nblists[nbl_ind_adress].nlist_sr;
+ }
+
+
+ vdwc = &nlist[eNL_VDWQQ];
+ vdw = &nlist[eNL_VDW];
+ coul = &nlist[eNL_QQ];
+
+ vdwc_adress = &nlist_adress[eNL_VDWQQ];
+ vdw_adress = &nlist_adress[eNL_VDW];
+ coul_adress = &nlist_adress[eNL_QQ];
+
+ /* We do not support solvent optimization with AdResS for now.
+ For this we would need hybrid solvent-other kernels */
+
+ /* no solvent as i charge group */
+ /* Loop over the atoms in the i charge group */
+ for (i = 0; i < nicg; i++)
+ {
+ i_atom = i0+i;
+ gid = GID(igid, jgid, ngid);
+ qi = charge[i_atom];
+
+ /* Create new i_atom for each energy group */
+ if (bDoVdW && bDoCoul)
+ {
+ new_i_nblist(vdwc, i_atom, shift, gid);
+ new_i_nblist(vdwc_adress, i_atom, shift, gid);
+
+ }
+ if (bDoVdW)
+ {
+ new_i_nblist(vdw, i_atom, shift, gid);
+ new_i_nblist(vdw_adress, i_atom, shift, gid);
+
+ }
+ if (bDoCoul)
+ {
+ new_i_nblist(coul, i_atom, shift, gid);
+ new_i_nblist(coul_adress, i_atom, shift, gid);
+ }
+ bDoVdW_i = (bDoVdW && bHaveVdW[type[i_atom]]);
+ bDoCoul_i = (bDoCoul && qi != 0);
+
+ /* Here we find out whether the energy groups interaction belong to a
+ * coarse-grained (vsite) or atomistic interaction. Note that, beacuse
+ * interactions between coarse-grained and other (atomistic) energygroups
+ * are excluded automatically by grompp, it is sufficient to check for
+ * the group id of atom i (igid) */
+ bEnergyGroupCG = !egp_explicit(fr, igid);
+
+ if (bDoVdW_i || bDoCoul_i)
+ {
+ /* Loop over the j charge groups */
+ for (j = 0; (j < nj); j++)
+ {
+ jcg = jjcg[j];
+
+ /* Check for large charge groups */
+ if (jcg == icg)
+ {
+ jj0 = i0 + i + 1;
+ }
+ else
+ {
+ jj0 = index[jcg];
+ }
+
+ jj1 = index[jcg+1];
+ /* Finally loop over the atoms in the j-charge group */
+ for (jj = jj0; jj < jj1; jj++)
+ {
+ bNotEx = NOTEXCL(bExcl, i, jj);
+
+ /* Now we have to exclude interactions which will be zero
+ * anyway due to the AdResS weights (in previous implementations
+ * this was done in the force kernel). This is necessary as
+ * pure interactions (those with b_hybrid=false, i.e. w_i*w_j==1 or 0)
+ * are put into neighbour lists which will be passed to the
+ * standard (optimized) kernels for speed. The interactions with
+ * b_hybrid=true are placed into the _adress neighbour lists and
+ * processed by the generic AdResS kernel.
+ */
+ if ( (bEnergyGroupCG &&
+ wf[i_atom] >= 1-GMX_REAL_EPS && wf[jj] >= 1-GMX_REAL_EPS ) ||
+ ( !bEnergyGroupCG && wf[jj] <= GMX_REAL_EPS ) )
+ {
+ continue;
+ }
+
+ b_hybrid = !((wf[i_atom] >= 1-GMX_REAL_EPS && wf[jj] >= 1-GMX_REAL_EPS) ||
+ (wf[i_atom] <= GMX_REAL_EPS && wf[jj] <= GMX_REAL_EPS));
+
+ if (bNotEx)
+ {
+ if (!bDoVdW_i)
+ {
+ if (charge[jj] != 0)
+ {
+ if (!b_hybrid)
+ {
+ add_j_to_nblist(coul, jj, bLR);
+ }
+ else
+ {
+ add_j_to_nblist(coul_adress, jj, bLR);
+ }
+ }
+ }
+ else if (!bDoCoul_i)
+ {
+ if (bHaveVdW[type[jj]])
+ {
+ if (!b_hybrid)
+ {
+ add_j_to_nblist(vdw, jj, bLR);
+ }
+ else
+ {
+ add_j_to_nblist(vdw_adress, jj, bLR);
+ }
+ }
+ }
+ else
+ {
+ if (bHaveVdW[type[jj]])
+ {
+ if (charge[jj] != 0)
+ {
+ if (!b_hybrid)
+ {
+ add_j_to_nblist(vdwc, jj, bLR);
+ }
+ else
+ {
+ add_j_to_nblist(vdwc_adress, jj, bLR);
+ }
+ }
+ else
+ {
+ if (!b_hybrid)
+ {
+ add_j_to_nblist(vdw, jj, bLR);
+ }
+ else
+ {
+ add_j_to_nblist(vdw_adress, jj, bLR);
+ }
+
+ }
+ }
+ else if (charge[jj] != 0)
+ {
+ if (!b_hybrid)
+ {
+ add_j_to_nblist(coul, jj, bLR);
+ }
+ else
+ {
+ add_j_to_nblist(coul_adress, jj, bLR);
+ }
+
+ }
+ }
+ }
+ }
+ }
+
+ close_i_nblist(vdw);
+ close_i_nblist(coul);
+ close_i_nblist(vdwc);
+ close_i_nblist(vdw_adress);
+ close_i_nblist(coul_adress);
+ close_i_nblist(vdwc_adress);
+ }
+ }
+}
+
+static void
+put_in_list_qmmm(gmx_bool gmx_unused bHaveVdW[],
+ int ngid,
+ t_mdatoms gmx_unused * md,
+ int icg,
+ int jgid,
+ int nj,
+ atom_id jjcg[],
+ atom_id index[],
+ t_excl bExcl[],
+ int shift,
+ t_forcerec * fr,
+ gmx_bool bLR,
+ gmx_bool gmx_unused bDoVdW,
+ gmx_bool gmx_unused bDoCoul,
+ int gmx_unused solvent_opt)
{
- t_nblist * coul;
- int i,j,jcg,igid,gid;
- atom_id jj,jj0,jj1,i_atom;
- int i0,nicg;
+ t_nblist * coul;
+ int i, j, jcg, igid, gid;
+ atom_id jj, jj0, jj1, i_atom;
+ int i0, nicg;
gmx_bool bNotEx;
-
+
/* Get atom range */
i0 = index[icg];
nicg = index[icg+1]-i0;
-
+
/* Get the i charge group info */
igid = GET_CGINFO_GID(fr->cginfo[icg]);
-
+
coul = &fr->QMMMlist;
-
+
/* Loop over atoms in the ith charge group */
- for (i=0;i<nicg;i++)
+ for (i = 0; i < nicg; i++)
{
i_atom = i0+i;
- gid = GID(igid,jgid,ngid);
+ gid = GID(igid, jgid, ngid);
/* Create new i_atom for each energy group */
- new_i_nblist(coul,bLR,i_atom,shift,gid);
-
+ new_i_nblist(coul, i_atom, shift, gid);
+
/* Loop over the j charge groups */
- for (j=0;j<nj;j++)
+ for (j = 0; j < nj; j++)
{
- jcg=jjcg[j];
-
+ jcg = jjcg[j];
+
/* Charge groups cannot have QM and MM atoms simultaneously */
- if (jcg!=icg)
+ if (jcg != icg)
{
jj0 = index[jcg];
jj1 = index[jcg+1];
/* Finally loop over the atoms in the j-charge group */
- for(jj=jj0; jj<jj1; jj++)
+ for (jj = jj0; jj < jj1; jj++)
{
- bNotEx = NOTEXCL(bExcl,i,jj);
- if(bNotEx)
- add_j_to_nblist(coul,jj,bLR);
+ bNotEx = NOTEXCL(bExcl, i, jj);
+ if (bNotEx)
+ {
+ add_j_to_nblist(coul, jj, bLR);
+ }
}
}
}
}
}
-static void
-put_in_list_cg(gmx_bool bHaveVdW[],
- int ngid,
- t_mdatoms * md,
- int icg,
- int jgid,
- int nj,
- atom_id jjcg[],
- atom_id index[],
- t_excl bExcl[],
- int shift,
- t_forcerec * fr,
- gmx_bool bLR,
- gmx_bool bDoVdW,
- gmx_bool bDoCoul,
- int solvent_opt)
+static void
+put_in_list_cg(gmx_bool gmx_unused bHaveVdW[],
+ int ngid,
+ t_mdatoms gmx_unused * md,
+ int icg,
+ int jgid,
+ int nj,
+ atom_id jjcg[],
+ atom_id index[],
+ t_excl bExcl[],
+ int shift,
+ t_forcerec * fr,
+ gmx_bool bLR,
+ gmx_bool gmx_unused bDoVdW,
+ gmx_bool gmx_unused bDoCoul,
+ int gmx_unused solvent_opt)
{
int cginfo;
- int igid,gid,nbl_ind;
+ int igid, gid, nbl_ind;
t_nblist * vdwc;
- int j,jcg;
+ int j, jcg;
cginfo = fr->cginfo[icg];
igid = GET_CGINFO_GID(cginfo);
- gid = GID(igid,jgid,ngid);
+ gid = GID(igid, jgid, ngid);
/* Unpack pointers to neighbourlist structs */
if (fr->nnblists == 1)
* If required, zero interactions could be removed here
* or in the force loop.
*/
- new_i_nblist(vdwc,bLR,index[icg],shift,gid);
+ new_i_nblist(vdwc, index[icg], shift, gid);
vdwc->iinr_end[vdwc->nri] = index[icg+1];
- for(j=0; (j<nj); j++)
+ for (j = 0; (j < nj); j++)
{
jcg = jjcg[j];
/* Skip the icg-icg pairs if all self interactions are excluded */
/* Here we add the j charge group jcg to the list,
* exclusions are also added to the list.
*/
- add_j_to_nblist_cg(vdwc,index[jcg],index[jcg+1],bExcl,icg==jcg,bLR);
+ add_j_to_nblist_cg(vdwc, index[jcg], index[jcg+1], bExcl, icg == jcg, bLR);
}
}
- close_i_nblist(vdwc);
+ close_i_nblist(vdwc);
}
-static void setexcl(atom_id start,atom_id end,t_blocka *excl,gmx_bool b,
+static void setexcl(atom_id start, atom_id end, t_blocka *excl, gmx_bool b,
t_excl bexcl[])
{
- atom_id i,k;
-
+ atom_id i, k;
+
if (b)
{
- for(i=start; i<end; i++)
+ for (i = start; i < end; i++)
{
- for(k=excl->index[i]; k<excl->index[i+1]; k++)
+ for (k = excl->index[i]; k < excl->index[i+1]; k++)
{
- SETEXCL(bexcl,i-start,excl->a[k]);
+ SETEXCL(bexcl, i-start, excl->a[k]);
}
}
}
else
{
- for(i=start; i<end; i++)
+ for (i = start; i < end; i++)
{
- for(k=excl->index[i]; k<excl->index[i+1]; k++)
+ for (k = excl->index[i]; k < excl->index[i+1]; k++)
{
- RMEXCL(bexcl,i-start,excl->a[k]);
+ RMEXCL(bexcl, i-start, excl->a[k]);
}
}
}
}
-int calc_naaj(int icg,int cgtot)
+int calc_naaj(int icg, int cgtot)
{
int naaj;
-
+
if ((cgtot % 2) == 1)
{
/* Odd number of charge groups, easy */
}
else if ((cgtot % 4) == 0)
{
- /* Multiple of four is hard */
+ /* Multiple of four is hard */
if (icg < cgtot/2)
{
if ((icg % 2) == 0)
{
- naaj=1+(cgtot/2);
+ naaj = 1+(cgtot/2);
}
else
{
- naaj=cgtot/2;
+ naaj = cgtot/2;
}
}
else
{
if ((icg % 2) == 1)
{
- naaj=1+(cgtot/2);
+ naaj = 1+(cgtot/2);
}
else
{
- naaj=cgtot/2;
+ naaj = cgtot/2;
}
}
}
/* cgtot/2 = odd */
if ((icg % 2) == 0)
{
- naaj=1+(cgtot/2);
+ naaj = 1+(cgtot/2);
}
else
{
- naaj=cgtot/2;
+ naaj = cgtot/2;
}
}
#ifdef DEBUG
- fprintf(log,"naaj=%d\n",naaj);
+ fprintf(log, "naaj=%d\n", naaj);
#endif
return naaj;
*
************************************************/
-static real calc_image_tric(rvec xi,rvec xj,matrix box,
- rvec b_inv,int *shift)
+static real calc_image_tric(rvec xi, rvec xj, matrix box,
+ rvec b_inv, int *shift)
{
/* This code assumes that the cut-off is smaller than
* a half times the smallest diagonal element of the box.
*/
- const real h25=2.5;
- real dx,dy,dz;
- real r2;
- int tx,ty,tz;
-
+ const real h25 = 2.5;
+ real dx, dy, dz;
+ real r2;
+ int tx, ty, tz;
+
/* Compute diff vector */
dz = xj[ZZ] - xi[ZZ];
dy = xj[YY] - xi[YY];
dx = xj[XX] - xi[XX];
-
- /* Perform NINT operation, using trunc operation, therefore
- * we first add 2.5 then subtract 2 again
- */
- tz = dz*b_inv[ZZ] + h25;
+
+ /* Perform NINT operation, using trunc operation, therefore
+ * we first add 2.5 then subtract 2 again
+ */
+ tz = dz*b_inv[ZZ] + h25;
tz -= 2;
dz -= tz*box[ZZ][ZZ];
dy -= tz*box[ZZ][YY];
dx -= tz*box[ZZ][XX];
- ty = dy*b_inv[YY] + h25;
+ ty = dy*b_inv[YY] + h25;
ty -= 2;
dy -= ty*box[YY][YY];
dx -= ty*box[YY][XX];
-
- tx = dx*b_inv[XX]+h25;
+
+ tx = dx*b_inv[XX]+h25;
tx -= 2;
dx -= tx*box[XX][XX];
-
+
/* Distance squared */
r2 = (dx*dx) + (dy*dy) + (dz*dz);
- *shift = XYZ2IS(tx,ty,tz);
+ *shift = XYZ2IS(tx, ty, tz);
return r2;
}
-static real calc_image_rect(rvec xi,rvec xj,rvec box_size,
- rvec b_inv,int *shift)
+static real calc_image_rect(rvec xi, rvec xj, rvec box_size,
+ rvec b_inv, int *shift)
{
- const real h15=1.5;
- real ddx,ddy,ddz;
- real dx,dy,dz;
- real r2;
- int tx,ty,tz;
-
+ const real h15 = 1.5;
+ real ddx, ddy, ddz;
+ real dx, dy, dz;
+ real r2;
+ int tx, ty, tz;
+
/* Compute diff vector */
dx = xj[XX] - xi[XX];
dy = xj[YY] - xi[YY];
dz = xj[ZZ] - xi[ZZ];
-
+
/* Perform NINT operation, using trunc operation, therefore
* we first add 1.5 then subtract 1 again
*/
tx--;
ty--;
tz--;
-
+
/* Correct diff vector for translation */
ddx = tx*box_size[XX] - dx;
ddy = ty*box_size[YY] - dy;
ddz = tz*box_size[ZZ] - dz;
-
+
/* Distance squared */
r2 = (ddx*ddx) + (ddy*ddy) + (ddz*ddz);
-
- *shift = XYZ2IS(tx,ty,tz);
-
+
+ *shift = XYZ2IS(tx, ty, tz);
+
return r2;
}
-static void add_simple(t_ns_buf *nsbuf,int nrj,atom_id cg_j,
- gmx_bool bHaveVdW[],int ngid,t_mdatoms *md,
- int icg,int jgid,t_block *cgs,t_excl bexcl[],
- int shift,t_forcerec *fr,put_in_list_t *put_in_list)
+static void add_simple(t_ns_buf *nsbuf, int nrj, atom_id cg_j,
+ gmx_bool bHaveVdW[], int ngid, t_mdatoms *md,
+ int icg, int jgid, t_block *cgs, t_excl bexcl[],
+ int shift, t_forcerec *fr, put_in_list_t *put_in_list)
{
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,fr->solvent_opt);
+ put_in_list(bHaveVdW, ngid, md, icg, jgid, nsbuf->ncg, nsbuf->jcg,
+ cgs->index, bexcl, shift, fr, FALSE, TRUE, TRUE, fr->solvent_opt);
/* Reset buffer contents */
nsbuf->ncg = nsbuf->nj = 0;
}
nsbuf->jcg[nsbuf->ncg++] = cg_j;
- nsbuf->nj += nrj;
+ nsbuf->nj += nrj;
}
-static void ns_inner_tric(rvec x[],int icg,int *i_egp_flags,
- int njcg,atom_id jcg[],
- matrix box,rvec b_inv,real rcut2,
- t_block *cgs,t_ns_buf **ns_buf,
- gmx_bool bHaveVdW[],int ngid,t_mdatoms *md,
- t_excl bexcl[],t_forcerec *fr,
+static void ns_inner_tric(rvec x[], int icg, int *i_egp_flags,
+ int njcg, atom_id jcg[],
+ matrix box, rvec b_inv, real rcut2,
+ t_block *cgs, t_ns_buf **ns_buf,
+ gmx_bool bHaveVdW[], int ngid, t_mdatoms *md,
+ t_excl bexcl[], t_forcerec *fr,
put_in_list_t *put_in_list)
{
- int shift;
- int j,nrj,jgid;
- int *cginfo=fr->cginfo;
- atom_id cg_j,*cgindex;
+ int shift;
+ int j, nrj, jgid;
+ int *cginfo = fr->cginfo;
+ atom_id cg_j, *cgindex;
t_ns_buf *nsbuf;
-
+
cgindex = cgs->index;
shift = CENTRAL;
- for(j=0; (j<njcg); j++)
+ for (j = 0; (j < njcg); j++)
{
cg_j = jcg[j];
nrj = cgindex[cg_j+1]-cgindex[cg_j];
- if (calc_image_tric(x[icg],x[cg_j],box,b_inv,&shift) < rcut2)
+ if (calc_image_tric(x[icg], x[cg_j], box, b_inv, &shift) < rcut2)
{
jgid = GET_CGINFO_GID(cginfo[cg_j]);
if (!(i_egp_flags[jgid] & EGP_EXCL))
{
- add_simple(&ns_buf[jgid][shift],nrj,cg_j,
- bHaveVdW,ngid,md,icg,jgid,cgs,bexcl,shift,fr,
+ add_simple(&ns_buf[jgid][shift], nrj, cg_j,
+ bHaveVdW, ngid, md, icg, jgid, cgs, bexcl, shift, fr,
put_in_list);
}
}
}
}
-static void ns_inner_rect(rvec x[],int icg,int *i_egp_flags,
- int njcg,atom_id jcg[],
- gmx_bool bBox,rvec box_size,rvec b_inv,real rcut2,
- t_block *cgs,t_ns_buf **ns_buf,
- gmx_bool bHaveVdW[],int ngid,t_mdatoms *md,
- t_excl bexcl[],t_forcerec *fr,
+static void ns_inner_rect(rvec x[], int icg, int *i_egp_flags,
+ int njcg, atom_id jcg[],
+ gmx_bool bBox, rvec box_size, rvec b_inv, real rcut2,
+ t_block *cgs, t_ns_buf **ns_buf,
+ gmx_bool bHaveVdW[], int ngid, t_mdatoms *md,
+ t_excl bexcl[], t_forcerec *fr,
put_in_list_t *put_in_list)
{
- int shift;
- int j,nrj,jgid;
- int *cginfo=fr->cginfo;
- atom_id cg_j,*cgindex;
+ int shift;
+ int j, nrj, jgid;
+ int *cginfo = fr->cginfo;
+ atom_id cg_j, *cgindex;
t_ns_buf *nsbuf;
cgindex = cgs->index;
if (bBox)
{
shift = CENTRAL;
- for(j=0; (j<njcg); j++)
+ for (j = 0; (j < njcg); j++)
{
cg_j = jcg[j];
nrj = cgindex[cg_j+1]-cgindex[cg_j];
- if (calc_image_rect(x[icg],x[cg_j],box_size,b_inv,&shift) < rcut2)
+ if (calc_image_rect(x[icg], x[cg_j], box_size, b_inv, &shift) < rcut2)
{
jgid = GET_CGINFO_GID(cginfo[cg_j]);
if (!(i_egp_flags[jgid] & EGP_EXCL))
{
- add_simple(&ns_buf[jgid][shift],nrj,cg_j,
- bHaveVdW,ngid,md,icg,jgid,cgs,bexcl,shift,fr,
+ add_simple(&ns_buf[jgid][shift], nrj, cg_j,
+ bHaveVdW, ngid, md, icg, jgid, cgs, bexcl, shift, fr,
put_in_list);
}
}
}
else
{
- for(j=0; (j<njcg); j++)
+ for (j = 0; (j < njcg); j++)
{
cg_j = jcg[j];
nrj = cgindex[cg_j+1]-cgindex[cg_j];
- if ((rcut2 == 0) || (distance2(x[icg],x[cg_j]) < rcut2)) {
+ if ((rcut2 == 0) || (distance2(x[icg], x[cg_j]) < rcut2))
+ {
jgid = GET_CGINFO_GID(cginfo[cg_j]);
if (!(i_egp_flags[jgid] & EGP_EXCL))
{
- add_simple(&ns_buf[jgid][CENTRAL],nrj,cg_j,
- bHaveVdW,ngid,md,icg,jgid,cgs,bexcl,CENTRAL,fr,
+ add_simple(&ns_buf[jgid][CENTRAL], nrj, cg_j,
+ bHaveVdW, ngid, md, icg, jgid, cgs, bexcl, CENTRAL, fr,
put_in_list);
}
}
static int ns_simple_core(t_forcerec *fr,
gmx_localtop_t *top,
t_mdatoms *md,
- matrix box,rvec box_size,
- t_excl bexcl[],atom_id *aaj,
- int ngid,t_ns_buf **ns_buf,
- put_in_list_t *put_in_list,gmx_bool bHaveVdW[])
+ matrix box, rvec box_size,
+ t_excl bexcl[], atom_id *aaj,
+ int ngid, t_ns_buf **ns_buf,
+ put_in_list_t *put_in_list, gmx_bool bHaveVdW[])
{
- int naaj,k;
- real rlist2;
- int nsearch,icg,jcg,igid,i0,nri,nn;
- int *cginfo;
- t_ns_buf *nsbuf;
+ int naaj, k;
+ real rlist2;
+ int nsearch, icg, jcg, igid, i0, nri, nn;
+ int *cginfo;
+ t_ns_buf *nsbuf;
/* atom_id *i_atoms; */
- t_block *cgs=&(top->cgs);
- t_blocka *excl=&(top->excls);
- rvec b_inv;
- int m;
- gmx_bool bBox,bTriclinic;
- int *i_egp_flags;
-
+ t_block *cgs = &(top->cgs);
+ t_blocka *excl = &(top->excls);
+ rvec b_inv;
+ int m;
+ gmx_bool bBox, bTriclinic;
+ int *i_egp_flags;
+
rlist2 = sqr(fr->rlist);
-
+
bBox = (fr->ePBC != epbcNONE);
if (bBox)
{
- for(m=0; (m<DIM); m++)
+ for (m = 0; (m < DIM); m++)
{
- b_inv[m] = divide_err(1.0,box_size[m]);
+ b_inv[m] = divide_err(1.0, box_size[m]);
}
bTriclinic = TRICLINIC(box);
}
{
bTriclinic = FALSE;
}
-
+
cginfo = fr->cginfo;
-
- nsearch=0;
- for (icg=fr->cg0; (icg<fr->hcg); icg++)
+
+ nsearch = 0;
+ for (icg = fr->cg0; (icg < fr->hcg); icg++)
{
/*
- i0 = cgs->index[icg];
- nri = cgs->index[icg+1]-i0;
- i_atoms = &(cgs->a[i0]);
- i_eg_excl = fr->eg_excl + ngid*md->cENER[*i_atoms];
- setexcl(nri,i_atoms,excl,TRUE,bexcl);
- */
- igid = GET_CGINFO_GID(cginfo[icg]);
+ i0 = cgs->index[icg];
+ nri = cgs->index[icg+1]-i0;
+ i_atoms = &(cgs->a[i0]);
+ i_eg_excl = fr->eg_excl + ngid*md->cENER[*i_atoms];
+ setexcl(nri,i_atoms,excl,TRUE,bexcl);
+ */
+ igid = GET_CGINFO_GID(cginfo[icg]);
i_egp_flags = fr->egp_flags + ngid*igid;
- setexcl(cgs->index[icg],cgs->index[icg+1],excl,TRUE,bexcl);
-
- naaj=calc_naaj(icg,cgs->nr);
+ setexcl(cgs->index[icg], cgs->index[icg+1], excl, TRUE, bexcl);
+
+ naaj = calc_naaj(icg, cgs->nr);
if (bTriclinic)
{
- ns_inner_tric(fr->cg_cm,icg,i_egp_flags,naaj,&(aaj[icg]),
- box,b_inv,rlist2,cgs,ns_buf,
- bHaveVdW,ngid,md,bexcl,fr,put_in_list);
+ ns_inner_tric(fr->cg_cm, icg, i_egp_flags, naaj, &(aaj[icg]),
+ box, b_inv, rlist2, cgs, ns_buf,
+ bHaveVdW, ngid, md, bexcl, fr, put_in_list);
}
else
{
- ns_inner_rect(fr->cg_cm,icg,i_egp_flags,naaj,&(aaj[icg]),
- bBox,box_size,b_inv,rlist2,cgs,ns_buf,
- bHaveVdW,ngid,md,bexcl,fr,put_in_list);
+ ns_inner_rect(fr->cg_cm, icg, i_egp_flags, naaj, &(aaj[icg]),
+ bBox, box_size, b_inv, rlist2, cgs, ns_buf,
+ bHaveVdW, ngid, md, bexcl, fr, put_in_list);
}
nsearch += naaj;
-
- for(nn=0; (nn<ngid); nn++)
+
+ for (nn = 0; (nn < ngid); nn++)
{
- for(k=0; (k<SHIFTS); k++)
+ for (k = 0; (k < SHIFTS); k++)
{
nsbuf = &(ns_buf[nn][k]);
if (nsbuf->ncg > 0)
{
- put_in_list(bHaveVdW,ngid,md,icg,nn,nsbuf->ncg,nsbuf->jcg,
- cgs->index,bexcl,k,fr,FALSE,TRUE,TRUE,fr->solvent_opt);
- nsbuf->ncg=nsbuf->nj=0;
+ put_in_list(bHaveVdW, ngid, md, icg, nn, nsbuf->ncg, nsbuf->jcg,
+ cgs->index, bexcl, k, fr, FALSE, TRUE, TRUE, fr->solvent_opt);
+ nsbuf->ncg = nsbuf->nj = 0;
}
}
}
/* setexcl(nri,i_atoms,excl,FALSE,bexcl); */
- setexcl(cgs->index[icg],cgs->index[icg+1],excl,FALSE,bexcl);
+ setexcl(cgs->index[icg], cgs->index[icg+1], excl, FALSE, bexcl);
}
- close_neighbor_lists(fr,FALSE);
-
+ close_neighbor_lists(fr, FALSE);
+
return nsearch;
}
*
************************************************/
-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;
-
+ real dcx, tmp;
+ int xgi0, xgi1, i;
+
if (xgi < 0)
{
*dx0 = 0;
else
{
dcx2[xgi] = 0;
- *dx0 = xgi;
- xgi0 = xgi-1;
- *dx1 = xgi;
- xgi1 = xgi+1;
+ *dx0 = xgi;
+ xgi0 = xgi-1;
+ *dx1 = xgi;
+ xgi1 = xgi+1;
}
-
- for(i=xgi0; i>=0; i--)
+
+ for (i = xgi0; i >= 0; i--)
{
dcx = (i+1)*gridx-x;
tmp = dcx*dcx;
if (tmp >= rc2)
+ {
break;
- *dx0 = i;
+ }
+ *dx0 = i;
dcx2[i] = tmp;
}
- for(i=xgi1; i<Nx; i++)
+ for (i = xgi1; i < Nx; i++)
{
dcx = i*gridx-x;
tmp = dcx*dcx;
{
break;
}
- *dx1 = i;
+ *dx1 = i;
dcx2[i] = tmp;
}
}
-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;
-
+ real dcx, tmp;
+ int g_min, g_max, shift_home;
+
if (xgi < 0)
{
g_min = 0;
}
else
{
- *g0 = xgi;
- *g1 = xgi;
+ *g0 = xgi;
+ *g1 = xgi;
dcx2[xgi] = 0;
}
}
-
+
while (*g0 > g_min)
{
/* Check one grid cell down */
(*g0)--;
dcx2[*g0] = tmp;
}
-
+
while (*g1 < g_max)
{
/* Check one grid cell up */
#define sqr(x) ((x)*(x))
-#define calc_dx2(XI,YI,ZI,y) (sqr(XI-y[XX]) + sqr(YI-y[YY]) + sqr(ZI-y[ZZ]))
-#define calc_cyl_dx2(XI,YI,y) (sqr(XI-y[XX]) + sqr(YI-y[YY]))
+#define calc_dx2(XI, YI, ZI, y) (sqr(XI-y[XX]) + sqr(YI-y[YY]) + sqr(ZI-y[ZZ]))
+#define calc_cyl_dx2(XI, YI, y) (sqr(XI-y[XX]) + sqr(YI-y[YY]))
/****************************************************
*
* F A S T N E I G H B O R S E A R C H I N G
*
- * Optimized neighboursearching routine using grid
+ * Optimized neighboursearching routine using grid
* at least 1x1x1, see GROMACS manual
*
****************************************************/
-static void get_cutoff2(t_forcerec *fr,gmx_bool bDoLongRange,
- real *rvdw2,real *rcoul2,
- real *rs2,real *rm2,real *rl2)
+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,
+ /* 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
*rvdw2 = *rs2;
*rcoul2 = *rs2;
}
- *rm2 = min(*rvdw2,*rcoul2);
- *rl2 = max(*rvdw2,*rcoul2);
+ *rm2 = min(*rvdw2, *rcoul2);
+ *rl2 = max(*rvdw2, *rcoul2);
}
-static void init_nsgrid_lists(t_forcerec *fr,int ngid,gmx_ns_t *ns)
+static void init_nsgrid_lists(t_forcerec *fr, int ngid, gmx_ns_t *ns)
{
- real rvdw2,rcoul2,rs2,rm2,rl2;
- int j;
+ real rvdw2, rcoul2, rs2, rm2, rl2;
+ int j;
- get_cutoff2(fr,TRUE,&rvdw2,&rcoul2,&rs2,&rm2,&rl2);
+ get_cutoff2(fr, TRUE, &rvdw2, &rcoul2, &rs2, &rm2, &rl2);
/* Short range buffers */
- snew(ns->nl_sr,ngid);
+ snew(ns->nl_sr, ngid);
/* Counters */
- snew(ns->nsr,ngid);
- snew(ns->nlr_ljc,ngid);
- snew(ns->nlr_one,ngid);
-
+ snew(ns->nsr, ngid);
+ snew(ns->nlr_ljc, ngid);
+ snew(ns->nlr_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);
+ snew(ns->nl_lr_ljc, ngid);
/* Long range VdW or Coul only buffers */
- snew(ns->nl_lr_one,ngid);
+ snew(ns->nl_lr_one, ngid);
- for(j=0; (j<ngid); j++) {
- snew(ns->nl_sr[j],MAX_CG);
- snew(ns->nl_lr_ljc[j],MAX_CG);
- snew(ns->nl_lr_one[j],MAX_CG);
+ for (j = 0; (j < ngid); j++)
+ {
+ snew(ns->nl_sr[j], MAX_CG);
+ snew(ns->nl_lr_ljc[j], MAX_CG);
+ snew(ns->nl_lr_one[j], MAX_CG);
}
if (debug)
{
fprintf(debug,
"ns5_core: rs2 = %g, rm2 = %g, rl2 = %g (nm^2)\n",
- rs2,rm2,rl2);
+ rs2, rm2, rl2);
}
}
-static int nsgrid_core(FILE *log,t_commrec *cr,t_forcerec *fr,
- matrix box,rvec box_size,int ngid,
+static int nsgrid_core(t_commrec *cr, t_forcerec *fr,
+ matrix box, int ngid,
gmx_localtop_t *top,
- t_grid *grid,rvec x[],
- t_excl bexcl[],gmx_bool *bExcludeAlleg,
- t_nrnb *nrnb,t_mdatoms *md,
- real *lambda,real *dvdlambda,
- gmx_grppairener_t *grppener,
+ t_grid *grid,
+ t_excl bexcl[], gmx_bool *bExcludeAlleg,
+ t_mdatoms *md,
put_in_list_t *put_in_list,
gmx_bool bHaveVdW[],
- gmx_bool bDoLongRange,gmx_bool bMakeQMMMnblist)
+ gmx_bool bDoLongRange, gmx_bool bMakeQMMMnblist)
{
- gmx_ns_t *ns;
- atom_id **nl_lr_ljc,**nl_lr_one,**nl_sr;
- int *nlr_ljc,*nlr_one,*nsr;
- gmx_domdec_t *dd=NULL;
- t_block *cgs=&(top->cgs);
- int *cginfo=fr->cginfo;
+ gmx_ns_t *ns;
+ atom_id **nl_lr_ljc, **nl_lr_one, **nl_sr;
+ int *nlr_ljc, *nlr_one, *nsr;
+ gmx_domdec_t *dd;
+ t_block *cgs = &(top->cgs);
+ int *cginfo = fr->cginfo;
/* atom_id *i_atoms,*cgsindex=cgs->index; */
- ivec sh0,sh1,shp;
- int cell_x,cell_y,cell_z;
- int d,tx,ty,tz,dx,dy,dz,cj;
+ ivec sh0, sh1, shp;
+ int cell_x, cell_y, cell_z;
+ int d, tx, ty, tz, dx, dy, dz, cj;
#ifdef ALLOW_OFFDIAG_LT_HALFDIAG
- int zsh_ty,zsh_tx,ysh_tx;
+ int zsh_ty, zsh_tx, ysh_tx;
#endif
- int dx0,dx1,dy0,dy1,dz0,dz1;
- int Nx,Ny,Nz,shift=-1,j,nrj,nns,nn=-1;
- real gridx,gridy,gridz,grid_x,grid_y,grid_z;
- real *dcx2,*dcy2,*dcz2;
- int zgi,ygi,xgi;
- int cg0,cg1,icg=-1,cgsnr,i0,igid,nri,naaj,max_jcg;
- int jcg0,jcg1,jjcg,cgj0,jgid;
- int *grida,*gridnra,*gridind;
- gmx_bool rvdw_lt_rcoul,rcoul_lt_rvdw;
- rvec xi,*cgcm,grid_offset;
- real r2,rs2,rvdw2,rcoul2,rm2,rl2,XI,YI,ZI,dcx,dcy,dcz,tmp1,tmp2;
- int *i_egp_flags;
- gmx_bool bDomDec,bTriclinicX,bTriclinicY;
- ivec ncpddc;
-
+ int dx0, dx1, dy0, dy1, dz0, dz1;
+ int Nx, Ny, Nz, shift = -1, j, nrj, nns, nn = -1;
+ real gridx, gridy, gridz, grid_x, grid_y, grid_z;
+ real *dcx2, *dcy2, *dcz2;
+ int zgi, ygi, xgi;
+ int cg0, cg1, icg = -1, cgsnr, i0, igid, nri, naaj, max_jcg;
+ int jcg0, jcg1, jjcg, cgj0, jgid;
+ int *grida, *gridnra, *gridind;
+ gmx_bool rvdw_lt_rcoul, rcoul_lt_rvdw;
+ rvec xi, *cgcm, grid_offset;
+ real r2, rs2, rvdw2, rcoul2, rm2, rl2, XI, YI, ZI, dcx, dcy, dcz, tmp1, tmp2;
+ int *i_egp_flags;
+ gmx_bool bDomDec, bTriclinicX, bTriclinicY;
+ ivec ncpddc;
+
ns = &fr->ns;
-
+
bDomDec = DOMAINDECOMP(cr);
- if (bDomDec)
- {
- dd = cr->dd;
- }
-
+ dd = cr->dd;
+
bTriclinicX = ((YY < grid->npbcdim &&
- (!bDomDec || dd->nc[YY]==1) && box[YY][XX] != 0) ||
+ (!bDomDec || dd->nc[YY] == 1) && box[YY][XX] != 0) ||
(ZZ < grid->npbcdim &&
- (!bDomDec || dd->nc[ZZ]==1) && box[ZZ][XX] != 0));
+ (!bDomDec || dd->nc[ZZ] == 1) && box[ZZ][XX] != 0));
bTriclinicY = (ZZ < grid->npbcdim &&
- (!bDomDec || dd->nc[ZZ]==1) && box[ZZ][YY] != 0);
-
+ (!bDomDec || dd->nc[ZZ] == 1) && box[ZZ][YY] != 0);
+
cgsnr = cgs->nr;
- get_cutoff2(fr,bDoLongRange,&rvdw2,&rcoul2,&rs2,&rm2,&rl2);
+ get_cutoff2(fr, bDoLongRange, &rvdw2, &rcoul2, &rs2, &rm2, &rl2);
rvdw_lt_rcoul = (rvdw2 >= rcoul2);
rcoul_lt_rvdw = (rcoul2 >= rvdw2);
-
+
if (bMakeQMMMnblist)
{
rm2 = rl2;
nl_lr_one = ns->nl_lr_one;
nlr_ljc = ns->nlr_ljc;
nlr_one = ns->nlr_one;
-
+
/* Unpack arrays */
cgcm = fr->cg_cm;
Nx = grid->n[XX];
gridind = grid->index;
gridnra = grid->nra;
nns = 0;
-
+
gridx = grid->cell_size[XX];
gridy = grid->cell_size[YY];
gridz = grid->cell_size[ZZ];
grid_x = 1/gridx;
grid_y = 1/gridy;
grid_z = 1/gridz;
- copy_rvec(grid->cell_offset,grid_offset);
- copy_ivec(grid->ncpddc,ncpddc);
+ copy_rvec(grid->cell_offset, grid_offset);
+ copy_ivec(grid->ncpddc, ncpddc);
dcx2 = grid->dcx2;
dcy2 = grid->dcy2;
dcz2 = grid->dcz2;
-
+
#ifdef ALLOW_OFFDIAG_LT_HALFDIAG
zsh_ty = floor(-box[ZZ][YY]/box[YY][YY]+0.5);
zsh_tx = floor(-box[ZZ][XX]/box[XX][XX]+0.5);
ysh_tx = floor(-box[YY][XX]/box[XX][XX]+0.5);
- if (zsh_tx!=0 && ysh_tx!=0)
+ if (zsh_tx != 0 && ysh_tx != 0)
{
/* This could happen due to rounding, when both ratios are 0.5 */
ysh_tx = 0;
}
#endif
-
+
debug_gmx();
if (fr->n_tpi)
cg1 = grid->icg1;
/* Set the shift range */
- for(d=0; d<DIM; d++)
+ for (d = 0; d < DIM; d++)
{
sh0[d] = -1;
sh1[d] = 1;
}
}
}
-
+
/* Loop over charge groups */
- for(icg=cg0; (icg < cg1); icg++)
+ for (icg = cg0; (icg < cg1); icg++)
{
igid = GET_CGINFO_GID(cginfo[icg]);
/* Skip this charge group if all energy groups are excluded! */
{
continue;
}
-
+
i0 = cgs->index[icg];
-
+
if (bMakeQMMMnblist)
- {
+ {
/* Skip this charge group if it is not a QM atom while making a
* QM/MM neighbourlist
*/
- if (md->bQM[i0]==FALSE)
+ if (md->bQM[i0] == FALSE)
{
- continue; /* MM particle, go to next particle */
+ continue; /* MM particle, go to next particle */
}
-
+
/* Compute the number of charge groups that fall within the control
* of this one (icg)
*/
- naaj = calc_naaj(icg,cgsnr);
+ naaj = calc_naaj(icg, cgsnr);
jcg0 = icg;
jcg1 = icg + naaj;
- max_jcg = cgsnr;
- }
+ max_jcg = cgsnr;
+ }
else
- {
+ {
/* make a normal neighbourlist */
-
+
if (bDomDec)
{
/* Get the j charge-group and dd cell shift ranges */
- dd_get_ns_ranges(cr->dd,icg,&jcg0,&jcg1,sh0,sh1);
+ dd_get_ns_ranges(cr->dd, icg, &jcg0, &jcg1, sh0, sh1);
max_jcg = 0;
}
else
/* Compute the number of charge groups that fall within the control
* of this one (icg)
*/
- naaj = calc_naaj(icg,cgsnr);
+ naaj = calc_naaj(icg, cgsnr);
jcg0 = icg;
jcg1 = icg + naaj;
-
+
if (fr->n_tpi)
{
/* The i-particle is awlways the test particle,
}
}
}
-
+
i_egp_flags = fr->egp_flags + igid*ngid;
-
+
/* Set the exclusions for the atoms in charge group icg using a bitmask */
- setexcl(i0,cgs->index[icg+1],&top->excls,TRUE,bexcl);
-
- ci2xyz(grid,icg,&cell_x,&cell_y,&cell_z);
-
- /* Changed iicg to icg, DvdS 990115
- * (but see consistency check above, DvdS 990330)
+ setexcl(i0, cgs->index[icg+1], &top->excls, TRUE, bexcl);
+
+ ci2xyz(grid, icg, &cell_x, &cell_y, &cell_z);
+
+ /* Changed iicg to icg, DvdS 990115
+ * (but see consistency check above, DvdS 990330)
*/
#ifdef NS5DB
- fprintf(log,"icg=%5d, naaj=%5d, cell %d %d %d\n",
- icg,naaj,cell_x,cell_y,cell_z);
+ fprintf(log, "icg=%5d, naaj=%5d, cell %d %d %d\n",
+ icg, naaj, cell_x, cell_y, cell_z);
#endif
/* Loop over shift vectors in three dimensions */
- for (tz=-shp[ZZ]; tz<=shp[ZZ]; tz++)
+ for (tz = -shp[ZZ]; tz <= shp[ZZ]; tz++)
{
ZI = cgcm[icg][ZZ]+tz*box[ZZ][ZZ];
/* Calculate range of cells in Z direction that have the shift tz */
zgi = cell_z + tz*Nz;
#define FAST_DD_NS
#ifndef FAST_DD_NS
- get_dx(Nz,gridz,rl2,zgi,ZI,&dz0,&dz1,dcz2);
+ get_dx(Nz, gridz, rl2, zgi, ZI, &dz0, &dz1, dcz2);
#else
- get_dx_dd(Nz,gridz,rl2,zgi,ZI-grid_offset[ZZ],
- ncpddc[ZZ],sh0[ZZ],sh1[ZZ],&dz0,&dz1,dcz2);
+ get_dx_dd(Nz, gridz, rl2, zgi, ZI-grid_offset[ZZ],
+ ncpddc[ZZ], sh0[ZZ], sh1[ZZ], &dz0, &dz1, dcz2);
#endif
if (dz0 > dz1)
{
continue;
}
- for (ty=-shp[YY]; ty<=shp[YY]; ty++)
+ for (ty = -shp[YY]; ty <= shp[YY]; ty++)
{
YI = cgcm[icg][YY]+ty*box[YY][YY]+tz*box[ZZ][YY];
/* Calculate range of cells in Y direction that have the shift ty */
ygi = cell_y + ty*Ny;
}
#ifndef FAST_DD_NS
- get_dx(Ny,gridy,rl2,ygi,YI,&dy0,&dy1,dcy2);
+ get_dx(Ny, gridy, rl2, ygi, YI, &dy0, &dy1, dcy2);
#else
- get_dx_dd(Ny,gridy,rl2,ygi,YI-grid_offset[YY],
- ncpddc[YY],sh0[YY],sh1[YY],&dy0,&dy1,dcy2);
+ get_dx_dd(Ny, gridy, rl2, ygi, YI-grid_offset[YY],
+ ncpddc[YY], sh0[YY], sh1[YY], &dy0, &dy1, dcy2);
#endif
if (dy0 > dy1)
{
continue;
}
- for (tx=-shp[XX]; tx<=shp[XX]; tx++)
+ for (tx = -shp[XX]; tx <= shp[XX]; tx++)
{
XI = cgcm[icg][XX]+tx*box[XX][XX]+ty*box[YY][XX]+tz*box[ZZ][XX];
/* Calculate range of cells in X direction that have the shift tx */
xgi = cell_x + tx*Nx;
}
#ifndef FAST_DD_NS
- get_dx(Nx,gridx,rl2,xgi*Nx,XI,&dx0,&dx1,dcx2);
+ get_dx(Nx, gridx, rl2, xgi*Nx, XI, &dx0, &dx1, dcx2);
#else
- get_dx_dd(Nx,gridx,rl2,xgi,XI-grid_offset[XX],
- ncpddc[XX],sh0[XX],sh1[XX],&dx0,&dx1,dcx2);
+ get_dx_dd(Nx, gridx, rl2, xgi, XI-grid_offset[XX],
+ ncpddc[XX], sh0[XX], sh1[XX], &dx0, &dx1, dcx2);
#endif
if (dx0 > dx1)
{
}
/* Adress: an explicit cg that has a weigthing function of 0 is excluded
* from the neigbour list as it will not interact */
- if (fr->adress_type != eAdressOff){
- if (md->wf[cgs->index[icg]]==0 && egp_explicit(fr, igid)){
+ if (fr->adress_type != eAdressOff)
+ {
+ if (md->wf[cgs->index[icg]] <= GMX_REAL_EPS && egp_explicit(fr, igid))
+ {
continue;
}
}
- /* Get shift vector */
- shift=XYZ2IS(tx,ty,tz);
+ /* Get shift vector */
+ shift = XYZ2IS(tx, ty, tz);
#ifdef NS5DB
- range_check(shift,0,SHIFTS);
+ range_check(shift, 0, SHIFTS);
#endif
- for(nn=0; (nn<ngid); nn++)
+ for (nn = 0; (nn < ngid); nn++)
{
nsr[nn] = 0;
nlr_ljc[nn] = 0;
- nlr_one[nn] = 0;
+ nlr_one[nn] = 0;
}
#ifdef NS5DB
- fprintf(log,"shift: %2d, dx0,1: %2d,%2d, dy0,1: %2d,%2d, dz0,1: %2d,%2d\n",
- shift,dx0,dx1,dy0,dy1,dz0,dz1);
- fprintf(log,"cgcm: %8.3f %8.3f %8.3f\n",cgcm[icg][XX],
- cgcm[icg][YY],cgcm[icg][ZZ]);
- fprintf(log,"xi: %8.3f %8.3f %8.3f\n",XI,YI,ZI);
+ fprintf(log, "shift: %2d, dx0,1: %2d,%2d, dy0,1: %2d,%2d, dz0,1: %2d,%2d\n",
+ shift, dx0, dx1, dy0, dy1, dz0, dz1);
+ fprintf(log, "cgcm: %8.3f %8.3f %8.3f\n", cgcm[icg][XX],
+ cgcm[icg][YY], cgcm[icg][ZZ]);
+ fprintf(log, "xi: %8.3f %8.3f %8.3f\n", XI, YI, ZI);
#endif
- for (dx=dx0; (dx<=dx1); dx++)
+ for (dx = dx0; (dx <= dx1); dx++)
{
tmp1 = rl2 - dcx2[dx];
- for (dy=dy0; (dy<=dy1); dy++)
+ for (dy = dy0; (dy <= dy1); dy++)
{
tmp2 = tmp1 - dcy2[dy];
if (tmp2 > 0)
{
- for (dz=dz0; (dz<=dz1); dz++) {
- if (tmp2 > dcz2[dz]) {
+ for (dz = dz0; (dz <= dz1); dz++)
+ {
+ if (tmp2 > dcz2[dz])
+ {
/* Find grid-cell cj in which possible neighbours are */
- cj = xyz2ci(Ny,Nz,dx,dy,dz);
-
+ cj = xyz2ci(Ny, Nz, dx, dy, dz);
+
/* Check out how many cgs (nrj) there in this cell */
nrj = gridnra[cj];
-
+
/* Find the offset in the cg list */
cgj0 = gridind[cj];
-
+
/* Check if all j's are out of range so we
* can skip the whole cell.
* Should save some time, especially with DD.
{
continue;
}
-
+
/* Loop over cgs */
- for (j=0; (j<nrj); j++)
+ for (j = 0; (j < nrj); j++)
{
jjcg = grida[cgj0+j];
-
+
/* check whether this guy is in range! */
if ((jjcg >= jcg0 && jjcg < jcg1) ||
(jjcg < max_jcg))
{
- r2=calc_dx2(XI,YI,ZI,cgcm[jjcg]);
- if (r2 < rl2) {
+ r2 = calc_dx2(XI, YI, ZI, cgcm[jjcg]);
+ if (r2 < rl2)
+ {
/* jgid = gid[cgsatoms[cgsindex[jjcg]]]; */
jgid = GET_CGINFO_GID(cginfo[jjcg]);
/* check energy group exclusions */
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,fr->solvent_opt);
- nsr[jgid]=0;
+ put_in_list(bHaveVdW, ngid, md, icg, jgid,
+ nsr[jgid], nl_sr[jgid],
+ cgs->index, /* cgsatoms, */ bexcl,
+ shift, fr, FALSE, TRUE, TRUE, fr->solvent_opt);
+ nsr[jgid] = 0;
}
- nl_sr[jgid][nsr[jgid]++]=jjcg;
- }
+ nl_sr[jgid][nsr[jgid]++] = jjcg;
+ }
else if (r2 < rm2)
{
if (nlr_ljc[jgid] >= MAX_CG)
{
/* 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;
+ 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;
+ nl_lr_ljc[jgid][nlr_ljc[jgid]++] = jjcg;
}
else
{
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;
+ 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;
+ nl_lr_one[jgid][nlr_one[jgid]++] = jjcg;
}
}
}
}
}
/* CHECK whether there is anything left in the buffers */
- for(nn=0; (nn<ngid); nn++)
+ for (nn = 0; (nn < ngid); nn++)
{
if (nsr[nn] > 0)
{
- put_in_list(bHaveVdW,ngid,md,icg,nn,nsr[nn],nl_sr[nn],
+ put_in_list(bHaveVdW, ngid, md, icg, nn, nsr[nn], nl_sr[nn],
cgs->index, /* cgsatoms, */ bexcl,
- shift,fr,FALSE,TRUE,TRUE,fr->solvent_opt);
+ shift, fr, FALSE, TRUE, TRUE, fr->solvent_opt);
}
-
+
if (nlr_ljc[nn] > 0)
{
- 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);
+ 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)
{
- 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);
+ 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);
}
}
}
}
}
/* setexcl(nri,i_atoms,&top->atoms.excl,FALSE,bexcl); */
- setexcl(cgs->index[icg],cgs->index[icg+1],&top->excls,FALSE,bexcl);
+ setexcl(cgs->index[icg], cgs->index[icg+1], &top->excls, FALSE, bexcl);
}
/* 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 neighbourlists */
- close_neighbor_lists(fr,bMakeQMMMnblist);
-
+ close_neighbor_lists(fr, bMakeQMMMnblist);
+
return nns;
}
-void ns_realloc_natoms(gmx_ns_t *ns,int natoms)
+void ns_realloc_natoms(gmx_ns_t *ns, int natoms)
{
int i;
-
+
if (natoms > ns->nra_alloc)
{
ns->nra_alloc = over_alloc_dd(natoms);
- srenew(ns->bexcl,ns->nra_alloc);
- for(i=0; i<ns->nra_alloc; i++)
+ srenew(ns->bexcl, ns->nra_alloc);
+ for (i = 0; i < ns->nra_alloc; i++)
{
ns->bexcl[i] = 0;
}
}
}
-void init_ns(FILE *fplog,const t_commrec *cr,
- gmx_ns_t *ns,t_forcerec *fr,
- const gmx_mtop_t *mtop,
- matrix box)
+void init_ns(FILE *fplog, const t_commrec *cr,
+ gmx_ns_t *ns, t_forcerec *fr,
+ const gmx_mtop_t *mtop)
{
- int mt,icg,nr_in_cg,maxcg,i,j,jcg,ngid,ncg;
+ int mt, icg, nr_in_cg, maxcg, i, j, jcg, ngid, ncg;
t_block *cgs;
char *ptr;
-
+
/* Compute largest charge groups size (# atoms) */
- nr_in_cg=1;
- for(mt=0; mt<mtop->nmoltype; mt++) {
+ nr_in_cg = 1;
+ for (mt = 0; mt < mtop->nmoltype; mt++)
+ {
cgs = &mtop->moltype[mt].cgs;
- for (icg=0; (icg < cgs->nr); icg++)
+ for (icg = 0; (icg < cgs->nr); icg++)
{
- nr_in_cg=max(nr_in_cg,(int)(cgs->index[icg+1]-cgs->index[icg]));
+ nr_in_cg = max(nr_in_cg, (int)(cgs->index[icg+1]-cgs->index[icg]));
}
}
/* Verify whether largest charge group is <= max cg.
- * This is determined by the type of the local exclusion type
+ * This is determined by the type of the local exclusion type
* Exclusions are stored in bits. (If the type is not large
* enough, enlarge it, unsigned char -> unsigned short -> unsigned long)
*/
maxcg = sizeof(t_excl)*8;
if (nr_in_cg > maxcg)
{
- gmx_fatal(FARGS,"Max #atoms in a charge group: %d > %d\n",
- nr_in_cg,maxcg);
+ gmx_fatal(FARGS, "Max #atoms in a charge group: %d > %d\n",
+ nr_in_cg, maxcg);
}
-
+
ngid = mtop->groups.grps[egcENER].nr;
- snew(ns->bExcludeAlleg,ngid);
- for(i=0; i<ngid; i++) {
+ snew(ns->bExcludeAlleg, ngid);
+ for (i = 0; i < ngid; i++)
+ {
ns->bExcludeAlleg[i] = TRUE;
- for(j=0; j<ngid; j++)
+ for (j = 0; j < ngid; j++)
{
if (!(fr->egp_flags[i*ngid+j] & EGP_EXCL))
{
}
}
}
-
- if (fr->bGrid) {
+
+ if (fr->bGrid)
+ {
/* Grid search */
- ns->grid = init_grid(fplog,fr);
- init_nsgrid_lists(fr,ngid,ns);
+ ns->grid = init_grid(fplog, fr);
+ init_nsgrid_lists(fr, ngid, ns);
}
else
{
/* Simple search */
- snew(ns->ns_buf,ngid);
- for(i=0; (i<ngid); i++)
+ snew(ns->ns_buf, ngid);
+ for (i = 0; (i < ngid); i++)
{
- snew(ns->ns_buf[i],SHIFTS);
+ snew(ns->ns_buf[i], SHIFTS);
}
ncg = ncg_mtop(mtop);
- snew(ns->simple_aaj,2*ncg);
- for(jcg=0; (jcg<ncg); jcg++)
+ snew(ns->simple_aaj, 2*ncg);
+ for (jcg = 0; (jcg < ncg); jcg++)
{
ns->simple_aaj[jcg] = jcg;
ns->simple_aaj[jcg+ncg] = jcg;
}
}
-
+
/* Create array that determines whether or not atoms have VdW */
- snew(ns->bHaveVdW,fr->ntype);
- for(i=0; (i<fr->ntype); i++)
+ snew(ns->bHaveVdW, fr->ntype);
+ for (i = 0; (i < fr->ntype); i++)
{
- for(j=0; (j<fr->ntype); j++)
+ for (j = 0; (j < fr->ntype); j++)
{
- ns->bHaveVdW[i] = (ns->bHaveVdW[i] ||
- (fr->bBHAM ?
- ((BHAMA(fr->nbfp,fr->ntype,i,j) != 0) ||
- (BHAMB(fr->nbfp,fr->ntype,i,j) != 0) ||
- (BHAMC(fr->nbfp,fr->ntype,i,j) != 0)) :
- ((C6(fr->nbfp,fr->ntype,i,j) != 0) ||
- (C12(fr->nbfp,fr->ntype,i,j) != 0))));
+ ns->bHaveVdW[i] = (ns->bHaveVdW[i] ||
+ (fr->bBHAM ?
+ ((BHAMA(fr->nbfp, fr->ntype, i, j) != 0) ||
+ (BHAMB(fr->nbfp, fr->ntype, i, j) != 0) ||
+ (BHAMC(fr->nbfp, fr->ntype, i, j) != 0)) :
+ ((C6(fr->nbfp, fr->ntype, i, j) != 0) ||
+ (C12(fr->nbfp, fr->ntype, i, j) != 0))));
}
}
- if (debug)
- pr_bvec(debug,0,"bHaveVdW",ns->bHaveVdW,fr->ntype,TRUE);
-
+ if (debug)
+ {
+ pr_bvec(debug, 0, "bHaveVdW", ns->bHaveVdW, fr->ntype, TRUE);
+ }
+
ns->nra_alloc = 0;
- ns->bexcl = NULL;
+ ns->bexcl = NULL;
if (!DOMAINDECOMP(cr))
{
- /* This could be reduced with particle decomposition */
- ns_realloc_natoms(ns,mtop->natoms);
+ ns_realloc_natoms(ns, mtop->natoms);
}
- ns->nblist_initialized=FALSE;
+ ns->nblist_initialized = FALSE;
/* nbr list debug dump */
{
- char *ptr=getenv("GMX_DUMP_NL");
+ char *ptr = getenv("GMX_DUMP_NL");
if (ptr)
{
- ns->dump_nl=strtol(ptr,NULL,10);
+ ns->dump_nl = strtol(ptr, NULL, 10);
if (fplog)
{
fprintf(fplog, "GMX_DUMP_NL = %d", ns->dump_nl);
}
else
{
- ns->dump_nl=0;
+ ns->dump_nl = 0;
}
}
}
-
-int search_neighbours(FILE *log,t_forcerec *fr,
- rvec x[],matrix box,
+
+int search_neighbours(FILE *log, t_forcerec *fr,
+ matrix box,
gmx_localtop_t *top,
gmx_groups_t *groups,
t_commrec *cr,
- t_nrnb *nrnb,t_mdatoms *md,
- real *lambda,real *dvdlambda,
- gmx_grppairener_t *grppener,
+ t_nrnb *nrnb, t_mdatoms *md,
gmx_bool bFillGrid,
gmx_bool bDoLongRangeNS)
{
- t_block *cgs=&(top->cgs);
- rvec box_size,grid_x0,grid_x1;
- int i,j,m,ngid;
- real min_size,grid_dens;
+ t_block *cgs = &(top->cgs);
+ rvec box_size, grid_x0, grid_x1;
+ int i, j, m, ngid;
+ real min_size, grid_dens;
int nsearch;
gmx_bool bGrid;
char *ptr;
gmx_bool *i_egp_flags;
- int cg_start,cg_end,start,end;
+ int cg_start, cg_end, start, end;
gmx_ns_t *ns;
t_grid *grid;
gmx_domdec_zones_t *dd_zones;
/* Set some local variables */
bGrid = fr->bGrid;
- ngid = groups->grps[egcENER].nr;
-
- for(m=0; (m<DIM); m++)
+ ngid = groups->grps[egcENER].nr;
+
+ for (m = 0; (m < DIM); m++)
{
box_size[m] = box[m][m];
}
-
+
if (fr->ePBC != epbcNONE)
{
- if (sqr(fr->rlistlong) >= max_cutoff2(fr->ePBC,box))
+ if (sqr(fr->rlistlong) >= max_cutoff2(fr->ePBC, box))
{
- gmx_fatal(FARGS,"One of the box vectors has become shorter than twice the cut-off length or box_yy-|box_zy| or box_zz has become smaller than the cut-off.");
+ gmx_fatal(FARGS, "One of the box vectors has become shorter than twice the cut-off length or box_yy-|box_zy| or box_zz has become smaller than the cut-off.");
}
if (!bGrid)
{
- min_size = min(box_size[XX],min(box_size[YY],box_size[ZZ]));
+ min_size = min(box_size[XX], min(box_size[YY], box_size[ZZ]));
if (2*fr->rlistlong >= min_size)
- gmx_fatal(FARGS,"One of the box diagonal elements has become smaller than twice the cut-off length.");
+ {
+ gmx_fatal(FARGS, "One of the box diagonal elements has become smaller than twice the cut-off length.");
+ }
}
}
-
+
if (DOMAINDECOMP(cr))
{
- ns_realloc_natoms(ns,cgs->index[cgs->nr]);
+ ns_realloc_natoms(ns, cgs->index[cgs->nr]);
}
debug_gmx();
-
+
/* Reset the neighbourlists */
- reset_neighbor_lists(fr,TRUE,TRUE);
-
+ reset_neighbor_lists(fr, TRUE, TRUE);
+
if (bGrid && bFillGrid)
{
-
+
grid = ns->grid;
if (DOMAINDECOMP(cr))
{
{
dd_zones = NULL;
- get_nsgrid_boundaries(grid->nboundeddim,box,NULL,NULL,NULL,NULL,
- cgs->nr,fr->cg_cm,grid_x0,grid_x1,&grid_dens);
+ get_nsgrid_boundaries(grid->nboundeddim, box, NULL, NULL, NULL, NULL,
+ cgs->nr, fr->cg_cm, grid_x0, grid_x1, &grid_dens);
- grid_first(log,grid,NULL,NULL,fr->ePBC,box,grid_x0,grid_x1,
- fr->rlistlong,grid_dens);
+ grid_first(log, grid, NULL, NULL, box, grid_x0, grid_x1,
+ fr->rlistlong, grid_dens);
}
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))
{
end = cgs->nr;
- fill_grid(log,dd_zones,grid,end,-1,end,fr->cg_cm);
+ fill_grid(dd_zones, grid, end, -1, end, fr->cg_cm);
grid->icg0 = 0;
grid->icg1 = dd_zones->izone[dd_zones->nizone-1].cg1;
}
else
{
- fill_grid(log,NULL,grid,cgs->nr,fr->cg0,fr->hcg,fr->cg_cm);
+ fill_grid(NULL, grid, cgs->nr, fr->cg0, fr->hcg, fr->cg_cm);
grid->icg0 = fr->cg0;
grid->icg1 = fr->hcg;
debug_gmx();
-
- if (PARTDECOMP(cr))
- mv_grid(cr,grid);
- debug_gmx();
}
-
- calc_elemnr(log,grid,start,end,cgs->nr);
+
+ calc_elemnr(grid, start, end, cgs->nr);
calc_ptrs(grid);
- grid_last(log,grid,start,end,cgs->nr);
-
+ grid_last(grid, start, end, cgs->nr);
+
if (gmx_debug_at)
{
- check_grid(debug,grid);
- print_grid(debug,grid);
+ check_grid(grid);
+ print_grid(debug, grid);
}
}
else if (fr->n_tpi)
/* Set the grid cell index for the test particle only.
* The cell to cg index is not corrected, but that does not matter.
*/
- fill_grid(log,NULL,ns->grid,fr->hcg,fr->hcg-1,fr->hcg,fr->cg_cm);
+ fill_grid(NULL, ns->grid, fr->hcg, fr->hcg-1, fr->hcg, fr->cg_cm);
}
debug_gmx();
-
- if (!fr->ns.bCGlist)
+
+ if (fr->adress_type == eAdressOff)
{
- put_in_list = put_in_list_at;
+ if (!fr->ns.bCGlist)
+ {
+ put_in_list = put_in_list_at;
+ }
+ else
+ {
+ put_in_list = put_in_list_cg;
+ }
}
else
{
- put_in_list = put_in_list_cg;
+ put_in_list = put_in_list_adress;
}
/* Do the core! */
if (bGrid)
{
- grid = ns->grid;
- nsearch = nsgrid_core(log,cr,fr,box,box_size,ngid,top,
- grid,x,ns->bexcl,ns->bExcludeAlleg,
- nrnb,md,lambda,dvdlambda,grppener,
- put_in_list,ns->bHaveVdW,
- bDoLongRangeNS,FALSE);
-
+ grid = ns->grid;
+ nsearch = nsgrid_core(cr, fr, box, ngid, top,
+ grid, ns->bexcl, ns->bExcludeAlleg,
+ md, put_in_list, ns->bHaveVdW,
+ bDoLongRangeNS, FALSE);
+
/* neighbour searching withouth QMMM! QM atoms have zero charge in
* the classical calculation. The charge-charge interaction
* between QM and MM atoms is handled in the QMMM core calculation
* and the QM MM atom pairs have just been put in the
* corresponding neighbourlists. in case of QMMM we still need to
* fill a special QMMM neighbourlist that contains all neighbours
- * of the QM atoms. If bQMMM is true, this list will now be made:
+ * of the QM atoms. If bQMMM is true, this list will now be made:
*/
- if (fr->bQMMM && fr->qr->QMMMscheme!=eQMMMschemeoniom)
+ if (fr->bQMMM && fr->qr->QMMMscheme != eQMMMschemeoniom)
{
- nsearch += nsgrid_core(log,cr,fr,box,box_size,ngid,top,
- grid,x,ns->bexcl,ns->bExcludeAlleg,
- nrnb,md,lambda,dvdlambda,grppener,
- put_in_list_qmmm,ns->bHaveVdW,
- bDoLongRangeNS,TRUE);
+ nsearch += nsgrid_core(cr, fr, box, ngid, top,
+ grid, ns->bexcl, ns->bExcludeAlleg,
+ md, put_in_list_qmmm, ns->bHaveVdW,
+ bDoLongRangeNS, TRUE);
}
}
- else
+ else
{
- nsearch = ns_simple_core(fr,top,md,box,box_size,
- ns->bexcl,ns->simple_aaj,
- ngid,ns->ns_buf,put_in_list,ns->bHaveVdW);
+ nsearch = ns_simple_core(fr, top, md, box, box_size,
+ ns->bexcl, ns->simple_aaj,
+ ngid, ns->ns_buf, put_in_list, ns->bHaveVdW);
}
debug_gmx();
#ifdef DEBUG
pr_nsblock(log);
#endif
-
- inc_nrnb(nrnb,eNR_NS,nsearch);
+
+ inc_nrnb(nrnb, eNR_NS, nsearch);
/* inc_nrnb(nrnb,eNR_LR,fr->nlr); */
-
+
return nsearch;
}
-int natoms_beyond_ns_buffer(t_inputrec *ir,t_forcerec *fr,t_block *cgs,
- matrix scale_tot,rvec *x)
+int natoms_beyond_ns_buffer(t_inputrec *ir, t_forcerec *fr, t_block *cgs,
+ matrix scale_tot, rvec *x)
{
- int cg0,cg1,cg,a0,a1,a,i,j;
- real rint,hbuf2,scale;
- rvec *cg_cm,cgsc;
+ int cg0, cg1, cg, a0, a1, a, i, j;
+ real rint, hbuf2, scale;
+ rvec *cg_cm, cgsc;
gmx_bool bIsotropic;
int nBeyond;
-
+
nBeyond = 0;
-
- rint = max(ir->rcoulomb,ir->rvdw);
+
+ rint = max(ir->rcoulomb, ir->rvdw);
if (ir->rlist < rint)
{
- gmx_fatal(FARGS,"The neighbor search buffer has negative size: %f nm",
+ gmx_fatal(FARGS, "The neighbor search buffer has negative size: %f nm",
ir->rlist - rint);
}
cg_cm = fr->cg_cm;
-
+
cg0 = fr->cg0;
cg1 = fr->hcg;
-
+
if (!EI_DYNAMICS(ir->eI) || !DYNAMIC_BOX(*ir))
{
hbuf2 = sqr(0.5*(ir->rlist - rint));
- for(cg=cg0; cg<cg1; cg++)
+ for (cg = cg0; cg < cg1; cg++)
{
a0 = cgs->index[cg];
a1 = cgs->index[cg+1];
- for(a=a0; a<a1; a++)
+ for (a = a0; a < a1; a++)
{
- if (distance2(cg_cm[cg],x[a]) > hbuf2)
+ if (distance2(cg_cm[cg], x[a]) > hbuf2)
{
nBeyond++;
}
else
{
bIsotropic = TRUE;
- scale = scale_tot[0][0];
- for(i=1; i<DIM; i++)
+ scale = scale_tot[0][0];
+ for (i = 1; i < DIM; i++)
{
/* With anisotropic scaling, the original spherical ns volumes become
* ellipsoids. To avoid costly transformations we use the minimum
* eigenvalue of the scaling matrix for determining the buffer size.
* Since the lower half is 0, the eigenvalues are the diagonal elements.
*/
- scale = min(scale,scale_tot[i][i]);
+ scale = min(scale, scale_tot[i][i]);
if (scale_tot[i][i] != scale_tot[i-1][i-1])
{
bIsotropic = FALSE;
}
- for(j=0; j<i; j++)
+ for (j = 0; j < i; j++)
{
if (scale_tot[i][j] != 0)
{
hbuf2 = sqr(0.5*(scale*ir->rlist - rint));
if (bIsotropic)
{
- for(cg=cg0; cg<cg1; cg++)
+ for (cg = cg0; cg < cg1; cg++)
{
- svmul(scale,cg_cm[cg],cgsc);
+ svmul(scale, cg_cm[cg], cgsc);
a0 = cgs->index[cg];
a1 = cgs->index[cg+1];
- for(a=a0; a<a1; a++)
+ for (a = a0; a < a1; a++)
{
- if (distance2(cgsc,x[a]) > hbuf2)
- {
+ if (distance2(cgsc, x[a]) > hbuf2)
+ {
nBeyond++;
}
}
else
{
/* Anistropic scaling */
- for(cg=cg0; cg<cg1; cg++)
+ for (cg = cg0; cg < cg1; cg++)
{
/* Since scale_tot contains the transpose of the scaling matrix,
* we need to multiply with the transpose.
*/
- tmvmul_ur0(scale_tot,cg_cm[cg],cgsc);
+ tmvmul_ur0(scale_tot, cg_cm[cg], cgsc);
a0 = cgs->index[cg];
a1 = cgs->index[cg+1];
- for(a=a0; a<a1; a++)
+ for (a = a0; a < a1; a++)
{
- if (distance2(cgsc,x[a]) > hbuf2)
+ if (distance2(cgsc, x[a]) > hbuf2)
{
nBeyond++;
}
}
}
}
-
+
return nBeyond;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff