diff --git a/PG-PuReMD/src/bonds.c b/PG-PuReMD/src/bonds.c
index fc125aaba649cd0109aec8e6a8d14b432e7120b5..27b7654e60a7ae2a8a33ad8cab9cc8f9c50e763d 100644
--- a/PG-PuReMD/src/bonds.c
+++ b/PG-PuReMD/src/bonds.c
@@ -51,7 +51,7 @@ void Bonds( reax_system * const system, control_params * const control,
     int start_i, end_i;
     int type_i, type_j;
     real ebond, pow_BOs_be2, exp_be12, CEbo;
-    real gp3, gp4, gp7, gp10, gp37;
+    real gp3, gp4, gp7, gp10;
     real exphu, exphua1, exphub1, exphuov, hulpov, estriph;
     real decobdbo, decobdboua, decobdboub;
     single_body_parameters *sbp_i, *sbp_j;
@@ -64,7 +64,6 @@ void Bonds( reax_system * const system, control_params * const control,
     gp4 = system->reax_param.gp.l[4];
     gp7 = system->reax_param.gp.l[7];
     gp10 = system->reax_param.gp.l[10];
-    gp37 = (int) system->reax_param.gp.l[37];
 
     for ( i = 0; i < system->n; ++i )
     {
@@ -120,11 +119,10 @@ void Bonds( reax_system * const system, control_params * const control,
                 /* Stabilisation terminal triple bond in C-O */
                 if ( bo_ij->BO >= 1.00 )
                 {
-                    if ( gp37 == 2 &&
-                            ( (strncmp( sbp_i->name, "C", sizeof(sbp_i->name) ) == 0
+                    if ( (strncmp( sbp_i->name, "C", sizeof(sbp_i->name) ) == 0
                                 && strncmp( sbp_j->name, "O", sizeof(sbp_j->name) ) == 0)
                             || (strncmp( sbp_i->name, "O", sizeof(sbp_i->name) ) == 0
-                                && strncmp( sbp_j->name, "C", sizeof(sbp_j->name) ) == 0) ) )
+                                && strncmp( sbp_j->name, "C", sizeof(sbp_j->name) ) == 0) )
                     {
                         //ba = SQR( bo_ij->BO - 2.5 );
                         exphu = EXP( -gp7 * SQR(bo_ij->BO - 2.5) );
diff --git a/PG-PuReMD/src/cuda/cuda_bonds.cu b/PG-PuReMD/src/cuda/cuda_bonds.cu
index d5ff0cac852532213c11f392f30a8c02c2dab7dd..a52df4c8b999cab1b27d8e96869562d123efbbfb 100644
--- a/PG-PuReMD/src/cuda/cuda_bonds.cu
+++ b/PG-PuReMD/src/cuda/cuda_bonds.cu
@@ -38,7 +38,7 @@ CUDA_GLOBAL void k_bonds( reax_atom *my_atoms, global_parameters gp,
     int start_i, end_i;
     int type_i, type_j;
     real pow_BOs_be2, exp_be12, CEbo, e_bond_l;
-    real gp3, gp4, gp7, gp10, gp37;
+    real gp3, gp4, gp7, gp10;
     real exphu, exphua1, exphub1, exphuov, hulpov;
     real decobdbo, decobdboua, decobdboub;
     single_body_parameters *sbp_i, *sbp_j;
@@ -60,7 +60,6 @@ CUDA_GLOBAL void k_bonds( reax_atom *my_atoms, global_parameters gp,
     gp4 = gp.l[4];
     gp7 = gp.l[7];
     gp10 = gp.l[10];
-    gp37 = (int) gp.l[37];
     e_bond_l = 0.0;
 
     start_i = Start_Index( i, bond_list );
@@ -97,11 +96,10 @@ CUDA_GLOBAL void k_bonds( reax_atom *my_atoms, global_parameters gp,
             /* Stabilisation terminal triple bond */
             if ( bo_ij->BO >= 1.00 )
             {
-                if ( gp37 == 2
-                        || ( (Cuda_strncmp( sbp_i->name, "C", sizeof(sbp_i->name) ) == 0
+                if ( (Cuda_strncmp( sbp_i->name, "C", sizeof(sbp_i->name) ) == 0
                             && Cuda_strncmp( sbp_j->name, "O", sizeof(sbp_j->name) ) == 0)
                         || (Cuda_strncmp( sbp_i->name, "O", sizeof(sbp_i->name) ) == 0
-                            && Cuda_strncmp( sbp_j->name, "C", sizeof(sbp_j->name) ) == 0) ) )
+                            && Cuda_strncmp( sbp_j->name, "C", sizeof(sbp_j->name) ) == 0) )
                 {
                     //ba = SQR( bo_ij->BO - 2.5 );
                     exphu = EXP( -gp7 * SQR(bo_ij->BO - 2.5) );
diff --git a/PG-PuReMD/src/reax_types.h b/PG-PuReMD/src/reax_types.h
index b2936f8104d3d6b9e1db3331f7c95d8f79295adf..619a6c71bad6b867b0c72ba4399918f4be8756fb 100644
--- a/PG-PuReMD/src/reax_types.h
+++ b/PG-PuReMD/src/reax_types.h
@@ -874,7 +874,7 @@ struct mpi_datatypes
  * l[34] = b_s_acks2 (ACKS2 bond softness)
  * l[35] = N/A
  * l[36] = N/A
- * l[37] = version number
+ * l[37] = N/A
  * l[38] = p_coa3
  * */
 struct global_parameters
diff --git a/PuReMD/src/bonds.c b/PuReMD/src/bonds.c
index f5ab71b1dcfea3fdd78277490519807d2e343227..629c7f20e7b1f5e2ed1d1b06049abdcf61919ba4 100644
--- a/PuReMD/src/bonds.c
+++ b/PuReMD/src/bonds.c
@@ -43,7 +43,7 @@ void Bonds( reax_system *system, control_params *control,
     int start_i, end_i;
     int type_i, type_j;
     real ebond, pow_BOs_be2, exp_be12, CEbo;
-    real gp3, gp4, gp7, gp10, gp37;
+    real gp3, gp4, gp7, gp10;
     real exphu, exphua1, exphub1, exphuov, hulpov, estriph;
     real decobdbo, decobdboua, decobdboub;
     single_body_parameters *sbp_i, *sbp_j;
@@ -56,7 +56,6 @@ void Bonds( reax_system *system, control_params *control,
     gp4 = system->reax_param.gp.l[4];
     gp7 = system->reax_param.gp.l[7];
     gp10 = system->reax_param.gp.l[10];
-    gp37 = (int) system->reax_param.gp.l[37];
     natoms = system->n;
 
     for ( i = 0; i < natoms; ++i )
@@ -113,11 +112,10 @@ void Bonds( reax_system *system, control_params *control,
                 /* Stabilisation terminal triple bond in C-O */
                 if ( bo_ij->BO >= 1.00 )
                 {
-                    if ( gp37 == 2 &&
-                            ( (strncmp( sbp_i->name, "C", sizeof(sbp_i->name) ) == 0
+                    if ( (strncmp( sbp_i->name, "C", sizeof(sbp_i->name) ) == 0
                                 && strncmp( sbp_j->name, "O", sizeof(sbp_j->name) ) == 0)
                             || (strncmp( sbp_i->name, "O", sizeof(sbp_i->name) ) == 0
-                                && strncmp( sbp_j->name, "C", sizeof(sbp_j->name) ) == 0) ) )
+                                && strncmp( sbp_j->name, "C", sizeof(sbp_j->name) ) == 0) )
                     {
                         //ba = SQR( bo_ij->BO - 2.5 );
                         exphu = EXP( -gp7 * SQR(bo_ij->BO - 2.5) );
diff --git a/PuReMD/src/reax_types.h b/PuReMD/src/reax_types.h
index 93d7ac3cb1d869882aa8a913c3a020189387f5de..f2cb2ead66a0d32606f92c6cc7a828c9668a74b9 100644
--- a/PuReMD/src/reax_types.h
+++ b/PuReMD/src/reax_types.h
@@ -591,7 +591,7 @@ struct mpi_datatypes
  * l[34] = b_s_acks2 (ACKS2 bond softness)
  * l[35] = N/A
  * l[36] = N/A
- * l[37] = version number
+ * l[37] = N/A
  * l[38] = p_coa3 */
 struct global_parameters
 {
diff --git a/README.md b/README.md
index 33712d54553fac5ae69a1522bd89536f5a1e9f18..bc289cc0762fa3f949f41c9c193d6d25e4252942 100644
--- a/README.md
+++ b/README.md
@@ -6,64 +6,84 @@
 
 # Introduction
 
-This repository contains the development version of the
-[Purdue Reactive Molecular Dynamics](https://www.cs.purdue.edu/puremd) (PuReMD) project.
+This repository contains the development version of the **Pu**ReMD **Re**active
+**M**olecular **D**ynamics (PuReMD) project.  PuReMD is an open-source highly
+performant range-limited atomic-level molecular dynamics code which implements
+the the reactive force field (ReaxFF) method coupled with a global atomic
+charge model.  Supported charge models include charge equilibration,
+electronegativity equilization, and atom-condensed Kohn-Sham approximated to
+second order.
 
 # Build Instructions
 
-## Developer
+## User
 
 To build, the following versions of software are required:
 
-- git
-- Autoconf v2.69+
-- Automake v1.15+
-- OpenMP v4.0+ compliant compiler (OpenMP versions only)
-- MPI v2+ compliant library (MPI versions only)
-- CUDA v6.0+ (CUDA versions only)
+- GNU make
+- C compiler with support for the c11 standard or newer and optionally OpenMP v4.0+ (shared-memory code)
+- C++ compiler with support for the c++14 standard or newer (CUDA, MPI+CUDA versions)
+- MPI v2+ compliant library (MPI, MPI+CUDA versions)
+- CUDA v6.0+ (CUDA, MPI+CUDA versions)
+- zlib v1.2.x or newer
 
-Instructions:
 ```bash
-git clone https://gitlab.msu.edu/SParTA/PuReMD.git
-cd PuReMD
-git submodule init
-git submodule update
-autoreconf -ivf
+# Download release tarball
+tar -xvf puremd-1.0.tar.gz
+cd puremd-1.0
 ./configure
 make
 ```
 
-To build tarball releases after configuring a specific build target, run the following:
+By default, the shared memory version with OpenMP support will be built. For other build targets,
+run ./configure --help and consult the documentation. An example of building the MPI+CUDA version
+is given below.
 
 ```bash
-make dist
+./configure --enable-openmp=no --enable-mpi-gpu=yes
 ```
 
-## User
+## Developer
 
+To build, the following versions of software are required:
+
+- git
+- Autoconf v2.69 or newer
+- Automake v1.15 or newer
+- libtool v2.2 or newer
+- GNU make
+- C compiler with support for the c11 standard or newer and optionally OpenMP v4.0+ (shared-memory code)
+- C++ compiler with support for the c++14 standard or newer (CUDA, MPI+CUDA versions)
+- MPI v2+ compliant library (MPI, MPI+CUDA versions)
+- CUDA v6.0+ (CUDA, MPI+CUDA versions)
+- zlib v1.2.x or newer
+
+Instructions:
 ```bash
-# Download release tarball
-tar -xvf puremd-1.0.tar.gz
-cd puremd-1.0
+git clone https://gitlab.msu.edu/SParTA/PuReMD.git
+cd PuReMD
+git submodule init
+git submodule update
+autoreconf -ivf
 ./configure
 make
 ```
 
-By default, the shared memory version with OpenMP support will be built. For other build targets,
-run ./configure --help and consult the documentation. An example of building the MPI+CUDA version
-is given below.
+To build tarball releases after configuring a specific build target, run the following:
 
 ```bash
-./configure --enable-openmp=no --enable-mpi-gpu=yes
+make dist
 ```
 
 # References
 
-Shared Memory:
+Shared-Memory Versions:
 - [Serial](https://www.cs.purdue.edu/puremd/docs/80859.pdf)
 - [CUDA (single GPU)](http://dx.doi.org/10.1016/j.jcp.2014.04.035)
-- [Charge Method Optimizations with OpenMP](https://doi.org/10.1109/ScalA.2016.006)
+- [Charge Model Optimizations with OpenMP](https://doi.org/10.1137/18M1224684)
 
-Distributed Memory:
+Distributed-Memory Versions:
 - [MPI (message passing interface)](https://www.cs.purdue.edu/puremd/docs/Parallel-Reactive-Molecular-Dynamics.pdf)
+- [Hybrid MPI+OpenMP optimization](https://doi.org/10.1177/1094342017746221)
+- [Charge Model Optimizations with MPI](https://doi.org/10.1145/3330345.3330359)
 - [CUDA+MPI (multi-GPU)](https://www.cs.purdue.edu/puremd/docs/pgpuremd.pdf)
diff --git a/sPuReMD/src/bonds.c b/sPuReMD/src/bonds.c
index 2e61e2617b30b007518d376c50808a69e7ac73dd..79546a42459f793b8d28c61c67f882f423456b8c 100644
--- a/sPuReMD/src/bonds.c
+++ b/sPuReMD/src/bonds.c
@@ -30,7 +30,7 @@ void Bonds( reax_system *system, control_params *control,
         reax_list **lists, output_controls *out_control )
 {
     int i;
-    real gp3, gp4, gp7, gp10, gp37, ebond_total;
+    real gp3, gp4, gp7, gp10, ebond_total;
     reax_list *bonds;
 
     bonds = lists[BONDS];
@@ -38,7 +38,6 @@ void Bonds( reax_system *system, control_params *control,
     gp4 = system->reax_param.gp.l[4];
     gp7 = system->reax_param.gp.l[7];
     gp10 = system->reax_param.gp.l[10];
-    gp37 = (int) system->reax_param.gp.l[37];
     ebond_total = 0.0;
 
 #if defined(_OPENMP)
@@ -119,11 +118,10 @@ void Bonds( reax_system *system, control_params *control,
                     /* Stabilisation terminal triple bond in C-O */
                     if ( bo_ij->BO >= 1.00 )
                     {
-                        if ( gp37 == 2 &&
-                                ( (strncmp( sbp_i->name, "C", sizeof(sbp_i->name) ) == 0
+                        if ( (strncmp( sbp_i->name, "C", sizeof(sbp_i->name) ) == 0
                                     && strncmp( sbp_j->name, "O", sizeof(sbp_j->name) ) == 0)
                                 || (strncmp( sbp_i->name, "O", sizeof(sbp_i->name) ) == 0
-                                    && strncmp( sbp_j->name, "C", sizeof(sbp_j->name) ) == 0) ) )
+                                    && strncmp( sbp_j->name, "C", sizeof(sbp_j->name) ) == 0) )
                         {
                             //ba = SQR( bo_ij->BO - 2.5 );
                             exphu = EXP( -gp7 * SQR(bo_ij->BO - 2.5) );
diff --git a/sPuReMD/src/reax_types.h b/sPuReMD/src/reax_types.h
index 5a57b49a36ffff1f3a2d5ecac58439717199a311..87cae33eacc6bd83ce1e53300b0c48bdc4aa77dd 100644
--- a/sPuReMD/src/reax_types.h
+++ b/sPuReMD/src/reax_types.h
@@ -503,7 +503,7 @@ typedef int (*write_function)( FILE *, const char *, ... );
  * l[34] = b_s_acks2 (ACKS2 bond softness)
  * l[35] = N/A
  * l[36] = N/A
- * l[37] = version number
+ * l[37] = N/A
  * l[38] = p_coa3 */
 struct global_parameters
 {