sPuReMD: add ACKS2 geometry-dependent forces to tabulated computation of van...

sPuReMD: add ACKS2 geometry-dependent forces to tabulated computation of van der Waals and Coulomb terms. Replace double precision equality check with difference comparison.

sPuReMD/src/two_body_interactions.c

@@ -113,8 +113,8 @@ void Bond_Energy( reax_system *system, control_params *control,
                     if ( bo_ij->BO >= 1.00 )
                     {
                         if ( gp37 == 2
-                                || (sbp_i->mass == 12.0000 && sbp_j->mass == 15.9990)
-                                || (sbp_j->mass == 12.0000 && sbp_i->mass == 15.9990) )
+                                || (FABS(sbp_i->mass - 12.0) < 0.0001 && FABS(sbp_j->mass - 15.9990) < 0.0001)
+                                || (FABS(sbp_j->mass - 12.0) < 0.0001 && FABS(sbp_i->mass - 15.9990) < 0.0001) )
                         {
                             //ba = SQR( bo_ij->BO - 2.5 );
                             exphu = EXP( -gp7 * SQR(bo_ij->BO - 2.5) );
@@ -617,6 +617,7 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
         int start_i, end_i;
         real r_ij, self_coef, base, dif;
         real e_vdW, e_ele;
+        real d, xcut, bond_softness, d_bond_softness, effpot_diff;
         real CEvd, CEclmb;
         rvec temp, ext_press;
         far_neighbor_data *nbr_pj;
@@ -737,6 +738,79 @@ void Tabulated_vdW_Coulomb_Energy( reax_system *system, control_params *control,
                 }
             }
         }
+
+        //TODO: better integration AND tabulation of ACKS2 code below with above code for performance
+#ifdef _OPENMP
+        #pragma omp barrier
+#endif
+
+        /* contribution to energy and gradients (atoms and cell)
+         * due to geometry-dependent terms in the ACKS2
+         * kinetic energy */
+        if ( control->charge_method == ACKS2_CM )
+        {
+#ifdef _OPENMP
+            #pragma omp for schedule(guided)
+#endif
+            for ( i = 0; i < system->N; ++i )
+            {
+                for ( pj = Start_Index(i, far_nbrs); pj < End_Index(i, far_nbrs); ++pj )
+                {
+                    nbr_pj = &far_nbrs->far_nbr_list[pj];
+                    j = nbr_pj->nbr;
+
+                    /* kinetic energy terms */
+                    xcut = 0.5 * ( system->reaxprm.sbp[ system->atoms[i].type ].b_s_acks2
+                            + system->reaxprm.sbp[ system->atoms[j].type ].b_s_acks2 );
+
+                    if ( far_nbrs->far_nbr_list[pj].d < xcut )
+                    {
+                        d = far_nbrs->far_nbr_list[pj].d / xcut;
+                        bond_softness = system->reaxprm.gp.l[34] * POW( d, 3.0 )
+                            * POW( 1.0 - d, 6.0 );
+
+                        if ( bond_softness > 0.0 )
+                        {
+                            /* Coulombic energy contribution */
+                            effpot_diff = workspace->s[0][system->N + i]
+                                - workspace->s[0][system->N + j];
+                            e_ele = -0.5 * KCALpMOL_to_EV * bond_softness
+                                * SQR( effpot_diff );
+                            e_ele_total += e_ele;
+
+                            /* forces contribution */
+                            d_bond_softness = system->reaxprm.gp.l[34]
+                                * 3.0 / xcut * POW( d, 2.0 )
+                                * POW( 1.0 - d, 5.0 ) * (1.0 - 3.0 * d);
+                            d_bond_softness = -0.5 * d_bond_softness
+                                * SQR( effpot_diff );
+                            d_bond_softness = KCALpMOL_to_EV * d_bond_softness
+                                / far_nbrs->far_nbr_list[pj].d;
+
+#if defined(DEBUG_FOCUS)
+                            fprintf( stderr, "%6d%6d%12.5f%12.5f%12.5f\n",
+                                    j + 1, i + 1,
+                                    d_bond_softness * nbr_pj->dvec[0],
+                                    d_bond_softness * nbr_pj->dvec[1],
+                                    d_bond_softness * nbr_pj->dvec[2] ); fflush( stderr );
+#endif
+
+#ifndef _OPENMP
+                            rvec_ScaledAdd( system->atoms[i].f,
+                                    -d_bond_softness, nbr_pj->dvec );
+                            rvec_ScaledAdd( system->atoms[j].f,
+                                    d_bond_softness, nbr_pj->dvec );
+#else
+                            rvec_ScaledAdd( workspace->f_local[tid * system->N + i],
+                                    -d_bond_softness, nbr_pj->dvec );
+                            rvec_ScaledAdd( workspace->f_local[tid * system->N + j],
+                                    d_bond_softness, nbr_pj->dvec );
+#endif
+                        }
+                    }
+                }
+            }
+        }
     }
 
     data->E_vdW += e_vdW_total;