Use cubic spline interpolation to speed up volmap

doc/cpptraj.lyx

@@ -35435,7 +35435,7 @@ volmap [out <filename>] <mask> [radscale <factor>] [stepfac <fac>]
 \end_layout
 
 \begin_layout LyX-Code
-   	[sphere] [radii {vdw | element}]
+   	[sphere] [radii {vdw | element}] [splinedx <spacing>]
 \end_layout
 
 \begin_layout LyX-Code
@@ -35528,6 +35528,14 @@ radii
  radii.
 \end_layout
 
+\begin_layout Description
+splinedx
+\begin_inset space ~
+\end_inset
+
+<spacing> Spacing to use for cubic spline interpolation (default 0.01 Ang.).
+\end_layout
+
 \begin_layout Description
 calcpeaks If specified, peaks in the grid density will be calculated and
  saved to set <setname> with aspect 
@@ -35737,6 +35745,25 @@ buffer
 .
 \end_layout
 
+\begin_layout Standard
+The calculation is sped up by using cubic splines to interpolate the exponential
+ function when calculating the Gaussians.
+ The details are in Roe & Brooks, 
+\begin_inset Quotes eld
+\end_inset
+
+Improving the Speed of Volumetric Density Map Generation via Cubic Spline
+ Interpolation
+\begin_inset Quotes erd
+\end_inset
+
+, J.
+ Mol.
+ Graph.
+ Model.
+ (2021).
+\end_layout
+
 \begin_layout Subsection
 volume
 \end_layout

src/Action_Volmap.cpp

@@ -15,6 +15,11 @@
 #ifdef _OPENMP
 # include <omp.h>
 #endif
+//#if defined(VOLMAP_USEFASTEXPS) || defined(VOLMAP_USEFASTEXP64)
+//# incl ude "FastExp_Schraudolph.h"
+//#elif defined(VOLMAP_USEFASTEXPIEEE)
+//# incl ude "FastExp_fastexp.h"
+//#endif
 
 const double Action_Volmap::sqrt_8_pi_cubed_ = sqrt(8.0*Constants::PI*Constants::PI*Constants::PI);
 
@@ -39,7 +44,8 @@ Action_Volmap::Action_Volmap() :
   peakcut_(0.05),
   buffer_(3.0),
   radscale_(1.0),
