From 61b1487cbd9364423129de4b506edd413c2524bd Mon Sep 17 00:00:00 2001
From: Axel Kohlmeyer <akohlmey@gmail.com>
Date: Thu, 20 Jul 2017 18:17:19 -0400
Subject: [PATCH] avoid division by zero in reaxff bond interaction
computations in very rare cases
this addresses the issue reported by stan and ishan
---
src/KOKKOS/pair_reaxc_kokkos.cpp | 3 ++-
src/USER-REAXC/reaxc_bonds.cpp | 3 ++-
2 files changed, 4 insertions(+), 2 deletions(-)
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 );
--
GitLab