From: Erik Lindahl Date: Fri, 27 Jun 2014 09:04:29 +0000 (+0200) Subject: Fixed free energy with LJ-PME X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=fc819eee87572d5608778642bcf22badaf3b4361;p=alexxy%2Fgromacs.git Fixed free energy with LJ-PME First, Only the normal LJ part of the potential was shifted for LJ-PME with potential-shift modifier in the free energy kernel, and not the Ewald component. This had an effect of roughly 0.05% on the total energy. Second, Berk found that the verlet nbnxn SIMD LJ-PME kernel had an issue that caused the potential and force to be slightly off. With this patch, both potential and force match the non-free-energy kernel perfectly with at lambda=0. Change-Id: I4bc53ce1bdb96ba06656b86470dbf8dbd1e81972 --- diff --git a/src/gromacs/gmxlib/ewald_util.c b/src/gromacs/gmxlib/ewald_util.c index b6af0d607b..d3e5f23748 100644 --- a/src/gromacs/gmxlib/ewald_util.c +++ b/src/gromacs/gmxlib/ewald_util.c @@ -226,7 +226,7 @@ void ewald_LRcorrection(int start, int end, { clear_mat(dxdf_lj); } - if ((calc_excl_corr || dipole_coeff != 0 || EVDW_PME(fr->vdwtype)) && !bFreeEnergy) + if ((calc_excl_corr || dipole_coeff != 0) && !bFreeEnergy) { for (i = start; (i < end); i++) { @@ -387,7 +387,7 @@ void ewald_LRcorrection(int start, int end, } } } - else if (calc_excl_corr || dipole_coeff != 0 || EVDW_PME(fr->vdwtype)) + else if (calc_excl_corr || dipole_coeff != 0) { for (i = start; (i < end); i++) { diff --git a/src/gromacs/gmxlib/nonbonded/nb_free_energy.c b/src/gromacs/gmxlib/nonbonded/nb_free_energy.c index 058b2a87eb..2bf8f90e3f 100644 --- a/src/gromacs/gmxlib/nonbonded/nb_free_energy.c +++ b/src/gromacs/gmxlib/nonbonded/nb_free_energy.c @@ -106,6 +106,8 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, const real * chargeB; real sigma6_min, sigma6_def, lam_power, sc_power, sc_r_power; real alpha_coul, alpha_vdw, lambda_coul, lambda_vdw, ewc_lj; + real ewcljrsq, ewclj, ewclj2, exponent, poly, vvdw_disp, vvdw_rep, sh_lj_ewald; + real ewclj6; const real * nbfp, *nbfp_grid; real * dvdl; real * Vv; @@ -178,6 +180,10 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, sh_ewald = fr->ic->sh_ewald; rvdw = fr->rvdw; sh_invrc6 = fr->ic->sh_invrc6; + sh_lj_ewald = fr->ic->sh_lj_ewald; + ewclj = fr->ewaldcoeff_lj; + ewclj2 = ewclj*ewclj; + ewclj6 = ewclj2*ewclj2*ewclj2; if (fr->coulomb_modifier == eintmodPOTSWITCH) { @@ -273,6 +279,14 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, */ bConvertLJEwaldToLJ6 = (bEwaldLJ && (fr->vdw_modifier != eintmodPOTSWITCH)); + /* We currently don't implement exclusion correction, needed with the Verlet cut-off scheme, without conversion */ + if (fr->cutoff_scheme == ecutsVERLET && + ((bEwald && !bConvertEwaldToCoulomb) || + (bEwaldLJ && !bConvertLJEwaldToLJ6))) + { + gmx_incons("Unimplemented non-bonded setup"); + } + /* fix compiler warnings */ nj1 = 0; n1C = n1V = 0; @@ -514,7 +528,7 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, Vcoul[i] = qq[i]*rinvC; FscalC[i] = Vcoul[i]; /* The shift for the Coulomb potential is stored in - * the RF parameter c_rf, which is 0 without shift + * the RF parameter c_rf, which is 0 without shift. */ Vcoul[i] -= qq[i]*fr->ic->c_rf; break; @@ -602,7 +616,6 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, switch (ivdw) { case GMX_NBKERNEL_VDW_LENNARDJONES: - case GMX_NBKERNEL_VDW_LJEWALD: /* cutoff LJ */ if (sc_r_power == 6.0) { @@ -610,19 +623,14 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, } else { - rinv6 = pow(rinvV, 6.0); + rinv6 = rinvV*rinvV; + rinv6 = rinv6*rinv6*rinv6; } Vvdw6 = c6[i]*rinv6; Vvdw12 = c12[i]*rinv6*rinv6; - if (fr->vdw_modifier == eintmodPOTSHIFT) - { - Vvdw[i] = ( (Vvdw12-c12[i]*sh_invrc6*sh_invrc6)*(1.0/12.0) - -(Vvdw6-c6[i]*sh_invrc6)*(1.0/6.0)); - } - else - { - Vvdw[i] = Vvdw12*(1.0/12.0) - Vvdw6*(1.0/6.0); - } + + Vvdw[i] = ( (Vvdw12 - c12[i]*sh_invrc6*sh_invrc6)*(1.0/12.0) + - (Vvdw6 - c6[i]*sh_invrc6)*(1.0/6.0)); FscalV[i] = Vvdw12 - Vvdw6; break; @@ -656,6 +664,41 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, FscalV[i] -= c12[i]*tabscale*FF*rV; break; + case GMX_NBKERNEL_VDW_LJEWALD: + if (sc_r_power == 6.0) + { + rinv6 = rpinvV; + } + else + { + rinv6 = rinvV*rinvV; + rinv6 = rinv6*rinv6*rinv6; + } + c6grid = nbfp_grid[tj[i]]; + + if (bConvertLJEwaldToLJ6) + { + /* cutoff LJ */ + Vvdw6 = c6[i]*rinv6; + Vvdw12 = c12[i]*rinv6*rinv6; + + Vvdw[i] = ( (Vvdw12 - c12[i]*sh_invrc6*sh_invrc6)*(1.0/12.0) + - (Vvdw6 - c6[i]*sh_invrc6 - c6grid*sh_lj_ewald)*(1.0/6.0)); + FscalV[i] = Vvdw12 - Vvdw6; + } + else + { + /* Normal LJ-PME */ + ewcljrsq = ewclj2*rV*rV; + exponent = exp(-ewcljrsq); + poly = exponent*(1.0 + ewcljrsq + ewcljrsq*ewcljrsq*0.5); + vvdw_disp = (c6[i]-c6grid*(1.0-poly))*rinv6; + vvdw_rep = c12[i]*rinv6*rinv6; + FscalV[i] = vvdw_rep - vvdw_disp - c6grid*(1.0/6.0)*exponent*ewclj6; + Vvdw[i] = (vvdw_rep - c12[i]*sh_invrc6*sh_invrc6)/12.0 - (vvdw_disp - c6[i]*sh_invrc6 - c6grid*sh_lj_ewald)/6.0; + } + break; + case GMX_NBKERNEL_VDW_NONE: Vvdw[i] = 0.0; FscalV[i] = 0.0; @@ -748,6 +791,10 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, v_lr = (ewtab[ewitab+2]-ewtabhalfspace*eweps*(ewtab[ewitab]+f_lr)); f_lr *= rinv; + /* Note that any possible Ewald shift has already been applied in + * the normal interaction part above. + */ + if (ii == jnr) { /* If we get here, the i particle (ii) has itself (jnr) @@ -776,6 +823,10 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, * the softcore to the entire VdW interaction, * including the reciprocal-space component. */ + /* We could also use the analytical form here + * iso a table, but that can cause issues for + * r close to 0 for non-interacting pairs. + */ real rs, frac, f_lr; int ri; @@ -783,8 +834,11 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, ri = (int)rs; frac = rs - ri; f_lr = (1 - frac)*tab_ewald_F_lj[ri] + frac*tab_ewald_F_lj[ri+1]; - FF = f_lr*rinv; - VV = tab_ewald_V_lj[ri] - ewtabhalfspace*frac*(tab_ewald_F_lj[ri] + f_lr); + /* TODO: Currently the Ewald LJ table does not contain + * the factor 1/6, we should add this. + */ + FF = f_lr*rinv/6.0; + VV = (tab_ewald_V_lj[ri] - ewtabhalfspace*frac*(tab_ewald_F_lj[ri] + f_lr))/6.0; if (ii == jnr) { @@ -799,11 +853,10 @@ gmx_nb_free_energy_kernel(const t_nblist * gmx_restrict nlist, for (i = 0; i < NSTATES; i++) { c6grid = nbfp_grid[tj[i]]; - vvtot += LFV[i]*c6grid*VV*(1.0/6.0); - Fscal += LFV[i]*c6grid*FF*(1.0/6.0); - dvdl_vdw += (DLF[i]*c6grid)*VV*(1.0/6.0); + vvtot += LFV[i]*c6grid*VV; + Fscal += LFV[i]*c6grid*FF; + dvdl_vdw += (DLF[i]*c6grid)*VV; } - } if (bDoForces) diff --git a/src/gromacs/mdlib/nbnxn_atomdata.c b/src/gromacs/mdlib/nbnxn_atomdata.c index 6836541212..d271a8ce7a 100644 --- a/src/gromacs/mdlib/nbnxn_atomdata.c +++ b/src/gromacs/mdlib/nbnxn_atomdata.c @@ -851,7 +851,7 @@ static void copy_lj_to_nbat_lj_comb_x8(const real *ljparam_type, } } -/* Sets the atom type and LJ data in nbnxn_atomdata_t */ +/* Sets the atom type in nbnxn_atomdata_t */ static void nbnxn_atomdata_set_atomtypes(nbnxn_atomdata_t *nbat, int ngrid, const nbnxn_search_t nbs, @@ -872,9 +872,30 @@ static void nbnxn_atomdata_set_atomtypes(nbnxn_atomdata_t *nbat, copy_int_to_nbat_int(nbs->a+ash, grid->cxy_na[i], ncz*grid->na_sc, type, nbat->ntype-1, nbat->type+ash); + } + } +} + +/* Sets the LJ combination rule parameters in nbnxn_atomdata_t */ +static void nbnxn_atomdata_set_ljcombparams(nbnxn_atomdata_t *nbat, + int ngrid, + const nbnxn_search_t nbs) +{ + int g, i, ncz, ash; + const nbnxn_grid_t *grid; - if (nbat->comb_rule != ljcrNONE) + if (nbat->comb_rule != ljcrNONE) + { + for (g = 0; g < ngrid; g++) + { + grid = &nbs->grid[g]; + + /* Loop over all columns and copy and fill */ + for (i = 0; i < grid->ncx*grid->ncy; i++) { + ncz = grid->cxy_ind[i+1] - grid->cxy_ind[i]; + ash = (grid->cell0 + grid->cxy_ind[i])*grid->na_sc; + if (nbat->XFormat == nbatX4) { copy_lj_to_nbat_lj_comb_x4(nbat->nbfp_comb, @@ -1096,6 +1117,9 @@ void nbnxn_atomdata_set(nbnxn_atomdata_t *nbat, nbnxn_atomdata_mask_fep(nbat, ngrid, nbs); } + /* This must be done after masking types for FEP */ + nbnxn_atomdata_set_ljcombparams(nbat, ngrid, nbs); + nbnxn_atomdata_set_energygroups(nbat, ngrid, nbs, atinfo); } diff --git a/src/gromacs/mdlib/nbnxn_search.c b/src/gromacs/mdlib/nbnxn_search.c index c06a5d10f1..ae8b82caf7 100644 --- a/src/gromacs/mdlib/nbnxn_search.c +++ b/src/gromacs/mdlib/nbnxn_search.c @@ -3440,7 +3440,6 @@ static void make_fep_list(const nbnxn_search_t nbs, * Note that the charge has been set to zero, * but we need to avoid 0/0, as perturbed atoms * can be on top of each other. - * (and the LJ parameters have not been zeroed) */ nbl->cj[cj_ind].excl &= ~(1U << (i*nbl->na_cj + j)); }