Version 4.6 introduced an error in the reaction-field correction
force term for perturbed interactions. A factor r_softcore^2 was
missing. The force calculation code is now slightly reorganized
and comments added to avoid such issues in the future.
Fixes #1318
Change-Id: I9105139f8975495c323008ce202cde517a69281a
case GMX_NBKERNEL_ELEC_EWALD:
/* simple cutoff (yes, ewald is done all on direct space for free energy) */
Vcoul[i] = qq[i]*rinvC;
case GMX_NBKERNEL_ELEC_EWALD:
/* simple cutoff (yes, ewald is done all on direct space for free energy) */
Vcoul[i] = qq[i]*rinvC;
- FscalC[i] = Vcoul[i]*rpinvC;
break;
case GMX_NBKERNEL_ELEC_REACTIONFIELD:
/* reaction-field */
break;
case GMX_NBKERNEL_ELEC_REACTIONFIELD:
/* reaction-field */
- Vcoul[i] = qq[i]*(rinvC+krf*rC*rC-crf);
- FscalC[i] = qq[i]*(rinvC*rpinvC-2.0*krf);
+ Vcoul[i] = qq[i]*(rinvC + krf*rC*rC-crf);
+ FscalC[i] = qq[i]*(rinvC - 2.0*krf*rC*rC);
break;
case GMX_NBKERNEL_ELEC_CUBICSPLINETABLE:
break;
case GMX_NBKERNEL_ELEC_CUBICSPLINETABLE:
VV = Y+epsC*Fp;
FF = Fp+Geps+2.0*Heps2;
Vcoul[i] = qq[i]*VV;
VV = Y+epsC*Fp;
FF = Fp+Geps+2.0*Heps2;
Vcoul[i] = qq[i]*VV;
- FscalC[i] = -qq[i]*tabscale*FF*rC*rpinvC;
+ FscalC[i] = -qq[i]*tabscale*FF*rC;
break;
case GMX_NBKERNEL_ELEC_GENERALIZEDBORN:
break;
case GMX_NBKERNEL_ELEC_GENERALIZEDBORN:
- Vvdw[i] = Vvdw12*(1.0/12.0)-Vvdw6*(1.0/6.0);
+ Vvdw[i] = Vvdw12*(1.0/12.0) - Vvdw6*(1.0/6.0);
- FscalV[i] = (Vvdw12-Vvdw6)*rpinvV;
+ FscalV[i] = Vvdw12 - Vvdw6;
break;
case GMX_NBKERNEL_VDW_BUCKINGHAM:
break;
case GMX_NBKERNEL_VDW_BUCKINGHAM:
VV = Y+epsV*Fp;
FF = Fp+Geps+2.0*Heps2;
Vvdw[i] += c6[i]*VV;
VV = Y+epsV*Fp;
FF = Fp+Geps+2.0*Heps2;
Vvdw[i] += c6[i]*VV;
- FscalV[i] -= c6[i]*tabscale*FF*rV*rpinvV;
+ FscalV[i] -= c6[i]*tabscale*FF*rV;
/* repulsion */
Y = VFtab[nnn+4];
/* repulsion */
Y = VFtab[nnn+4];
VV = Y+epsV*Fp;
FF = Fp+Geps+2.0*Heps2;
Vvdw[i] += c12[i]*VV;
VV = Y+epsV*Fp;
FF = Fp+Geps+2.0*Heps2;
Vvdw[i] += c12[i]*VV;
- FscalV[i] -= c12[i]*tabscale*FF*rV*rpinvV;
+ FscalV[i] -= c12[i]*tabscale*FF*rV;
break;
case GMX_NBKERNEL_VDW_NONE:
break;
case GMX_NBKERNEL_VDW_NONE:
Vvdw[i] = (rV < rvdw) ? Vvdw[i] : 0.0;
}
}
Vvdw[i] = (rV < rvdw) ? Vvdw[i] : 0.0;
}
}
+
+ /* FscalC (and FscalV) now contain: dV/drC * rC
+ * Now we multiply by rC^-p, so it will be: dV/drC * rC^1-p
+ * Further down we first multiply by r^p-2 and then by
+ * the vector r, which in total gives: dV/drC * (r/rC)^1-p
+ */
+ FscalC[i] *= rpinvC;
+ FscalV[i] *= rpinvV;