QEq refactor (4f8e6274) · Commits · 郑智淋 / lammps

src/KOKKOS/fix_qeq_reax_kokkos.cpp

+15 −16

Original line number	Diff line number	Diff line
		@@ -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);

src/USER-OMP/fix_qeq_reax_omp.cpp

+29 −115

Original line number	Diff line number	Diff line
		@@ -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

src/USER-REAXC/fix_qeq_reax.cpp

+62 −116

Original line number	Diff line number	Diff line
		@@ -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],"imax") == 0) {
		if (iarg+1 > narg-1)
		error->all(FLERR,"Illegal fix qeq/reax command");
		imax = force->numeric(FLERR,arg[iarg+1]);
		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];
		}
		}

src/USER-REAXC/fix_qeq_reax.h

+3 −1

Original line number	Diff line number	Diff line
		@@ -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;

Admin message