Commit 61b1487c authored by Axel Kohlmeyer's avatar Axel Kohlmeyer
Browse files

avoid division by zero in reaxff bond interaction computations in very rare cases 

this addresses the issue reported by stan and ishan
parent e53583d9

src/KOKKOS/pair_reaxc_kokkos.cpp

+2 −1
Original line number Diff line number Diff line
@@ -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

src/USER-REAXC/reaxc_bonds.cpp

+2 −1
Original line number Diff line number Diff line
@@ -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 );

CRA Git | Maintained and supported by SUSTech CRA and CCSE