bugfix for fix qeq/reax to make it usable without pair reax/c

This commit is contained in:
Axel Kohlmeyer
2019-02-05 11:55:02 +01:00
parent b417cfda9b
commit 0c4e76ce84

View File

@ -126,14 +126,12 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) :
reaxc = NULL;
reaxc = (PairReaxC *) force->pair_match("reax/c",0);
if (reaxc) {
s_hist = t_hist = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
for (int i = 0; i < atom->nmax; i++)
for (int j = 0; j < nprev; ++j)
s_hist[i][j] = t_hist[i][j] = 0;
}
s_hist = t_hist = NULL;
grow_arrays(atom->nmax);
atom->add_callback(0);
for (int i = 0; i < atom->nmax; i++)
for (int j = 0; j < nprev; ++j)
s_hist[i][j] = t_hist[i][j] = 0;
}
/* ---------------------------------------------------------------------- */
@ -557,17 +555,11 @@ void FixQEqReax::init_matvec()
b_s[i] = -chi[ atom->type[i] ];
b_t[i] = -1.0;
/* linear extrapolation for s & t from previous solutions */
//s[i] = 2 * s_hist[i][0] - s_hist[i][1];
//t[i] = 2 * t_hist[i][0] - t_hist[i][1];
/* quadratic extrapolation for s & t from previous solutions */
//s[i] = s_hist[i][2] + 3 * ( s_hist[i][0] - s_hist[i][1] );
t[i] = t_hist[i][2] + 3 * ( t_hist[i][0] - t_hist[i][1]);
/* cubic extrapolation for s & t from previous solutions */
s[i] = 4*(s_hist[i][0]+s_hist[i][2])-(6*s_hist[i][1]+s_hist[i][3]);
//t[i] = 4*(t_hist[i][0]+t_hist[i][2])-(6*t_hist[i][1]+t_hist[i][3]);
}
}
@ -615,6 +607,7 @@ void FixQEqReax::compute_H()
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
dx = x[j][0] - x[i][0];
dy = x[j][1] - x[i][1];