From 41f518ede7be331f766895e91fee6c09ca9d98e5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 9 Apr 2024 00:57:15 -0400 Subject: [PATCH] repulsive is r^12 not r^6 --- src/EXTRA-PAIR/pair_pedone.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/EXTRA-PAIR/pair_pedone.cpp b/src/EXTRA-PAIR/pair_pedone.cpp index 4e634b3dfc..59b883b2c5 100644 --- a/src/EXTRA-PAIR/pair_pedone.cpp +++ b/src/EXTRA-PAIR/pair_pedone.cpp @@ -107,7 +107,7 @@ void PairPedone::compute(int eflag, int vflag) r = sqrt(rsq); dr = r - r0[itype][jtype]; dexp = exp(-alpha[itype][jtype] * dr); - fpair = pedone1[itype][jtype] * (dexp * dexp - dexp) / r - pedone2[itype][jtype] * r6inv * r2inv; + fpair = pedone1[itype][jtype] * (dexp * dexp - dexp) / r - pedone2[itype][jtype] * r6inv * r6inv * r2inv; fpair *= factor_lj; f[i][0] += delx * fpair; @@ -120,7 +120,7 @@ void PairPedone::compute(int eflag, int vflag) } if (eflag) { - evdwl = d0[itype][jtype] * (dexp * dexp - 2.0 * dexp) - c0[itype][jtype] * r6inv - + evdwl = d0[itype][jtype] * (dexp * dexp - 2.0 * dexp) - c0[itype][jtype] * r6inv * r6inv - offset[itype][jtype]; evdwl *= factor_lj; } @@ -224,12 +224,12 @@ double PairPedone::init_one(int i, int j) if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set"); pedone1[i][j] = 2.0 * d0[i][j] * alpha[i][j]; - pedone2[i][j] = 6.0 * c0[i][j]; + pedone2[i][j] = 12.0 * c0[i][j]; if (offset_flag) { double alpha_dr = -alpha[i][j] * (cut[i][j] - r0[i][j]); offset[i][j] = - d0[i][j] * (exp(2.0 * alpha_dr) - 2.0 * exp(alpha_dr)) - c0[i][j] / pow(cut[i][j], 6.0); + d0[i][j] * (exp(2.0 * alpha_dr) - 2.0 * exp(alpha_dr)) - c0[i][j] / pow(cut[i][j], 12.0); } else offset[i][j] = 0.0; @@ -362,10 +362,10 @@ double PairPedone::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq dexp = exp(-alpha[itype][jtype] * dr); r2inv = 1.0 / rsq; r6inv = r2inv * r2inv * r2inv; - fforce = pedone1[itype][jtype] * (dexp * dexp - dexp) / r - pedone2[itype][jtype] * r6inv * r2inv; + fforce = pedone1[itype][jtype] * (dexp * dexp - dexp) / r - pedone2[itype][jtype] * r6inv * r6inv * r2inv; fforce *= factor_lj; - phi = d0[itype][jtype] * (dexp * dexp - 2.0 * dexp) - c0[itype][jtype] * r6inv - + phi = d0[itype][jtype] * (dexp * dexp - 2.0 * dexp) - c0[itype][jtype] * r6inv * r6inv - offset[itype][jtype]; return factor_lj * phi; }