This also required removing the generic group kernel code.
Removed domdec group scheme sorting code.
Also removed the nbnxm grid size call from the domdec module off it
turned out the results were not used.
Change-Id: I1a49967bbeca33dcd62834e8a244aaf02e80a95e
domdec -> mdlib
domdec -> pulling
fileio -> gmxlib
-gmxlib -> listed_forces
mdlib -> essentialdynamics
mdlib -> imd
mdlib -> ewald
}
}
-//! Comparator for sorting charge groups.
-static bool comp_cgsort(const gmx_cgsort_t &a, const gmx_cgsort_t &b)
-{
-
- if (a.nsc == b.nsc)
- {
- return a.ind_gl < b.ind_gl;
- }
- return a.nsc < b.nsc;
-}
-
/*! \brief Order data in \p dataToSort according to \p sort
*
* Note: both buffers should have at least \p sort.size() elements.
}
}
-//! Returns whether a < b */
-static bool compareCgsort(const gmx_cgsort_t &a,
- const gmx_cgsort_t &b)
-{
- return (a.nsc < b.nsc ||
- (a.nsc == b.nsc && a.ind_gl < b.ind_gl));
-}
-
-//! Do sorting of charge groups.
-static void orderedSort(gmx::ArrayRef<const gmx_cgsort_t> stationary,
- gmx::ArrayRef<gmx_cgsort_t> moved,
- std::vector<gmx_cgsort_t> *sort1)
-{
- /* The new indices are not very ordered, so we qsort them */
- std::sort(moved.begin(), moved.end(), comp_cgsort);
-
- /* stationary is already ordered, so now we can merge the two arrays */
- sort1->resize(stationary.size() + moved.size());
- std::merge(stationary.begin(), stationary.end(),
- moved.begin(), moved.end(),
- sort1->begin(),
- compareCgsort);
-}
-
-/*! \brief Set the sorting order for systems with charge groups, returned in sort->sort.
- *
- * The order is according to the global charge group index.
- * This adds and removes charge groups that moved between domains.
- */
-static void dd_sort_order(const gmx_domdec_t *dd,
- const t_forcerec *fr,
- int ncg_home_old,
- gmx_domdec_sort_t *sort)
-{
- const int *a = fr->ns->grid->cell_index;
-
- const int movedValue = NSGRID_SIGNAL_MOVED_FAC*fr->ns->grid->ncells;
-
- if (ncg_home_old >= 0 && !sort->sorted.empty())
- {
- GMX_RELEASE_ASSERT(sort->sorted.size() == static_cast<size_t>(ncg_home_old),
- "The sorting buffer should contain the old home charge group indices");
-
- std::vector<gmx_cgsort_t> &stationary = sort->stationary;
- std::vector<gmx_cgsort_t> &moved = sort->moved;
-
- /* The charge groups that remained in the same ns grid cell
- * are completely ordered. So we can sort efficiently by sorting
- * the charge groups that did move into the stationary list.
- * Note: push_back() seems to be slightly slower than direct access.
- */
- stationary.clear();
- moved.clear();
- for (int i = 0; i < dd->ncg_home; i++)
- {
- /* Check if this cg did not move to another node */
- if (a[i] < movedValue)
- {
- gmx_cgsort_t entry;
- entry.nsc = a[i];
- entry.ind_gl = dd->globalAtomGroupIndices[i];
- entry.ind = i;
-
- if (i >= ncg_home_old || a[i] != sort->sorted[i].nsc)
- {
- /* This cg is new on this node or moved ns grid cell */
- moved.push_back(entry);
- }
- else
- {
- /* This cg did not move */
- stationary.push_back(entry);
- }
- }
- }
-
- if (debug)
- {
- fprintf(debug, "ordered sort cgs: stationary %zu moved %zu\n",
- stationary.size(), moved.size());
- }
- /* Sort efficiently */
- orderedSort(stationary, moved, &sort->sorted);
- }
- else
- {
- std::vector<gmx_cgsort_t> &cgsort = sort->sorted;
- cgsort.clear();
- cgsort.reserve(dd->ncg_home);
- int numCGNew = 0;
- for (int i = 0; i < dd->ncg_home; i++)
- {
- /* Sort on the ns grid cell indices
- * and the global topology index
- */
- gmx_cgsort_t entry;
- entry.nsc = a[i];
- entry.ind_gl = dd->globalAtomGroupIndices[i];
- entry.ind = i;
- cgsort.push_back(entry);
- if (cgsort[i].nsc < movedValue)
- {
- numCGNew++;
- }
- }
- if (debug)
- {
- fprintf(debug, "qsort cgs: %d new home %d\n", dd->ncg_home, numCGNew);
- }
- /* Determine the order of the charge groups using qsort */
- std::sort(cgsort.begin(), cgsort.end(), comp_cgsort);
-
- /* Remove the charge groups which are no longer at home here */
- cgsort.resize(numCGNew);
- }
-}
-
//! Returns the sorting order for atoms based on the nbnxn grid order in sort
static void dd_sort_order_nbnxn(const t_forcerec *fr,
std::vector<gmx_cgsort_t> *sort)
}
//! Returns the sorting state for DD.
-static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state,
- int ncg_home_old)
+static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state *state)
{
gmx_domdec_sort_t *sort = dd->comm->sort.get();
- switch (fr->cutoff_scheme)
- {
- case ecutsGROUP:
- dd_sort_order(dd, fr, ncg_home_old, sort);
- break;
- case ecutsVERLET:
- dd_sort_order_nbnxn(fr, &sort->sorted);
- break;
- default:
- gmx_incons("unimplemented");
- }
+ dd_sort_order_nbnxn(fr, &sort->sorted);
const gmx::RangePartitioning &atomGrouping = dd->atomGrouping();
/* Set the home atom number */
dd->comm->atomRanges.setEnd(DDAtomRanges::Type::Home, dd->atomGrouping().fullRange().end());
- if (fr->cutoff_scheme == ecutsVERLET)
- {
- /* The atoms are now exactly in grid order, update the grid order */
- fr->nbv->setLocalAtomOrder();
- }
- else
- {
- /* Copy the sorted ns cell indices back to the ns grid struct */
- for (gmx::index i = 0; i < cgsort.ssize(); i++)
- {
- fr->ns->grid->cell_index[i] = cgsort[i].nsc;
- }
- fr->ns->grid->nr = cgsort.size();
- }
+ /* The atoms are now exactly in grid order, update the grid order */
+ fr->nbv->setLocalAtomOrder();
}
//! Accumulates load statistics.
t_block *cgs_gl;
int64_t step_pcoupl;
rvec cell_ns_x0, cell_ns_x1;
- int ncgindex_set, ncg_home_old = -1, ncg_moved, nat_f_novirsum;
+ int ncgindex_set, ncg_moved, nat_f_novirsum;
gmx_bool bBoxChanged, bNStGlobalComm, bDoDLB, bCheckWhetherToTurnDlbOn, bLogLoad;
- gmx_bool bRedist, bResortAll;
- ivec ncells_old = {0, 0, 0}, ncells_new = {0, 0, 0}, np;
+ gmx_bool bRedist;
+ ivec np;
real grid_density;
char sbuf[22];
/* Check if we should sort the charge groups */
const bool bSortCG = (bMasterState || bRedist);
- ncg_home_old = dd->ncg_home;
-
/* When repartitioning we mark atom groups that will move to neighboring
* DD cells, but we do not move them right away for performance reasons.
* Thus we need to keep track of how many charge groups will move for
comm_dd_ns_cell_sizes(dd, &ddbox, cell_ns_x0, cell_ns_x1, step);
}
- switch (fr->cutoff_scheme)
- {
- case ecutsGROUP:
- copy_ivec(fr->ns->grid->n, ncells_old);
- grid_first(fplog, fr->ns->grid, dd, &ddbox,
- state_local->box, cell_ns_x0, cell_ns_x1,
- fr->rlist, grid_density);
- break;
- case ecutsVERLET:
- fr->nbv->getLocalNumCells(&ncells_old[XX], &ncells_old[YY]);
- break;
- default:
- gmx_incons("unimplemented");
- }
/* We need to store tric_dir for dd_get_ns_ranges called from ns.c */
copy_ivec(ddbox.tric_dir, comm->tric_dir);
*/
set_zones_ncg_home(dd);
- switch (fr->cutoff_scheme)
- {
- case ecutsVERLET:
- set_zones_size(dd, state_local->box, &ddbox, 0, 1, ncg_moved);
-
- nbnxn_put_on_grid(fr->nbv.get(), state_local->box,
- 0,
- comm->zones.size[0].bb_x0,
- comm->zones.size[0].bb_x1,
- comm->updateGroupsCog.get(),
- 0, dd->ncg_home,
- comm->zones.dens_zone0,
- fr->cginfo,
- state_local->x,
- ncg_moved, bRedist ? comm->movedBuffer.data() : nullptr);
-
- fr->nbv->getLocalNumCells(&ncells_new[XX], &ncells_new[YY]);
- break;
- case ecutsGROUP:
- fill_grid(&comm->zones, fr->ns->grid, dd->ncg_home,
- 0, dd->ncg_home, fr->cg_cm);
+ set_zones_size(dd, state_local->box, &ddbox, 0, 1, ncg_moved);
- copy_ivec(fr->ns->grid->n, ncells_new);
- break;
- default:
- gmx_incons("unimplemented");
- }
-
- bResortAll = bMasterState;
-
- /* Check if we can user the old order and ns grid cell indices
- * of the charge groups to sort the charge groups efficiently.
- */
- if (ncells_new[XX] != ncells_old[XX] ||
- ncells_new[YY] != ncells_old[YY] ||
- ncells_new[ZZ] != ncells_old[ZZ])
- {
- bResortAll = TRUE;
- }
+ nbnxn_put_on_grid(fr->nbv.get(), state_local->box,
+ 0,
+ comm->zones.size[0].bb_x0,
+ comm->zones.size[0].bb_x1,
+ comm->updateGroupsCog.get(),
+ 0, dd->ncg_home,
+ comm->zones.dens_zone0,
+ fr->cginfo,
+ state_local->x,
+ ncg_moved, bRedist ? comm->movedBuffer.data() : nullptr);
if (debug)
{
fprintf(debug, "Step %s, sorting the %d home charge groups\n",
gmx_step_str(step, sbuf), dd->ncg_home);
}
- dd_sort_state(dd, fr->cg_cm, fr, state_local,
- bResortAll ? -1 : ncg_home_old);
+ dd_sort_state(dd, fr->cg_cm, fr, state_local);
/* After sorting and compacting we set the correct size */
dd_resize_state(state_local, f, comm->atomRanges.numHomeAtoms());
comm);
}
- int *moved;
- if (fr->cutoff_scheme == ecutsVERLET)
- {
- moved = getMovedBuffer(comm, 0, dd->ncg_home);
- }
- else
- {
- moved = fr->ns->grid->cell_index;
- }
+ int *moved = getMovedBuffer(comm, 0, dd->ncg_home);
clear_and_mark_ind(move,
dd->globalAtomGroupIndices, dd->atomGrouping(), dd->globalAtomIndices,
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012,2014,2015,2017,2018, 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.
- *
- * 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 research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "nb_generic.h"
-
-#include <cmath>
-
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/gmxlib/nonbonded/nb_kernel.h"
-#include "gromacs/gmxlib/nonbonded/nonbonded.h"
-#include "gromacs/math/functions.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/utility/fatalerror.h"
-
-void
-gmx_nb_generic_kernel(t_nblist * nlist,
- rvec * xx,
- rvec * ff,
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- nb_kernel_data_t * kernel_data,
- t_nrnb * nrnb)
-{
- int ntype, table_nelements, ielec, ivdw;
- real facel;
- int n, ii, is3, ii3, k, nj0, nj1, jnr, j3, ggid, nnn, n0;
- real shX, shY, shZ;
- real fscal, felec, fvdw, velec, vvdw, tx, ty, tz;
- real rinvsq;
- real iq;
- real qq, vctot;
- int nti, nvdwparam;
- int tj;
- real rt, r, eps, eps2, Y, F, Geps, Heps2, VV, FF, Fp, fijD, fijR;
- real rinvsix;
- real vvdwtot;
- real vvdw_rep, vvdw_disp;
- real ix, iy, iz, fix, fiy, fiz;
- real jx, jy, jz;
- real dx, dy, dz, rsq, rinv;
- real c6, c12, c6grid, cexp1, cexp2, br;
- real * charge;
- real * shiftvec;
- real * vdwparam, *vdwgridparam;
- int * type;
- real * fshift;
- real * velecgrp;
- real * vvdwgrp;
- real tabscale;
- real * VFtab;
- real * x;
- real * f;
- int ewitab;
- real ewtabscale, eweps, ewrt, ewtabhalfspace;
- real * ewtab;
- real rcoulomb2, rvdw, rvdw2, sh_dispersion, sh_repulsion;
- real rcutoff, rcutoff2;
- real d, d2, sw, dsw, rinvcorr;
- real elec_swV3, elec_swV4, elec_swV5, elec_swF2, elec_swF3, elec_swF4;
- real vdw_swV3, vdw_swV4, vdw_swV5, vdw_swF2, vdw_swF3, vdw_swF4;
- real ewclj, ewclj2, ewclj6, ewcljrsq, poly, exponent, sh_lj_ewald;
- gmx_bool bExactElecCutoff, bExactVdwCutoff, bExactCutoff;
- gmx_bool do_tab;
-
- x = xx[0];
- f = ff[0];
- ielec = nlist->ielec;
- ivdw = nlist->ivdw;
-
- fshift = fr->fshift[0];
- velecgrp = kernel_data->energygrp_elec;
- vvdwgrp = kernel_data->energygrp_vdw;
-
- do_tab = (ielec == GMX_NBKERNEL_ELEC_CUBICSPLINETABLE ||
- ivdw == GMX_NBKERNEL_VDW_CUBICSPLINETABLE);
- if (do_tab)
- {
- tabscale = kernel_data->table_elec_vdw->scale;
- VFtab = kernel_data->table_elec_vdw->data;
- }
- else
- {
- tabscale = 0;
- VFtab = nullptr;
- }
-
- const interaction_const_t *ic = fr->ic;
-
- if (ielec == GMX_NBKERNEL_ELEC_EWALD)
- {
- ewtab = ic->tabq_coul_FDV0;
- ewtabscale = ic->tabq_scale;
- ewtabhalfspace = 0.5/ewtabscale;
- }
- else
- {
- ewtab = nullptr;
- ewtabhalfspace = ewtabscale = 0;
- }
-
- rcoulomb2 = ic->rcoulomb*ic->rcoulomb;
- rvdw = ic->rvdw;
- rvdw2 = rvdw*rvdw;
- sh_dispersion = ic->dispersion_shift.cpot;
- sh_repulsion = ic->repulsion_shift.cpot;
- sh_lj_ewald = ic->sh_lj_ewald;
-
- ewclj = ic->ewaldcoeff_lj;
- ewclj2 = ewclj*ewclj;
- ewclj6 = ewclj2*ewclj2*ewclj2;
-
- if (ic->coulomb_modifier == eintmodPOTSWITCH)
- {
- d = ic->rcoulomb - ic->rcoulomb_switch;
- elec_swV3 = -10.0/(d*d*d);
- elec_swV4 = 15.0/(d*d*d*d);
- elec_swV5 = -6.0/(d*d*d*d*d);
- elec_swF2 = -30.0/(d*d*d);
- elec_swF3 = 60.0/(d*d*d*d);
- elec_swF4 = -30.0/(d*d*d*d*d);
- }
- else
- {
- /* Avoid warnings from stupid compilers (looking at you, Clang!) */
- elec_swV3 = elec_swV4 = elec_swV5 = elec_swF2 = elec_swF3 = elec_swF4 = 0.0;
- }
- if (ic->vdw_modifier == eintmodPOTSWITCH)
- {
- d = ic->rvdw - ic->rvdw_switch;
- vdw_swV3 = -10.0/(d*d*d);
- vdw_swV4 = 15.0/(d*d*d*d);
- vdw_swV5 = -6.0/(d*d*d*d*d);
- vdw_swF2 = -30.0/(d*d*d);
- vdw_swF3 = 60.0/(d*d*d*d);
- vdw_swF4 = -30.0/(d*d*d*d*d);
- }
- else
- {
- /* Avoid warnings from stupid compilers (looking at you, Clang!) */
- vdw_swV3 = vdw_swV4 = vdw_swV5 = vdw_swF2 = vdw_swF3 = vdw_swF4 = 0.0;
- }
-
- bExactElecCutoff = (ic->coulomb_modifier != eintmodNONE) || ic->eeltype == eelRF_ZERO;
- bExactVdwCutoff = (ic->vdw_modifier != eintmodNONE);
- bExactCutoff = bExactElecCutoff && bExactVdwCutoff;
-
- if (bExactCutoff)
- {
- rcutoff = ( ic->rcoulomb > ic->rvdw ) ? ic->rcoulomb : ic->rvdw;
- rcutoff2 = rcutoff*rcutoff;
- }
- else
- {
- /* Fix warnings for stupid compilers */
- rcutoff2 = 1e30;
- }
-
- /* avoid compiler warnings for cases that cannot happen */
- nnn = 0;
- eps = 0.0;
- eps2 = 0.0;
-
- /* 3 VdW parameters for Buckingham, otherwise 2 */
- nvdwparam = (ivdw == GMX_NBKERNEL_VDW_BUCKINGHAM) ? 3 : 2;
- table_nelements = 12;
-
- charge = mdatoms->chargeA;
- type = mdatoms->typeA;
- facel = fr->ic->epsfac;
- shiftvec = fr->shift_vec[0];
- vdwparam = fr->nbfp;
- ntype = fr->ntype;
- vdwgridparam = fr->ljpme_c6grid;
-
- for (n = 0; (n < nlist->nri); n++)
- {
- is3 = 3*nlist->shift[n];
- shX = shiftvec[is3];
- shY = shiftvec[is3+1];
- shZ = shiftvec[is3+2];
- nj0 = nlist->jindex[n];
- nj1 = nlist->jindex[n+1];
- ii = nlist->iinr[n];
- ii3 = 3*ii;
- ix = shX + x[ii3+0];
- iy = shY + x[ii3+1];
- iz = shZ + x[ii3+2];
- iq = facel*charge[ii];
- nti = nvdwparam*ntype*type[ii];
- vctot = 0;
- vvdwtot = 0;
- fix = 0;
- fiy = 0;
- fiz = 0;
-
- for (k = nj0; (k < nj1); k++)
- {
- jnr = nlist->jjnr[k];
- j3 = 3*jnr;
- jx = x[j3+0];
- jy = x[j3+1];
- jz = x[j3+2];
- dx = ix - jx;
- dy = iy - jy;
- dz = iz - jz;
- rsq = dx*dx+dy*dy+dz*dz;
- rinv = gmx::invsqrt(rsq);
- rinvsq = rinv*rinv;
- felec = 0;
- fvdw = 0;
- velec = 0;
- vvdw = 0;
-
- if (bExactCutoff && rsq >= rcutoff2)
- {
- continue;
- }
-
- if (ielec == GMX_NBKERNEL_ELEC_CUBICSPLINETABLE || ivdw == GMX_NBKERNEL_VDW_CUBICSPLINETABLE)
- {
- r = rsq*rinv;
- rt = r*tabscale;
- n0 = rt;
- eps = rt-n0;
- eps2 = eps*eps;
- nnn = table_nelements*n0;
- }
-
- /* Coulomb interaction. ielec==0 means no interaction */
- if (ielec != GMX_NBKERNEL_ELEC_NONE)
- {
- qq = iq*charge[jnr];
-
- switch (ielec)
- {
- case GMX_NBKERNEL_ELEC_NONE:
- break;
-
- case GMX_NBKERNEL_ELEC_COULOMB:
- /* Vanilla cutoff coulomb */
- velec = qq*rinv;
- felec = velec*rinvsq;
- /* The shift for the Coulomb potential is stored in
- * the RF parameter c_rf, which is 0 without shift
- */
- velec -= qq*ic->c_rf;
- break;
-
- case GMX_NBKERNEL_ELEC_REACTIONFIELD:
- /* Reaction-field */
- velec = qq*(rinv + ic->k_rf*rsq-ic->c_rf);
- felec = qq*(rinv*rinvsq - 2.0*ic->k_rf);
- break;
-
- case GMX_NBKERNEL_ELEC_CUBICSPLINETABLE:
- /* Tabulated coulomb */
- Y = VFtab[nnn];
- F = VFtab[nnn+1];
- Geps = eps*VFtab[nnn+2];
- Heps2 = eps2*VFtab[nnn+3];
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+2.0*Heps2;
- velec = qq*VV;
- felec = -qq*FF*tabscale*rinv;
- break;
-
- case GMX_NBKERNEL_ELEC_EWALD:
- ewrt = rsq*rinv*ewtabscale;
- ewitab = ewrt;
- eweps = ewrt-ewitab;
- ewitab = 4*ewitab;
- felec = ewtab[ewitab]+eweps*ewtab[ewitab+1];
- rinvcorr = (ic->coulomb_modifier == eintmodPOTSHIFT) ? rinv - ic->sh_ewald : rinv;
- velec = qq*(rinvcorr-(ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+felec)));
- felec = qq*rinv*(rinvsq-felec);
- break;
-
- default:
- gmx_fatal(FARGS, "Death & horror! No generic coulomb interaction for ielec=%d.\n", ielec);
- }
- if (ic->coulomb_modifier == eintmodPOTSWITCH)
- {
- d = rsq*rinv - ic->rcoulomb_switch;
- d = (d > 0.0) ? d : 0.0;
- d2 = d*d;
- sw = 1.0+d2*d*(elec_swV3+d*(elec_swV4+d*elec_swV5));
- dsw = d2*(elec_swF2+d*(elec_swF3+d*elec_swF4));
- /* Apply switch function. Note that felec=f/r since it will be multiplied
- * by the i-j displacement vector. This means felec'=f'/r=-(v*sw)'/r=
- * -(v'*sw+v*dsw)/r=-v'*sw/r-v*dsw/r=felec*sw-v*dsw/r
- */
- felec = felec*sw - rinv*velec*dsw;
- /* Once we have used velec to update felec we can modify velec too */
- velec *= sw;
- }
- if (bExactElecCutoff)
- {
- felec = (rsq < rcoulomb2) ? felec : 0.0;
- velec = (rsq < rcoulomb2) ? velec : 0.0;
- }
- vctot += velec;
- } /* End of coulomb interactions */
-
-
- /* VdW interaction. ivdw==0 means no interaction */
- if (ivdw != GMX_NBKERNEL_VDW_NONE)
- {
- tj = nti+nvdwparam*type[jnr];
-
- switch (ivdw)
- {
- case GMX_NBKERNEL_VDW_NONE:
- break;
-
- case GMX_NBKERNEL_VDW_LENNARDJONES:
- /* Vanilla Lennard-Jones cutoff */
- c6 = vdwparam[tj];
- c12 = vdwparam[tj+1];
- rinvsix = rinvsq*rinvsq*rinvsq;
- vvdw_disp = c6*rinvsix;
- vvdw_rep = c12*rinvsix*rinvsix;
- fvdw = (vvdw_rep-vvdw_disp)*rinvsq;
- if (ic->vdw_modifier == eintmodPOTSHIFT)
- {
- vvdw = (vvdw_rep + c12*sh_repulsion)/12.0 - (vvdw_disp + c6*sh_dispersion)/6.0;
- }
- else
- {
- vvdw = vvdw_rep/12.0-vvdw_disp/6.0;
- }
- break;
-
- case GMX_NBKERNEL_VDW_BUCKINGHAM:
- /* Buckingham */
- c6 = vdwparam[tj];
- cexp1 = vdwparam[tj+1];
- cexp2 = vdwparam[tj+2];
-
- rinvsix = rinvsq*rinvsq*rinvsq;
- vvdw_disp = c6*rinvsix;
- br = cexp2*rsq*rinv;
- vvdw_rep = cexp1*std::exp(-br);
- fvdw = (br*vvdw_rep-vvdw_disp)*rinvsq;
- if (ic->vdw_modifier == eintmodPOTSHIFT)
- {
- vvdw = (vvdw_rep-cexp1*std::exp(-cexp2*rvdw))-(vvdw_disp + c6*sh_dispersion)/6.0;
- }
- else
- {
- vvdw = vvdw_rep-vvdw_disp/6.0;
- }
- break;
-
- case GMX_NBKERNEL_VDW_CUBICSPLINETABLE:
- /* Tabulated VdW */
- c6 = vdwparam[tj];
- c12 = vdwparam[tj+1];
- Y = VFtab[nnn+4];
- F = VFtab[nnn+5];
- Geps = eps*VFtab[nnn+6];
- Heps2 = eps2*VFtab[nnn+7];
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+2.0*Heps2;
- vvdw_disp = c6*VV;
- fijD = c6*FF;
- Y = VFtab[nnn+8];
- F = VFtab[nnn+9];
- Geps = eps*VFtab[nnn+10];
- Heps2 = eps2*VFtab[nnn+11];
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+2.0*Heps2;
- vvdw_rep = c12*VV;
- fijR = c12*FF;
- fvdw = -(fijD+fijR)*tabscale*rinv;
- vvdw = vvdw_disp + vvdw_rep;
- break;
-
-
- case GMX_NBKERNEL_VDW_LJEWALD:
- /* LJ-PME */
- rinvsix = rinvsq*rinvsq*rinvsq;
- ewcljrsq = ewclj2*rsq;
- exponent = std::exp(-ewcljrsq);
- poly = exponent*(1.0 + ewcljrsq + ewcljrsq*ewcljrsq*0.5);
- c6 = vdwparam[tj];
- c12 = vdwparam[tj+1];
- c6grid = vdwgridparam[tj];
- vvdw_disp = (c6-c6grid*(1.0-poly))*rinvsix;
- vvdw_rep = c12*rinvsix*rinvsix;
- fvdw = (vvdw_rep - vvdw_disp - c6grid*(1.0/6.0)*exponent*ewclj6)*rinvsq;
- if (ic->vdw_modifier == eintmodPOTSHIFT)
- {
- vvdw = (vvdw_rep + c12*sh_repulsion)/12.0 - (vvdw_disp + c6*sh_dispersion - c6grid*sh_lj_ewald)/6.0;
- }
- else
- {
- vvdw = vvdw_rep/12.0-vvdw_disp/6.0;
- }
- break;
-
- default:
- gmx_fatal(FARGS, "Death & horror! No generic VdW interaction for ivdw=%d.\n", ivdw);
- }
- if (ic->vdw_modifier == eintmodPOTSWITCH)
- {
- d = rsq*rinv - ic->rvdw_switch;
- d = (d > 0.0) ? d : 0.0;
- d2 = d*d;
- sw = 1.0+d2*d*(vdw_swV3+d*(vdw_swV4+d*vdw_swV5));
- dsw = d2*(vdw_swF2+d*(vdw_swF3+d*vdw_swF4));
- /* See coulomb interaction for the force-switch formula */
- fvdw = fvdw*sw - rinv*vvdw*dsw;
- vvdw *= sw;
- }
- if (bExactVdwCutoff)
- {
- fvdw = (rsq < rvdw2) ? fvdw : 0.0;
- vvdw = (rsq < rvdw2) ? vvdw : 0.0;
- }
- vvdwtot += vvdw;
- } /* end VdW interactions */
-
- fscal = felec+fvdw;
-
- tx = fscal*dx;
- ty = fscal*dy;
- tz = fscal*dz;
- fix = fix + tx;
- fiy = fiy + ty;
- fiz = fiz + tz;
- f[j3+0] = f[j3+0] - tx;
- f[j3+1] = f[j3+1] - ty;
- f[j3+2] = f[j3+2] - tz;
- }
-
- f[ii3+0] = f[ii3+0] + fix;
- f[ii3+1] = f[ii3+1] + fiy;
- f[ii3+2] = f[ii3+2] + fiz;
- fshift[is3] = fshift[is3]+fix;
- fshift[is3+1] = fshift[is3+1]+fiy;
- fshift[is3+2] = fshift[is3+2]+fiz;
- ggid = nlist->gid[n];
- velecgrp[ggid] += vctot;
- vvdwgrp[ggid] += vvdwtot;
- }
- /* Estimate flops, average for generic kernel:
- * 12 flops per outer iteration
- * 50 flops per inner iteration
- */
- inc_nrnb(nrnb, eNR_NBKERNEL_GENERIC, nlist->nri*12 + nlist->jindex[n]*50);
-}
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2012,2014,2015, 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.
- *
- * 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 research papers on the package. Check out http://www.gromacs.org.
- */
-
-#ifndef _nb_generic_h_
-#define _nb_generic_h_
-
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/gmxlib/nonbonded/nb_kernel.h"
-#include "gromacs/math/vectypes.h"
-#include "gromacs/mdtypes/mdatom.h"
-#include "gromacs/mdtypes/nblist.h"
-
-void
-gmx_nb_generic_kernel(t_nblist * nlist,
- rvec * x,
- rvec * f,
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- nb_kernel_data_t * kernel_data,
- t_nrnb * nrnb);
-
-
-#endif
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018, 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.
- *
- * 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 research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "nb_generic_cg.h"
-
-#include <cmath>
-
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/gmxlib/nonbonded/nb_kernel.h"
-#include "gromacs/gmxlib/nonbonded/nonbonded.h"
-#include "gromacs/math/functions.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/utility/fatalerror.h"
-
-void
-gmx_nb_generic_cg_kernel(t_nblist * nlist,
- rvec * xx,
- rvec * ff,
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- nb_kernel_data_t * kernel_data,
- t_nrnb * nrnb)
-{
- int ntype, table_nelements, ielec, ivdw;
- real facel;
- int n, is3, i3, k, nj0, nj1, j3, ggid, nnn, n0;
- int ai0, ai1, ai, aj0, aj1, aj;
- real shX, shY, shZ;
- real fscal, tx, ty, tz;
- real rinvsq;
- real iq;
- real qq, vcoul, krsq, vctot;
- int nti, nvdwparam;
- int tj;
- real rt, r, eps, eps2, Y, F, Geps, Heps2, VV, FF, Fp, fijD, fijR;
- real rinvsix;
- real Vvdwtot;
- real Vvdw_rep, Vvdw_disp;
- real ix, iy, iz, fix, fiy, fiz;
- real jx, jy, jz;
- real dx, dy, dz, rsq, rinv;
- real c6, c12, cexp1, cexp2, br;
- real * charge;
- real * shiftvec;
- real * vdwparam;
- int * type;
- t_excl * excl;
- real * fshift;
- real * Vc;
- real * Vvdw;
- real tabscale;
- real * VFtab;
- real * x;
- real * f;
- gmx_bool do_tab;
-
- x = xx[0];
- f = ff[0];
- ielec = nlist->ielec;
- ivdw = nlist->ivdw;
-
- fshift = fr->fshift[0];
- Vc = kernel_data->energygrp_elec;
- Vvdw = kernel_data->energygrp_vdw;
-
- do_tab = (ielec == GMX_NBKERNEL_ELEC_CUBICSPLINETABLE ||
- ivdw == GMX_NBKERNEL_VDW_CUBICSPLINETABLE);
- if (do_tab)
- {
- tabscale = kernel_data->table_elec_vdw->scale;
- VFtab = kernel_data->table_elec_vdw->data;
- }
- else
- {
- tabscale = 0;
- VFtab = nullptr;
- }
-
- /* avoid compiler warnings for cases that cannot happen */
- nnn = 0;
- eps = 0.0;
- eps2 = 0.0;
-
- /* 3 VdW parameters for buckingham, otherwise 2 */
- nvdwparam = (nlist->ivdw == 2) ? 3 : 2;
- table_nelements = (ielec == 3) ? 4 : 0;
- table_nelements += (ivdw == 3) ? 8 : 0;
-
- charge = mdatoms->chargeA;
- type = mdatoms->typeA;
- facel = fr->ic->epsfac;
- shiftvec = fr->shift_vec[0];
- vdwparam = fr->nbfp;
- ntype = fr->ntype;
-
- for (n = 0; (n < nlist->nri); n++)
- {
- is3 = 3*nlist->shift[n];
- shX = shiftvec[is3];
- shY = shiftvec[is3+1];
- shZ = shiftvec[is3+2];
- nj0 = nlist->jindex[n];
- nj1 = nlist->jindex[n+1];
- ai0 = nlist->iinr[n];
- ai1 = nlist->iinr_end[n];
- vctot = 0;
- Vvdwtot = 0;
- fix = 0;
- fiy = 0;
- fiz = 0;
-
- for (k = nj0; (k < nj1); k++)
- {
- aj0 = nlist->jjnr[k];
- aj1 = nlist->jjnr_end[k];
- excl = &nlist->excl[k*MAX_CGCGSIZE];
-
- for (ai = ai0; (ai < ai1); ai++)
- {
- i3 = ai*3;
- ix = shX + x[i3+0];
- iy = shY + x[i3+1];
- iz = shZ + x[i3+2];
- iq = facel*charge[ai];
- nti = nvdwparam*ntype*type[ai];
-
- /* Note that this code currently calculates
- * all LJ and Coulomb interactions,
- * even if the LJ parameters or charges are zero.
- * If required, this can be optimized.
- */
-
- for (aj = aj0; (aj < aj1); aj++)
- {
- /* Check if this interaction is excluded */
- if (excl[aj-aj0] & (1<<(ai-ai0)))
- {
- continue;
- }
-
- j3 = aj*3;
- jx = x[j3+0];
- jy = x[j3+1];
- jz = x[j3+2];
- dx = ix - jx;
- dy = iy - jy;
- dz = iz - jz;
- rsq = dx*dx+dy*dy+dz*dz;
- rinv = gmx::invsqrt(rsq);
- rinvsq = rinv*rinv;
- fscal = 0;
-
- if (ielec == 3 || ivdw == 3)
- {
- r = rsq*rinv;
- rt = r*tabscale;
- n0 = rt;
- eps = rt-n0;
- eps2 = eps*eps;
- nnn = table_nelements*n0;
- }
-
- /* Coulomb interaction. ielec==0 means no interaction */
- if (ielec > 0)
- {
- qq = iq*charge[aj];
-
- switch (ielec)
- {
- case 1:
- /* Vanilla cutoff coulomb */
- vcoul = qq*rinv;
- fscal = vcoul*rinvsq;
- break;
-
- case 2:
- /* Reaction-field */
- krsq = fr->ic->k_rf*rsq;
- vcoul = qq*(rinv + krsq - fr->ic->c_rf);
- fscal = qq*(rinv-2.0*krsq)*rinvsq;
- break;
-
- case 3:
- /* Tabulated coulomb */
- Y = VFtab[nnn];
- F = VFtab[nnn+1];
- Geps = eps*VFtab[nnn+2];
- Heps2 = eps2*VFtab[nnn+3];
- nnn += 4;
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+2.0*Heps2;
- vcoul = qq*VV;
- fscal = -qq*FF*tabscale*rinv;
- break;
-
- case 4:
- /* GB */
- gmx_fatal(FARGS, "Death & horror! GB generic interaction not implemented.\n");
- default:
- gmx_fatal(FARGS, "Death & horror! No generic coulomb interaction for ielec=%d.\n", ielec);
- }
- vctot = vctot+vcoul;
- } /* End of coulomb interactions */
-
-
- /* VdW interaction. ivdw==0 means no interaction */
- if (ivdw > 0)
- {
- tj = nti+nvdwparam*type[aj];
-
- switch (ivdw)
- {
- case 1:
- /* Vanilla Lennard-Jones cutoff */
- c6 = vdwparam[tj];
- c12 = vdwparam[tj+1];
-
- rinvsix = rinvsq*rinvsq*rinvsq;
- Vvdw_disp = c6*rinvsix;
- Vvdw_rep = c12*rinvsix*rinvsix;
- fscal += (12.0*Vvdw_rep-6.0*Vvdw_disp)*rinvsq;
- Vvdwtot = Vvdwtot+Vvdw_rep-Vvdw_disp;
- break;
-
- case 2:
- /* Buckingham */
- c6 = vdwparam[tj];
- cexp1 = vdwparam[tj+1];
- cexp2 = vdwparam[tj+2];
-
- rinvsix = rinvsq*rinvsq*rinvsq;
- Vvdw_disp = c6*rinvsix;
- br = cexp2*rsq*rinv;
- Vvdw_rep = cexp1*std::exp(-br);
- fscal += (br*Vvdw_rep-6.0*Vvdw_disp)*rinvsq;
- Vvdwtot = Vvdwtot+Vvdw_rep-Vvdw_disp;
- break;
-
- case 3:
- /* Tabulated VdW */
- c6 = vdwparam[tj];
- c12 = vdwparam[tj+1];
-
- Y = VFtab[nnn];
- F = VFtab[nnn+1];
- Geps = eps*VFtab[nnn+2];
- Heps2 = eps2*VFtab[nnn+3];
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+2.0*Heps2;
- Vvdw_disp = c6*VV;
- fijD = c6*FF;
- nnn += 4;
- Y = VFtab[nnn];
- F = VFtab[nnn+1];
- Geps = eps*VFtab[nnn+2];
- Heps2 = eps2*VFtab[nnn+3];
- Fp = F+Geps+Heps2;
- VV = Y+eps*Fp;
- FF = Fp+Geps+2.0*Heps2;
- Vvdw_rep = c12*VV;
- fijR = c12*FF;
- fscal += -(fijD+fijR)*tabscale*rinv;
- Vvdwtot = Vvdwtot + Vvdw_disp + Vvdw_rep;
- break;
-
- default:
- gmx_fatal(FARGS, "Death & horror! No generic VdW interaction for ivdw=%d.\n", ivdw);
- }
- } /* end VdW interactions */
-
-
- tx = fscal*dx;
- ty = fscal*dy;
- tz = fscal*dz;
- f[i3+0] += tx;
- f[i3+1] += ty;
- f[i3+2] += tz;
- f[j3+0] -= tx;
- f[j3+1] -= ty;
- f[j3+2] -= tz;
- fix += tx;
- fiy += ty;
- fiz += tz;
- }
- }
- }
-
- fshift[is3] += fix;
- fshift[is3+1] += fiy;
- fshift[is3+2] += fiz;
- ggid = nlist->gid[n];
- Vc[ggid] += vctot;
- Vvdw[ggid] += Vvdwtot;
- }
- /* Estimate flops, average for generic cg kernel:
- * 12 flops per outer iteration
- * 100 flops per inner iteration
- */
- inc_nrnb(nrnb, eNR_NBKERNEL_GENERIC_CG, nlist->nri*12 + nlist->jindex[n]*100);
-}
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2012,2014,2015, 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.
- *
- * 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 research papers on the package. Check out http://www.gromacs.org.
- */
-
-#ifndef _nb_generic_cg_h_
-#define _nb_generic_cg_h_
-
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/gmxlib/nonbonded/nb_kernel.h"
-#include "gromacs/math/vectypes.h"
-#include "gromacs/mdtypes/mdatom.h"
-#include "gromacs/mdtypes/nblist.h"
-
-void
-gmx_nb_generic_cg_kernel(t_nblist * nlist,
- rvec * x,
- rvec * f,
- t_forcerec * fr,
- t_mdatoms * mdatoms,
- nb_kernel_data_t * kernel_data,
- t_nrnb * nrnb);
-
-#endif
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018,2019, 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.
- *
- * 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 research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "nonbonded.h"
-
-#include <cassert>
-#include <cstdio>
-#include <cstdlib>
-
-#include "thread_mpi/threads.h"
-
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/gmxlib/nonbonded/nb_free_energy.h"
-#include "gromacs/gmxlib/nonbonded/nb_generic.h"
-#include "gromacs/gmxlib/nonbonded/nb_generic_cg.h"
-#include "gromacs/gmxlib/nonbonded/nb_kernel.h"
-#include "gromacs/listed_forces/bonded.h"
-#include "gromacs/math/utilities.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/mdtypes/enerdata.h"
-#include "gromacs/mdtypes/forcerec.h"
-#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/mdtypes/mdatom.h"
-#include "gromacs/mdtypes/nblist.h"
-#include "gromacs/pbcutil/ishift.h"
-#include "gromacs/pbcutil/mshift.h"
-#include "gromacs/pbcutil/pbc.h"
-#include "gromacs/tables/forcetable.h"
-#include "gromacs/utility/arraysize.h"
-#include "gromacs/utility/basedefinitions.h"
-#include "gromacs/utility/cstringutil.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/smalloc.h"
-
-void
-gmx_nonbonded_set_kernel_pointers(t_nblist *nl)
-{
- const char * elec;
- const char * elec_mod;
- const char * vdw;
- const char * vdw_mod;
- const char * geom;
-
- nl->kernelptr_vf = nullptr;
- nl->kernelptr_v = nullptr;
- nl->kernelptr_f = nullptr;
-
- elec = gmx_nbkernel_elec_names[nl->ielec];
- elec_mod = eintmod_names[nl->ielecmod];
- vdw = gmx_nbkernel_vdw_names[nl->ivdw];
- vdw_mod = eintmod_names[nl->ivdwmod];
- geom = gmx_nblist_geometry_names[nl->igeometry];
-
- if (nl->type == GMX_NBLIST_INTERACTION_FREE_ENERGY)
- {
- nl->kernelptr_vf = reinterpret_cast<void *>(gmx_nb_free_energy_kernel);
- nl->kernelptr_f = reinterpret_cast<void *>(gmx_nb_free_energy_kernel);
- nl->simd_padding_width = 1;
- }
- else if (!gmx_strcasecmp_min(geom, "CG-CG"))
- {
- nl->kernelptr_vf = reinterpret_cast<void *>(gmx_nb_generic_cg_kernel);
- nl->kernelptr_f = reinterpret_cast<void *>(gmx_nb_generic_cg_kernel);
- nl->simd_padding_width = 1;
- }
- else
- {
- /* "Pick" the only remaining kernel, the generic one.
- * We only do this for particle-particle kernels; by leaving the water-optimized kernel
- * pointers to NULL, the water optimization will automatically be disabled for this interaction.
- */
- if (!gmx_strcasecmp_min(geom, "Particle-Particle"))
- {
- nl->kernelptr_vf = reinterpret_cast<void *>(gmx_nb_generic_kernel);
- nl->kernelptr_f = reinterpret_cast<void *>(gmx_nb_generic_kernel);
- nl->simd_padding_width = 1;
- if (debug)
- {
- fprintf(debug,
- "WARNING - Slow generic NB kernel used for neighborlist with\n"
- " Elec: '%s', Modifier: '%s'\n"
- " Vdw: '%s', Modifier: '%s'\n",
- elec, elec_mod, vdw, vdw_mod);
- }
- }
- }
-}
-
-void do_nonbonded(const t_forcerec *fr,
- rvec x[],
- rvec f_shortrange[],
- const t_mdatoms *mdatoms,
- const t_blocka *excl,
- gmx_grppairener_t *grppener,
- t_nrnb *nrnb,
- real *lambda,
- real *dvdl,
- int nls,
- int eNL,
- int flags)
-{
- t_nblist * nlist;
- int n, n0, n1, i, i0, i1;
- t_nblists * nblists;
- nb_kernel_data_t kernel_data;
- nb_kernel_t * kernelptr = nullptr;
- rvec * f;
-
- kernel_data.flags = flags;
- kernel_data.exclusions = excl;
- kernel_data.lambda = lambda;
- kernel_data.dvdl = dvdl;
-
- if (eNL >= 0)
- {
- i0 = eNL;
- i1 = i0+1;
- }
- else
- {
- i0 = 0;
- i1 = eNL_NR;
- }
-
- if (nls >= 0)
- {
- n0 = nls;
- n1 = nls+1;
- }
- else
- {
- n0 = 0;
- n1 = fr->nnblists;
- }
-
- for (n = n0; (n < n1); n++)
- {
- nblists = &fr->nblists[n];
-
- /* Tabulated kernels hard-code a lot of assumptions about the
- * structure of these tables, but that's not worth fixing with
- * the group scheme due for removal soon. As a token
- * improvement, this assertion will stop code segfaulting if
- * someone assumes that extending the group-scheme table-type
- * enumeration is something that GROMACS supports. */
- static_assert(etiNR == 3, "");
-
- kernel_data.table_elec = nblists->table_elec;
- kernel_data.table_vdw = nblists->table_vdw;
- kernel_data.table_elec_vdw = nblists->table_elec_vdw;
-
- {
- {
- /* Short-range */
- if (!(flags & GMX_NONBONDED_DO_SR))
- {
- continue;
- }
- kernel_data.energygrp_elec = grppener->ener[egCOULSR];
- kernel_data.energygrp_vdw = grppener->ener[fr->bBHAM ? egBHAMSR : egLJSR];
- nlist = nblists->nlist_sr;
- f = f_shortrange;
- }
-
- for (i = i0; (i < i1); i++)
- {
- if (nlist[i].nri > 0)
- {
- if (flags & GMX_NONBONDED_DO_POTENTIAL)
- {
- /* Potential and force */
- kernelptr = reinterpret_cast<nb_kernel_t *>(nlist[i].kernelptr_vf);
- }
- else
- {
- /* Force only, no potential */
- kernelptr = reinterpret_cast<nb_kernel_t *>(nlist[i].kernelptr_f);
- }
-
- if (nlist[i].type != GMX_NBLIST_INTERACTION_FREE_ENERGY && (flags & GMX_NONBONDED_DO_FOREIGNLAMBDA))
- {
- /* We don't need the non-perturbed interactions */
- continue;
- }
- /* Neighborlists whose kernelptr==NULL will always be empty */
- if (kernelptr != nullptr)
- {
- (*kernelptr)(&(nlist[i]), x, f, const_cast<t_forcerec*>(fr),
- const_cast<t_mdatoms*>(mdatoms), &kernel_data, nrnb);
- }
- else
- {
- gmx_fatal(FARGS, "Non-empty neighborlist does not have any kernel pointer assigned.");
- }
- }
- }
- }
- }
-}
#ifndef GMX_GMXLIB_NONBONDED_NONBONDED_H
#define GMX_GMXLIB_NONBONDED_NONBONDED_H
-#include <stdio.h>
-
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/math/vectypes.h"
-#include "gromacs/mdtypes/mdatom.h"
-#include "gromacs/mdtypes/nblist.h"
-#include "gromacs/topology/block.h"
-#include "gromacs/utility/basedefinitions.h"
-
-struct gmx_grppairener_t;
-struct t_forcerec;
-
-void
-gmx_nonbonded_set_kernel_pointers(t_nblist * nl);
-
-
#define GMX_NONBONDED_DO_FORCE (1<<1)
#define GMX_NONBONDED_DO_SHIFTFORCE (1<<2)
#define GMX_NONBONDED_DO_FOREIGNLAMBDA (1<<3)
#define GMX_NONBONDED_DO_POTENTIAL (1<<4)
#define GMX_NONBONDED_DO_SR (1<<5)
-void
-do_nonbonded(const t_forcerec *fr,
- rvec x[],
- rvec f_shortrange[],
- const t_mdatoms *md,
- const t_blocka *excl,
- gmx_grppairener_t *grppener,
- t_nrnb *nrnb,
- real *lambda,
- real dvdlambda[],
- int nls,
- int eNL,
- int flags);
-
#endif
#include "gromacs/mdlib/forcerec_threading.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
#include "gromacs/mdlib/md_support.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdlib/rf_util.h"
#include "gromacs/mdlib/wall.h"
if (!needGroupSchemeTables)
{
bSomeNormalNbListsAreInUse = TRUE;
- fr->nnblists = 1;
}
else
{
}
}
}
- if (bSomeNormalNbListsAreInUse)
- {
- fr->nnblists = negptable + 1;
- }
- else
- {
- fr->nnblists = negptable;
- }
- if (fr->nnblists > 1)
- {
- snew(fr->gid2nblists, ir->opts.ngener*ir->opts.ngener);
- }
}
- snew(fr->nblists, fr->nnblists);
-
/* This code automatically gives table length tabext without cut-off's,
* in that case grompp should already have checked that we do not need
* normal tables and we only generate tables for 1-4 interactions.
/* make tables for ordinary interactions */
if (bSomeNormalNbListsAreInUse)
{
- make_nbf_tables(fp, ic, rtab, tabfn, nullptr, nullptr, &fr->nblists[0]);
+ make_nbf_tables(fp, ic, rtab, tabfn, nullptr, nullptr, nullptr);
m = 1;
}
else
egp_flags = ir->opts.egp_flags[GID(egi, egj, ir->opts.ngener)];
if ((egp_flags & EGP_TABLE) && !(egp_flags & EGP_EXCL))
{
- if (fr->nnblists > 1)
- {
- fr->gid2nblists[GID(egi, egj, ir->opts.ngener)] = m;
- }
/* Read the table file with the two energy groups names appended */
make_nbf_tables(fp, ic, rtab, tabfn,
*mtop->groups.groupNames[nm_ind[egi]],
*mtop->groups.groupNames[nm_ind[egj]],
- &fr->nblists[m]);
+ nullptr);
m++;
}
- else if (fr->nnblists > 1)
- {
- fr->gid2nblists[GID(egi, egj, ir->opts.ngener)] = 0;
- }
}
}
}
fr->print_force = print_force;
-
- /* Initialize neighbor search */
- snew(fr->ns, 1);
- init_ns(fp, cr, fr->ns, fr, mtop);
-
/* Initialize the thread working data for bonded interactions */
init_bonded_threading(fp, mtop->groups.groups[SimulationAtomGroupType::EnergyOutput].nr,
&fr->bondedThreading);
}
}
-void done_forcerec(t_forcerec *fr, int numMolBlocks, int numEnergyGroups)
+void done_forcerec(t_forcerec *fr, int numMolBlocks)
{
if (fr == nullptr)
{
done_interaction_const(fr->ic);
sfree(fr->shift_vec);
sfree(fr->fshift);
- sfree(fr->nblists);
- done_ns(fr->ns, numEnergyGroups);
sfree(fr->ewc_t);
tear_down_bonded_threading(fr->bondedThreading);
GMX_RELEASE_ASSERT(fr->gpuBonded == nullptr, "Should have been deleted earlier, when used");
}
//! Destroy a forcerec.
-void done_forcerec(t_forcerec *fr, int numMolBlocks, int numEnergyGroups);
+void done_forcerec(t_forcerec *fr, int numMolBlocks);
/*! \brief Print the contents of the forcerec to a file
*
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018,2019, 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.
- *
- * 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 research papers on the package. Check out http://www.gromacs.org.
- */
-#include "gmxpre.h"
-
-#include "ns.h"
-
-#include <cmath>
-#include <cstdlib>
-#include <cstring>
-
-#include <algorithm>
-
-#include "gromacs/domdec/domdec.h"
-#include "gromacs/domdec/domdec_struct.h"
-#include "gromacs/gmxlib/network.h"
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/gmxlib/nonbonded/nonbonded.h"
-#include "gromacs/math/functions.h"
-#include "gromacs/math/utilities.h"
-#include "gromacs/math/vec.h"
-#include "gromacs/math/vecdump.h"
-#include "gromacs/mdlib/nsgrid.h"
-#include "gromacs/mdlib/qmmm.h"
-#include "gromacs/mdtypes/commrec.h"
-#include "gromacs/mdtypes/forcerec.h"
-#include "gromacs/mdtypes/group.h"
-#include "gromacs/mdtypes/inputrec.h"
-#include "gromacs/mdtypes/md_enums.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"
-
-/*
- * 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[], int i, int j)
-{
- e[j] = e[j] | (1<<i);
-}
-static void RMEXCL_(t_excl e[], int i, int j)
-{
- e[j] = e[j] & ~(1<<i);
-}
-static gmx_bool ISEXCL_(t_excl e[], int i, int j)
-{
- return (gmx_bool)(e[j] & (1<<i));
-}
-static gmx_bool NOTEXCL_(t_excl e[], int i, int j)
-{
- return !(ISEXCL(e, i, j));
-}
-#else
-#define SETEXCL(e, i, j) (e)[int(j)] |= (1<<(int(i)))
-#define RMEXCL(e, i, j) (e)[int(j)] &= (~(1<<(int(i))))
-#define ISEXCL(e, i, j) static_cast<gmx_bool>((e)[(int(j))] & (1<<(int(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;
-
- 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
- *
- ************************************************/
-
-void reallocate_nblist(t_nblist *nl)
-{
- if (gmx_debug_at)
- {
- 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->gid, nl->maxnri);
- srenew(nl->shift, nl->maxnri);
- srenew(nl->jindex, nl->maxnri+1);
-}
-
-
-static void init_nblist(t_nblist *nl_sr,
- int maxsr,
- int ivdw, int ivdwmod,
- int ielec, int ielecmod,
- int igeometry, int type)
-{
- t_nblist *nl;
- int homenr;
-
- {
- nl = nl_sr;
- homenr = maxsr;
-
- if (nl == nullptr)
- {
- return;
- }
-
-
- /* Set coul/vdw in neighborlist, and for the normal loops we determine
- * an index of which one to call.
- */
- nl->ivdw = ivdw;
- nl->ivdwmod = ivdwmod;
- nl->ielec = ielec;
- nl->ielecmod = ielecmod;
- nl->type = type;
- nl->igeometry = igeometry;
-
- if (nl->type == GMX_NBLIST_INTERACTION_FREE_ENERGY)
- {
- nl->igeometry = GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE;
- }
-
- gmx_nonbonded_set_kernel_pointers(nl);
-
- /* 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.
- * 4 seems to be a reasonable factor, which only causes reallocation
- * during runs with tiny and many energygroups.
- */
- nl->maxnri = homenr*4;
- nl->maxnrj = 0;
- nl->nri = -1;
- nl->nrj = 0;
- nl->iinr = nullptr;
- nl->gid = nullptr;
- nl->shift = nullptr;
- nl->jindex = nullptr;
- nl->jjnr = nullptr;
- nl->excl_fep = nullptr;
- reallocate_nblist(nl);
- nl->jindex[0] = 0;
-
- if (debug)
- {
- fprintf(debug, "Initiating neighbourlist (ielec=%d, ivdw=%d, type=%d) for %s interactions,\nwith %d SR atoms.\n",
- nl->ielec, nl->ivdw, nl->type, gmx_nblist_geometry_names[nl->igeometry], maxsr);
- }
- }
-}
-
-void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr)
-{
- int maxsr, maxsr_wat;
- int ielec, ivdw, ielecmod, ivdwmod, type;
- 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 = std::min(fr->nWatMol, (homenr+2)/3);
-
- /* 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;
-
- fr->ns->bCGlist = (getenv("GMX_NBLISTCG") != nullptr);
- if (!fr->ns->bCGlist)
- {
- igeometry_def = GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE;
- }
- else
- {
- igeometry_def = GMX_NBLIST_GEOMETRY_CG_CG;
- if (log != nullptr)
- {
- 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(&nbl->nlist_sr[eNL_VDWQQ],
- maxsr, ivdw, ivdwmod, ielec, ielecmod, igeometry_def, type);
- init_nblist(&nbl->nlist_sr[eNL_VDW],
- maxsr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, igeometry_def, type);
- init_nblist(&nbl->nlist_sr[eNL_QQ],
- maxsr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_def, type);
- init_nblist(&nbl->nlist_sr[eNL_VDWQQ_WATER],
- maxsr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_w, type);
- init_nblist(&nbl->nlist_sr[eNL_QQ_WATER],
- maxsr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_w, type);
- init_nblist(&nbl->nlist_sr[eNL_VDWQQ_WATERWATER],
- maxsr_wat, ivdw, ivdwmod, ielec, ielecmod, igeometry_ww, type);
- init_nblist(&nbl->nlist_sr[eNL_QQ_WATERWATER],
- maxsr_wat, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, igeometry_ww, type);
-
- /* Did we get the solvent loops so we can use optimized water kernels? */
- if (nbl->nlist_sr[eNL_VDWQQ_WATER].kernelptr_vf == nullptr
- || nbl->nlist_sr[eNL_QQ_WATER].kernelptr_vf == nullptr
- || nbl->nlist_sr[eNL_VDWQQ_WATERWATER].kernelptr_vf == nullptr
- || nbl->nlist_sr[eNL_QQ_WATERWATER].kernelptr_vf == nullptr)
- {
- fr->solvent_opt = esolNO;
- if (log != nullptr)
- {
- fprintf(log, "Note: The available nonbonded kernels do not support water optimization - disabling.\n");
- }
- }
-
- if (fr->efep != efepNO)
- {
- init_nblist(&nbl->nlist_sr[eNL_VDWQQ_FREE],
- maxsr, ivdw, ivdwmod, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
- init_nblist(&nbl->nlist_sr[eNL_VDW_FREE],
- maxsr, ivdw, ivdwmod, GMX_NBKERNEL_ELEC_NONE, eintmodNONE, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
- init_nblist(&nbl->nlist_sr[eNL_QQ_FREE],
- maxsr, GMX_NBKERNEL_VDW_NONE, eintmodNONE, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_FREE_ENERGY);
- }
- }
- /* QMMM MM list */
- if (fr->bQMMM && fr->qr->QMMMscheme != eQMMMschemeoniom)
- {
- if (nullptr == fr->QMMMlist)
- {
- snew(fr->QMMMlist, 1);
- }
- init_nblist(fr->QMMMlist,
- maxsr, 0, 0, ielec, ielecmod, GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE, GMX_NBLIST_INTERACTION_STANDARD);
- }
-
- if (log != nullptr)
- {
- fprintf(log, "\n");
- }
-
- fr->ns->nblist_initialized = TRUE;
-}
-
-static void reset_nblist(t_nblist *nl)
-{
- GMX_RELEASE_ASSERT(nl, "Should only reset valid nblists");
-
- nl->nri = -1;
- nl->nrj = 0;
- if (nl->jindex)
- {
- nl->jindex[0] = 0;
- }
-}
-
-static void reset_neighbor_lists(t_forcerec *fr)
-{
- int n, i;
-
- if (fr->bQMMM && fr->qr->QMMMscheme != eQMMMschemeoniom)
- {
- /* only reset the short-range nblist */
- reset_nblist(fr->QMMMlist);
- }
-
- for (n = 0; n < fr->nnblists; n++)
- {
- for (i = 0; i < eNL_NR; i++)
- {
- reset_nblist( &(fr->nblists[n].nlist_sr[i]) );
- }
- }
-}
-
-
-
-
-static inline void new_i_nblist(t_nblist *nlist, int i_atom, int shift, int gid)
-{
- int 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->gid[nri] != gid))
- {
- /* 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->gid[nri] != -1)))
- {
- /* If so increase the counter */
- nlist->nri++;
- nri++;
- if (nlist->nri >= nlist->maxnri)
- {
- nlist->maxnri += over_alloc_large(nlist->nri);
- reallocate_nblist(nlist);
- }
- }
- /* Set the number of neighbours and the atom number */
- nlist->jindex[nri+1] = nlist->jindex[nri];
- nlist->iinr[nri] = i_atom;
- 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)
-{
- int nri = nlist->nri;
- int len;
-
- if (nri >= 0)
- {
- /* 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.
- */
- while ((nlist->nrj % nlist->simd_padding_width) != 0)
- {
- /* 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];
- /* 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 == 0) && (nlist->nri > 0))
- {
- nlist->nri--;
- }
- }
-}
-
-static 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.
- */
- if (nlist->nrj == 0)
- {
- /* Some assembly kernels do not support empty lists,
- * make sure here that we don't generate any empty lists.
- * With the current ns code this branch is taken in two cases:
- * No i-particles at all: nri=-1 here
- * There are i-particles, but no j-particles; nri=0 here
- */
- nlist->nri = 0;
- }
- else
- {
- /* Close list number nri by incrementing the count */
- nlist->nri++;
- }
-}
-
-static inline void close_neighbor_lists(t_forcerec *fr, gmx_bool bMakeQMMMnblist)
-{
- int n, i;
-
- if (bMakeQMMMnblist)
- {
- close_nblist(fr->QMMMlist);
- }
-
- for (n = 0; n < fr->nnblists; n++)
- {
- for (i = 0; (i < eNL_NR); i++)
- {
- close_nblist(&(fr->nblists[n].nlist_sr[i]));
- }
- }
-}
-
-
-static inline void add_j_to_nblist(t_nblist *nlist, int j_atom)
-{
- int nrj = nlist->nrj;
-
- if (nlist->nrj >= nlist->maxnrj)
- {
- nlist->maxnrj = round_up_to_simd_width(over_alloc_small(nlist->nrj + 1), nlist->simd_padding_width);
-
- if (gmx_debug_at)
- {
- fprintf(debug, "Increasing SR nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
- nlist->ielec, nlist->ivdw, nlist->type, nlist->igeometry, nlist->maxnrj);
- }
-
- srenew(nlist->jjnr, nlist->maxnrj);
- }
-
- nlist->jjnr[nrj] = j_atom;
- nlist->nrj++;
-}
-
-static inline void add_j_to_nblist_cg(t_nblist *nlist,
- int j_start, int j_end,
- const t_excl *bexcl, gmx_bool i_is_j)
-{
- 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 SR nblist (ielec=%d,ivdw=%d,type=%d,igeometry=%d) j size to %d\n",
- 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;
- nlist->jjnr_end[nrj] = j_end;
-
- 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);
- }
-
- /* Set the exclusions */
- 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++)
- {
- nlist->excl[nrj*MAX_CGCGSIZE + j] |= (1<<j) - 1;
- }
- }
-
- nlist->nrj++;
-}
-
-typedef void
- put_in_list_t (const gmx_bool bHaveVdW[],
- int ngid,
- const t_mdatoms *md,
- int icg,
- int jgid,
- int nj,
- const int jjcg[],
- const int index[],
- const t_excl bExcl[],
- int shift,
- t_forcerec * fr,
- gmx_bool bDoVdW,
- gmx_bool bDoCoul,
- int solvent_opt);
-
-static void
-put_in_list_at(const gmx_bool bHaveVdW[],
- int ngid,
- const t_mdatoms *md,
- int icg,
- int jgid,
- int nj,
- const int jjcg[],
- const int index[],
- const t_excl bExcl[],
- int shift,
- t_forcerec * fr,
- 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 = nullptr;
- t_nblist * vdw_free = nullptr;
- t_nblist * coul_free = nullptr;
- t_nblist * vdwc_ww = nullptr;
- t_nblist * coul_ww = nullptr;
-
- int i, j, jcg, igid, gid, nbl_ind;
- int jj, jj0, jj1, i_atom;
- int i0, nicg;
-
- int *cginfo;
- int *type, *typeB;
- real *charge, *chargeB;
- real qi, qiB;
- 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;
- chargeB = md->chargeB;
- 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;
-
- bFreeEnergy = FALSE;
- if (md->nPerturbed)
- {
- /* 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++)
- {
- bFreeEnergy |= bPert[i0+i];
- }
- /* Loop over the j charge groups */
- 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++)
- {
- bFreeEnergy |= bPert[jj];
- }
- }
- }
-
- /* Unpack pointers to neighbourlist structs */
- if (fr->nnblists == 1)
- {
- nbl_ind = 0;
- }
- else
- {
- nbl_ind = fr->gid2nblists[GID(igid, jgid, ngid)];
- }
- nlist = fr->nblists[nbl_ind].nlist_sr;
-
- if (iwater != esolNO)
- {
- 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];
- }
- else
- {
- vdwc = &nlist[eNL_VDWQQ];
- vdw = &nlist[eNL_VDW];
- coul = &nlist[eNL_QQ];
- }
-
- if (!bFreeEnergy)
- {
- if (iwater != esolNO)
- {
- /* Loop over the atoms in the i charge group */
- i_atom = i0;
- gid = GID(igid, jgid, ngid);
- /* Create new i_atom for each energy group */
- if (bDoCoul && bDoVdW)
- {
- new_i_nblist(vdwc, i_atom, shift, gid);
- new_i_nblist(vdwc_ww, 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(coul_ww, i_atom, shift, gid);
- }
- /* Loop over the j charge groups */
- for (j = 0; (j < nj); j++)
- {
- jcg = jjcg[j];
-
- if (jcg == icg)
- {
- continue;
- }
-
- 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);
- }
- else
- {
- /* One entry for the entire water-water interaction */
- if (!bDoVdW)
- {
- add_j_to_nblist(coul_ww, jj0);
- }
- else
- {
- add_j_to_nblist(vdwc_ww, jj0);
- }
- }
- }
- 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);
- }
- else
- {
- /* One entry for the entire water-water interaction */
- if (!bDoVdW)
- {
- add_j_to_nblist(coul_ww, jj0);
- }
- else
- {
- add_j_to_nblist(vdwc_ww, jj0);
- }
- }
- }
- 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)
- {
- for (jj = jj0; (jj < jj1); jj++)
- {
- if (charge[jj] != 0)
- {
- add_j_to_nblist(coul, jj);
- }
- }
- }
- else if (!bDoCoul)
- {
- for (jj = jj0; (jj < jj1); jj++)
- {
- if (bHaveVdW[type[jj]])
- {
- add_j_to_nblist(vdw, jj);
- }
- }
- }
- else
- {
- /* _charge_ _groups_ interact with both coulomb and LJ */
- /* Check which atoms we should add to the lists! */
- for (jj = jj0; (jj < jj1); jj++)
- {
- if (bHaveVdW[type[jj]])
- {
- if (charge[jj] != 0)
- {
- add_j_to_nblist(vdwc, jj);
- }
- else
- {
- add_j_to_nblist(vdw, jj);
- }
- }
- else if (charge[jj] != 0)
- {
- add_j_to_nblist(coul, jj);
- }
- }
- }
- }
- }
- close_i_nblist(vdw);
- close_i_nblist(coul);
- close_i_nblist(vdwc);
- close_i_nblist(coul_ww);
- 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++)
- {
- 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);
- }
- if (bDoVdW)
- {
- new_i_nblist(vdw, i_atom, shift, gid);
- }
- if (bDoCoul)
- {
- 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)
- {
- /* 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);
-
- if (bNotEx)
- {
- if (!bDoVdW_i)
- {
- if (charge[jj] != 0)
- {
- add_j_to_nblist(coul, jj);
- }
- }
- else if (!bDoCoul_i)
- {
- if (bHaveVdW[type[jj]])
- {
- add_j_to_nblist(vdw, jj);
- }
- }
- else
- {
- if (bHaveVdW[type[jj]])
- {
- if (charge[jj] != 0)
- {
- add_j_to_nblist(vdwc, jj);
- }
- else
- {
- add_j_to_nblist(vdw, jj);
- }
- }
- else if (charge[jj] != 0)
- {
- add_j_to_nblist(coul, jj);
- }
- }
- }
- }
- }
- }
- close_i_nblist(vdw);
- close_i_nblist(coul);
- close_i_nblist(vdwc);
- }
- }
- }
- else
- {
- /* we are doing free energy */
- 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++)
- {
- i_atom = i0+i;
- 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, 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));
- /* 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
- * for a TIP4p i charge group.
- * For SPC type water the first atom has LJ and charge,
- * so there is no such problem.
- */
- if (iwater == esolNO)
- {
- bDoCoul_i_sol = bDoCoul_i;
- }
- else
- {
- bDoCoul_i_sol = bDoCoul;
- }
-
- if (bDoVdW_i || bDoCoul_i_sol)
- {
- /* 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 */
- bFree = bPert[i_atom];
- 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)
- {
- bNotEx = NOTEXCL(bExcl, i, jj);
-
- if (bNotEx)
- {
- if (bFreeJ)
- {
- if (!bDoVdW_i)
- {
- if (charge[jj] != 0 || chargeB[jj] != 0)
- {
- add_j_to_nblist(coul_free, jj);
- }
- }
- else if (!bDoCoul_i)
- {
- if (bHaveVdW[type[jj]] || bHaveVdW[typeB[jj]])
- {
- add_j_to_nblist(vdw_free, jj);
- }
- }
- else
- {
- if (bHaveVdW[type[jj]] || bHaveVdW[typeB[jj]])
- {
- if (charge[jj] != 0 || chargeB[jj] != 0)
- {
- add_j_to_nblist(vdwc_free, jj);
- }
- else
- {
- add_j_to_nblist(vdw_free, jj);
- }
- }
- else if (charge[jj] != 0 || chargeB[jj] != 0)
- {
- add_j_to_nblist(coul_free, jj);
- }
- }
- }
- else if (!bDoVdW_i)
- {
- /* This is done whether or not bWater is set */
- if (charge[jj] != 0)
- {
- add_j_to_nblist(coul, jj);
- }
- }
- else if (!bDoCoul_i_sol)
- {
- if (bHaveVdW[type[jj]])
- {
- add_j_to_nblist(vdw, jj);
- }
- }
- else
- {
- if (bHaveVdW[type[jj]])
- {
- if (charge[jj] != 0)
- {
- add_j_to_nblist(vdwc, jj);
- }
- else
- {
- add_j_to_nblist(vdw, jj);
- }
- }
- else if (charge[jj] != 0)
- {
- add_j_to_nblist(coul, jj);
- }
- }
- }
- }
- }
- }
- }
- close_i_nblist(vdw);
- close_i_nblist(coul);
- close_i_nblist(vdwc);
- close_i_nblist(vdw_free);
- close_i_nblist(coul_free);
- close_i_nblist(vdwc_free);
- }
- }
-}
-
-static void
-put_in_list_qmmm(const gmx_bool gmx_unused bHaveVdW[],
- int ngid,
- const t_mdatoms * /* md */,
- int icg,
- int jgid,
- int nj,
- const int jjcg[],
- const int index[],
- const t_excl bExcl[],
- int shift,
- t_forcerec * fr,
- gmx_bool gmx_unused bDoVdW,
- gmx_bool gmx_unused bDoCoul,
- int gmx_unused solvent_opt)
-{
- t_nblist * coul;
- int i, j, jcg, igid, gid;
- int 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++)
- {
- i_atom = i0+i;
- gid = GID(igid, jgid, ngid);
- /* Create new i_atom for each energy group */
- new_i_nblist(coul, i_atom, shift, gid);
-
- /* Loop over the j charge groups */
- for (j = 0; j < nj; j++)
- {
- jcg = jjcg[j];
-
- /* Charge groups cannot have QM and MM atoms simultaneously */
- 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++)
- {
- bNotEx = NOTEXCL(bExcl, i, jj);
- if (bNotEx)
- {
- add_j_to_nblist(coul, jj);
- }
- }
- }
- }
- close_i_nblist(coul);
- }
-}
-
-static void
-put_in_list_cg(const gmx_bool gmx_unused bHaveVdW[],
- int ngid,
- const t_mdatoms * /* md */,
- int icg,
- int jgid,
- int nj,
- const int jjcg[],
- const int index[],
- const t_excl bExcl[],
- int shift,
- t_forcerec * fr,
- gmx_bool gmx_unused bDoVdW,
- gmx_bool gmx_unused bDoCoul,
- int gmx_unused solvent_opt)
-{
- int cginfo;
- int igid, gid, nbl_ind;
- t_nblist * vdwc;
- int j, jcg;
-
- cginfo = fr->cginfo[icg];
-
- igid = GET_CGINFO_GID(cginfo);
- gid = GID(igid, jgid, ngid);
-
- /* Unpack pointers to neighbourlist structs */
- if (fr->nnblists == 1)
- {
- nbl_ind = 0;
- }
- else
- {
- nbl_ind = fr->gid2nblists[gid];
- }
- vdwc = &fr->nblists[nbl_ind].nlist_sr[eNL_VDWQQ];
-
- /* Make a new neighbor list for charge group icg.
- * Currently simply one neighbor list is made with LJ and Coulomb.
- * If required, zero interactions could be removed here
- * or in the force loop.
- */
- new_i_nblist(vdwc, index[icg], shift, gid);
- vdwc->iinr_end[vdwc->nri] = index[icg+1];
-
- for (j = 0; (j < nj); j++)
- {
- jcg = jjcg[j];
- /* Skip the icg-icg pairs if all self interactions are excluded */
- if (!(jcg == icg && GET_CGINFO_EXCL_INTRA(cginfo)))
- {
- /* 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);
- }
- }
-
- close_i_nblist(vdwc);
-}
-
-static void setexcl(int start, int end, t_blocka *excl, gmx_bool b,
- t_excl bexcl[])
-{
- int i, k;
-
- if (b)
- {
- for (i = start; i < end; i++)
- {
- for (k = excl->index[i]; k < excl->index[i+1]; k++)
- {
- SETEXCL(bexcl, i-start, excl->a[k]);
- }
- }
- }
- else
- {
- for (i = start; i < end; i++)
- {
- for (k = excl->index[i]; k < excl->index[i+1]; k++)
- {
- RMEXCL(bexcl, i-start, excl->a[k]);
- }
- }
- }
-}
-
-int calc_naaj(int icg, int cgtot)
-{
- int naaj;
-
- if ((cgtot % 2) == 1)
- {
- /* Odd number of charge groups, easy */
- naaj = 1 + (cgtot/2);
- }
- else if ((cgtot % 4) == 0)
- {
- /* Multiple of four is hard */
- if (icg < cgtot/2)
- {
- if ((icg % 2) == 0)
- {
- naaj = 1+(cgtot/2);
- }
- else
- {
- naaj = cgtot/2;
- }
- }
- else
- {
- if ((icg % 2) == 1)
- {
- naaj = 1+(cgtot/2);
- }
- else
- {
- naaj = cgtot/2;
- }
- }
- }
- else
- {
- /* cgtot/2 = odd */
- if ((icg % 2) == 0)
- {
- naaj = 1+(cgtot/2);
- }
- else
- {
- naaj = cgtot/2;
- }
- }
-#ifdef DEBUG
- fprintf(log, "naaj=%d\n", naaj);
-#endif
-
- return naaj;
-}
-
-/************************************************
- *
- * S I M P L E C O R E S T U F F
- *
- ************************************************/
-
-static real calc_image_tric(const rvec xi, const rvec xj, matrix box,
- const 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;
-
- /* 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 = static_cast<int>(dz*b_inv[ZZ] + h25);
- tz -= 2;
- dz -= tz*box[ZZ][ZZ];
- dy -= tz*box[ZZ][YY];
- dx -= tz*box[ZZ][XX];
-
- ty = static_cast<int>(dy*b_inv[YY] + h25);
- ty -= 2;
- dy -= ty*box[YY][YY];
- dx -= ty*box[YY][XX];
-
- tx = static_cast<int>(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);
-
- return r2;
-}
-
-static real calc_image_rect(const rvec xi, const rvec xj, const rvec box_size,
- const rvec b_inv, int *shift)
-{
- 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 = static_cast<int>(dx*b_inv[XX] + h15);
- ty = static_cast<int>(dy*b_inv[YY] + h15);
- tz = static_cast<int>(dz*b_inv[ZZ] + h15);
- 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);
-
- return r2;
-}
-
-static void add_simple(t_ns_buf * nsbuf,
- int nrj,
- int cg_j,
- gmx_bool bHaveVdW[],
- int ngid,
- const 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, TRUE, TRUE, fr->solvent_opt);
- /* Reset buffer contents */
- nsbuf->ncg = nsbuf->nj = 0;
- }
- nsbuf->jcg[nsbuf->ncg++] = cg_j;
- nsbuf->nj += nrj;
-}
-
-static void ns_inner_tric(rvec x[],
- int icg,
- const int *i_egp_flags,
- int njcg,
- const int jcg[],
- matrix box,
- rvec b_inv,
- real rcut2,
- t_block *cgs,
- t_ns_buf **ns_buf,
- gmx_bool bHaveVdW[],
- int ngid,
- const 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;
- int cg_j, *cgindex;
-
- cgindex = cgs->index;
- shift = CENTRAL;
- 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)
- {
- 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,
- put_in_list);
- }
- }
- }
-}
-
-static void ns_inner_rect(rvec x[],
- int icg,
- const int *i_egp_flags,
- int njcg,
- const int 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,
- const 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;
- int cg_j, *cgindex;
-
- cgindex = cgs->index;
- if (bBox)
- {
- shift = CENTRAL;
- 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)
- {
- 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,
- put_in_list);
- }
- }
- }
- }
- else
- {
- 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))
- {
- 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,
- put_in_list);
- }
- }
- }
- }
-}
-
-/* ns_simple_core needs to be adapted for QMMM still 2005 */
-
-static int ns_simple_core(t_forcerec *fr,
- gmx_localtop_t *top,
- const t_mdatoms *md,
- matrix box,
- rvec box_size,
- t_excl bexcl[],
- int *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, igid, nn;
- int *cginfo;
- t_ns_buf *nsbuf;
- /* int *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;
-
- rlist2 = gmx::square(fr->rlist);
-
- bBox = (fr->ePBC != epbcNONE);
- if (bBox)
- {
- for (m = 0; (m < DIM); m++)
- {
- if (gmx_numzero(box_size[m]))
- {
- gmx_fatal(FARGS, "Dividing by zero box size!");
- }
- b_inv[m] = 1.0/box_size[m];
- }
- bTriclinic = TRICLINIC(box);
- }
- else
- {
- bTriclinic = FALSE;
- }
-
- cginfo = fr->cginfo;
-
- 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]);
- 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);
- 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);
- }
- 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);
- }
- nsearch += naaj;
-
- for (nn = 0; (nn < ngid); nn++)
- {
- 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, 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);
- }
- close_neighbor_lists(fr, FALSE);
-
- return nsearch;
-}
-
-/************************************************
- *
- * N S 5 G R I D S T U F F
- *
- ************************************************/
-
-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)
-{
- real dcx, tmp;
- int g_min, g_max, shift_home;
-
- if (xgi < 0)
- {
- g_min = 0;
- g_max = Nx - 1;
- *g0 = 0;
- *g1 = -1;
- }
- else if (xgi >= Nx)
- {
- g_min = 0;
- g_max = Nx - 1;
- *g0 = Nx;
- *g1 = Nx - 1;
- }
- else
- {
- if (ncpddc == 0)
- {
- g_min = 0;
- g_max = Nx - 1;
- }
- else
- {
- if (xgi < ncpddc)
- {
- shift_home = 0;
- }
- else
- {
- shift_home = -1;
- }
- g_min = (shift_min == shift_home ? 0 : ncpddc);
- g_max = (shift_max == shift_home ? ncpddc - 1 : Nx - 1);
- }
- if (shift_min > 0)
- {
- *g0 = g_min;
- *g1 = g_min - 1;
- }
- else if (shift_max < 0)
- {
- *g0 = g_max + 1;
- *g1 = g_max;
- }
- else
- {
- *g0 = xgi;
- *g1 = xgi;
- dcx2[xgi] = 0;
- }
- }
-
- while (*g0 > g_min)
- {
- /* Check one grid cell down */
- dcx = ((*g0 - 1) + 1)*gridx - x;
- tmp = dcx*dcx;
- if (tmp >= rc2)
- {
- break;
- }
- (*g0)--;
- dcx2[*g0] = tmp;
- }
-
- while (*g1 < g_max)
- {
- /* Check one grid cell up */
- dcx = (*g1 + 1)*gridx - x;
- tmp = dcx*dcx;
- if (tmp >= rc2)
- {
- break;
- }
- (*g1)++;
- dcx2[*g1] = tmp;
- }
-}
-
-
-#define calc_dx2(XI, YI, ZI, y) (gmx::square((XI)-(y)[XX]) + gmx::square((YI)-(y)[YY]) + gmx::square((ZI)-(y)[ZZ]))
-#define calc_cyl_dx2(XI, YI, y) (gmx::square((XI)-(y)[XX]) + gmx::square((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
- * at least 1x1x1, see GROMACS manual
- *
- ****************************************************/
-
-
-static void get_cutoff2(t_forcerec *fr, real *rs2)
-{
- *rs2 = gmx::square(fr->rlist);
-}
-
-static void init_nsgrid_lists(t_forcerec *fr, int ngid, gmx_ns_t *ns)
-{
- real rs2;
- int j;
-
- get_cutoff2(fr, &rs2);
-
- /* Short range buffers */
- snew(ns->nl_sr, ngid);
- /* Counters */
- snew(ns->nsr, ngid);
-
- for (j = 0; (j < ngid); j++)
- {
- snew(ns->nl_sr[j], MAX_CG);
- }
- if (debug)
- {
- fprintf(debug,
- "ns5_core: rs2 = %g (nm^2)\n",
- rs2);
- }
-}
-
-static int nsgrid_core(const t_commrec *cr,
- t_forcerec *fr,
- matrix box,
- int ngid,
- gmx_localtop_t *top,
- t_grid *grid,
- t_excl bexcl[],
- const gmx_bool *bExcludeAlleg,
- const t_mdatoms *md,
- put_in_list_t *put_in_list,
- gmx_bool bHaveVdW[],
- gmx_bool bMakeQMMMnblist)
-{
- gmx_ns_t *ns;
- int **nl_sr;
- int *nsr;
- gmx_domdec_t *dd;
- const t_block *cgs = &(top->cgs);
- int *cginfo = fr->cginfo;
- /* int *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;
-#ifdef ALLOW_OFFDIAG_LT_HALFDIAG
- 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;
- real *dcx2, *dcy2, *dcz2;
- int zgi, ygi, xgi;
- int cg0, cg1, icg = -1, cgsnr, i0, igid, naaj, max_jcg;
- int jcg0, jcg1, jjcg, cgj0, jgid;
- int *grida, *gridnra, *gridind;
- rvec *cgcm, grid_offset;
- real r2, rs2, XI, YI, ZI, tmp1, tmp2;
- int *i_egp_flags;
- gmx_bool bDomDec, bTriclinicX, bTriclinicY;
- ivec ncpddc;
-
- ns = fr->ns;
-
- bDomDec = DOMAINDECOMP(cr);
- dd = cr->dd;
-
- bTriclinicX = ((YY < grid->npbcdim &&
- (!bDomDec || dd->nc[YY] == 1) && box[YY][XX] != 0) ||
- (ZZ < grid->npbcdim &&
- (!bDomDec || dd->nc[ZZ] == 1) && box[ZZ][XX] != 0));
- bTriclinicY = (ZZ < grid->npbcdim &&
- (!bDomDec || dd->nc[ZZ] == 1) && box[ZZ][YY] != 0);
-
- cgsnr = cgs->nr;
-
- get_cutoff2(fr, &rs2);
-
- nl_sr = ns->nl_sr;
- nsr = ns->nsr;
-
- /* Unpack arrays */
- cgcm = fr->cg_cm;
- Nx = grid->n[XX];
- Ny = grid->n[YY];
- Nz = grid->n[ZZ];
- grida = grid->a;
- 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;
- 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)
- {
- /* This could happen due to rounding, when both ratios are 0.5 */
- ysh_tx = 0;
- }
-#endif
-
- if (fr->n_tpi)
- {
- /* We only want a list for the test particle */
- cg0 = cgsnr - 1;
- }
- else
- {
- cg0 = grid->icg0;
- }
- cg1 = grid->icg1;
-
- /* Set the shift range */
- for (d = 0; d < DIM; d++)
- {
- sh0[d] = -1;
- sh1[d] = 1;
- /* Check if we need periodicity shifts.
- * Without PBC or with domain decomposition we don't need them.
- */
- if (d >= ePBC2npbcdim(fr->ePBC) || (bDomDec && dd->nc[d] > 1))
- {
- shp[d] = 0;
- }
- else
- {
- if (d == XX &&
- box[XX][XX] - std::abs(box[YY][XX]) - std::abs(box[ZZ][XX]) < std::sqrt(rs2))
- {
- shp[d] = 2;
- }
- else
- {
- shp[d] = 1;
- }
- }
- }
-
- /* Loop over charge groups */
- for (icg = cg0; (icg < cg1); icg++)
- {
- igid = GET_CGINFO_GID(cginfo[icg]);
- /* Skip this charge group if all energy groups are excluded! */
- if (bExcludeAlleg[igid])
- {
- 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])
- {
- 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);
- jcg0 = icg;
- jcg1 = icg + naaj;
- 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);
- max_jcg = 0;
- }
- else
- {
- /* Compute the number of charge groups that fall within the control
- * of this one (icg)
- */
- naaj = calc_naaj(icg, cgsnr);
- jcg0 = icg;
- jcg1 = icg + naaj;
-
- if (fr->n_tpi)
- {
- /* The i-particle is awlways the test particle,
- * so we want all j-particles
- */
- max_jcg = cgsnr - 1;
- }
- else
- {
- max_jcg = jcg1 - cgsnr;
- }
- }
- }
-
- 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)
- */
-#ifdef NS5DB
- 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++)
- {
- 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;
- get_dx_dd(Nz, gridz, rs2, zgi, ZI-grid_offset[ZZ],
- ncpddc[ZZ], sh0[ZZ], sh1[ZZ], &dz0, &dz1, dcz2);
- if (dz0 > dz1)
- {
- continue;
- }
- 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 */
- if (bTriclinicY)
- {
- ygi = static_cast<int>(Ny + (YI - grid_offset[YY])*grid_y) - Ny;
- }
- else
- {
- ygi = cell_y + ty*Ny;
- }
- get_dx_dd(Ny, gridy, rs2, ygi, YI-grid_offset[YY],
- ncpddc[YY], sh0[YY], sh1[YY], &dy0, &dy1, dcy2);
- if (dy0 > dy1)
- {
- continue;
- }
- 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 */
- if (bTriclinicX)
- {
- xgi = static_cast<int>(Nx + (XI - grid_offset[XX])*grid_x) - Nx;
- }
- else
- {
- xgi = cell_x + tx*Nx;
- }
- get_dx_dd(Nx, gridx, rs2, xgi, XI-grid_offset[XX],
- ncpddc[XX], sh0[XX], sh1[XX], &dx0, &dx1, dcx2);
- if (dx0 > dx1)
- {
- continue;
- }
- /* Get shift vector */
- shift = XYZ2IS(tx, ty, tz);
-#ifdef NS5DB
- range_check(shift, 0, SHIFTS);
-#endif
- for (nn = 0; (nn < ngid); nn++)
- {
- nsr[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);
-#endif
- for (dx = dx0; (dx <= dx1); dx++)
- {
- tmp1 = rs2 - dcx2[dx];
- for (dy = dy0; (dy <= dy1); dy++)
- {
- tmp2 = tmp1 - dcy2[dy];
- if (tmp2 > 0)
- {
- 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);
-
- /* 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.
- */
- if (nrj == 0 ||
- (grida[cgj0] >= max_jcg &&
- (grida[cgj0] >= jcg1 || grida[cgj0+nrj-1] < jcg0)))
- {
- continue;
- }
-
- /* Loop over cgs */
- 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 < rs2)
- {
- /* jgid = gid[cgsatoms[cgsindex[jjcg]]]; */
- jgid = GET_CGINFO_GID(cginfo[jjcg]);
- /* check energy group exclusions */
- if (!(i_egp_flags[jgid] & EGP_EXCL))
- {
- 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, TRUE, TRUE, fr->solvent_opt);
- nsr[jgid] = 0;
- }
- nl_sr[jgid][nsr[jgid]++] = jjcg;
- }
- }
- nns++;
- }
- }
- }
- }
- }
- }
- }
- /* CHECK whether there is anything left in the buffers */
- for (nn = 0; (nn < ngid); nn++)
- {
- if (nsr[nn] > 0)
- {
- put_in_list(bHaveVdW, ngid, md, icg, nn, nsr[nn], nl_sr[nn],
- cgs->index, /* cgsatoms, */ bexcl,
- shift, fr, TRUE, TRUE, fr->solvent_opt);
- }
- }
- }
- }
- }
- 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.
- */
-
- /* Close neighbourlists */
- close_neighbor_lists(fr, bMakeQMMMnblist);
-
- return nns;
-}
-
-static 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++)
- {
- 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)
-{
- int icg, nr_in_cg, maxcg, i, j, jcg, ngid, ncg;
-
- /* Compute largest charge groups size (# atoms) */
- nr_in_cg = 1;
- for (const gmx_moltype_t &molt : mtop->moltype)
- {
- const t_block *cgs = &molt.cgs;
- for (icg = 0; (icg < cgs->nr); icg++)
- {
- nr_in_cg = std::max(nr_in_cg, (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
- * 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);
- }
-
- ngid = mtop->groups.groups[SimulationAtomGroupType::EnergyOutput].nr;
- snew(ns->bExcludeAlleg, ngid);
- for (i = 0; i < ngid; i++)
- {
- ns->bExcludeAlleg[i] = TRUE;
- for (j = 0; j < ngid; j++)
- {
- if (!(fr->egp_flags[i*ngid+j] & EGP_EXCL))
- {
- ns->bExcludeAlleg[i] = FALSE;
- }
- }
- }
-
- if (fr->bGrid)
- {
- /* Grid search */
- 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[i], SHIFTS);
- }
- ncg = ncg_mtop(mtop);
- 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++)
- {
- 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))));
- }
- }
- if (debug)
- {
- pr_bvec(debug, 0, "bHaveVdW", ns->bHaveVdW, fr->ntype, TRUE);
- }
-
- ns->nra_alloc = 0;
- ns->bexcl = nullptr;
- if (!DOMAINDECOMP(cr))
- {
- ns_realloc_natoms(ns, mtop->natoms);
- }
-
- ns->nblist_initialized = FALSE;
-
- /* nbr list debug dump */
- {
- char *ptr = getenv("GMX_DUMP_NL");
- if (ptr)
- {
- ns->dump_nl = strtol(ptr, nullptr, 10);
- if (fplog)
- {
- fprintf(fplog, "GMX_DUMP_NL = %d", ns->dump_nl);
- }
- }
- else
- {
- ns->dump_nl = 0;
- }
- }
-}
-
-void done_ns(gmx_ns_t *ns, int numEnergyGroups)
-{
- if (ns->bexcl != nullptr)
- {
- sfree(ns->bexcl);
- }
- sfree(ns->bExcludeAlleg);
- if (ns->ns_buf)
- {
- for (int i = 0; i < numEnergyGroups; i++)
- {
- sfree(ns->ns_buf[i]);
- }
- sfree(ns->ns_buf);
- }
- if (ns->nl_sr != nullptr)
- {
- for (int i = 0; i < numEnergyGroups; i++)
- {
- if (ns->nl_sr[i] != nullptr)
- {
- sfree(ns->nl_sr[i]);
- }
- }
- sfree(ns->nl_sr);
- }
- sfree(ns->simple_aaj);
- sfree(ns->bHaveVdW);
- sfree(ns->nsr);
- done_grid(ns->grid);
- sfree(ns);
-}
-
-int search_neighbours(FILE *log,
- t_forcerec *fr,
- matrix box,
- gmx_localtop_t *top,
- const SimulationGroups *groups,
- const t_commrec *cr,
- t_nrnb *nrnb,
- const t_mdatoms *md,
- gmx_bool bFillGrid)
-{
- const t_block *cgs = &(top->cgs);
- rvec box_size, grid_x0, grid_x1;
- int m, ngid;
- real min_size, grid_dens;
- int nsearch;
- gmx_bool bGrid;
- int start, end;
- gmx_ns_t *ns;
- t_grid *grid;
- gmx_domdec_zones_t *dd_zones;
- put_in_list_t *put_in_list;
-
- ns = fr->ns;
-
- /* Set some local variables */
- bGrid = fr->bGrid;
- ngid = groups->groups[SimulationAtomGroupType::EnergyOutput].nr;
-
- for (m = 0; (m < DIM); m++)
- {
- box_size[m] = box[m][m];
- }
-
- if (fr->ePBC != epbcNONE)
- {
- if (gmx::square(fr->rlist) >= 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.");
- }
- if (!bGrid)
- {
- min_size = std::min(box_size[XX], std::min(box_size[YY], box_size[ZZ]));
- if (2*fr->rlist >= min_size)
- {
- 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]);
- }
-
- /* Reset the neighbourlists */
- reset_neighbor_lists(fr);
-
- if (bGrid && bFillGrid)
- {
-
- grid = ns->grid;
- if (DOMAINDECOMP(cr))
- {
- dd_zones = domdec_zones(cr->dd);
- }
- else
- {
- dd_zones = nullptr;
-
- get_nsgrid_boundaries(grid->nboundeddim, box, nullptr, nullptr, nullptr, nullptr,
- cgs->nr, fr->cg_cm, grid_x0, grid_x1, &grid_dens);
-
- grid_first(log, grid, nullptr, nullptr, box, grid_x0, grid_x1,
- fr->rlist, grid_dens);
- }
-
- start = 0;
- end = cgs->nr;
-
- if (DOMAINDECOMP(cr))
- {
- end = cgs->nr;
- 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(nullptr, grid, cgs->nr, fr->cg0, fr->hcg, fr->cg_cm);
- grid->icg0 = fr->cg0;
- grid->icg1 = fr->hcg;
- }
-
- calc_elemnr(grid, start, end, cgs->nr);
- calc_ptrs(grid);
- grid_last(grid, start, end, cgs->nr);
-
- if (gmx_debug_at)
- {
- 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(nullptr, ns->grid, fr->hcg, fr->hcg-1, fr->hcg, fr->cg_cm);
- }
-
- if (!fr->ns->bCGlist)
- {
- put_in_list = put_in_list_at;
- }
- else
- {
- put_in_list = put_in_list_cg;
- }
-
- /* Do the core! */
- if (bGrid)
- {
- grid = ns->grid;
- nsearch = nsgrid_core(cr, fr, box, ngid, top,
- grid, ns->bexcl, ns->bExcludeAlleg,
- md, put_in_list, ns->bHaveVdW,
- 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
- * (see QMMM.c). The VDW however, we'd like to compute classically
- * 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:
- */
- if (fr->bQMMM && fr->qr->QMMMscheme != eQMMMschemeoniom)
- {
- nsearch += nsgrid_core(cr, fr, box, ngid, top,
- grid, ns->bexcl, ns->bExcludeAlleg,
- md, put_in_list_qmmm, ns->bHaveVdW,
- TRUE);
- }
- }
- 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);
- }
-
- inc_nrnb(nrnb, eNR_NS, nsearch);
-
- return nsearch;
-}
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2018,2019, 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.
- *
- * 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 research papers on the package. Check out http://www.gromacs.org.
- */
-#ifndef GMX_MDLIB_NS_H
-#define GMX_MDLIB_NS_H
-
-#include <stdio.h>
-
-#include "gromacs/math/vectypes.h"
-#include "gromacs/utility/basedefinitions.h"
-
-struct SimulationGroups;
-struct gmx_localtop_t;
-struct gmx_mtop_t;
-struct gmx_ns_t;
-struct t_commrec;
-struct t_forcerec;
-struct t_mdatoms;
-struct t_nblist;
-struct t_nrnb;
-
-/****************************************************
- *
- * U T I L I T I E S May be found in ns.c
- *
- ****************************************************/
-
-void init_neighbor_list(FILE *log, t_forcerec *fr, int homenr);
-/*
- * nn is the number of energy terms in the energy matrix
- * (ngener*(ngener-1))/2
- * start is the first atom on this processor
- * homenr is the number of atoms on this processor
- */
-
-int calc_naaj(int icg, int cgtot);
-/* Calculate the number of charge groups to interact with for icg */
-
-/****************************************************
- *
- * N E I G H B O R S E A R C H I N G
- *
- * Calls either ns5_core (when grid selected in .mdp file)
- * or ns_simple_core (when simple selected in .mdp file)
- *
- * Return total number of pairs searched
- *
- ****************************************************/
-void init_ns(FILE *fplog, const t_commrec *cr,
- gmx_ns_t *ns, t_forcerec *fr,
- const gmx_mtop_t *mtop);
-
-//! Destructor.
-void done_ns(gmx_ns_t *ns, int numEnergyGroups);
-
-int search_neighbours(FILE *log,
- t_forcerec *fr,
- matrix box,
- gmx_localtop_t *top,
- const SimulationGroups *groups,
- const t_commrec *cr,
- t_nrnb *nrnb,
- const t_mdatoms *md,
- gmx_bool bFillGrid);
-
-
-/* Debugging routines from wnblist.c */
-void dump_nblist(FILE *out, const t_commrec *cr, t_forcerec *fr, int nDNL);
-
-int read_nblist(FILE *in, FILE *out, int **mat, int natoms, gmx_bool bSymm);
-/* Returns total number of neighbors. If bSymm the matrix is symmetrized. */
-
-void reallocate_nblist(t_nblist *nl);
-/* List reallocation, only exported for Verlet scheme use with FEP */
-
-#endif
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2017,2018,2019, 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.
#include "gromacs/gmxlib/nrnb.h"
#include "gromacs/math/units.h"
#include "gromacs/math/vec.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/forcerec.h"
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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.
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/forcerec.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/utility/cstringutil.h"
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2017,2018,2019, 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.
#include "gromacs/gmxlib/nrnb.h"
#include "gromacs/math/units.h"
#include "gromacs/math/vec.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/utility/fatalerror.h"
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2017,2018,2019, 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.
#include "gromacs/gmxlib/nrnb.h"
#include "gromacs/math/units.h"
#include "gromacs/math/vec.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/qmmm.h"
#include "gromacs/mdtypes/forcerec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/math/functions.h"
#include "gromacs/math/units.h"
#include "gromacs/math/vec.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/qm_gamess.h"
#include "gromacs/mdlib/qm_gaussian.h"
#include "gromacs/mdlib/qm_mopac.h"
#include <gtest/gtest.h>
#include "gromacs/fileio/gmxfio.h"
+#include "gromacs/gmxlib/nrnb.h"
#include "gromacs/gmxlib/nonbonded/nonbonded.h"
#include "gromacs/math/paddedvector.h"
#include "gromacs/math/vec.h"
+++ /dev/null
-/*
- * This file is part of the GROMACS molecular simulation package.
- *
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2018, 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.
- *
- * 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 research papers on the package. Check out http://www.gromacs.org.
- */
-/* This file is completely threadsafe - keep it that way! */
-#include "gmxpre.h"
-
-#include <cstdio>
-#include <cstring>
-
-#include <algorithm>
-
-#include "gromacs/domdec/domdec.h"
-#include "gromacs/domdec/domdec_struct.h"
-#include "gromacs/fileio/gmxfio.h"
-#include "gromacs/gmxlib/nrnb.h"
-#include "gromacs/mdlib/ns.h"
-#include "gromacs/mdtypes/commrec.h"
-#include "gromacs/mdtypes/forcerec.h"
-#include "gromacs/mdtypes/md_enums.h"
-#include "gromacs/mdtypes/nblist.h"
-#include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/futil.h"
-#include "gromacs/utility/smalloc.h"
-
-static void write_nblist(FILE *out, gmx_domdec_t *dd, t_nblist *nblist, int nDNL)
-{
- int i, nii, ii, j, zi, zj0, zj1, aj, zj, nj;
- int ca1[DD_MAXZONE], np[DD_MAXZONE];
- gmx_domdec_zones_t *dd_zones;
-
- if (nblist->nri > 0)
- {
- fprintf(out, "elec: %s, vdw: %s, type: %s, Solvent opt: %s\n",
- gmx_nbkernel_elec_names[nblist->ielec],
- gmx_nbkernel_vdw_names[nblist->ivdw],
- gmx_nblist_interaction_names[nblist->type],
- gmx_nblist_geometry_names[nblist->igeometry]);
- fprintf(out, "nri: %d npair: %d\n", nblist->nri, nblist->nrj);
- if (dd)
- {
- dd_zones = domdec_zones(dd);
-
- for (zi = 0; zi < dd_zones->n; zi++)
- {
- ca1[zi] = dd->atomGrouping().block(dd_zones->cg_range[zi + 1]).begin();
- }
- i = 0;
- for (zi = 0; zi < dd_zones->nizone && zi < dd_zones->n; zi++)
- {
- zj0 = dd_zones->izone[zi].j0;
- zj1 = dd_zones->izone[zi].j1;
- for (zj = zj0; zj < zj1; zj++)
- {
- np[zj] = 0;
- }
- while (i < nblist->nri && nblist->iinr[i] < ca1[zi])
- {
- for (j = nblist->jindex[i]; (j < nblist->jindex[i+1]); j++)
- {
- aj = nblist->jjnr[j];
- zj = zj0;
- while (aj >= ca1[zj])
- {
- zj++;
- }
- np[zj]++;
- }
- i++;
- }
- fprintf(out, "DD zone %d:", zi);
- for (zj = zj0; zj < zj1; zj++)
- {
- fprintf(out, " %d %d", zj, np[zj]);
- }
- fprintf(out, "\n");
- }
- }
- if (nDNL >= 2)
- {
- for (i = 0; i < nblist->nri; i++)
- {
- nii = 1;
- if (nDNL >= 3 && nblist->igeometry != GMX_NBLIST_GEOMETRY_PARTICLE_PARTICLE)
- {
- nii = 3;
- }
- nj = nblist->jindex[i+1] - nblist->jindex[i];
- fprintf(out, "i: %d shift: %d gid: %d nj: %d\n",
- ddglatnr(dd, nblist->iinr[i]),
- nblist->shift[i], nblist->gid[i], nj);
- for (ii = 0; ii < nii; ii++)
- {
- for (j = nblist->jindex[i]; (j < nblist->jindex[i+1]); j++)
- {
- fprintf(out, " i: %5d j: %5d\n",
- ddglatnr(dd, nblist->iinr[i]+ii),
- ddglatnr(dd, nblist->jjnr[j]));
- }
- }
- }
- }
- fflush(out);
- }
-}
-
-
-
-void dump_nblist(FILE *out, const t_commrec *cr, t_forcerec *fr, int nDNL)
-{
- int n, i;
-
- fprintf(out, "Neighborlist:\n");
-
- for (n = 0; (n < fr->nnblists); n++)
- {
- for (i = 0; (i < eNL_NR); i++)
- {
- write_nblist(out, cr->dd, &fr->nblists[n].nlist_sr[i], nDNL);
- }
- }
-
-}
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/mdoutf.h"
#include "gromacs/mdlib/membed.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/resethandler.h"
#include "gromacs/mdlib/sighandler.h"
#include "gromacs/mdlib/simulationsignal.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/mdoutf.h"
#include "gromacs/mdlib/membed.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/resethandler.h"
#include "gromacs/mdlib/sighandler.h"
#include "gromacs/mdlib/simulationsignal.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/mdatoms.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/stat.h"
#include "gromacs/mdlib/tgroup.h"
#include "gromacs/mdlib/trajectory_writing.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/mdoutf.h"
#include "gromacs/mdlib/membed.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/resethandler.h"
#include "gromacs/mdlib/sighandler.h"
#include "gromacs/mdlib/simulationsignal.h"
free_gpu_resources(fr, physicalNodeComm);
free_gpu(nonbondedDeviceInfo);
free_gpu(pmeDeviceInfo);
- done_forcerec(fr, mtop.molblock.size(), mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].nr);
+ done_forcerec(fr, mtop.molblock.size());
sfree(fcd);
if (doMembed)
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/force_flags.h"
#include "gromacs/mdlib/mdatoms.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdlib/tgroup.h"
#include "gromacs/mdlib/update.h"
#include "gromacs/mdlib/vsite.h"
struct nonbonded_verlet_t;
struct bonded_threading_t;
struct t_forcetable;
-struct t_nblist;
-struct t_nblists;
struct t_QMMMrec;
namespace gmx
* The maximum cg size in cginfo is 63
* because we only have space for 6 bits in cginfo,
* this cg size entry is actually only read with domain decomposition.
- * But there is a smaller limit due to the t_excl data structure
- * which is defined in nblist.h.
*/
#define SET_CGINFO_GID(cgi, gid) (cgi) = (((cgi) & ~255) | (gid))
#define GET_CGINFO_GID(cgi) ( (cgi) & 255)
/* Table stuff */
gmx_bool bcoultab = FALSE;
gmx_bool bvdwtab = FALSE;
- /* The normal tables are in the nblists struct(s) below */
struct t_forcetable *pairsTable; /* for 1-4 interactions, [pairs] and [pairs_nb] */
int cg_nalloc = 0;
rvec *shift_vec = nullptr;
- /* The neighborlists including tables */
- int nnblists = 0;
- int *gid2nblists = nullptr;
- struct t_nblists *nblists = nullptr;
-
- int cutoff_scheme = 0; /* group- or Verlet-style cutoff */
- gmx_bool bNonbonded = FALSE; /* true if nonbonded calculations are *not* turned off */
+ int cutoff_scheme = 0; /* group- or Verlet-style cutoff */
+ gmx_bool bNonbonded = FALSE; /* true if nonbonded calculations are *not* turned off */
/* The Nbnxm Verlet non-bonded machinery */
std::unique_ptr<nonbonded_verlet_t> nbv;
*/
int n_tpi = 0;
- /* Neighbor searching stuff */
- struct gmx_ns_t *ns = nullptr;
-
/* QMMM stuff */
gmx_bool bQMMM = FALSE;
struct t_QMMMrec *qr = nullptr;
return domainSetup_;
}
- //! Returns the number of cells along x and y for the local grid
- void getLocalNumCells(int *numCellsX,
- int *numCellsY) const
- {
- *numCellsX = grids_[0].dimensions().numCells[XX];
- *numCellsY = grids_[0].dimensions().numCells[YY];
- }
-
//! Returns the total number of atoms in the grid set, including padding
int numGridAtomsTotal() const
{
wallcycle_stop(wcycle, ewcNB_XF_BUF_OPS);
}
-void nonbonded_verlet_t::getLocalNumCells(int *numCellsX,
- int *numCellsY) const
-{
- pairSearch_->gridSet().getLocalNumCells(numCellsX, numCellsY);
-}
-
gmx::ArrayRef<const int> nonbonded_verlet_t::getGridIndices() const
{
return pairSearch_->gridSet().cells();
return kernelSetup_;
}
- /*! \brief Returns the number of x and y cells in the local grid */
- void getLocalNumCells(int *numCellsX,
- int *numCellsY) const;
-
-
//! Returns the outer radius for the pair list
real pairlistInnerRadius() const;
#include "gromacs/math/utilities.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/gmx_omp_nthreads.h"
-#include "gromacs/mdlib/ns.h"
#include "gromacs/mdtypes/group.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/nbnxm/atomdata.h"
nlist->jindex[nlist->nri] = nlist->nrj;
}
+/* Rellocate FEP list for size nl->maxnri, TODO: replace by C++ */
+static void reallocate_nblist(t_nblist *nl)
+{
+ if (gmx_debug_at)
+ {
+ 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);
+ srenew(nl->gid, nl->maxnri);
+ srenew(nl->shift, nl->maxnri);
+ srenew(nl->jindex, nl->maxnri+1);
+}
+
/* For load balancing of the free-energy lists over threads, we set
* the maximum nrj size of an i-entry to 40. This leads to good
* load balancing in the worst case scenario of a single perturbed
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff