fix subtle bug in stress tally for dihedral style table

This commit is contained in:
Axel Kohlmeyer
2022-12-26 16:47:29 -05:00
parent e82eb9ecc4
commit 0df140876c
2 changed files with 9 additions and 10 deletions

View File

@ -13,7 +13,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Andrew Jewett (jewett.aij g m ail)
Contributing author: Andrew Jewett (jewett.aij at gmail)
The cyclic tridiagonal matrix solver was borrowed from
the "tridiag.c" written by Gerard Jungman for GSL
------------------------------------------------------------------------- */
@ -375,7 +375,7 @@ static double cyc_splintD(double const *xa, double const *ya, double const *y2a,
// --------------------------------------------
// ------- Calculate the dihedral angle -------
// --------------------------------------------
static const int g_dim=3;
static constexpr int g_dim=3;
static double Phi(double const *x1, //array holding x,y,z coords atom 1
double const *x2, // : : : : 2
@ -692,10 +692,9 @@ void DihedralTable::compute(int eflag, int vflag)
}
if (evflag)
ev_tally(i1,i2,i3,i4,
nlocal,newton_bond,edihedral,
ev_tally(i1,i2,i3,i4,nlocal,newton_bond,edihedral,
f1,f3,f4,
vb12[0],vb12[1],vb12[2],
-vb12[0],-vb12[1],-vb12[2],
vb23[0],vb23[1],vb23[2],
vb34[0],vb34[1],vb34[2]);
}

View File

@ -16,20 +16,20 @@
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "omp_compat.h"
#include "dihedral_table_omp.h"
#include <cmath>
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "neighbor.h"
#include "force.h"
#include "math_const.h"
#include "math_extra.h"
#include "suffix.h"
#include <cmath>
#include "omp_compat.h"
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace MathExtra;
@ -380,7 +380,7 @@ void DihedralTableOMP::eval(int nfrom, int nto, ThrData * const thr)
if (EVFLAG)
ev_tally_thr(this,i1,i2,i3,i4,nlocal,NEWTON_BOND,edihedral,f1,f3,f4,
vb12[0],vb12[1],vb12[2],vb23[0],vb23[1],vb23[2],vb34[0],
-vb12[0],-vb12[1],-vb12[2],vb23[0],vb23[1],vb23[2],vb34[0],
vb34[1],vb34[2],thr);
}
}