Precompute terms in overlap integral

src/REAXFF/fix_qeq_rel_reaxff.cpp

@@ -50,12 +50,12 @@ void FixQEqRelReaxFF::calc_chi_eff()
   const auto x = (const double *const *) atom->x;
   const int *type = atom->type;
 
-  double dx, dy, dz, dist_sq, overlap, sum_n, sum_d, expa, expb, chia, phia, phib, p, m;
+  double dx, dy, dz, dist_sq, overlap, sum_n, sum_d, chia, phia, phib;
   int i, j;
 
   // check ghost atoms are stored up to the distance cutoff for overlap integrals
   const double comm_cutoff = MAX(neighbor->cutneighmax, comm->cutghostuser);
-  if(comm_cutoff*comm_cutoff < dist_cutoff_sq) {
+  if (comm_cutoff*comm_cutoff < dist_cutoff_sq) {
     error->all(FLERR, Error::NOLASTLINE,
                "Comm cutoff {} is smaller than distance cutoff {} for overlap integrals in fix {}. "
                "Increase accordingly using comm_modify cutoff",
@@ -71,7 +71,6 @@ void FixQEqRelReaxFF::calc_chi_eff()
 
     // compute chi_eff for each local atom
     for (i = 0; i < nn; i++) {
-      expa = gauss_exp[type[i]];
       chia = chi[type[i]];
       if (efield->varflag != FixEfield::ATOM) {
         phia = -factor * (x[i][0] * efield->ex + x[i][1] * efield->ey + x[i][2] * efield->ez);
@@ -86,15 +85,12 @@ void FixQEqRelReaxFF::calc_chi_eff()
           dx = x[i][0] - x[j][0];
           dy = x[i][1] - x[j][1];
           dz = x[i][2] - x[j][2];
-	  dist_sq = (dx*dx + dy*dy + dz*dz);
+          dist_sq = (dx*dx + dy*dy + dz*dz);
 
         if (dist_sq < dist_cutoff_sq) {
-          expb = gauss_exp[type[j]];
 
           // overlap integral of two normalised 1s Gaussian type orbitals
-          p = expa + expb;
-          m = expa * expb / p;
-          overlap = pow((4.0 * m / p), 0.75) * exp(-m * dist_sq);
+          overlap = prefactor[type[i]][type[j]] * exp(-expfactor[type[i]][type[j]] * dist_sq);
 
           if (efield->varflag != FixEfield::ATOM) {
             phib = -factor * (x[j][0] * efield->ex + x[j][1] * efield->ey + x[j][2] * efield->ez);

src/REAXFF/fix_qtpie_reaxff.cpp

@@ -111,6 +111,8 @@ FixQtpieReaxFF::FixQtpieReaxFF(LAMMPS *lmp, int narg, char **arg) :
     iarg++;
   }
   shld = nullptr;
+  prefactor = nullptr;
+  expfactor = nullptr;
 
   nn = nt = n_cap = 0;
   nmax = 0;
@@ -172,8 +174,10 @@ FixQtpieReaxFF::~FixQtpieReaxFF()
   FixQtpieReaxFF::deallocate_matrix();
 
   memory->destroy(shld);
-
   memory->destroy(gauss_exp);
+  memory->destroy(prefactor);
+  memory->destroy(expfactor);
+
   if (!reaxflag) {
     memory->destroy(chi);
     memory->destroy(eta);
@@ -502,6 +506,7 @@ void FixQtpieReaxFF::init()
 
   init_shielding();
   init_taper();
+  init_olap();
 
   if (utils::strmatch(update->integrate_style,"^respa"))
     nlevels_respa = (dynamic_cast<Respa *>(update->integrate))->nlevels;
@@ -569,6 +574,31 @@ void FixQtpieReaxFF::init_taper()
 
 /* ---------------------------------------------------------------------- */
 
+void FixQtpieReaxFF::init_olap()
+{
+  int i,j;
+  int ntypes;
+  double expa,expb,expsum,expnorm;
+
+  ntypes = atom->ntypes;
+  if (prefactor == nullptr)
+    memory->create(prefactor,ntypes+1,ntypes+1,"qtpie:overlap_prefactor");
+  if (expfactor == nullptr)
+    memory->create(expfactor,ntypes+1,ntypes+1,"qtpie:overlap_expfactor");
+
+  for (i = 1; i <= ntypes; ++i)
+    for (j = 1; j <= ntypes; ++j) {
+      expa = gauss_exp[i];
+      expb = gauss_exp[j];
+      expsum = expa + expb;
+      expnorm = expa * expb / expsum;
+      prefactor[i][j] = pow((4.0 * expnorm / expsum), 0.75);
+      expfactor[i][j] = expnorm;
+    }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixQtpieReaxFF::setup_pre_force(int vflag)
 {
   if (reaxff) {
@@ -1131,7 +1161,7 @@ void FixQtpieReaxFF::calc_chi_eff()
   const auto x = (const double * const *)atom->x;
   const int *type = atom->type;
 
-  double dx,dy,dz,dist_sq,overlap,sum_n,sum_d,expa,expb,chia,chib,phia,phib,p,m;
+  double dx,dy,dz,dist_sq,overlap,sum_n,sum_d,chia,chib,phia,phib;
   int i,j;
 
   // check ghost atoms are stored up to the distance cutoff for overlap integrals
@@ -1154,7 +1184,6 @@ void FixQtpieReaxFF::calc_chi_eff()
 
   // compute chi_eff for each local atom
   for (i = 0; i < nn; i++) {
-    expa = gauss_exp[type[i]];
     chia = chi[type[i]];
     if (efield) {
       if (efield->varflag != FixEfield::ATOM) {
@@ -1174,13 +1203,10 @@ void FixQtpieReaxFF::calc_chi_eff()
       dist_sq = (dx*dx + dy*dy + dz*dz);
 
       if (dist_sq < dist_cutoff_sq) {
-        expb = gauss_exp[type[j]];
         chib = chi[type[j]];
 
         // overlap integral of two normalised 1s Gaussian type orbitals
-        p = expa + expb;
-        m = expa * expb / p;
-        overlap = pow((4.0 * m / p), 0.75) * exp(-m * dist_sq);
+        overlap = prefactor[type[i]][type[j]] * exp(-expfactor[type[i]][type[j]] * dist_sq);
 
         if (efield) {
           if (efield->varflag != FixEfield::ATOM) {

src/REAXFF/fix_qtpie_reaxff.h

@@ -90,12 +90,15 @@ class FixQtpieReaxFF : public Fix {
   char *pertype_option;    // argument to determine how per-type info is obtained
   char *gauss_file;        // input file for gaussian orbital exponents
   double *gauss_exp;       // array of gaussian orbital exponents for each atom type
+  double **prefactor;      // factor used in computation of overlap integrals
+  double **expfactor;      // factor used in exponential term of overlap integrals
   double dist_cutoff_sq;   // separation distance squared beyond which overlap integrals are neglected
   double scale;            // scaling factor for electric polarization effects
 
   void pertype_parameters(char *);
   void init_shielding();
   void init_taper();
+  void init_olap();
   void allocate_storage();
   void deallocate_storage();
   void reallocate_storage();