* 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 <stdio.h>
#include <math.h>
+#include <stdio.h>
+
+#include "gromacs/legacyheaders/coulomb.h"
+#include "gromacs/legacyheaders/macros.h"
+#include "gromacs/legacyheaders/names.h"
+#include "gromacs/legacyheaders/txtdump.h"
+#include "gromacs/legacyheaders/typedefs.h"
+#include "gromacs/legacyheaders/types/commrec.h"
+#include "gromacs/math/units.h"
#include "gromacs/math/utilities.h"
-#include "typedefs.h"
-#include "vec.h"
-#include "coulomb.h"
-#include "smalloc.h"
-#include "physics.h"
-#include "txtdump.h"
-#include "gromacs/fileio/futil.h"
-#include "names.h"
-#include "writeps.h"
-#include "macros.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/smalloc.h"
real calc_ewaldcoeff_q(real rc, real dtol)
{
double Vexcl_q, dvdl_excl_q, dvdl_excl_lj; /* Necessary for precision */
double Vexcl_lj;
real one_4pi_eps;
- real v, vc, qiA, qiB, dr, dr2, rinv, fscal, enercorr;
+ real v, vc, qiA, qiB, dr2, rinv, enercorr;
real Vself_q[2], Vself_lj[2], Vdipole[2], rinv2, ewc_q = fr->ewaldcoeff_q, ewcdr;
- real ewc_lj = fr->ewaldcoeff_lj;
- real c6Ai = 0, c6Bi = 0, c6A = 0, c6B = 0, ewcdr2, ewcdr4, ewcdr5, c6L = 0, rinv6;
+ real ewc_lj = fr->ewaldcoeff_lj, ewc_lj2 = ewc_lj * ewc_lj;
+ real c6Ai = 0, c6Bi = 0, c6A = 0, c6B = 0, ewcdr2, ewcdr4, c6L = 0, rinv6;
rvec df, dx, mutot[2], dipcorrA, dipcorrB;
tensor dxdf_q, dxdf_lj;
real vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
vr0_q = ewc_q*M_2_SQRTPI;
if (EVDW_PME(fr->vdwtype))
{
- vr0_lj = -pow(ewc_lj, 6)/6.0;
+ vr0_lj = -pow(ewc_lj, 6)/6.0;
}
AA = excl->a;
{
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++)
{
}
if (qqA != 0.0 || c6A != 0.0)
{
- rvec_sub(x[i], x[k], dx);
+ real fscal;
+
fscal = 0;
+ rvec_sub(x[i], x[k], dx);
if (bMolPBC)
{
/* Cheap pbc_dx, assume excluded pairs are at short distance. */
{
rinv = gmx_invsqrt(dr2);
rinv2 = rinv*rinv;
- dr = 1.0/rinv;
if (qqA != 0.0)
{
+ real dr;
+
+ dr = 1.0/rinv;
ewcdr = ewc_q*dr;
vc = qqA*gmx_erf(ewcdr)*rinv;
Vexcl_q += vc;
/* Relative accuracy at R_ERF_R_INACC of 2e-5 */
#define R_ERF_R_INACC 0.1
#endif
+ /* fscal is the scalar force pre-multiplied by rinv,
+ * to normalise the relative position vector dx */
if (ewcdr > R_ERF_R_INACC)
{
fscal = rinv2*(vc - qqA*ewc_q*M_2_SQRTPI*exp(-ewcdr*ewcdr));
}
/* The force vector is obtained by multiplication with
- * the distance vector
+ * the relative position vector
*/
svmul(fscal, dx, df);
rvec_inc(f[k], df);
if (c6A != 0.0)
{
rinv6 = rinv2*rinv2*rinv2;
- ewcdr = ewc_lj*dr;
- ewcdr2 = ewcdr*ewcdr;
+ ewcdr2 = ewc_lj2*dr2;
ewcdr4 = ewcdr2*ewcdr2;
- ewcdr5 = ewcdr4*ewcdr;
/* We get the excluded long-range contribution from -C6*(1-g(r))
* g(r) is also defined in the manual under LJ-PME
*/
vc = -c6A*rinv6*(1.0 - exp(-ewcdr2)*(1 + ewcdr2 + 0.5*ewcdr4));
Vexcl_lj += vc;
- /* The force is the derivative of the potential vc */
- fscal = 6.0*vc*rinv + c6A*rinv6*exp(-ewcdr2)*ewc_lj*ewcdr5;
+ /* The force is the derivative of the potential vc.
+ * fscal is the scalar force pre-multiplied by rinv,
+ * to normalise the relative position vector dx */
+ fscal = 6.0*vc*rinv2 + c6A*rinv6*exp(-ewcdr2)*ewc_lj2*ewcdr4;
- /* The force vector is obtained by
- * multiplication with the distance
- * vector
+ /* The force vector is obtained by multiplication with
+ * the relative position vector
*/
svmul(fscal, dx, df);
rvec_inc(f[k], df);
}
}
}
- 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++)
{
}
if (qqA != 0.0 || qqB != 0.0 || c6A != 0.0 || c6B != 0.0)
{
+ real fscal;
+
fscal = 0;
qqL = L1_q*qqA + lambda_q*qqB;
if (EVDW_PME(fr->vdwtype))
{
rinv = gmx_invsqrt(dr2);
rinv2 = rinv*rinv;
- dr = 1.0/rinv;
if (qqA != 0.0 || qqB != 0.0)
{
+ real dr;
+
+ dr = 1.0/rinv;
v = gmx_erf(ewc_q*dr)*rinv;
vc = qqL*v;
Vexcl_q += vc;
+ /* fscal is the scalar force pre-multiplied by rinv,
+ * to normalise the relative position vector dx */
fscal = rinv2*(vc-qqL*ewc_q*M_2_SQRTPI*exp(-ewc_q*ewc_q*dr2));
dvdl_excl_q += (qqB - qqA)*v;
+ /* The force vector is obtained by multiplication with
+ * the relative position vector
+ */
svmul(fscal, dx, df);
rvec_inc(f[k], df);
rvec_dec(f[i], df);
if ((c6A != 0.0 || c6B != 0.0) && EVDW_PME(fr->vdwtype))
{
rinv6 = rinv2*rinv2*rinv2;
- ewcdr = ewc_lj*dr;
- ewcdr2 = ewcdr*ewcdr;
+ ewcdr2 = ewc_lj2*dr2;
ewcdr4 = ewcdr2*ewcdr2;
- ewcdr5 = ewcdr4*ewcdr;
v = -rinv6*(1.0 - exp(-ewcdr2)*(1 + ewcdr2 + 0.5*ewcdr4));
vc = c6L*v;
Vexcl_lj += vc;
- fscal = 6.0*vc*rinv + c6L*rinv6*exp(-ewcdr2)*ewc_lj*ewcdr5;
+ /* fscal is the scalar force pre-multiplied by rinv,
+ * to normalise the relative position vector dx */
+ fscal = 6.0*vc*rinv2 + c6L*rinv6*exp(-ewcdr2)*ewc_lj2*ewcdr4;
dvdl_excl_lj += (c6B - c6A)*v;
+
+ /* The force vector is obtained by multiplication with
+ * the relative position vector
+ */
svmul(fscal, dx, df);
rvec_inc(f[k], df);
rvec_dec(f[i], df);
Vself_q[q] = ewc_q*one_4pi_eps*fr->q2sum[q]*M_1_SQRTPI;
if (EVDW_PME(fr->vdwtype))
{
- Vself_lj[q] = -pow(ewc_lj, 6)*fr->c6sum[q];
+ Vself_lj[q] = fr->c6sum[q]*0.5*vr0_lj;
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff