diff --git a/include/utils/PolyatomicDisplacementFunctionBase.h b/include/utils/PolyatomicDisplacementFunctionBase.h
index 413015d2..bb840ac7 100644
--- a/include/utils/PolyatomicDisplacementFunctionBase.h
+++ b/include/utils/PolyatomicDisplacementFunctionBase.h
@@ -51,6 +51,11 @@ class PolyatomicDisplacementFunctionBase
   /// computes the integral of the displacement function
   void computeDisplacementFunctionIntegral();
 
+  /// obtain quadrature rule from GSL
+  static void gslQuadRule(unsigned int quad_order,
+                          std::vector<Real> & quad_points,
+                          std::vector<Real> & quad_weights);
+
   ///@{ some getters needed for accessing this pointer in odeRHS
   unsigned int nSpecies() const { return _n_species; }
   unsigned int problemSize() const { return _problem_size; }
diff --git a/src/utils/PolyatomicDisplacementFunctionBase.C b/src/utils/PolyatomicDisplacementFunctionBase.C
index c108fb33..3deadf7e 100644
--- a/src/utils/PolyatomicDisplacementFunctionBase.C
+++ b/src/utils/PolyatomicDisplacementFunctionBase.C
@@ -72,23 +72,13 @@ PolyatomicDisplacementFunctionBase::PolyatomicDisplacementFunctionBase(
   else
     _Ecap = Ecap;
 
-  // set up integration rule
-  auto * qp_table = gsl_integration_glfixed_table_alloc(_quad_order);
-  _quad_points.resize(_quad_order);
-  _quad_weights.resize(_quad_order);
-  for (std::size_t j = 0; j < _quad_order; ++j)
-  {
-    double point, weight;
-    gsl_integration_glfixed_point(-1.0, 1.0, j, &point, &weight, qp_table);
-    _quad_points[j] = point;
-    _quad_weights[j] = weight;
-  }
-  gsl_integration_glfixed_table_free(qp_table);
-
   // set up the gsl ODE machinery
   auto func = &PolyatomicDisplacementFunction::gslInterface;
   _sys = {func, NULL, _problem_size, this};
   _ode_driver = gsl_odeiv2_driver_alloc_y_new(&_sys, gsl_odeiv2_step_rk4, 10.0, 1e-2, 1.0e-3);
+
+  // set up integration rule
+  gslQuadRule(_quad_order, _quad_points, _quad_weights);
 }
 
 PolyatomicDisplacementFunctionBase::~PolyatomicDisplacementFunctionBase()
@@ -163,6 +153,25 @@ PolyatomicDisplacementFunctionBase::computeDisplacementFunctionIntegral()
   }
 }
 
+void
+PolyatomicDisplacementFunctionBase::gslQuadRule(unsigned int quad_order,
+                                                std::vector<Real> & quad_points,
+                                                std::vector<Real> & quad_weights)
+{
+  // set up integration rule
+  auto * qp_table = gsl_integration_glfixed_table_alloc(quad_order);
+  quad_points.resize(quad_order);
+  quad_weights.resize(quad_order);
+  for (std::size_t j = 0; j < quad_order; ++j)
+  {
+    Real point, weight;
+    gsl_integration_glfixed_point(-1.0, 1.0, j, &point, &weight, qp_table);
+    quad_points[j] = point;
+    quad_weights[j] = weight;
+  }
+  gsl_integration_glfixed_table_free(qp_table);
+}
+
 Real
 PolyatomicDisplacementFunctionBase::stoppingPower(unsigned int species, Real energy)
 {