From 4f8e6274790fcc5eb66ed3cb289b9b75bf3daac5 Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 8 Sep 2020 13:16:45 -0600 Subject: [PATCH 1/6] QEq refactor --- src/KOKKOS/fix_qeq_reax_kokkos.cpp | 31 +++-- src/USER-OMP/fix_qeq_reax_omp.cpp | 144 +++++------------------ src/USER-REAXC/fix_qeq_reax.cpp | 178 ++++++++++------------------- src/USER-REAXC/fix_qeq_reax.h | 4 +- 4 files changed, 109 insertions(+), 248 deletions(-) diff --git a/src/KOKKOS/fix_qeq_reax_kokkos.cpp b/src/KOKKOS/fix_qeq_reax_kokkos.cpp index eac044da28..a71b52e0f1 100644 --- a/src/KOKKOS/fix_qeq_reax_kokkos.cpp +++ b/src/KOKKOS/fix_qeq_reax_kokkos.cpp @@ -112,7 +112,7 @@ void FixQEqReaxKokkos::init() ("FixQEqReax::params",ntypes+1); params = k_params.template view(); - 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::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(); 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(); 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(); 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(); h_d = k_d.h_view; } // init_storage + FixQEqReaxKokkosZeroFunctor zero_functor(this); Kokkos::parallel_for(ignum,zero_functor); @@ -779,8 +780,7 @@ void FixQEqReaxKokkos::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::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::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::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); diff --git a/src/USER-OMP/fix_qeq_reax_omp.cpp b/src/USER-OMP/fix_qeq_reax_omp.cpp index ef0f918cbf..92c1b73e13 100644 --- a/src/USER-OMP/fix_qeq_reax_omp.cpp +++ b/src/USER-OMP/fix_qeq_reax_omp.cpp @@ -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 diff --git a/src/USER-REAXC/fix_qeq_reax.cpp b/src/USER-REAXC/fix_qeq_reax.cpp index 4a61cfa03b..28270b562a 100644 --- a/src/USER-REAXC/fix_qeq_reax.cpp +++ b/src/USER-REAXC/fix_qeq_reax.cpp @@ -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]; } } + diff --git a/src/USER-REAXC/fix_qeq_reax.h b/src/USER-REAXC/fix_qeq_reax.h index 96a174b908..5658167aca 100644 --- a/src/USER-REAXC/fix_qeq_reax.h +++ b/src/USER-REAXC/fix_qeq_reax.h @@ -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; From 1273179d03941f43505ead10354518edb12f852e Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Tue, 8 Sep 2020 13:42:25 -0600 Subject: [PATCH 2/6] Fix compile error --- src/USER-REAXC/fix_qeq_reax.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/USER-REAXC/fix_qeq_reax.cpp b/src/USER-REAXC/fix_qeq_reax.cpp index 28270b562a..8e6e114861 100644 --- a/src/USER-REAXC/fix_qeq_reax.cpp +++ b/src/USER-REAXC/fix_qeq_reax.cpp @@ -85,7 +85,7 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) : 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]); + imax = utils::numeric(FLERR,arg[iarg+1],false,lmp); iarg++; } else error->all(FLERR,"Illegal fix qeq/reax command"); iarg++; From bd72ef7996c345098a1327c5a28c4538ae04d049 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 9 Sep 2020 22:03:11 -0400 Subject: [PATCH 3/6] add API to library/python interface to extract the MPI communicator --- doc/src/pg_lib_properties.rst | 5 +++++ python/lammps.py | 24 ++++++++++++++++++++++++ src/library.cpp | 34 ++++++++++++++++++++++++++++++++++ src/library.h | 1 + 4 files changed, 64 insertions(+) diff --git a/doc/src/pg_lib_properties.rst b/doc/src/pg_lib_properties.rst index 789c7e2441..cdc7617d91 100644 --- a/doc/src/pg_lib_properties.rst +++ b/doc/src/pg_lib_properties.rst @@ -21,6 +21,11 @@ event as atoms are migrating between sub-domains. ----------------------- +.. doxygenfunction:: lammps_get_mpi_comm + :project: progguide + +----------------------- + .. doxygenfunction:: lammps_get_natoms :project: progguide diff --git a/python/lammps.py b/python/lammps.py index b76608af7d..6e910b49bf 100644 --- a/python/lammps.py +++ b/python/lammps.py @@ -332,6 +332,8 @@ class lammps(object): self.lib.lammps_version.argtypes = [c_void_p] + self.lib.lammps_get_mpi_comm.argtypes = [c_void_p] + self.lib.lammps_decode_image_flags.argtypes = [self.c_imageint, POINTER(c_int*3)] self.lib.lammps_extract_atom.argtypes = [c_void_p, c_char_p] @@ -601,6 +603,28 @@ class lammps(object): # ------------------------------------------------------------------------- + def get_mpi_comm(self): + """Get the MPI communicator in use by the current LAMMPS instance + + This is a wrapper around the :cpp:func:`lammps_get_mpi_comm` function + of the C-library interface. It will return ``None`` if either the + LAMMPS library was compiled without MPI support or the mpi4py + Python module is not available. + + :return: MPI communicator + :rtype: MPI_Comm + """ + + if self.has_mpi4py and self.has_mpi_support: + from mpi4py import MPI + f_comm = self.lib.lammps_get_mpi_comm(self.lmp) + c_comm = MPI.Comm.f2py(f_comm) + return c_comm + else: + return None + + # ------------------------------------------------------------------------- + def file(self, path): """Read LAMMPS commands from a file. diff --git a/src/library.cpp b/src/library.cpp index 72105bdd15..0e0c4070ea 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -532,6 +532,40 @@ int lammps_version(void *handle) return atoi(lmp->universe->num_ver); } + +/* ---------------------------------------------------------------------- */ + +/** Return current LAMMPS world communicator as integer + * +\verbatim embed:rst + +This will take the LAMMPS "world" communicator and convert it to an +integer using ``MPI_Comm_c2f()``, so it is equivalent to the +corresponding MPI communicator in Fortran. This way it can be safely +passed around between different programming languages. To convert it +to the C language representation use ``MPI_Comm_f2c()``. + +If LAMMPS was compiled with MPI_STUBS, this function returns -1. + +.. versionadded:: 15Sep2020 + +\endverbatim + * \sa lammps_open_fortran + * + * \param handle pointer to a previously created LAMMPS instance + * \return Fortran representation of the LAMMPS world communicator */ + +int lammps_get_mpi_comm(void *handle) +{ +#ifdef MPI_STUBS + return -1; +#else + LAMMPS *lmp = (LAMMPS *) handle; + MPI_Fint f_comm = MPI_Comm_c2f(lmp->world); + return f_comm; +#endif +} + /* ---------------------------------------------------------------------- */ /** Return the total number of atoms in the system. diff --git a/src/library.h b/src/library.h index 0ffd111d7b..f5a02502e7 100644 --- a/src/library.h +++ b/src/library.h @@ -98,6 +98,7 @@ void lammps_commands_string(void *handle, const char *str); * ----------------------------------------------------------------------- */ int lammps_version(void *handle); +int lammps_get_mpi_comm(void* handle); double lammps_get_natoms(void *handle); double lammps_get_thermo(void *handle, char *keyword); From e813e2d30ae8983c03e1c0e11d04252903e6dd86 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 9 Sep 2020 22:12:47 -0400 Subject: [PATCH 4/6] add minimal unit test for lammps_get_mpi_comm() API --- unittest/c-library/test_library_properties.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/unittest/c-library/test_library_properties.cpp b/unittest/c-library/test_library_properties.cpp index 87b4d54f86..5e8c410cce 100644 --- a/unittest/c-library/test_library_properties.cpp +++ b/unittest/c-library/test_library_properties.cpp @@ -41,5 +41,13 @@ protected: } }; +TEST_F(LAMMPS_properties, get_mpi_comm) { + int f_comm = lammps_get_mpi_comm(lmp); + if (lammps_config_has_mpi_support()) + EXPECT_GE(f_comm,0); + else + EXPECT_EQ(f_comm,-1); +}; + TEST_F(LAMMPS_properties, box) { }; From 4a8d6016e412049e012359310c7d7da37478547c Mon Sep 17 00:00:00 2001 From: Stan Moore Date: Fri, 11 Sep 2020 22:49:44 -0600 Subject: [PATCH 5/6] Update docs --- doc/src/fix_qeq_reax.rst | 8 ++++++-- src/USER-REAXC/fix_qeq_reax.cpp | 2 +- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/doc/src/fix_qeq_reax.rst b/doc/src/fix_qeq_reax.rst index 47dca86401..512c596b4d 100644 --- a/doc/src/fix_qeq_reax.rst +++ b/doc/src/fix_qeq_reax.rst @@ -21,6 +21,7 @@ Syntax * tolerance = precision to which charges will be equilibrated * params = reax/c or a filename * args = *dual* (optional) +* args = *maxiter value* (optional) Examples """""""" @@ -28,7 +29,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 +70,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 +106,7 @@ Related commands Default """"""" -none +maxiter 200 ---------- diff --git a/src/USER-REAXC/fix_qeq_reax.cpp b/src/USER-REAXC/fix_qeq_reax.cpp index 8e6e114861..fea8b39f54 100644 --- a/src/USER-REAXC/fix_qeq_reax.cpp +++ b/src/USER-REAXC/fix_qeq_reax.cpp @@ -82,7 +82,7 @@ FixQEqReax::FixQEqReax(LAMMPS *lmp, int narg, char **arg) : int iarg = 8; while (iarg < narg) { if (strcmp(arg[iarg],"dual") == 0) dual_enabled = 1; - else if (strcmp(arg[iarg],"imax") == 0) { + 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); From da0cdb0de42bcd06c6367a719402d7a4303b49fc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 12 Sep 2020 19:06:30 -0400 Subject: [PATCH 6/6] update formatting of keyword summary --- doc/src/fix_qeq_reax.rst | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/doc/src/fix_qeq_reax.rst b/doc/src/fix_qeq_reax.rst index 512c596b4d..8752888c4c 100644 --- a/doc/src/fix_qeq_reax.rst +++ b/doc/src/fix_qeq_reax.rst @@ -20,8 +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) -* args = *maxiter value* (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 """"""""