diff --git a/src/USER-REAXC/reaxc_allocate.cpp b/src/USER-REAXC/reaxc_allocate.cpp
index 1851ac5bd9..de2f117bbc 100644
--- a/src/USER-REAXC/reaxc_allocate.cpp
+++ b/src/USER-REAXC/reaxc_allocate.cpp
@@ -176,6 +176,7 @@ void DeAllocate_Workspace( control_params *control, storage *workspace )
   sfree( workspace->dDelta_lp_temp, "dDelta_lp_temp" );
   sfree( workspace->Delta_e, "Delta_e" );
   sfree( workspace->Delta_boc, "Delta_boc" );
+  sfree( workspace->Delta_val, "Delta_val" );
   sfree( workspace->nlp, "nlp" );
   sfree( workspace->nlp_temp, "nlp_temp" );
   sfree( workspace->Clp, "Clp" );
@@ -311,6 +312,7 @@ int Allocate_Workspace( reax_system *system, control_params *control,
     smalloc( total_real, "dDelta_lp_temp", comm );
   workspace->Delta_e = (real*) smalloc( total_real, "Delta_e", comm );
   workspace->Delta_boc = (real*) smalloc( total_real, "Delta_boc", comm );
+  workspace->Delta_val = (real*) smalloc( total_real, "Delta_val", comm );
   workspace->nlp = (real*) smalloc( total_real, "nlp", comm );
   workspace->nlp_temp = (real*) smalloc( total_real, "nlp_temp", comm );
   workspace->Clp = (real*) smalloc( total_real, "Clp", comm );
diff --git a/src/USER-REAXC/reaxc_bond_orders.cpp b/src/USER-REAXC/reaxc_bond_orders.cpp
index 02cee90d16..26c1648a77 100644
--- a/src/USER-REAXC/reaxc_bond_orders.cpp
+++ b/src/USER-REAXC/reaxc_bond_orders.cpp
@@ -1137,6 +1137,8 @@ void BO( reax_system *system, control_params *control, simulation_data *data,
     workspace->Delta_e[j] = workspace->total_bond_order[j] - sbp_j->valency_e;
     workspace->Delta_boc[j] = workspace->total_bond_order[j] -
       sbp_j->valency_boc;
+    workspace->Delta_val[j] = workspace->total_bond_order[j] -
+      sbp_j->valency_val;
 
     workspace->vlpex[j] = workspace->Delta_e[j] - 
       2.0 * (int)(workspace->Delta_e[j]/2.0);    
diff --git a/src/USER-REAXC/reaxc_nonbonded.cpp b/src/USER-REAXC/reaxc_nonbonded.cpp
index 419779a68f..e03c9ef102 100644
--- a/src/USER-REAXC/reaxc_nonbonded.cpp
+++ b/src/USER-REAXC/reaxc_nonbonded.cpp
@@ -116,9 +116,9 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
 	
 	e_vdW = twbp->D * (exp1 - 2.0 * exp2);
 	data->my_en.e_vdW += Tap * e_vdW;
-	
+
 	CEvd = dTap * e_vdW - 
-	  Tap * twbp->D * (twbp->alpha / twbp->r_vdW) * (exp1 - exp2);
+	  Tap * twbp->D * (twbp->alpha / twbp->r_vdW) * (exp1 - exp2) / r_ij;
       }
 
       if(system->reax_param.gp.vdw_type==2 || system->reax_param.gp.vdw_type==3)
@@ -127,7 +127,7 @@ void vdW_Coulomb_Energy( reax_system *system, control_params *control,
 	  data->my_en.e_vdW += Tap * e_core;
 
 	  de_core = -(twbp->acore/twbp->rcore) * e_core;
-	  CEvd += dTap * e_core + Tap * de_core;
+	  CEvd += dTap * e_core + Tap * de_core / r_ij;
 	}
 
       /*Coulomb Calculations*/
diff --git a/src/USER-REAXC/reaxc_types.h b/src/USER-REAXC/reaxc_types.h
index 7c4dfc3c5b..bf54758d92 100644
--- a/src/USER-REAXC/reaxc_types.h
+++ b/src/USER-REAXC/reaxc_types.h
@@ -748,7 +748,7 @@ typedef struct
   /* bond order related storage */
   real *total_bond_order;
   real *Deltap, *Deltap_boc;
-  real *Delta, *Delta_lp, *Delta_lp_temp, *Delta_e, *Delta_boc;
+  real *Delta, *Delta_lp, *Delta_lp_temp, *Delta_e, *Delta_boc, *Delta_val;
   real *dDelta_lp, *dDelta_lp_temp;
   real *nlp, *nlp_temp, *Clp, *vlpex;
   rvec *dDeltap_self;
diff --git a/src/USER-REAXC/reaxc_valence_angles.cpp b/src/USER-REAXC/reaxc_valence_angles.cpp
index 84b41ba5fc..8ffd103a9d 100644
--- a/src/USER-REAXC/reaxc_valence_angles.cpp
+++ b/src/USER-REAXC/reaxc_valence_angles.cpp
@@ -235,7 +235,7 @@ void Valence_Angles( reax_system *system, control_params *control,
 	  
 	  
 	  if( (j < system->n) && (BOA_jk > 0.0) && 
-	      (bo_ij->BO * bo_jk->BO > SQR(control->thb_cut)/*0*/) ) {
+              (bo_ij->BO * bo_jk->BO > SQR(control->thb_cut)/*0*/) ) {
 	    r_jk = pbond_jk->d;		      
 	    thbh = &( system->reax_param.thbp[ type_i ][ type_j ][ type_k ] );
 	    
@@ -343,7 +343,7 @@ void Valence_Angles( reax_system *system, control_params *control,
 		p_coa3 = system->reax_param.gp.l[38];
 		p_coa4 = system->reax_param.gp.l[30];
 		
-		exp_coa2 = exp( p_coa2 * workspace->Delta_boc[j] );
+		exp_coa2 = exp( p_coa2 * workspace->Delta_val[j] );
 		data->my_en.e_coa += e_coa = 
 		  p_coa1 / (1. + exp_coa2) *
 		  exp( -p_coa3 * SQR(workspace->total_bond_order[i]-BOA_ij) ) *