ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

correctly handle the case of n=1. clean up some ugliness

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-02-24 18:03:47 -05:00
parent d81ca27e96
commit 7efacdc911
2 changed files with 13 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

15
src/USER-MISC/dihedral_nharmonic.cpp
View File

@ -181,15 +181,14 @@ void DihedralNHarmonic::compute(int eflag, int vflag)
// p = sum (i=1,n) a_i * c**(i-1) // p = sum (i=1,n) a_i * c**(i-1)
// pd = dp/dc // pd = dp/dc
c_ = c; c_ = 1.0;
p = this->a[type][0]; p = a[type][0];
pd = this->a[type][1]; pd = 0.0;
for (int i = 1; i < nterms[type]-1; i++) { for (int i = 1; i < nterms[type]; i++) {
p += c_ * this->a[type][i]; pd += c_ * i * a[type][i];
pd += c_ * static_cast<double>(i+1) * this->a[type][i+1];
c_ *= c; c_ *= c;
p += c_ * a[type][i];
} }
p += c_ * this->a[type][nterms[type]-1];
if (eflag) edihedral = p; if (eflag) edihedral = p;
@ -275,7 +274,7 @@ void DihedralNHarmonic::allocate()
void DihedralNHarmonic::coeff(int narg, char **arg) void DihedralNHarmonic::coeff(int narg, char **arg)
{ {
if (narg < 4 ) error->all(FLERR,"Incorrect args for dihedral coefficients"); if (narg < 3) error->all(FLERR,"Incorrect args for dihedral coefficients");
int n = utils::inumeric(FLERR,arg[1],false,lmp); int n = utils::inumeric(FLERR,arg[1],false,lmp);
if (narg != n + 2) if (narg != n + 2)

12
src/USER-OMP/dihedral_nharmonic_omp.cpp
View File

@ -203,15 +203,15 @@ void DihedralNHarmonicOMP::eval(int nfrom, int nto, ThrData * const thr)
// force & energy // force & energy
// p = sum (i=1,n) a_i * c**(i-1) // p = sum (i=1,n) a_i * c**(i-1)
// pd = dp/dc // pd = dp/dc
c_ = c;
c_ = 1.0;
p = a[type][0]; p = a[type][0];
pd = a[type][1]; pd = 0.0;
for (int i = 1; i < nterms[type]-1; i++) { for (int i = 1; i < nterms[type]; i++) {
p += c_ * a[type][i]; pd += c_ * i * a[type][i];
pd += c_ * static_cast<double>(i+1) * a[type][i+1];
c_ *= c; c_ *= c;
p += c_ * a[type][i];
} }
p += c_ * a[type][nterms[type]-1];
if (EFLAG) edihedral = p; if (EFLAG) edihedral = p;