-  stepfac_(4.1)
+  stepfac_(4.1),
+  splineDx_(0.01) // Recommendation from Roe & Brooks JMGM 2021
 {}
 
 void Action_Volmap::Help() const {
@@ -55,7 +61,7 @@ void Action_Volmap::Help() const {
 
 void Action_Volmap::RawHelp() const {
   mprintf("\t[out <filename>] <mask> [radscale <factor>] [stepfac <fac>]\n"
-          "\t[sphere] [radii {vdw | element}]\n"
+          "\t[sphere] [radii {vdw | element}] [splinedx <spacing>]\n"
           "\t[calcpeaks] [peakcut <cutoff>] [peakfile <xyzfile>]\n"
           "\t{ data <existing set> |\n"
           "\t  name <setname> <dx> [<dy> <dz>]\n"
@@ -82,6 +88,8 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d
     radscale_ = 0.5;
     stepfac_ = 1.0;
   }
+  //splineDx_ = actionArgs.getKeyDouble("splinedx", 1.0/5000.0);
+  splineDx_ = actionArgs.getKeyDouble("splinedx", 0.01);
   radscale_ = 1.0 / actionArgs.getKeyDouble("radscale", radscale_);
   stepfac_ = actionArgs.getKeyDouble("stepfac", stepfac_);
   std::string radarg = actionArgs.GetStringKey("radii");
@@ -214,7 +222,11 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d
     xmin_ = oxyz[0];
     ymin_ = oxyz[1];
     zmin_ = oxyz[2];
-  } 
+  }
+  // Warn for small grid spacing
+  if (dx_ < 0.4 || dy_ < 0.4 || dz_ < 0.4)
+    mprintf("Warning: Grid spacings smaller than 0.4 Ang. may be very slow;\n"
+            "Warning:  consider using larger grid spacings.\n");
   //std::string density = actionArgs.GetStringKey("density"); // FIXME obsolete?
   // Get the required mask
   std::string reqmask = actionArgs.GetMaskNext();
@@ -295,7 +307,26 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d
     mprintf("\tWhen smearing Gaussian, voxels farther than radii/2 will be skipped.\n");
   mprintf("\tDividing radii by %f\n", 1.0/radscale_);
   mprintf("\tFactor for determining number of bins to smear Gaussian is %f\n", stepfac_);
-
+# if defined(VOLMAP_USEEXP)
+  mprintf("\tUsing system exp() function for evaluating Gaussians.\n");
+//# elif defined(VOLMAP_USEFASTEXPS)
+//  mprintf("\tUsing exp() from N. Schraudolph, Neural Computation 11, 853–862 (1999).\n");
+//# elif defined(VOLMAP_USEFASTEXP64)
+//  mprintf("\tUsing 64 bit version of exp() from N. Schraudolph, Neural Computation 11, 853–862 (1999).\n");
+//# elif defined(VOLMAP_USEFASTEXPIEEE)
+//  mprintf("\tUsing exp() from Schraudolph & Malossi et al.\n");
+# else /* VOLMAP_USEEXP */
+  mprintf("\tExponential for Gaussians will be approximated using cubic splines\n"
+          "\t  with a spacing of %g Ang.\n", splineDx_);
+  mprintf("# Citation: Roe, D. R.; Brooks, B. R.; \"Improving the Speed of Volumetric\n"
+          "#           Density Map Generation via Cubic Spline Interpolation\".\n"
+          "#           Journal of Molecular Graphics and Modelling (2021).\n");
+# if defined(VOLMAP_USEACCURATE)
+  mprintf("\tSplines using more accurate but slower table lookup.\n");
+# elif defined(VOLMAP_USEXTABLE)
+  mprintf("\tSplines using less accurate lookup with tabled X values.\n");
+# endif
+# endif /* VOLMAP_USEEXP */
   if (outfile != 0)
     mprintf("\tDensity will wrtten to '%s'\n", outfile->DataFilename().full());
   mprintf("\tGrid dataset name is '%s'\n", grid_->legend());
@@ -381,8 +412,14 @@ Action::RetType Action_Volmap::Setup(ActionSetup& setup) {
   //mprintf("DEBUG: nx= %i  ny= %i  nz= %i\n", nxstep, nystep, nzstep);
   //mprintf("DEBUG: %g %g %g %g\n", maxx, maxy, maxz, maxDist);
   maxDist *= (-1.0 / (2.0 * maxRad * maxRad));
-  if (debug_ > 1)
-    mprintf("DEBUG: maxDist= %g\n", maxDist);
+  //mprintf("DEBUG: max= %g\n", maxDist);
+# ifndef VOLMAP_USEEXP
+  // Set up the interpolation table
+  if (table_.FillTable( exp, splineDx_, maxDist, 1.0 )) return Action::ERR;
+  table_.PrintMemUsage("\t");
+  if (debug_ > 0)
+    table_.PrintTableInfo("DEBUG: ");
+# endif
 
   if ((int)Atoms_.size() < densitymask_.Nselected())
     mprintf("Warning: %i atoms have 0.0 radii and will be skipped.\n",
@@ -511,9 +548,39 @@ Action::RetType Action_Volmap::DoAction(int frameNum, ActionFrame& frm) {
                 if (dist2 < rcut2) {
                   //mprintf("DEBUG: rhalf= %g  dist2= %g  exfac= %g  exp= %g\n", rhalf, dist2, exfac, exfac*dist2);
 #                 ifdef _OPENMP
+                  // NOTE: It is OK to call table_.Yval() here because in OpenMP
+                  //       local variables of called functions are private.
+#                 if defined(VOLMAP_USEEXP)
                   GRID_THREAD_[mythread].incrementBy(xval, yval, zval, norm * exp(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXPS)
+//                  GRID_THREAD[mythread].incrementBy(xval, yval, zval, norm * FASTEXPS(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXP64)
+//                  GRID_THREAD[mythread].incrementBy(xval, yval, zval, norm * fast_exps_64(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXPIEEE)
+//                  GRID_THREAD[mythread].incrementBy(xval, yval, zval, norm * fastexp::exp<double, fastexp::IEEE, 4>(exfac * dist2));
+#                 elif defined(VOLMAP_USEACCURATE)
+                  GRID_THREAD_[mythread].incrementBy(xval, yval, zval, norm * table_.Yval_accurate(exfac * dist2));
+#                 elif defined(VOLMAP_USEXTABLE)
+                  GRID_THREAD_[mythread].incrementBy(xval, yval, xval, norm * table_.Yval_xtable(exfac * dist2));
+#                 else
+                  GRID_THREAD_[mythread].incrementBy(xval, yval, zval, norm * table_.Yval(exfac * dist2));
+#                 endif
 #                 else /* OPENMP */
+#                 if defined(VOLMAP_USEEXP)
                   grid_->Increment(xval, yval, zval, norm * exp(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXPS)
+//                  grid_->Increment(xval, yval, zval, norm * FASTEXPS(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXP64)
+//                  grid_->Increment(xval, yval, zval, norm * fast_exps_64(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXPIEEE)
+//                  grid_->Increment(xval, yval, zval, norm * fastexp::exp<double, fastexp::IEEE, 4>(exfac * dist2));
+#                 elif defined(VOLMAP_USEACCURATE)
+                  grid_->Increment(xval, yval, zval, norm * table_.Yval_accurate(exfac * dist2));
+#                 elif defined(VOLMAP_USEXTABLE)
+                  grid_->Increment(xval, yval, zval, norm * table_.Yval_xtable(exfac * dist2));
+#                 else
+                  grid_->Increment(xval, yval, zval, norm * table_.Yval(exfac * dist2));
+#                 endif
 #                 endif /* OPENMP */
                 }
               } // END loop over zval
@@ -531,6 +598,7 @@ Action::RetType Action_Volmap::DoAction(int frameNum, ActionFrame& frm) {
 }
 
 #ifdef MPI
+/** Sync grid to the master process. */
 int Action_Volmap::SyncAction() {
 # ifdef _OPENMP
   CombineGridThreads();

src/Action_Volmap.h

@@ -2,6 +2,7 @@
 #define INC_ACTION_VOLMAP_H
 #include "Action.h"
 #include "Grid.h"
+#include "SplineFxnTable.h"
 #ifdef VOLMAP_DOUBLE
 # define VOLMAP_DS_T DataSet_GridDbl
 # define VOLMAP_T double
@@ -10,6 +11,14 @@
 # define VOLMAP_T float
 #endif
 class VOLMAP_DS_T;
+/// Calculate atomic volumetric density maps from trajectory data.
+/** By default the grid type used is single-precision, mostly to save space.
+  * A double-precision grid can be used by compiling with the
+  * VOLMAP_DOUBLE define.
+  * Also by default the exp() function will be approximated with cubic spline
+  * interpolation. To use the system exp() function, compile with the
+  * VOLMAP_USEEXP define.
+  */
 class Action_Volmap : public Action {
   public:
     Action_Volmap();
@@ -51,6 +60,8 @@ class Action_Volmap : public Action {
     double radscale_;       ///< The scaling factor to divide all radii by
     double stepfac_;        ///< Factor for determining how many steps to smear Gaussian
     static const double sqrt_8_pi_cubed_;
+    SplineFxnTable table_;
+    double splineDx_;
 #   ifdef _OPENMP
     typedef std::vector< Grid<VOLMAP_T> > Garray;
     Garray GRID_THREAD_;

src/Ewald.cpp

@@ -13,18 +13,17 @@
 
 /// CONSTRUCTOR
 Ewald::Ewald() :
-  sumq_(0.0),
-  sumq2_(0.0),
   ew_coeff_(0.0),
   lw_coeff_(0.0),
   switch_width_(0.0),
   cutoff_(0.0),
   cut2_(0.0),
   cut2_0_(0.0),
   dsumTol_(0.0),
-  erfcTableDx_(0.0),
-  one_over_Dx_(0.0),
-  debug_(0)
+  debug_(0),
+  sumq_(0.0),
+  sumq2_(0.0),
+  Vdw_Recip_term_(0)
 {
 # ifdef DEBUG_EWALD
   // Save fractional translations for 1 cell in each direction (and primary cell).
@@ -91,44 +90,9 @@ double Ewald::erfc_func(double xIn) {
   return erfc;
 }
 
-// Ewald::FillErfcTable()
-void Ewald::FillErfcTable(double cutoffIn, double dxdr) {
-  one_over_Dx_ = 1.0 / erfcTableDx_;
-  unsigned int erfcTableSize = (unsigned int)(dxdr * one_over_Dx_ * cutoffIn * 1.5);
-  Darray erfc_X, erfc_Y;
-  erfc_X.reserve( erfcTableSize );
-  erfc_Y.reserve( erfcTableSize );
-  // Save X and Y values so we can calc the spline coefficients
-  double xval = 0.0;
-  for (unsigned int i = 0; i != erfcTableSize; i++) {
-    double yval = erfc_func( xval );
-    erfc_X.push_back( xval );
-    erfc_Y.push_back( yval );
-    xval += erfcTableDx_;
-  }
-  Spline cspline;
-  cspline.CubicSpline_Coeff(erfc_X, erfc_Y);
-  erfc_X.clear();
-  // Store values in Spline table
-  erfc_table_.reserve( erfcTableSize * 4 ); // Y B C D
-  for (unsigned int i = 0; i != erfcTableSize; i++) {
-    erfc_table_.push_back( erfc_Y[i] );
-    erfc_table_.push_back( cspline.B_coeff()[i] );
-    erfc_table_.push_back( cspline.C_coeff()[i] );
-    erfc_table_.push_back( cspline.D_coeff()[i] );
-  }
-  // Memory saved Y values plus spline B, C, and D coefficient arrays.
-  mprintf("\tMemory used by Erfc table and splines: %s\n",
-          ByteString(erfc_table_.size() * sizeof(double), BYTE_DECIMAL).c_str());
-}
-
 // Ewald::ERFC()
 double Ewald::ERFC(double xIn) const {
-  int xidx = ((int)(one_over_Dx_ * xIn));
-  double dx = xIn - ((double)xidx * erfcTableDx_);
-  xidx *= 4;
-  return erfc_table_[xidx] + 
-         dx*(erfc_table_[xidx+1] + dx*(erfc_table_[xidx+2] + dx*erfc_table_[xidx+3]));
+  return table_.Yval( xIn);
 }
 
 /** Determine Ewald coefficient from cutoff and direct sum tolerance.
@@ -200,7 +164,7 @@ void Ewald::CalculateC6params(Topology const& topIn, AtomMask const& maskIn) {
 }
 
 /** Set up exclusion lists for selected atoms. */
-void Ewald::SetupExcluded(Topology const& topIn, AtomMask const& maskIn)
+void Ewald::SetupExclusionList(Topology const& topIn, AtomMask const& maskIn)
 {
   // Use distance of 4 (up to dihedrals)
   if (Excluded_.SetupExcluded(topIn.Atoms(), maskIn, 4,
@@ -223,7 +187,7 @@ int Ewald::CheckInput(Box const& boxIn, int debugIn, double cutoffIn, double dsu
   ew_coeff_ = ew_coeffIn;
   lw_coeff_ = lw_coeffIn;
   switch_width_ = switch_widthIn;
-  erfcTableDx_ = erfcTableDxIn;
+  double erfcTableDx = erfcTableDxIn;
   // Check input
   if (cutoff_ < Constants::SMALL) {
     mprinterr("Error: Direct space cutoff (%g) is too small.\n", cutoff_);
@@ -252,9 +216,14 @@ int Ewald::CheckInput(Box const& boxIn, int debugIn, double cutoffIn, double dsu
     dsumTol_ = 1E-5;
   if (DABS(ew_coeff_) < Constants::SMALL)
     ew_coeff_ = FindEwaldCoefficient( cutoff_, dsumTol_ );
-  if (erfcTableDx_ <= 0.0) erfcTableDx_ = 1.0 / 5000;
+  if (erfcTableDx <= 0.0) erfcTableDx = 1.0 / 5000;
   // TODO make this optional
-  FillErfcTable( cutoff_, ew_coeff_ ); 
+  if (table_.FillTable( erfc_func, erfcTableDx, 0.0, cutoff_*ew_coeff_*1.5 )) {
+    mprinterr("Error: Could not set up spline table for ERFC\n");
+    return 1;
+  }
+  table_.PrintMemUsage("\t");
+  table_.PrintTableInfo("\t");
   // TODO do for C6 as well
   // TODO for C6 correction term
   if (lw_coeff_ < 0.0)