* To help us fund GROMACS development, we humbly ask that you cite
* the research papers on the package. Check out http://www.gromacs.org.
*/
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "gmxpre.h"
+
+#include "gromacs/legacyheaders/nonbonded.h"
+
+#include "config.h"
#include <stdio.h>
#include <stdlib.h>
#include "thread_mpi/threads.h"
-#include "typedefs.h"
-#include "txtdump.h"
-#include "ns.h"
-#include "gromacs/math/vec.h"
+#include "gromacs/gmxlib/nonbonded/nb_free_energy.h"
+#include "gromacs/gmxlib/nonbonded/nb_generic.h"
+#include "gromacs/gmxlib/nonbonded/nb_generic_adress.h"
+#include "gromacs/gmxlib/nonbonded/nb_generic_cg.h"
+#include "gromacs/gmxlib/nonbonded/nb_kernel.h"
+#include "gromacs/legacyheaders/force.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/nrnb.h"
+#include "gromacs/legacyheaders/ns.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/listed-forces/bonded.h"
#include "gromacs/math/utilities.h"
-#include "macros.h"
-#include "gromacs/utility/cstringutil.h"
-#include "force.h"
-#include "names.h"
-#include "force.h"
-#include "bondf.h"
-#include "nrnb.h"
-#include "nonbonded.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/ishift.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
#include "gromacs/simd/simd.h"
-
-#include "mshift.h"
+#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
-#include "nb_kernel.h"
-#include "nb_free_energy.h"
-#include "nb_generic.h"
-#include "nb_generic_cg.h"
-#include "nb_generic_adress.h"
-
/* Different default (c) and SIMD instructions interaction-specific kernels */
-#include "nb_kernel_c/nb_kernel_c.h"
+#include "gromacs/gmxlib/nonbonded/nb_kernel_c/nb_kernel_c.h"
#if (defined GMX_SIMD_X86_SSE2) && !(defined GMX_DOUBLE)
-# include "nb_kernel_sse2_single/nb_kernel_sse2_single.h"
+# include "gromacs/gmxlib/nonbonded/nb_kernel_sse2_single/nb_kernel_sse2_single.h"
#endif
#if (defined GMX_SIMD_X86_SSE4_1) && !(defined GMX_DOUBLE)
-# include "nb_kernel_sse4_1_single/nb_kernel_sse4_1_single.h"
+# include "gromacs/gmxlib/nonbonded/nb_kernel_sse4_1_single/nb_kernel_sse4_1_single.h"
#endif
#if (defined GMX_SIMD_X86_AVX_128_FMA) && !(defined GMX_DOUBLE)
-# include "nb_kernel_avx_128_fma_single/nb_kernel_avx_128_fma_single.h"
+# include "gromacs/gmxlib/nonbonded/nb_kernel_avx_128_fma_single/nb_kernel_avx_128_fma_single.h"
#endif
#if (defined GMX_SIMD_X86_AVX_256_OR_HIGHER) && !(defined GMX_DOUBLE)
-# include "nb_kernel_avx_256_single/nb_kernel_avx_256_single.h"
+# include "gromacs/gmxlib/nonbonded/nb_kernel_avx_256_single/nb_kernel_avx_256_single.h"
#endif
#if (defined GMX_SIMD_X86_SSE2 && defined GMX_DOUBLE)
-# include "nb_kernel_sse2_double/nb_kernel_sse2_double.h"
+# include "gromacs/gmxlib/nonbonded/nb_kernel_sse2_double/nb_kernel_sse2_double.h"
#endif
#if (defined GMX_SIMD_X86_SSE4_1 && defined GMX_DOUBLE)
-# include "nb_kernel_sse4_1_double/nb_kernel_sse4_1_double.h"
+# include "gromacs/gmxlib/nonbonded/nb_kernel_sse4_1_double/nb_kernel_sse4_1_double.h"
#endif
#if (defined GMX_SIMD_X86_AVX_128_FMA && defined GMX_DOUBLE)
-# include "nb_kernel_avx_128_fma_double/nb_kernel_avx_128_fma_double.h"
+# include "gromacs/gmxlib/nonbonded/nb_kernel_avx_128_fma_double/nb_kernel_avx_128_fma_double.h"
#endif
#if (defined GMX_SIMD_X86_AVX_256_OR_HIGHER && defined GMX_DOUBLE)
-# include "nb_kernel_avx_256_double/nb_kernel_avx_256_double.h"
+# include "gromacs/gmxlib/nonbonded/nb_kernel_avx_256_double/nb_kernel_avx_256_double.h"
#endif
#if (defined GMX_SIMD_SPARC64_HPC_ACE && defined GMX_DOUBLE)
-# include "nb_kernel_sparc64_hpc_ace_double/nb_kernel_sparc64_hpc_ace_double.h"
+# include "gromacs/gmxlib/nonbonded/nb_kernel_sparc64_hpc_ace_double/nb_kernel_sparc64_hpc_ace_double.h"
#endif
void
-gmx_nonbonded_set_kernel_pointers(FILE *log, t_nblist *nl)
+gmx_nonbonded_set_kernel_pointers(FILE *log, t_nblist *nl, gmx_bool bElecAndVdwSwitchDiffers)
{
const char * elec;
const char * elec_mod;
}
}
- /* Give up. If this was a water kernel, leave the pointer as NULL, which
- * will disable water optimization in NS. If it is a particle kernel, set
- * the pointer to the generic NB kernel.
+ /* For now, the accelerated kernels cannot handle the combination of switch functions for both
+ * electrostatics and VdW that use different switch radius or switch cutoff distances
+ * (both of them enter in the switch function calculation). This would require
+ * us to evaluate two completely separate switch functions for every interaction.
+ * Instead, we disable such kernels by setting the pointer to NULL.
+ * This will cause the generic kernel (which can handle it) to be called instead.
+ *
+ * Note that we typically already enable tabulated coulomb interactions for this case,
+ * so this is mostly a safe-guard to make sure we call the generic kernel if the
+ * tables are disabled.
+ */
+ if ((nl->ielec != GMX_NBKERNEL_ELEC_NONE) && (nl->ielecmod == eintmodPOTSWITCH) &&
+ (nl->ivdw != GMX_NBKERNEL_VDW_NONE) && (nl->ivdwmod == eintmodPOTSWITCH) &&
+ bElecAndVdwSwitchDiffers)
+ {
+ nl->kernelptr_vf = NULL;
+ nl->kernelptr_f = NULL;
+ }
+
+ /* Give up, pick a generic one instead.
+ * 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 (nl->kernelptr_vf == NULL && !gmx_strcasecmp_min(geom, "Particle-Particle"))
{
fprintf(debug,
"WARNING - Slow generic NB kernel used for neighborlist with\n"
" Elec: '%s', Modifier: '%s'\n"
- " Vdw: '%s', Modifier: '%s'\n"
- " Geom: '%s', Other: '%s'\n\n",
- elec, elec_mod, vdw, vdw_mod, geom, other);
+ " Vdw: '%s', Modifier: '%s'\n",
+ elec, elec_mod, vdw, vdw_mod);
}
}
}
-
return;
}
{
(*kernelptr)(&(nlist[i]), x, f, fr, mdatoms, &kernel_data, nrnb);
}
+ else
+ {
+ gmx_fatal(FARGS, "Non-empty neighborlist does not have any kernel pointer assigned.");
+ }
}
}
}
}
}
-
-static void
-nb_listed_warning_rlimit(const rvec *x, int ai, int aj, int * global_atom_index, real r, real rlimit)
-{
- gmx_warning("Listed nonbonded interaction between particles %d and %d\n"
- "at distance %.3f which is larger than the table limit %.3f nm.\n\n"
- "This is likely either a 1,4 interaction, or a listed interaction inside\n"
- "a smaller molecule you are decoupling during a free energy calculation.\n"
- "Since interactions at distances beyond the table cannot be computed,\n"
- "they are skipped until they are inside the table limit again. You will\n"
- "only see this message once, even if it occurs for several interactions.\n\n"
- "IMPORTANT: This should not happen in a stable simulation, so there is\n"
- "probably something wrong with your system. Only change the table-extension\n"
- "distance in the mdp file if you are really sure that is the reason.\n",
- glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r, rlimit);
-
- if (debug)
- {
- fprintf(debug,
- "%8f %8f %8f\n%8f %8f %8f\n1-4 (%d,%d) interaction not within cut-off! r=%g. Ignored\n",
- x[ai][XX], x[ai][YY], x[ai][ZZ], x[aj][XX], x[aj][YY], x[aj][ZZ],
- glatnr(global_atom_index, ai), glatnr(global_atom_index, aj), r);
- }
-}
-
-
-
-/* This might logically belong better in the nb_generic.c module, but it is only
- * used in do_nonbonded_listed(), and we want it to be inlined there to avoid an
- * extra functional call for every single pair listed in the topology.
- */
-static real
-nb_evaluate_single(real r2, real tabscale, real *vftab,
- real qq, real c6, real c12, real *velec, real *vvdw)
-{
- real rinv, r, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VVe, FFe, VVd, FFd, VVr, FFr, fscal;
- int ntab;
-
- /* Do the tabulated interactions - first table lookup */
- rinv = gmx_invsqrt(r2);
- r = r2*rinv;
- rtab = r*tabscale;
- ntab = rtab;
- eps = rtab-ntab;
- eps2 = eps*eps;
- ntab = 12*ntab;
- /* Electrostatics */
- Y = vftab[ntab];
- F = vftab[ntab+1];
- Geps = eps*vftab[ntab+2];
- Heps2 = eps2*vftab[ntab+3];
- Fp = F+Geps+Heps2;
- VVe = Y+eps*Fp;
- FFe = Fp+Geps+2.0*Heps2;
- /* Dispersion */
- Y = vftab[ntab+4];
- F = vftab[ntab+5];
- Geps = eps*vftab[ntab+6];
- Heps2 = eps2*vftab[ntab+7];
- Fp = F+Geps+Heps2;
- VVd = Y+eps*Fp;
- FFd = Fp+Geps+2.0*Heps2;
- /* Repulsion */
- Y = vftab[ntab+8];
- F = vftab[ntab+9];
- Geps = eps*vftab[ntab+10];
- Heps2 = eps2*vftab[ntab+11];
- Fp = F+Geps+Heps2;
- VVr = Y+eps*Fp;
- FFr = Fp+Geps+2.0*Heps2;
-
- *velec = qq*VVe;
- *vvdw = c6*VVd+c12*VVr;
-
- fscal = -(qq*FFe+c6*FFd+c12*FFr)*tabscale*rinv;
-
- return fscal;
-}
-
-
-real
-do_nonbonded_listed(int ftype, int nbonds,
- const t_iatom iatoms[], const t_iparams iparams[],
- const rvec x[], rvec f[], rvec fshift[],
- const t_pbc *pbc, const t_graph *g,
- real *lambda, real *dvdl,
- const t_mdatoms *md,
- const t_forcerec *fr, gmx_grppairener_t *grppener,
- int *global_atom_index)
-{
- int ielec, ivdw;
- real qq, c6, c12;
- rvec dx;
- ivec dt;
- int i, j, itype, ai, aj, gid;
- int fshift_index;
- real r2, rinv;
- real fscal, velec, vvdw;
- real * energygrp_elec;
- real * energygrp_vdw;
- static gmx_bool warned_rlimit = FALSE;
- /* Free energy stuff */
- gmx_bool bFreeEnergy;
- real LFC[2], LFV[2], DLF[2], lfac_coul[2], lfac_vdw[2], dlfac_coul[2], dlfac_vdw[2];
- real qqB, c6B, c12B, sigma2_def, sigma2_min;
-
-
- switch (ftype)
- {
- case F_LJ14:
- case F_LJC14_Q:
- energygrp_elec = grppener->ener[egCOUL14];
- energygrp_vdw = grppener->ener[egLJ14];
- break;
- case F_LJC_PAIRS_NB:
- energygrp_elec = grppener->ener[egCOULSR];
- energygrp_vdw = grppener->ener[egLJSR];
- break;
- default:
- energygrp_elec = NULL; /* Keep compiler happy */
- energygrp_vdw = NULL; /* Keep compiler happy */
- gmx_fatal(FARGS, "Unknown function type %d in do_nonbonded14", ftype);
- break;
- }
-
- if (fr->efep != efepNO)
- {
- /* Lambda factor for state A=1-lambda and B=lambda */
- LFC[0] = 1.0 - lambda[efptCOUL];
- LFV[0] = 1.0 - lambda[efptVDW];
- LFC[1] = lambda[efptCOUL];
- LFV[1] = lambda[efptVDW];
-
- /*derivative of the lambda factor for state A and B */
- DLF[0] = -1;
- DLF[1] = 1;
-
- /* precalculate */
- sigma2_def = pow(fr->sc_sigma6_def, 1.0/3.0);
- sigma2_min = pow(fr->sc_sigma6_min, 1.0/3.0);
-
- for (i = 0; i < 2; i++)
- {
- lfac_coul[i] = (fr->sc_power == 2 ? (1-LFC[i])*(1-LFC[i]) : (1-LFC[i]));
- dlfac_coul[i] = DLF[i]*fr->sc_power/fr->sc_r_power*(fr->sc_power == 2 ? (1-LFC[i]) : 1);
- lfac_vdw[i] = (fr->sc_power == 2 ? (1-LFV[i])*(1-LFV[i]) : (1-LFV[i]));
- dlfac_vdw[i] = DLF[i]*fr->sc_power/fr->sc_r_power*(fr->sc_power == 2 ? (1-LFV[i]) : 1);
- }
- }
- else
- {
- sigma2_min = sigma2_def = 0;
- }
-
- bFreeEnergy = FALSE;
- for (i = 0; (i < nbonds); )
- {
- itype = iatoms[i++];
- ai = iatoms[i++];
- aj = iatoms[i++];
- gid = GID(md->cENER[ai], md->cENER[aj], md->nenergrp);
-
- /* Get parameters */
- switch (ftype)
- {
- case F_LJ14:
- bFreeEnergy =
- (fr->efep != efepNO &&
- ((md->nPerturbed && (md->bPerturbed[ai] || md->bPerturbed[aj])) ||
- iparams[itype].lj14.c6A != iparams[itype].lj14.c6B ||
- iparams[itype].lj14.c12A != iparams[itype].lj14.c12B));
- qq = md->chargeA[ai]*md->chargeA[aj]*fr->epsfac*fr->fudgeQQ;
- c6 = iparams[itype].lj14.c6A;
- c12 = iparams[itype].lj14.c12A;
- break;
- case F_LJC14_Q:
- qq = iparams[itype].ljc14.qi*iparams[itype].ljc14.qj*fr->epsfac*iparams[itype].ljc14.fqq;
- c6 = iparams[itype].ljc14.c6;
- c12 = iparams[itype].ljc14.c12;
- break;
- case F_LJC_PAIRS_NB:
- qq = iparams[itype].ljcnb.qi*iparams[itype].ljcnb.qj*fr->epsfac;
- c6 = iparams[itype].ljcnb.c6;
- c12 = iparams[itype].ljcnb.c12;
- break;
- default:
- /* Cannot happen since we called gmx_fatal() above in this case */
- qq = c6 = c12 = 0; /* Keep compiler happy */
- break;
- }
-
- /* To save flops in the optimized kernels, c6/c12 have 6.0/12.0 derivative prefactors
- * included in the general nfbp array now. This means the tables are scaled down by the
- * same factor, so when we use the original c6/c12 parameters from iparams[] they must
- * be scaled up.
- */
- c6 *= 6.0;
- c12 *= 12.0;
-
- /* Do we need to apply full periodic boundary conditions? */
- if (fr->bMolPBC == TRUE)
- {
- fshift_index = pbc_dx_aiuc(pbc, x[ai], x[aj], dx);
- }
- else
- {
- fshift_index = CENTRAL;
- rvec_sub(x[ai], x[aj], dx);
- }
- r2 = norm2(dx);
-
- if (r2 >= fr->tab14.r*fr->tab14.r)
- {
- if (warned_rlimit == FALSE)
- {
- nb_listed_warning_rlimit(x, ai, aj, global_atom_index, sqrt(r2), fr->tab14.r);
- warned_rlimit = TRUE;
- }
- continue;
- }
-
- if (bFreeEnergy)
- {
- /* Currently free energy is only supported for F_LJ14, so no need to check for that if we got here */
- qqB = md->chargeB[ai]*md->chargeB[aj]*fr->epsfac*fr->fudgeQQ;
- c6B = iparams[itype].lj14.c6B*6.0;
- c12B = iparams[itype].lj14.c12B*12.0;
-
- fscal = nb_free_energy_evaluate_single(r2, fr->sc_r_power, fr->sc_alphacoul, fr->sc_alphavdw,
- fr->tab14.scale, fr->tab14.data, qq, c6, c12, qqB, c6B, c12B,
- LFC, LFV, DLF, lfac_coul, lfac_vdw, dlfac_coul, dlfac_vdw,
- fr->sc_sigma6_def, fr->sc_sigma6_min, sigma2_def, sigma2_min, &velec, &vvdw, dvdl);
- }
- else
- {
- /* Evaluate tabulated interaction without free energy */
- fscal = nb_evaluate_single(r2, fr->tab14.scale, fr->tab14.data, qq, c6, c12, &velec, &vvdw);
- }
-
- energygrp_elec[gid] += velec;
- energygrp_vdw[gid] += vvdw;
- svmul(fscal, dx, dx);
-
- /* Add the forces */
- rvec_inc(f[ai], dx);
- rvec_dec(f[aj], dx);
-
- if (g)
- {
- /* Correct the shift forces using the graph */
- ivec_sub(SHIFT_IVEC(g, ai), SHIFT_IVEC(g, aj), dt);
- fshift_index = IVEC2IS(dt);
- }
- if (fshift_index != CENTRAL)
- {
- rvec_inc(fshift[fshift_index], dx);
- rvec_dec(fshift[CENTRAL], dx);
- }
- }
- return 0.0;
-}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff