Merge remote-tracking branch 'origin/master' into python_interface_guide

doc/src/fix_qeq_reax.rst

@@ -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
 
 ----------
 

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
 

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.
 

src/KOKKOS/fix_qeq_reax_kokkos.cpp

@@ -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 < imax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
-  for (loop = 1; (loop < loopmax) && (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 < imax) && (sqrt(sig_new)/b_norm > tolerance); loop++) {
-  for (loop = 1; (loop < loopmax) && (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

@@ -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;
+  int n = atom->nlocal;
-  N = atom->nlocal + atom->nghost;
+
+  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 i;
-  int *ilist;
-
-  if (reaxc) {
-    nn = reaxc->list->inum;
-    ilist = reaxc->list->ilist;
-  } else {
-    nn = list->inum;
-    ilist = list->ilist;
-  }
 
   // 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 ii;
-    int *ilist;
     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 ii;
-    int *ilist;
     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 ii;
-    int *ilist;
     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,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>11) error->all(FLERR,"Illegal fix qeq/reax command");
-
-  if (narg<8 || narg>9) 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) {
+  imax = 200;
-    if (strcmp(arg[8],"dual") == 0) dual_enabled = 1;
+
-    else error->all(FLERR,"Illegal fix qeq/reax command");
+  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;
+  nn = n_cap = 0;
-  N = nmax = 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 i,ii,n,m;
-  int *ilist, *numneigh;
 
   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;
+  int n = atom->nlocal;
-  N = atom->nlocal + atom->nghost;
+
+  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 ii, i;
-  int *ilist;
-
-  if (reaxc) {
-    nn = reaxc->list->inum;
-    ilist = reaxc->list->ilist;
-  } else {
-    nn = list->inum;
-    ilist = list->ilist;
-  }
 
   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 jj;
-  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);
@@ -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 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;
-  }
 
   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 ii;
-  int *ilist;
-
-  if (reaxc) {
-    nn = reaxc->list->inum;
-    ilist = reaxc->list->ilist;
-  } else {
-    nn = list->inum;
-    ilist = list->ilist;
-  }
 
   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

@@ -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;

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.

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);
 

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) {
 };