From a1fbbb07241495656429dbb560e93769368609a4 Mon Sep 17 00:00:00 2001 From: Mark Abraham Date: Tue, 23 Sep 2014 06:50:34 +0200 Subject: [PATCH] Unite code for handling listed pairs Some implementation details are closer to the other listed interactions than the group-kernel non-bondeds. This functionality will also survive the death of the group kernels, so may as well be separated now. This is pure code motion, to prepare for transition to the listed-forces module. Had to add an #include guard to a not-really-related file. Some not-quite related clean-up in nonbonded.h Added temporary cycle suppression until this moves into the listed-forces module Change-Id: Ie8865f5ca1a4512d75e70e4b7f25b90f6ed019ec --- docs/doxygen/cycle-suppressions.txt | 1 + src/gromacs/gmxlib/nonbonded/nb_free_energy.c | 195 ------- src/gromacs/gmxlib/nonbonded/nb_free_energy.h | 24 +- src/gromacs/gmxlib/nonbonded/nonbonded.c | 263 --------- src/gromacs/gmxlib/nonbonded/pairs.c | 507 ++++++++++++++++++ src/gromacs/gmxlib/nonbonded/pairs.h | 71 +++ src/gromacs/legacyheaders/nonbonded.h | 17 - src/gromacs/legacyheaders/types/group.h | 5 +- src/gromacs/listed-forces/bonded.cpp | 2 +- 9 files changed, 591 insertions(+), 494 deletions(-) create mode 100644 src/gromacs/gmxlib/nonbonded/pairs.c create mode 100644 src/gromacs/gmxlib/nonbonded/pairs.h diff --git a/docs/doxygen/cycle-suppressions.txt b/docs/doxygen/cycle-suppressions.txt index b87e5469ca..c49d9ed7ea 100644 --- a/docs/doxygen/cycle-suppressions.txt +++ b/docs/doxygen/cycle-suppressions.txt @@ -10,3 +10,4 @@ topology -> fileio topology -> legacyheaders pbcutil -> fileio pbcutil -> legacyheaders +listed-forces -> gmxlib diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.c b/src/gromacs/gmxlib/nonbonded/nb_free_energy.c index 6b5d6d89dc..b4a0200fd1 100644 --- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.c +++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.c @@ -935,198 +935,3 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, #pragma omp atomic inc_nrnb(nrnb, eNR_NBKERNEL_FREE_ENERGY, nlist->nri*12 + nlist->jindex[n]*150); } - -real -nb_free_energy_evaluate_single(real r2, real sc_r_power, real alpha_coul, real alpha_vdw, - real tabscale, real *vftab, - real qqA, real c6A, real c12A, real qqB, real c6B, real c12B, - real LFC[2], real LFV[2], real DLF[2], - real lfac_coul[2], real lfac_vdw[2], real dlfac_coul[2], real dlfac_vdw[2], - real sigma6_def, real sigma6_min, real sigma2_def, real sigma2_min, - real *velectot, real *vvdwtot, real *dvdl) -{ - real r, rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal; - real qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2], sigma_powm2[2]; - real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw; - real rpinv, r_coul, r_vdw, velecsum, vvdwsum; - real fscal_vdw[2], fscal_elec[2]; - real velec[2], vvdw[2]; - int i, ntab; - - const real half = 0.5; - const real one = 1.0; - const real two = 2.0; - const real six = 6.0; - const real fourtyeight = 48.0; - - qq[0] = qqA; - qq[1] = qqB; - c6[0] = c6A; - c6[1] = c6B; - c12[0] = c12A; - c12[1] = c12B; - - if (sc_r_power == six) - { - rpm2 = r2*r2; /* r4 */ - rp = rpm2*r2; /* r6 */ - } - else if (sc_r_power == fourtyeight) - { - rp = r2*r2*r2; /* r6 */ - rp = rp*rp; /* r12 */ - rp = rp*rp; /* r24 */ - rp = rp*rp; /* r48 */ - rpm2 = rp/r2; /* r46 */ - } - else - { - rp = pow(r2, half*sc_r_power); /* not currently supported as input, but can handle it */ - rpm2 = rp/r2; - } - - /* Loop over state A(0) and B(1) */ - for (i = 0; i < 2; i++) - { - if ((c6[i] > 0) && (c12[i] > 0)) - { - /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively. - * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5. - */ - sigma6[i] = half*c12[i]/c6[i]; - sigma2[i] = pow(half*c12[i]/c6[i], 1.0/3.0); - /* should be able to get rid of this ^^^ internal pow call eventually. Will require agreement on - what data to store externally. Can't be fixed without larger scale changes, so not 5.0 */ - if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */ - { - sigma6[i] = sigma6_min; - sigma2[i] = sigma2_min; - } - } - else - { - sigma6[i] = sigma6_def; - sigma2[i] = sigma2_def; - } - if (sc_r_power == six) - { - sigma_pow[i] = sigma6[i]; - sigma_powm2[i] = sigma6[i]/sigma2[i]; - } - else if (sc_r_power == fourtyeight) - { - sigma_pow[i] = sigma6[i]*sigma6[i]; /* sigma^12 */ - sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */ - sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */ - sigma_powm2[i] = sigma_pow[i]/sigma2[i]; - } - else - { /* not really supported as input, but in here for testing the general case*/ - sigma_pow[i] = pow(sigma2[i], sc_r_power/2); - sigma_powm2[i] = sigma_pow[i]/(sigma2[i]); - } - } - - /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/ - if ((c12[0] > 0) && (c12[1] > 0)) - { - alpha_vdw_eff = 0; - alpha_coul_eff = 0; - } - else - { - alpha_vdw_eff = alpha_vdw; - alpha_coul_eff = alpha_coul; - } - - /* Loop over A and B states again */ - for (i = 0; i < 2; i++) - { - fscal_elec[i] = 0; - fscal_vdw[i] = 0; - velec[i] = 0; - vvdw[i] = 0; - - /* Only spend time on A or B state if it is non-zero */ - if ( (qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0) ) - { - /* Coulomb */ - rpinv = one/(alpha_coul_eff*lfac_coul[i]*sigma_pow[i]+rp); - r_coul = pow(rpinv, -one/sc_r_power); - - /* Electrostatics table lookup data */ - rtab = r_coul*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; - VV = Y+eps*Fp; - FF = Fp+Geps+two*Heps2; - velec[i] = qq[i]*VV; - fscal_elec[i] = -qq[i]*FF*r_coul*rpinv*tabscale; - - /* Vdw */ - rpinv = one/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp); - r_vdw = pow(rpinv, -one/sc_r_power); - /* Vdw table lookup data */ - rtab = r_vdw*tabscale; - ntab = rtab; - eps = rtab-ntab; - eps2 = eps*eps; - ntab = 12*ntab; - /* Dispersion */ - Y = vftab[ntab+4]; - F = vftab[ntab+5]; - Geps = eps*vftab[ntab+6]; - Heps2 = eps2*vftab[ntab+7]; - Fp = F+Geps+Heps2; - VV = Y+eps*Fp; - FF = Fp+Geps+two*Heps2; - vvdw[i] = c6[i]*VV; - fscal_vdw[i] = -c6[i]*FF; - - /* Repulsion */ - Y = vftab[ntab+8]; - F = vftab[ntab+9]; - Geps = eps*vftab[ntab+10]; - Heps2 = eps2*vftab[ntab+11]; - Fp = F+Geps+Heps2; - VV = Y+eps*Fp; - FF = Fp+Geps+two*Heps2; - vvdw[i] += c12[i]*VV; - fscal_vdw[i] -= c12[i]*FF; - fscal_vdw[i] *= r_vdw*rpinv*tabscale; - } - } - /* Now we have velec[i], vvdw[i], and fscal[i] for both states */ - /* Assemble A and B states */ - velecsum = 0; - vvdwsum = 0; - dvdl_coul = 0; - dvdl_vdw = 0; - fscal = 0; - for (i = 0; i < 2; i++) - { - velecsum += LFC[i]*velec[i]; - vvdwsum += LFV[i]*vvdw[i]; - - fscal += (LFC[i]*fscal_elec[i]+LFV[i]*fscal_vdw[i])*rpm2; - - dvdl_coul += velec[i]*DLF[i] + LFC[i]*alpha_coul_eff*dlfac_coul[i]*fscal_elec[i]*sigma_pow[i]; - dvdl_vdw += vvdw[i]*DLF[i] + LFV[i]*alpha_vdw_eff*dlfac_vdw[i]*fscal_vdw[i]*sigma_pow[i]; - } - - dvdl[efptCOUL] += dvdl_coul; - dvdl[efptVDW] += dvdl_vdw; - - *velectot = velecsum; - *vvdwtot = vvdwsum; - - return fscal; -} diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.h b/src/gromacs/gmxlib/nonbonded/nb_free_energy.h index 8b76221810..c9366f377b 100644 --- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.h +++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.h @@ -42,22 +42,12 @@ #include "gromacs/legacyheaders/typedefs.h" void -gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, - rvec * gmx_restrict xx, - rvec * gmx_restrict ff, - t_forcerec * gmx_restrict fr, - const t_mdatoms * gmx_restrict mdatoms, - nb_kernel_data_t * gmx_restrict kernel_data, - t_nrnb * gmx_restrict nrnb); - -real - nb_free_energy_evaluate_single(real r2, real sc_r_power, real alpha_coul, - real alpha_vdw, real tabscale, real *vftab, - real qqA, real c6A, real c12A, real qqB, - real c6B, real c12B, real LFC[2], real LFV[2], real DLF[2], - real lfac_coul[2], real lfac_vdw[2], real dlfac_coul[2], - real dlfac_vdw[2], real sigma6_def, real sigma6_min, - real sigma2_def, real sigma2_min, - real *velectot, real *vvdwtot, real *dvdl); + gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, + rvec * gmx_restrict xx, + rvec * gmx_restrict ff, + t_forcerec * gmx_restrict fr, + const t_mdatoms * gmx_restrict mdatoms, + nb_kernel_data_t * gmx_restrict kernel_data, + t_nrnb * gmx_restrict nrnb); #endif diff --git a/src/gromacs/gmxlib/nonbonded/nonbonded.c b/src/gromacs/gmxlib/nonbonded/nonbonded.c index c21e3201b2..f66f73bb0f 100644 --- a/src/gromacs/gmxlib/nonbonded/nonbonded.c +++ b/src/gromacs/gmxlib/nonbonded/nonbonded.c @@ -443,266 +443,3 @@ void do_nonbonded(t_forcerec *fr, } } } - -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) - { - /* This check isn't race free. But it doesn't matter because if a race occurs the only - * disadvantage is that the warning is printed twice */ - 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; -} diff --git a/src/gromacs/gmxlib/nonbonded/pairs.c b/src/gromacs/gmxlib/nonbonded/pairs.c new file mode 100644 index 0000000000..d24d280cf3 --- /dev/null +++ b/src/gromacs/gmxlib/nonbonded/pairs.c @@ -0,0 +1,507 @@ +/* + * 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, 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 "pairs.h" + +#include + +#include "gromacs/legacyheaders/types/group.h" +#include "gromacs/listed-forces/bonded.h" +#include "gromacs/math/vec.h" +#include "gromacs/pbcutil/ishift.h" +#include "gromacs/pbcutil/mshift.h" +#include "gromacs/pbcutil/pbc.h" +#include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/fatalerror.h" + +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; +} + +static real +nb_free_energy_evaluate_single(real r2, real sc_r_power, real alpha_coul, real alpha_vdw, + real tabscale, real *vftab, + real qqA, real c6A, real c12A, real qqB, real c6B, real c12B, + real LFC[2], real LFV[2], real DLF[2], + real lfac_coul[2], real lfac_vdw[2], real dlfac_coul[2], real dlfac_vdw[2], + real sigma6_def, real sigma6_min, real sigma2_def, real sigma2_min, + real *velectot, real *vvdwtot, real *dvdl) +{ + real r, rp, rpm2, rtab, eps, eps2, Y, F, Geps, Heps2, Fp, VV, FF, fscal; + real qq[2], c6[2], c12[2], sigma6[2], sigma2[2], sigma_pow[2], sigma_powm2[2]; + real alpha_coul_eff, alpha_vdw_eff, dvdl_coul, dvdl_vdw; + real rpinv, r_coul, r_vdw, velecsum, vvdwsum; + real fscal_vdw[2], fscal_elec[2]; + real velec[2], vvdw[2]; + int i, ntab; + + const real half = 0.5; + const real one = 1.0; + const real two = 2.0; + const real six = 6.0; + const real fourtyeight = 48.0; + + qq[0] = qqA; + qq[1] = qqB; + c6[0] = c6A; + c6[1] = c6B; + c12[0] = c12A; + c12[1] = c12B; + + if (sc_r_power == six) + { + rpm2 = r2*r2; /* r4 */ + rp = rpm2*r2; /* r6 */ + } + else if (sc_r_power == fourtyeight) + { + rp = r2*r2*r2; /* r6 */ + rp = rp*rp; /* r12 */ + rp = rp*rp; /* r24 */ + rp = rp*rp; /* r48 */ + rpm2 = rp/r2; /* r46 */ + } + else + { + rp = pow(r2, half*sc_r_power); /* not currently supported as input, but can handle it */ + rpm2 = rp/r2; + } + + /* Loop over state A(0) and B(1) */ + for (i = 0; i < 2; i++) + { + if ((c6[i] > 0) && (c12[i] > 0)) + { + /* The c6 & c12 coefficients now contain the constants 6.0 and 12.0, respectively. + * Correct for this by multiplying with (1/12.0)/(1/6.0)=6.0/12.0=0.5. + */ + sigma6[i] = half*c12[i]/c6[i]; + sigma2[i] = pow(half*c12[i]/c6[i], 1.0/3.0); + /* should be able to get rid of this ^^^ internal pow call eventually. Will require agreement on + what data to store externally. Can't be fixed without larger scale changes, so not 5.0 */ + if (sigma6[i] < sigma6_min) /* for disappearing coul and vdw with soft core at the same time */ + { + sigma6[i] = sigma6_min; + sigma2[i] = sigma2_min; + } + } + else + { + sigma6[i] = sigma6_def; + sigma2[i] = sigma2_def; + } + if (sc_r_power == six) + { + sigma_pow[i] = sigma6[i]; + sigma_powm2[i] = sigma6[i]/sigma2[i]; + } + else if (sc_r_power == fourtyeight) + { + sigma_pow[i] = sigma6[i]*sigma6[i]; /* sigma^12 */ + sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^24 */ + sigma_pow[i] = sigma_pow[i]*sigma_pow[i]; /* sigma^48 */ + sigma_powm2[i] = sigma_pow[i]/sigma2[i]; + } + else + { /* not really supported as input, but in here for testing the general case*/ + sigma_pow[i] = pow(sigma2[i], sc_r_power/2); + sigma_powm2[i] = sigma_pow[i]/(sigma2[i]); + } + } + + /* only use softcore if one of the states has a zero endstate - softcore is for avoiding infinities!*/ + if ((c12[0] > 0) && (c12[1] > 0)) + { + alpha_vdw_eff = 0; + alpha_coul_eff = 0; + } + else + { + alpha_vdw_eff = alpha_vdw; + alpha_coul_eff = alpha_coul; + } + + /* Loop over A and B states again */ + for (i = 0; i < 2; i++) + { + fscal_elec[i] = 0; + fscal_vdw[i] = 0; + velec[i] = 0; + vvdw[i] = 0; + + /* Only spend time on A or B state if it is non-zero */ + if ( (qq[i] != 0) || (c6[i] != 0) || (c12[i] != 0) ) + { + /* Coulomb */ + rpinv = one/(alpha_coul_eff*lfac_coul[i]*sigma_pow[i]+rp); + r_coul = pow(rpinv, -one/sc_r_power); + + /* Electrostatics table lookup data */ + rtab = r_coul*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; + VV = Y+eps*Fp; + FF = Fp+Geps+two*Heps2; + velec[i] = qq[i]*VV; + fscal_elec[i] = -qq[i]*FF*r_coul*rpinv*tabscale; + + /* Vdw */ + rpinv = one/(alpha_vdw_eff*lfac_vdw[i]*sigma_pow[i]+rp); + r_vdw = pow(rpinv, -one/sc_r_power); + /* Vdw table lookup data */ + rtab = r_vdw*tabscale; + ntab = rtab; + eps = rtab-ntab; + eps2 = eps*eps; + ntab = 12*ntab; + /* Dispersion */ + Y = vftab[ntab+4]; + F = vftab[ntab+5]; + Geps = eps*vftab[ntab+6]; + Heps2 = eps2*vftab[ntab+7]; + Fp = F+Geps+Heps2; + VV = Y+eps*Fp; + FF = Fp+Geps+two*Heps2; + vvdw[i] = c6[i]*VV; + fscal_vdw[i] = -c6[i]*FF; + + /* Repulsion */ + Y = vftab[ntab+8]; + F = vftab[ntab+9]; + Geps = eps*vftab[ntab+10]; + Heps2 = eps2*vftab[ntab+11]; + Fp = F+Geps+Heps2; + VV = Y+eps*Fp; + FF = Fp+Geps+two*Heps2; + vvdw[i] += c12[i]*VV; + fscal_vdw[i] -= c12[i]*FF; + fscal_vdw[i] *= r_vdw*rpinv*tabscale; + } + } + /* Now we have velec[i], vvdw[i], and fscal[i] for both states */ + /* Assemble A and B states */ + velecsum = 0; + vvdwsum = 0; + dvdl_coul = 0; + dvdl_vdw = 0; + fscal = 0; + for (i = 0; i < 2; i++) + { + velecsum += LFC[i]*velec[i]; + vvdwsum += LFV[i]*vvdw[i]; + + fscal += (LFC[i]*fscal_elec[i]+LFV[i]*fscal_vdw[i])*rpm2; + + dvdl_coul += velec[i]*DLF[i] + LFC[i]*alpha_coul_eff*dlfac_coul[i]*fscal_elec[i]*sigma_pow[i]; + dvdl_vdw += vvdw[i]*DLF[i] + LFV[i]*alpha_vdw_eff*dlfac_vdw[i]*fscal_vdw[i]*sigma_pow[i]; + } + + dvdl[efptCOUL] += dvdl_coul; + dvdl[efptVDW] += dvdl_vdw; + + *velectot = velecsum; + *vvdwtot = vvdwsum; + + 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) + { + /* This check isn't race free. But it doesn't matter because if a race occurs the only + * disadvantage is that the warning is printed twice */ + 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; +} diff --git a/src/gromacs/gmxlib/nonbonded/pairs.h b/src/gromacs/gmxlib/nonbonded/pairs.h new file mode 100644 index 0000000000..4582b08aef --- /dev/null +++ b/src/gromacs/gmxlib/nonbonded/pairs.h @@ -0,0 +1,71 @@ +/* + * 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, 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 _pairs_h +#define _pairs_h + +#include "gromacs/legacyheaders/types/forcerec.h" +#include "gromacs/legacyheaders/types/mdatom.h" +#include "gromacs/math/vec.h" +#include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/real.h" + +#ifdef __cplusplus +extern "C" { +#endif +#if 0 +} /* fixes auto-indentation problems */ +#endif + +struct t_graph; +struct t_pbc; + +/* Calculate VdW/charge listed pair interactions (usually 1-4 interactions). + * global_atom_index is only passed for printing error messages. + */ +real +do_nonbonded_listed(int ftype, int nbonds, const t_iatom iatoms[], const t_iparams iparams[], + const rvec x[], rvec f[], rvec fshift[], + const struct t_pbc *pbc, const struct t_graph *g, + real *lambda, real *dvdl, const t_mdatoms *md, const t_forcerec *fr, + gmx_grppairener_t *grppener, int *global_atom_index); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/src/gromacs/legacyheaders/nonbonded.h b/src/gromacs/legacyheaders/nonbonded.h index c985396e59..3b7b75e949 100644 --- a/src/gromacs/legacyheaders/nonbonded.h +++ b/src/gromacs/legacyheaders/nonbonded.h @@ -38,9 +38,6 @@ #ifndef _nonbonded_h #define _nonbonded_h -#include "gromacs/legacyheaders/genborn.h" -#include "gromacs/legacyheaders/network.h" -#include "gromacs/legacyheaders/tgroup.h" #include "gromacs/legacyheaders/typedefs.h" #ifdef __cplusplus @@ -50,10 +47,6 @@ extern "C" { } /* fixes auto-indentation problems */ #endif -struct t_graph; -struct t_pbc; - - void gmx_nonbonded_setup(t_forcerec * fr, gmx_bool bGenericKernelOnly); @@ -83,16 +76,6 @@ do_nonbonded(t_forcerec *fr, t_nrnb *nrnb, real *lambda, real dvdlambda[], int nls, int eNL, int flags); -/* Calculate VdW/charge listed pair interactions (usually 1-4 interactions). - * global_atom_index is only passed for printing error messages. - */ -real -do_nonbonded_listed(int ftype, int nbonds, const t_iatom iatoms[], const t_iparams iparams[], - const rvec x[], rvec f[], rvec fshift[], - const struct t_pbc *pbc, const struct t_graph *g, - real *lambda, real *dvdl, const t_mdatoms *md, const t_forcerec *fr, - gmx_grppairener_t *grppener, int *global_atom_index); - #ifdef __cplusplus } #endif diff --git a/src/gromacs/legacyheaders/types/group.h b/src/gromacs/legacyheaders/types/group.h index 035b12b215..ea2fd293c7 100644 --- a/src/gromacs/legacyheaders/types/group.h +++ b/src/gromacs/legacyheaders/types/group.h @@ -34,7 +34,8 @@ * 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_LEGACYHEADERS_TYPES_GROUP_H +#define GMX_LEGACYHEADERS_TYPES_GROUP_H #include "gromacs/legacyheaders/types/simple.h" @@ -90,3 +91,5 @@ typedef struct { #ifdef __cplusplus } #endif + +#endif diff --git a/src/gromacs/listed-forces/bonded.cpp b/src/gromacs/listed-forces/bonded.cpp index dfb846d9e7..6afb559bba 100644 --- a/src/gromacs/listed-forces/bonded.cpp +++ b/src/gromacs/listed-forces/bonded.cpp @@ -55,11 +55,11 @@ #include +#include "gromacs/gmxlib/nonbonded/pairs.h" #include "gromacs/legacyheaders/disre.h" #include "gromacs/legacyheaders/force.h" #include "gromacs/legacyheaders/macros.h" #include "gromacs/legacyheaders/names.h" -#include "gromacs/legacyheaders/nonbonded.h" #include "gromacs/legacyheaders/ns.h" #include "gromacs/legacyheaders/orires.h" #include "gromacs/legacyheaders/txtdump.h" -- 2.22.0