use morse interaction energy from bond style morse not pair style morse

This commit is contained in:
Axel Kohlmeyer
2022-10-04 07:02:44 -04:00
parent 4a06559da5
commit f2b3d8e8c3
2 changed files with 7 additions and 6 deletions

View File

@ -88,7 +88,7 @@ The style *morse* computes pairwise interactions with the formula
.. math::
E = D_0 \left[ e^{- 2 \alpha (r - r_0)} - 2 e^{- \alpha (r - r_0)} \right] \qquad r < r_c
E = D_0 \left[ 1 - e^{-\alpha (r - r_0)} \right]^2 \qquad r < r_c
and the coefficients:
@ -171,8 +171,9 @@ Related commands
:doc:`pair_coeff <pair_coeff>`,
:doc:`pair_style hybrid/overlay <pair_hybrid>`,
:doc:`pair_style lj/cut <pair_lj>`,
:doc:`pair_style morse <pair_morse>`,
:doc:`bond_style morse <bond_morse>`,
:doc:`bond_style harmonic <bond_harmonic>`
:doc:`bond_style quartic <bond_quartic>`
Default
"""""""

View File

@ -168,11 +168,11 @@ void PairList::compute(int eflag, int vflag)
} else if (par.style == MORSE) {
const double r = sqrt(rsq);
const double dr = par.param.morse.r0 - r;
const double dexp = exp(par.param.morse.alpha * dr);
const double dr = r - par.param.morse.r0;
const double dexp = exp(-par.param.morse.alpha * dr);
fpair = 2.0 * par.param.morse.d0 * par.param.morse.alpha * (dexp * dexp - dexp) / r;
if (eflag_either) epair = par.param.morse.d0 * (dexp * dexp - 2.0 * dexp) - par.offset;
if (eflag_either) epair = par.param.morse.d0 * (dexp * dexp - 2.0 * dexp + 1.0) - par.offset;
} else if (par.style == LJ126) {
@ -385,7 +385,7 @@ void PairList::init_style()
} else if (par.style == MORSE) {
const double dr = par.param.morse.r0 - sqrt(par.cutsq);
const double dexp = exp(par.param.morse.alpha * dr);
par.offset = par.param.morse.d0 * (dexp * dexp - 2.0 * dexp);
par.offset = par.param.morse.d0 * (dexp * dexp - 2.0 * dexp - 1.0);
} else if (par.style == LJ126) {
const double r6inv = par.cutsq * par.cutsq * par.cutsq;