diff --git a/src/KOKKOS/pair_reaxc_kokkos.cpp b/src/KOKKOS/pair_reaxc_kokkos.cpp
index 6be09548da..841b7fbea9 100644
--- a/src/KOKKOS/pair_reaxc_kokkos.cpp
+++ b/src/KOKKOS/pair_reaxc_kokkos.cpp
@@ -3335,7 +3335,8 @@ void PairReaxCKokkos<DeviceType>::operator()(PairReaxComputeBond1<NEIGHFLAG,EVFL
     const F_FLOAT BO_pi_i = d_BO_pi(i,j_index);
     const F_FLOAT BO_pi2_i = d_BO_pi2(i,j_index);
 
-    pow_BOs_be2 = pow(BO_s_i,p_be2);
+    if (BO_s_i == 0.0) pow_BOs_be2 = 0.0;
+    else pow_BOs_be2 = pow(BO_s_i,p_be2);
     exp_be12 = exp(p_be1*(1.0-pow_BOs_be2));
     CEbo = -De_s*exp_be12*(1.0-p_be1*p_be2*pow_BOs_be2);
     ebond = -De_s*BO_s_i*exp_be12
diff --git a/src/USER-REAXC/reaxc_bonds.cpp b/src/USER-REAXC/reaxc_bonds.cpp
index a8a1298166..cf536fe606 100644
--- a/src/USER-REAXC/reaxc_bonds.cpp
+++ b/src/USER-REAXC/reaxc_bonds.cpp
@@ -82,7 +82,8 @@ void Bonds( reax_system *system, control_params *control,
       bo_ij = &( bonds->select.bond_list[pj].bo_data );
 
       /* calculate the constants */
-      pow_BOs_be2 = pow( bo_ij->BO_s, twbp->p_be2 );
+      if (bo_ij->BO_s == 0.0) pow_BOs_be2 = 0.0;
+      else pow_BOs_be2 = pow( bo_ij->BO_s, twbp->p_be2 );
       exp_be12 = exp( twbp->p_be1 * ( 1.0 - pow_BOs_be2 ) );
       CEbo = -twbp->De_s * exp_be12 *
 	( 1.0 - twbp->p_be1 * twbp->p_be2 * pow_BOs_be2 );