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

Merge pull request #2354 from stanmoore1/qeq_refactor

Browse Source 
small fix qeq/reax refactoring for base and accelerated variants
This commit is contained in:
Axel Kohlmeyer
2020-09-14 12:18:32 -04:00
committed by GitHub
parent 101d39142e da0cdb0de4
commit 50b80e078d
5 changed files with 122 additions and 251 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
doc/src/fix_qeq_reax.rst
View File

@ -20,7 +20,14 @@ Syntax
* cutlo,cuthi = lo and hi cutoff for Taper radius
* tolerance = precision to which charges will be equilibrated
* params = reax/c or a filename
* args = *dual* (optional)
* one or more keywords or keyword/value pairs may be appended
.. parsed-literal::
keyword = *dual* or *maxiter*
*dual* = process S and T matrix in parallel (only for qeq/reax/omp)
*maxiter* N = limit the number of iterations to *N*
Examples
""""""""
@ -28,7 +35,7 @@ Examples
.. code-block:: LAMMPS
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 reax/c
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 param.qeq
fix 1 all qeq/reax 1 0.0 10.0 1.0e-6 param.qeq maxiter 500
Description
"""""""""""
@ -69,6 +76,9 @@ The optional *dual* keyword allows to perform the optimization
of the S and T matrices in parallel. This is only supported for
the *qeq/reax/omp* style. Otherwise they are processed separately.
The optional *maxiter* keyword allows changing the max number
of iterations in the linear solver. The default value is 200.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
@ -102,7 +112,7 @@ Related commands
Default
"""""""
none
maxiter 200
----------

31
src/KOKKOS/fix_qeq_reax_kokkos.cpp
View File

@ -112,7 +112,7 @@ void FixQEqReaxKokkos<DeviceType>::init()
("FixQEqReax::params",ntypes+1);
params = k_params.template view<DeviceType>();
for (n = 1; n <= ntypes; n++) {
for (int n = 1; n <= ntypes; n++) {
k_params.h_view(n).chi = chi[n];
k_params.h_view(n).eta = eta[n];
k_params.h_view(n).gamma = gamma[n];
@ -351,34 +351,35 @@ void FixQEqReaxKokkos<DeviceType>::allocate_array()
if (atom->nmax > nmax) {
nmax = atom->nmax;
k_o = DAT::tdual_ffloat_1d("qeq/kk:h_o",nmax);
k_o = DAT::tdual_ffloat_1d("qeq/kk:o",nmax);
d_o = k_o.template view<DeviceType>();
h_o = k_o.h_view;
d_Hdia_inv = typename AT::t_ffloat_1d("qeq/kk:h_Hdia_inv",nmax);
d_Hdia_inv = typename AT::t_ffloat_1d("qeq/kk:Hdia_inv",nmax);
d_b_s = typename AT::t_ffloat_1d("qeq/kk:h_b_s",nmax);
d_b_s = typename AT::t_ffloat_1d("qeq/kk:b_s",nmax);
d_b_t = typename AT::t_ffloat_1d("qeq/kk:h_b_t",nmax);
d_b_t = typename AT::t_ffloat_1d("qeq/kk:b_t",nmax);
k_s = DAT::tdual_ffloat_1d("qeq/kk:h_s",nmax);
k_s = DAT::tdual_ffloat_1d("qeq/kk:s",nmax);
d_s = k_s.template view<DeviceType>();
h_s = k_s.h_view;
k_t = DAT::tdual_ffloat_1d("qeq/kk:h_t",nmax);
k_t = DAT::tdual_ffloat_1d("qeq/kk:t",nmax);
d_t = k_t.template view<DeviceType>();
h_t = k_t.h_view;
d_p = typename AT::t_ffloat_1d("qeq/kk:h_p",nmax);
d_p = typename AT::t_ffloat_1d("qeq/kk:p",nmax);
d_r = typename AT::t_ffloat_1d("qeq/kk:h_r",nmax);
d_r = typename AT::t_ffloat_1d("qeq/kk:r",nmax);
k_d = DAT::tdual_ffloat_1d("qeq/kk:h_d",nmax);
k_d = DAT::tdual_ffloat_1d("qeq/kk:d",nmax);
d_d = k_d.template view<DeviceType>();
h_d = k_d.h_view;
}
// init_storage
FixQEqReaxKokkosZeroFunctor<DeviceType> zero_functor(this);
Kokkos::parallel_for(ignum,zero_functor);
@ -779,8 +780,7 @@ void FixQEqReaxKokkos<DeviceType>::cg_solve1()
F_FLOAT sig_new = dot_sqr;
int loop;
const int loopmax = 200;
for (loop = 1; (loop < loopmax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
for (loop = 1; (loop < imax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
// comm->forward_comm_fix(this); //Dist_vector( d );
pack_flag = 1;
@ -848,7 +848,7 @@ void FixQEqReaxKokkos<DeviceType>::cg_solve1()
Kokkos::parallel_for(inum,vecsum2_functor);
}
if (loop >= loopmax && comm->me == 0) {
if (loop >= imax && comm->me == 0) {
char str[128];
sprintf(str,"Fix qeq/reax cg_solve1 convergence failed after %d iterations "
"at " BIGINT_FORMAT " step: %f",loop,update->ntimestep,sqrt(sig_new)/b_norm);
@ -918,8 +918,7 @@ void FixQEqReaxKokkos<DeviceType>::cg_solve2()
F_FLOAT sig_new = dot_sqr;
int loop;
const int loopmax = 200;
for (loop = 1; (loop < loopmax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
for (loop = 1; (loop < imax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
// comm->forward_comm_fix(this); //Dist_vector( d );
pack_flag = 1;
@ -987,7 +986,7 @@ void FixQEqReaxKokkos<DeviceType>::cg_solve2()
Kokkos::parallel_for(inum,vecsum2_functor);
}
if (loop >= loopmax && comm->me == 0) {
if (loop >= imax && comm->me == 0) {
char str[128];
sprintf(str,"Fix qeq/reax cg_solve2 convergence failed after %d iterations "
"at " BIGINT_FORMAT " step: %f",loop,update->ntimestep,sqrt(sig_new)/b_norm);

144
src/USER-OMP/fix_qeq_reax_omp.cpp
View File

@ -63,8 +63,6 @@ using namespace FixConst;
FixQEqReaxOMP::FixQEqReaxOMP(LAMMPS *lmp, int narg, char **arg) :
FixQEqReax(lmp, narg, arg)
{
if (narg<8 || narg>9) error->all(FLERR,"Illegal fix qeq/reax/omp command");
b_temp = NULL;
// ASPC: Kolafa, J. Comp. Chem., 25(3), 335 (2003)
@ -85,6 +83,11 @@ FixQEqReaxOMP::~FixQEqReaxOMP()
void FixQEqReaxOMP::post_constructor()
{
grow_arrays(atom->nmax);
for (int i = 0; i < atom->nmax; i++)
for (int j = 0; j < nprev; ++j)
s_hist[i][j] = t_hist[i][j] = 0;
pertype_parameters(pertype_option);
}
@ -148,7 +151,6 @@ void FixQEqReaxOMP::init()
void FixQEqReaxOMP::compute_H()
{
int inum, *ilist, *numneigh, **firstneigh;
double SMALL = 0.0001;
int *type = atom->type;
@ -156,17 +158,6 @@ void FixQEqReaxOMP::compute_H()
double **x = atom->x;
int *mask = atom->mask;
if (reaxc) {
inum = reaxc->list->inum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
firstneigh = reaxc->list->firstneigh;
} else {
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
}
int ai, num_nbrs;
// sumscan of the number of neighbors per atom to determine the offsets
@ -175,7 +166,7 @@ void FixQEqReaxOMP::compute_H()
num_nbrs = 0;
for (int itr_i = 0; itr_i < inum; ++itr_i) {
for (int itr_i = 0; itr_i < nn; ++itr_i) {
ai = ilist[itr_i];
H.firstnbr[ai] = num_nbrs;
num_nbrs += numneigh[ai];
@ -197,7 +188,7 @@ void FixQEqReaxOMP::compute_H()
#if defined(_OPENMP)
#pragma omp for schedule(guided)
#endif
for (int ii = 0; ii < inum; ii++) {
for (int ii = 0; ii < nn; ii++) {
int i = ilist[ii];
if (mask[i] & groupbit) {
jlist = firstneigh[i];
@ -214,7 +205,7 @@ void FixQEqReaxOMP::compute_H()
flag = 0;
if (r_sqr <= SQR(swb)) {
if (j < n) flag = 1;
if (j < atom->nlocal) flag = 1;
else if (tag[i] < tag[j]) flag = 1;
else if (tag[i] == tag[j]) {
if (dz > SMALL) flag = 1;
@ -251,11 +242,6 @@ void FixQEqReaxOMP::compute_H()
void FixQEqReaxOMP::init_storage()
{
int NN;
if (reaxc) NN = reaxc->list->inum + reaxc->list->gnum;
else NN = list->inum + list->gnum;
#if defined(_OPENMP)
#pragma omp parallel for schedule(static)
#endif
@ -284,8 +270,21 @@ void FixQEqReaxOMP::pre_force(int /* vflag */)
if (update->ntimestep % nevery) return;
if (comm->me == 0) t_start = MPI_Wtime();
n = atom->nlocal;
N = atom->nlocal + atom->nghost;
int n = atom->nlocal;
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
firstneigh = reaxc->list->firstneigh;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
}
// grow arrays if necessary
// need to be atom->nmax in length
@ -365,16 +364,7 @@ void FixQEqReaxOMP::init_matvec()
/* fill-in H matrix */
compute_H();
int nn,i;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
int i;
// Should really be more careful with initialization and first (aspc_order+2) MD steps
if (do_aspc) {
@ -450,24 +440,12 @@ void FixQEqReaxOMP::init_matvec()
int FixQEqReaxOMP::CG( double *b, double *x)
{
int i, imax;
int i;
double alpha, beta, b_norm;
double sig_old, sig_new;
double my_buf[2], buf[2];
int nn;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
imax = 200;
pack_flag = 1;
sparse_matvec( &H, x, q );
comm->reverse_comm_fix( this); //Coll_Vector( q );
@ -579,8 +557,7 @@ void FixQEqReaxOMP::sparse_matvec( sparse_matrix *A, double *x, double *b)
#endif
{
int i, j, itr_j;
int nn, NN, ii;
int *ilist;
int ii;
int nthreads = comm->nthreads;
#if defined(_OPENMP)
int tid = omp_get_thread_num();
@ -588,16 +565,6 @@ void FixQEqReaxOMP::sparse_matvec( sparse_matrix *A, double *x, double *b)
int tid = 0;
#endif
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
}
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif
@ -655,17 +622,6 @@ void FixQEqReaxOMP::calculate_Q()
int i;
double *q = atom->q;
int nn;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
double tmp1, tmp2;
tmp1 = tmp2 = 0.0;
#if defined(_OPENMP)
@ -718,10 +674,6 @@ void FixQEqReaxOMP::vector_sum( double* dest, double c, double* v,
double d, double* y, int k)
{
int i;
int *ilist;
if (reaxc) ilist = reaxc->list->ilist;
else ilist = list->ilist;
#if defined(_OPENMP)
#pragma omp parallel for schedule(static) private(i)
@ -737,10 +689,6 @@ void FixQEqReaxOMP::vector_sum( double* dest, double c, double* v,
void FixQEqReaxOMP::vector_add( double* dest, double c, double* v, int k)
{
int i;
int *ilist;
if (reaxc) ilist = reaxc->list->ilist;
else ilist = list->ilist;
#if defined(_OPENMP)
#pragma omp parallel for schedule(static) private(i)
@ -765,24 +713,12 @@ int FixQEqReaxOMP::dual_CG( double *b1, double *b2, double *x1, double *x2)
startTimeBase = MPI_Wtime();
#endif
int i, imax;
int i;
double alpha_s, alpha_t, beta_s, beta_t, b_norm_s, b_norm_t;
double sig_old_s, sig_old_t, sig_new_s, sig_new_t;
double my_buf[4], buf[4];
int nn;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
imax = 200;
pack_flag = 5; // forward 2x d and reverse 2x q
dual_sparse_matvec( &H, x1, x2, q );
comm->reverse_comm_fix(this); //Coll_Vector( q );
@ -975,8 +911,7 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x1, double *x2
#endif
{
int i, j, itr_j;
int nn, NN, ii;
int *ilist;
int ii;
int indxI, indxJ;
int nthreads = comm->nthreads;
@ -986,16 +921,6 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x1, double *x2
int tid = 0;
#endif
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
}
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif
@ -1077,8 +1002,7 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x, double *b )
#endif
{
int i, j, itr_j;
int nn, NN, ii;
int *ilist;
int ii;
int indxI, indxJ;
int nthreads = comm->nthreads;
@ -1088,16 +1012,6 @@ void FixQEqReaxOMP::dual_sparse_matvec( sparse_matrix *A, double *x, double *b )
int tid = 0;
#endif
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
}
#if defined(_OPENMP)
#pragma omp for schedule(dynamic,50)
#endif

178
src/USER-REAXC/fix_qeq_reax.cpp
View File

@ -62,9 +62,7 @@ static const char cite_fix_qeq_reax[] =
FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg), pertype_option(NULL)
{
if (lmp->citeme) lmp->citeme->add(cite_fix_qeq_reax);
if (narg<8 || narg>9) error->all(FLERR,"Illegal fix qeq/reax command");
if (narg<8 || narg>11) error->all(FLERR,"Illegal fix qeq/reax command");
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery <= 0) error->all(FLERR,"Illegal fix qeq/reax command");
@ -79,14 +77,23 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) :
// dual CG support only available for USER-OMP variant
// check for compatibility is in Fix::post_constructor()
dual_enabled = 0;
if (narg == 9) {
if (strcmp(arg[8],"dual") == 0) dual_enabled = 1;
else error->all(FLERR,"Illegal fix qeq/reax command");
imax = 200;
int iarg = 8;
while (iarg < narg) {
if (strcmp(arg[iarg],"dual") == 0) dual_enabled = 1;
else if (strcmp(arg[iarg],"maxiter") == 0) {
if (iarg+1 > narg-1)
error->all(FLERR,"Illegal fix qeq/reax command");
imax = utils::numeric(FLERR,arg[iarg+1],false,lmp);
iarg++;
} else error->all(FLERR,"Illegal fix qeq/reax command");
iarg++;
}
shld = NULL;
n = n_cap = 0;
N = nmax = 0;
nn = n_cap = 0;
NN = nmax = 0;
m_fill = m_cap = 0;
pack_flag = 0;
s = NULL;
@ -123,11 +130,7 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) :
reaxc = (PairReaxC *) force->pair_match("^reax/c",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;
}
/* ---------------------------------------------------------------------- */
@ -161,6 +164,13 @@ FixQEqReax::~FixQEqReax()
void FixQEqReax::post_constructor()
{
if (lmp->citeme) lmp->citeme->add(cite_fix_qeq_reax);
grow_arrays(atom->nmax);
for (int i = 0; i < atom->nmax; i++)
for (int j = 0; j < nprev; ++j)
s_hist[i][j] = t_hist[i][j] = 0;
pertype_parameters(pertype_option);
if (dual_enabled)
error->all(FLERR,"Dual keyword only supported with fix qeq/reax/omp");
@ -287,8 +297,7 @@ void FixQEqReax::reallocate_storage()
void FixQEqReax::allocate_matrix()
{
int i,ii,inum,m;
int *ilist, *numneigh;
int i,ii,n,m;
int mincap;
double safezone;
@ -306,18 +315,8 @@ void FixQEqReax::allocate_matrix()
// determine the total space for the H matrix
if (reaxc) {
inum = reaxc->list->inum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
} else {
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
}
m = 0;
for (ii = 0; ii < inum; ii++) {
for (ii = 0; ii < nn; ii++) {
i = ilist[ii];
m += numneigh[i];
}
@ -432,6 +431,20 @@ void FixQEqReax::init_taper()
void FixQEqReax::setup_pre_force(int vflag)
{
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
firstneigh = reaxc->list->firstneigh;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
}
deallocate_storage();
allocate_storage();
@ -495,8 +508,21 @@ void FixQEqReax::pre_force(int /*vflag*/)
if (update->ntimestep % nevery) return;
if (comm->me == 0) t_start = MPI_Wtime();
n = atom->nlocal;
N = atom->nlocal + atom->nghost;
int n = atom->nlocal;
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
firstneigh = reaxc->list->firstneigh;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
}
// grow arrays if necessary
// need to be atom->nmax in length
@ -540,16 +566,7 @@ void FixQEqReax::init_matvec()
/* fill-in H matrix */
compute_H();
int nn, ii, i;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
int ii, i;
for (ii = 0; ii < nn; ++ii) {
i = ilist[ii];
@ -578,7 +595,7 @@ void FixQEqReax::init_matvec()
void FixQEqReax::compute_H()
{
int inum, jnum, *ilist, *jlist, *numneigh, **firstneigh;
int jnum;
int i, j, ii, jj, flag;
double dx, dy, dz, r_sqr;
const double SMALL = 0.0001;
@ -588,22 +605,10 @@ void FixQEqReax::compute_H()
double **x = atom->x;
int *mask = atom->mask;
if (reaxc) {
inum = reaxc->list->inum;
ilist = reaxc->list->ilist;
numneigh = reaxc->list->numneigh;
firstneigh = reaxc->list->firstneigh;
} else {
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
}
// fill in the H matrix
m_fill = 0;
r_sqr = 0;
for (ii = 0; ii < inum; ii++) {
for (ii = 0; ii < nn; ii++) {
i = ilist[ii];
if (mask[i] & groupbit) {
jlist = firstneigh[i];
@ -621,7 +626,7 @@ void FixQEqReax::compute_H()
flag = 0;
if (r_sqr <= SQR(swb)) {
if (j < n) flag = 1;
if (j < atom->nlocal) flag = 1;
else if (tag[i] < tag[j]) flag = 1;
else if (tag[i] == tag[j]) {
if (dz > SMALL) flag = 1;
@ -676,21 +681,11 @@ double FixQEqReax::calculate_H( double r, double gamma)
int FixQEqReax::CG( double *b, double *x)
{
int i, j, imax;
int i, j;
double tmp, alpha, beta, b_norm;
double sig_old, sig_new;
int nn, jj;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
imax = 200;
int jj;
pack_flag = 1;
sparse_matvec( &H, x, q);
@ -748,18 +743,7 @@ int FixQEqReax::CG( double *b, double *x)
void FixQEqReax::sparse_matvec( sparse_matrix *A, double *x, double *b)
{
int i, j, itr_j;
int nn, NN, ii;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
NN = reaxc->list->inum + reaxc->list->gnum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
NN = list->inum + list->gnum;
ilist = list->ilist;
}
int ii;
for (ii = 0; ii < nn; ++ii) {
i = ilist[ii];
@ -794,16 +778,7 @@ void FixQEqReax::calculate_Q()
double u, s_sum, t_sum;
double *q = atom->q;
int nn, ii;
int *ilist;
if (reaxc) {
nn = reaxc->list->inum;
ilist = reaxc->list->ilist;
} else {
nn = list->inum;
ilist = list->ilist;
}
int ii;
s_sum = parallel_vector_acc( s, nn);
t_sum = parallel_vector_acc( t, nn);
@ -988,12 +963,6 @@ double FixQEqReax::parallel_norm( double *v, int n)
double my_sum, norm_sqr;
int ii;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
my_sum = 0.0;
norm_sqr = 0.0;
@ -1016,12 +985,6 @@ double FixQEqReax::parallel_dot( double *v1, double *v2, int n)
double my_dot, res;
int ii;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
my_dot = 0.0;
res = 0.0;
@ -1044,12 +1007,6 @@ double FixQEqReax::parallel_vector_acc( double *v, int n)
double my_acc, res;
int ii;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
my_acc = 0.0;
res = 0.0;
@ -1070,12 +1027,6 @@ void FixQEqReax::vector_sum( double* dest, double c, double* v,
double d, double* y, int k)
{
int kk;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
for (--k; k>=0; --k) {
kk = ilist[k];
@ -1089,12 +1040,6 @@ void FixQEqReax::vector_sum( double* dest, double c, double* v,
void FixQEqReax::vector_add( double* dest, double c, double* v, int k)
{
int kk;
int *ilist;
if (reaxc)
ilist = reaxc->list->ilist;
else
ilist = list->ilist;
for (--k; k>=0; --k) {
kk = ilist[k];
@ -1102,3 +1047,4 @@ void FixQEqReax::vector_add( double* dest, double c, double* v, int k)
dest[kk] += c * v[kk];
}
}

4
src/USER-REAXC/fix_qeq_reax.h
View File

@ -57,12 +57,13 @@ class FixQEqReax : public Fix {
protected:
int nevery,reaxflag;
int n, N, m_fill;
int nn, NN, m_fill;
int n_cap, nmax, m_cap;
int pack_flag;
int nlevels_respa;
class NeighList *list;
class PairReaxC *reaxc;
int *ilist, *jlist, *numneigh, **firstneigh;
double swa, swb; // lower/upper Taper cutoff radius
double Tap[8]; // Taper function
@ -94,6 +95,7 @@ class FixQEqReax : public Fix {
//CG storage
double *p, *q, *r, *d;
int imax;
//GMRES storage
//double *g,*y;