diff --git a/src/stan/callbacks/stream_logger.hpp b/src/stan/callbacks/stream_logger.hpp
index 6e6cca07e64..677add7c707 100644
--- a/src/stan/callbacks/stream_logger.hpp
+++ b/src/stan/callbacks/stream_logger.hpp
@@ -14,7 +14,7 @@ namespace callbacks {
  * <code>logger</code> that writes messages to separate
  * std::stringstream outputs.
  */
-class stream_logger : public logger {
+class stream_logger final : public logger {
  private:
   std::ostream& debug_;
   std::ostream& info_;
diff --git a/src/stan/callbacks/tee_writer.hpp b/src/stan/callbacks/tee_writer.hpp
index 8eca4af61a7..c491832e9ec 100644
--- a/src/stan/callbacks/tee_writer.hpp
+++ b/src/stan/callbacks/tee_writer.hpp
@@ -16,7 +16,7 @@ namespace callbacks {
  * For any call to this writer, it will tee the call to both writers
  * provided in the constructor.
  */
-class tee_writer : public writer {
+class tee_writer final : public writer {
  public:
   /**
    * Constructor accepting two writers.
diff --git a/src/stan/callbacks/unique_stream_writer.hpp b/src/stan/callbacks/unique_stream_writer.hpp
new file mode 100644
index 00000000000..b576a1d77a5
--- /dev/null
+++ b/src/stan/callbacks/unique_stream_writer.hpp
@@ -0,0 +1,127 @@
+#ifndef STAN_CALLBACKS_UNIQUE_STREAM_WRITER_HPP
+#define STAN_CALLBACKS_UNIQUE_STREAM_WRITER_HPP
+
+#include <stan/callbacks/writer.hpp>
+#include <ostream>
+#include <vector>
+#include <string>
+
+namespace stan {
+namespace callbacks {
+
+/**
+ * <code>unique_stream_writer</code> is an implementation
+ * of <code>writer</code> that holds a unique pointer to the stream it is
+ * writing to.
+ * @tparam Stream A type with with a valid `operator<<(std::string)`
+ */
+template <typename Stream>
+class unique_stream_writer final : public writer {
+ public:
+  /**
+   * Constructs a unique stream writer with an output stream
+   * and an optional prefix for comments.
+   *
+   * @param[in, out] A unique pointer to a type inheriting from `std::ostream`
+   * @param[in] comment_prefix string to stream before each comment line.
+   *  Default is "".
+   */
+  explicit unique_stream_writer(std::unique_ptr<Stream>&& output,
+                                const std::string& comment_prefix = "")
+      : output_(std::move(output)), comment_prefix_(comment_prefix) {}
+
+  unique_stream_writer();
+  unique_stream_writer(unique_stream_writer& other) = delete;
+  unique_stream_writer(unique_stream_writer&& other)
+      : output_(std::move(other.output_)),
+        comment_prefix_(std::move(other.comment_prefix_)) {}
+  /**
+   * Virtual destructor
+   */
+  virtual ~unique_stream_writer() {}
+
+  /**
+   * Writes a set of names on a single line in csv format followed
+   * by a newline.
+   *
+   * Note: the names are not escaped.
+   *
+   * @param[in] names Names in a std::vector
+   */
+  void operator()(const std::vector<std::string>& names) {
+    write_vector(names);
+  }
+  /**
+   * Get the underlying stream
+   */
+  auto& get_stream() { return *output_; }
+
+  /**
+   * Writes a set of values in csv format followed by a newline.
+   *
+   * Note: the precision of the output is determined by the settings
+   *  of the stream on construction.
+   *
+   * @param[in] state Values in a std::vector
+   */
+  void operator()(const std::vector<double>& state) { write_vector(state); }
+
+  /**
+   * Writes the comment_prefix to the stream followed by a newline.
+   */
+  void operator()() {
+    std::stringstream streamer;
+    streamer << comment_prefix_ << std::endl;
+    *output_ << streamer.str();
+  }
+
+  /**
+   * Writes the comment_prefix then the message followed by a newline.
+   *
+   * @param[in] message A string
+   */
+  void operator()(const std::string& message) {
+    std::stringstream streamer;
+    streamer << comment_prefix_ << message << std::endl;
+    *output_ << streamer.str();
+  }
+
+ private:
+  /**
+   * Output stream
+   */
+  std::unique_ptr<Stream> output_;
+
+  /**
+   * Comment prefix to use when printing comments: strings and blank lines
+   */
+  std::string comment_prefix_;
+
+  /**
+   * Writes a set of values in csv format followed by a newline.
+   *
+   * Note: the precision of the output is determined by the settings
+   *  of the stream on construction.
+   *
+   * @param[in] v Values in a std::vector
+   */
+  template <class T>
+  void write_vector(const std::vector<T>& v) {
+    if (v.empty())
+      return;
+    using const_iter = typename std::vector<T>::const_iterator;
+    const_iter last = v.end();
+    --last;
+    std::stringstream streamer;
+    for (const_iter it = v.begin(); it != last; ++it) {
+      streamer << *it << ",";
+    }
+    streamer << v.back() << std::endl;
+    *output_ << streamer.str();
+  }
+};
+
+}  // namespace callbacks
+}  // namespace stan
+
+#endif
diff --git a/src/stan/mcmc/hmc/hamiltonians/softabs_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/softabs_metric.hpp
index 10f24dffe48..da2e392c90b 100644
--- a/src/stan/mcmc/hmc/hamiltonians/softabs_metric.hpp
+++ b/src/stan/mcmc/hmc/hamiltonians/softabs_metric.hpp
@@ -174,26 +174,26 @@ class softabs_metric : public base_hamiltonian<Model, softabs_point, BaseRNG> {
 
   // Threshold below which a power series
   // approximation of the softabs function is used
-  static double lower_softabs_thresh;
+  static constexpr double lower_softabs_thresh = 1e-4;
 
   // Threshold above which an asymptotic
   // approximation of the softabs function is used
-  static double upper_softabs_thresh;
+  static constexpr double upper_softabs_thresh = 18;
 
   // Threshold below which an exact derivative is
   // used in the Jacobian calculation instead of
   // finite differencing
-  static double jacobian_thresh;
+  static constexpr double jacobian_thresh = 1e-10;
 };
 
 template <class Model, class BaseRNG>
-double softabs_metric<Model, BaseRNG>::lower_softabs_thresh = 1e-4;
+constexpr double softabs_metric<Model, BaseRNG>::lower_softabs_thresh;
 
 template <class Model, class BaseRNG>
-double softabs_metric<Model, BaseRNG>::upper_softabs_thresh = 18;
+constexpr double softabs_metric<Model, BaseRNG>::upper_softabs_thresh;
 
 template <class Model, class BaseRNG>
-double softabs_metric<Model, BaseRNG>::jacobian_thresh = 1e-10;
+constexpr double softabs_metric<Model, BaseRNG>::jacobian_thresh;
 }  // namespace mcmc
 }  // namespace stan
 #endif
diff --git a/src/stan/services/sample/hmc_nuts_dense_e_adapt.hpp b/src/stan/services/sample/hmc_nuts_dense_e_adapt.hpp
index cb644a346c6..2733cff1f40 100644
--- a/src/stan/services/sample/hmc_nuts_dense_e_adapt.hpp
+++ b/src/stan/services/sample/hmc_nuts_dense_e_adapt.hpp
@@ -66,7 +66,6 @@ int hmc_nuts_dense_e_adapt(
     callbacks::writer& sample_writer, callbacks::writer& diagnostic_writer) {
   boost::ecuyer1988 rng = util::create_rng(random_seed, chain);
 
-  std::vector<int> disc_vector;
   std::vector<double> cont_vector = util::initialize(
       model, init, rng, init_radius, true, logger, init_writer);
 
@@ -156,6 +155,209 @@ int hmc_nuts_dense_e_adapt(
       interrupt, logger, init_writer, sample_writer, diagnostic_writer);
 }
 
+/**
+ * Runs multiple chains of NUTS with adaptation using dense Euclidean metric
+ * with a pre-specified Euclidean metric.
+ *
+ * @tparam Model Model class
+ * @tparam InitContextPtr A pointer with underlying type derived from
+ `stan::io::var_context`
+ * @tparam InitInvContextPtr A pointer with underlying type derived from
+ `stan::io::var_context`
+ * @tparam SamplerWriter A type derived from `stan::callbacks::writer`
+ * @tparam DiagnosticWriter A type derived from `stan::callbacks::writer`
+ * @tparam InitWriter A type derived from `stan::callbacks::writer`
+ * @param[in] model Input model to test (with data already instantiated)
+ * @param[in] num_chains The number of chains to run in parallel. `init`,
+ * `init_inv_metric`, `init_writer`, `sample_writer`, and `diagnostic_writer`
+ must
+ * be the same length as this value.
+ * @param[in] init An std vector of init var contexts for initialization of each
+ * chain.
+ * @param[in] init_inv_metric An std vector of var contexts exposing an initial
+ * diagonal inverse Euclidean metric for each chain (must be positive definite)
+ * @param[in] random_seed random seed for the random number generator
+ * @param[in] init_chain_id first chain id. The pseudo random number generator
+ * will advance by for each chain by an integer sequence from `init_chain_id` to
+ * `init_chain_id+num_chains-1`
+ * @param[in] init_radius radius to initialize
+ * @param[in] num_warmup Number of warmup samples
+ * @param[in] num_samples Number of samples
+ * @param[in] num_thin Number to thin the samples
+ * @param[in] save_warmup Indicates whether to save the warmup iterations
+ * @param[in] refresh Controls the output
+ * @param[in] stepsize initial stepsize for discrete evolution
+ * @param[in] stepsize_jitter uniform random jitter of stepsize
+ * @param[in] max_depth Maximum tree depth
+ * @param[in] delta adaptation target acceptance statistic
+ * @param[in] gamma adaptation regularization scale
+ * @param[in] kappa adaptation relaxation exponent
+ * @param[in] t0 adaptation iteration offset
+ * @param[in] init_buffer width of initial fast adaptation interval
+ * @param[in] term_buffer width of final fast adaptation interval
+ * @param[in] window initial width of slow adaptation interval
+ * @param[in,out] interrupt Callback for interrupts
+ * @param[in,out] logger Logger for messages
+ * @param[in,out] init_writer std vector of Writer callbacks for unconstrained
+ inits of each chain.
+ * @param[in,out] sample_writer std vector of Writers for draws of each chain.
+ * @param[in,out] diagnostic_writer std vector of Writers for diagnostic
+ * information of each chain.
+ * @return error_codes::OK if successful
+ */
+template <class Model, typename InitContextPtr, typename InitInvContextPtr,
+          typename InitWriter, typename SampleWriter, typename DiagnosticWriter>
+int hmc_nuts_dense_e_adapt(
+    Model& model, size_t num_chains, const std::vector<InitContextPtr>& init,
+    const std::vector<InitInvContextPtr>& init_inv_metric,
+    unsigned int random_seed, unsigned int init_chain_id, double init_radius,
+    int num_warmup, int num_samples, int num_thin, bool save_warmup,
+    int refresh, double stepsize, double stepsize_jitter, int max_depth,
+    double delta, double gamma, double kappa, double t0,
+    unsigned int init_buffer, unsigned int term_buffer, unsigned int window,
+    callbacks::interrupt& interrupt, callbacks::logger& logger,
+    std::vector<InitWriter>& init_writer,
+    std::vector<SampleWriter>& sample_writer,
+    std::vector<DiagnosticWriter>& diagnostic_writer) {
+  if (num_chains == 1) {
+    return hmc_nuts_dense_e_adapt(
+        model, *init[0], *init_inv_metric[0], random_seed, init_chain_id,
+        init_radius, num_warmup, num_samples, num_thin, save_warmup, refresh,
+        stepsize, stepsize_jitter, max_depth, delta, gamma, kappa, t0,
+        init_buffer, term_buffer, window, interrupt, logger, init_writer[0],
+        sample_writer[0], diagnostic_writer[0]);
+  }
+  using sample_t = stan::mcmc::adapt_dense_e_nuts<Model, boost::ecuyer1988>;
+  std::vector<boost::ecuyer1988> rngs;
+  rngs.reserve(num_chains);
+  std::vector<std::vector<double>> cont_vectors;
+  cont_vectors.reserve(num_chains);
+  std::vector<sample_t> samplers;
+  samplers.reserve(num_chains);
+  try {
+    for (int i = 0; i < num_chains; ++i) {
+      rngs.emplace_back(util::create_rng(random_seed, init_chain_id + i));
+      cont_vectors.emplace_back(util::initialize(
+          model, *init[i], rngs[i], init_radius, true, logger, init_writer[i]));
+      Eigen::MatrixXd inv_metric = util::read_dense_inv_metric(
+          *init_inv_metric[i], model.num_params_r(), logger);
+      util::validate_dense_inv_metric(inv_metric, logger);
+
+      samplers.emplace_back(model, rngs[i]);
+      samplers[i].set_metric(inv_metric);
+      samplers[i].set_nominal_stepsize(stepsize);
+      samplers[i].set_stepsize_jitter(stepsize_jitter);
+      samplers[i].set_max_depth(max_depth);
+
+      samplers[i].get_stepsize_adaptation().set_mu(log(10 * stepsize));
+      samplers[i].get_stepsize_adaptation().set_delta(delta);
+      samplers[i].get_stepsize_adaptation().set_gamma(gamma);
+      samplers[i].get_stepsize_adaptation().set_kappa(kappa);
+      samplers[i].get_stepsize_adaptation().set_t0(t0);
+      samplers[i].set_window_params(num_warmup, init_buffer, term_buffer,
+                                    window, logger);
+    }
+  } catch (const std::domain_error& e) {
+    return error_codes::CONFIG;
+  }
+  tbb::parallel_for(tbb::blocked_range<size_t>(0, num_chains, 1),
+                    [num_warmup, num_samples, num_thin, refresh, save_warmup,
+                     num_chains, init_chain_id, &samplers, &model, &rngs,
+                     &interrupt, &logger, &sample_writer, &cont_vectors,
+                     &diagnostic_writer](const tbb::blocked_range<size_t>& r) {
+                      for (size_t i = r.begin(); i != r.end(); ++i) {
+                        util::run_adaptive_sampler(
+                            samplers[i], model, cont_vectors[i], num_warmup,
+                            num_samples, num_thin, refresh, save_warmup,
+                            rngs[i], interrupt, logger, sample_writer[i],
+                            diagnostic_writer[i], init_chain_id + i,
+                            num_chains);
+                      }
+                    },
+                    tbb::simple_partitioner());
+  return error_codes::OK;
+}
+
+/**
+ * Runs multiple chains of NUTS with adaptation using dense Euclidean metric,
+ * with identity matrix as initial inv_metric.
+ *
+ * @tparam Model Model class
+ * @tparam InitContextPtr A pointer with underlying type derived from
+ * `stan::io::var_context`
+ * @tparam InitWriter A type derived from `stan::callbacks::writer`
+ * @tparam SamplerWriter A type derived from `stan::callbacks::writer`
+ * @tparam DiagnosticWriter A type derived from `stan::callbacks::writer`
+ * @param[in] model Input model to test (with data already instantiated)
+ * @param[in] num_chains The number of chains to run in parallel. `init`,
+ * `init_writer`, `sample_writer`, and `diagnostic_writer` must be the same
+ * length as this value.
+ * @param[in] init An std vector of init var contexts for initialization of each
+ * chain.
+ * @param[in] random_seed random seed for the random number generator
+ * @param[in] init_chain_id first chain id. The pseudo random number generator
+ * will advance by for each chain by an integer sequence from `init_chain_id` to
+ * `init_chain_id+num_chains-1`
+ * @param[in] init_radius radius to initialize
+ * @param[in] num_warmup Number of warmup samples
+ * @param[in] num_samples Number of samples
+ * @param[in] num_thin Number to thin the samples
+ * @param[in] save_warmup Indicates whether to save the warmup iterations
+ * @param[in] refresh Controls the output
+ * @param[in] stepsize initial stepsize for discrete evolution
+ * @param[in] stepsize_jitter uniform random jitter of stepsize
+ * @param[in] max_depth Maximum tree depth
+ * @param[in] delta adaptation target acceptance statistic
+ * @param[in] gamma adaptation regularization scale
+ * @param[in] kappa adaptation relaxation exponent
+ * @param[in] t0 adaptation iteration offset
+ * @param[in] init_buffer width of initial fast adaptation interval
+ * @param[in] term_buffer width of final fast adaptation interval
+ * @param[in] window initial width of slow adaptation interval
+ * @param[in,out] interrupt Callback for interrupts
+ * @param[in,out] logger Logger for messages
+ * @param[in,out] init_writer std vector of Writer callbacks for unconstrained
+ * inits of each chain.
+ * @param[in,out] sample_writer std vector of Writers for draws of each chain.
+ * @param[in,out] diagnostic_writer std vector of Writers for diagnostic
+ * information of each chain.
+ * @return error_codes::OK if successful
+ */
+template <class Model, typename InitContextPtr, typename InitWriter,
+          typename SampleWriter, typename DiagnosticWriter>
+int hmc_nuts_dense_e_adapt(
+    Model& model, size_t num_chains, const std::vector<InitContextPtr>& init,
+    unsigned int random_seed, unsigned int init_chain_id, double init_radius,
+    int num_warmup, int num_samples, int num_thin, bool save_warmup,
+    int refresh, double stepsize, double stepsize_jitter, int max_depth,
+    double delta, double gamma, double kappa, double t0,
+    unsigned int init_buffer, unsigned int term_buffer, unsigned int window,
+    callbacks::interrupt& interrupt, callbacks::logger& logger,
+    std::vector<InitWriter>& init_writer,
+    std::vector<SampleWriter>& sample_writer,
+    std::vector<DiagnosticWriter>& diagnostic_writer) {
+  if (num_chains == 1) {
+    return hmc_nuts_dense_e_adapt(
+        model, *init[0], random_seed, init_chain_id, init_radius, num_warmup,
+        num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+        max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+        interrupt, logger, init_writer[0], sample_writer[0],
+        diagnostic_writer[0]);
+  }
+  std::vector<std::unique_ptr<stan::io::dump>> unit_e_metrics;
+  unit_e_metrics.reserve(num_chains);
+  for (size_t i = 0; i < num_chains; ++i) {
+    unit_e_metrics.emplace_back(std::make_unique<stan::io::dump>(
+        util::create_unit_e_dense_inv_metric(model.num_params_r())));
+  }
+  return hmc_nuts_dense_e_adapt(
+      model, num_chains, init, unit_e_metrics, random_seed, init_chain_id,
+      init_radius, num_warmup, num_samples, num_thin, save_warmup, refresh,
+      stepsize, stepsize_jitter, max_depth, delta, gamma, kappa, t0,
+      init_buffer, term_buffer, window, interrupt, logger, init_writer,
+      sample_writer, diagnostic_writer);
+}
+
 }  // namespace sample
 }  // namespace services
 }  // namespace stan
diff --git a/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp b/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp
index 8d03373e1a5..9689fa4a9a1 100644
--- a/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp
+++ b/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp
@@ -24,6 +24,11 @@ namespace sample {
  * with a pre-specified Euclidean metric.
  *
  * @tparam Model Model class
+ * @tparam InitContextPtr A type derived from `stan::io::var_context`
+ * @tparam InitMetricContext A type derived from `stan::io::var_context`
+ * @tparam SamplerWriter A type derived from `stan::callbacks::writer`
+ * @tparam DiagnosticWriter A type derived from `stan::callbacks::writer`
+ * @tparam InitWriter A type derived from `stan::callbacks::writer`
  * @param[in] model Input model to test (with data already instantiated)
  * @param[in] init var context for initialization
  * @param[in] init_inv_metric var context exposing an initial diagonal
@@ -53,7 +58,7 @@ namespace sample {
  * @param[in,out] diagnostic_writer Writer for diagnostic information
  * @return error_codes::OK if successful
  */
-template <class Model>
+template <typename Model>
 int hmc_nuts_diag_e_adapt(
     Model& model, const stan::io::var_context& init,
     const stan::io::var_context& init_inv_metric, unsigned int random_seed,
@@ -66,7 +71,6 @@ int hmc_nuts_diag_e_adapt(
     callbacks::writer& sample_writer, callbacks::writer& diagnostic_writer) {
   boost::ecuyer1988 rng = util::create_rng(random_seed, chain);
 
-  std::vector<int> disc_vector;
   std::vector<double> cont_vector = util::initialize(
       model, init, rng, init_radius, true, logger, init_writer);
 
@@ -133,7 +137,7 @@ int hmc_nuts_diag_e_adapt(
  * @param[in,out] diagnostic_writer Writer for diagnostic information
  * @return error_codes::OK if successful
  */
-template <class Model>
+template <typename Model>
 int hmc_nuts_diag_e_adapt(
     Model& model, const stan::io::var_context& init, unsigned int random_seed,
     unsigned int chain, double init_radius, int num_warmup, int num_samples,
@@ -143,10 +147,8 @@ int hmc_nuts_diag_e_adapt(
     unsigned int window, callbacks::interrupt& interrupt,
     callbacks::logger& logger, callbacks::writer& init_writer,
     callbacks::writer& sample_writer, callbacks::writer& diagnostic_writer) {
-  stan::io::dump dmp
+  stan::io::dump unit_e_metric
       = util::create_unit_e_diag_inv_metric(model.num_params_r());
-  stan::io::var_context& unit_e_metric = dmp;
-
   return hmc_nuts_diag_e_adapt(
       model, init, unit_e_metric, random_seed, chain, init_radius, num_warmup,
       num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
@@ -154,6 +156,209 @@ int hmc_nuts_diag_e_adapt(
       interrupt, logger, init_writer, sample_writer, diagnostic_writer);
 }
 
+/**
+ * Runs multiple chains of HMC with NUTS with adaptation using diagonal
+ * Euclidean metric with a pre-specified Euclidean metric.
+ *
+ * @tparam Model Model class
+ * @tparam InitContextPtr A pointer with underlying type derived from
+ `stan::io::var_context`
+ * @tparam InitInvContextPtr A pointer with underlying type derived from
+ `stan::io::var_context`
+ * @tparam SamplerWriter A type derived from `stan::callbacks::writer`
+ * @tparam DiagnosticWriter A type derived from `stan::callbacks::writer`
+ * @tparam InitWriter A type derived from `stan::callbacks::writer`
+ * @param[in] model Input model to test (with data already instantiated)
+ * @param[in] num_chains The number of chains to run in parallel. `init`,
+ * `init_inv_metric`, `init_writer`, `sample_writer`, and `diagnostic_writer`
+ must
+ * be the same length as this value.
+ * @param[in] init An std vector of init var contexts for initialization of each
+ * chain.
+ * @param[in] init_inv_metric An std vector of var contexts exposing an initial
+ diagonal inverse Euclidean metric for each chain (must be positive definite)
+ * @param[in] random_seed random seed for the random number generator
+ * @param[in] init_chain_id first chain id. The pseudo random number generator
+ * will advance for each chain by an integer sequence from `init_chain_id` to
+ * `init_chain_id + num_chains - 1`
+ * @param[in] init_radius radius to initialize
+ * @param[in] num_warmup Number of warmup samples
+ * @param[in] num_samples Number of samples
+ * @param[in] num_thin Number to thin the samples
+ * @param[in] save_warmup Indicates whether to save the warmup iterations
+ * @param[in] refresh Controls the output
+ * @param[in] stepsize initial stepsize for discrete evolution
+ * @param[in] stepsize_jitter uniform random jitter of stepsize
+ * @param[in] max_depth Maximum tree depth
+ * @param[in] delta adaptation target acceptance statistic
+ * @param[in] gamma adaptation regularization scale
+ * @param[in] kappa adaptation relaxation exponent
+ * @param[in] t0 adaptation iteration offset
+ * @param[in] init_buffer width of initial fast adaptation interval
+ * @param[in] term_buffer width of final fast adaptation interval
+ * @param[in] window initial width of slow adaptation interval
+ * @param[in,out] interrupt Callback for interrupts
+ * @param[in,out] logger Logger for messages
+ * @param[in,out] init_writer std vector of Writer callbacks for unconstrained
+ * inits of each chain.
+ * @param[in,out] sample_writer std vector of Writers for draws of each chain.
+ * @param[in,out] diagnostic_writer std vector of Writers for diagnostic
+ * information of each chain.
+ * @return error_codes::OK if successful
+ */
+template <class Model, typename InitContextPtr, typename InitInvContextPtr,
+          typename InitWriter, typename SampleWriter, typename DiagnosticWriter>
+int hmc_nuts_diag_e_adapt(
+    Model& model, size_t num_chains, const std::vector<InitContextPtr>& init,
+    const std::vector<InitInvContextPtr>& init_inv_metric,
+    unsigned int random_seed, unsigned int init_chain_id, double init_radius,
+    int num_warmup, int num_samples, int num_thin, bool save_warmup,
+    int refresh, double stepsize, double stepsize_jitter, int max_depth,
+    double delta, double gamma, double kappa, double t0,
+    unsigned int init_buffer, unsigned int term_buffer, unsigned int window,
+    callbacks::interrupt& interrupt, callbacks::logger& logger,
+    std::vector<InitWriter>& init_writer,
+    std::vector<SampleWriter>& sample_writer,
+    std::vector<DiagnosticWriter>& diagnostic_writer) {
+  if (num_chains == 1) {
+    return hmc_nuts_diag_e_adapt(
+        model, *init[0], *init_inv_metric[0], random_seed, init_chain_id,
+        init_radius, num_warmup, num_samples, num_thin, save_warmup, refresh,
+        stepsize, stepsize_jitter, max_depth, delta, gamma, kappa, t0,
+        init_buffer, term_buffer, window, interrupt, logger, init_writer[0],
+        sample_writer[0], diagnostic_writer[0]);
+  }
+  using sample_t = stan::mcmc::adapt_diag_e_nuts<Model, boost::ecuyer1988>;
+  std::vector<boost::ecuyer1988> rngs;
+  rngs.reserve(num_chains);
+  std::vector<std::vector<double>> cont_vectors;
+  cont_vectors.reserve(num_chains);
+  std::vector<sample_t> samplers;
+  samplers.reserve(num_chains);
+  try {
+    for (int i = 0; i < num_chains; ++i) {
+      rngs.emplace_back(util::create_rng(random_seed, init_chain_id + i));
+      cont_vectors.emplace_back(util::initialize(
+          model, *init[i], rngs[i], init_radius, true, logger, init_writer[i]));
+      samplers.emplace_back(model, rngs[i]);
+      Eigen::VectorXd inv_metric = util::read_diag_inv_metric(
+          *init_inv_metric[i], model.num_params_r(), logger);
+      util::validate_diag_inv_metric(inv_metric, logger);
+
+      samplers[i].set_metric(inv_metric);
+      samplers[i].set_nominal_stepsize(stepsize);
+      samplers[i].set_stepsize_jitter(stepsize_jitter);
+      samplers[i].set_max_depth(max_depth);
+
+      samplers[i].get_stepsize_adaptation().set_mu(log(10 * stepsize));
+      samplers[i].get_stepsize_adaptation().set_delta(delta);
+      samplers[i].get_stepsize_adaptation().set_gamma(gamma);
+      samplers[i].get_stepsize_adaptation().set_kappa(kappa);
+      samplers[i].get_stepsize_adaptation().set_t0(t0);
+      samplers[i].set_window_params(num_warmup, init_buffer, term_buffer,
+                                    window, logger);
+    }
+  } catch (const std::domain_error& e) {
+    return error_codes::CONFIG;
+  }
+  tbb::parallel_for(tbb::blocked_range<size_t>(0, num_chains, 1),
+                    [num_warmup, num_samples, num_thin, refresh, save_warmup,
+                     num_chains, init_chain_id, &samplers, &model, &rngs,
+                     &interrupt, &logger, &sample_writer, &cont_vectors,
+                     &diagnostic_writer](const tbb::blocked_range<size_t>& r) {
+                      for (size_t i = r.begin(); i != r.end(); ++i) {
+                        util::run_adaptive_sampler(
+                            samplers[i], model, cont_vectors[i], num_warmup,
+                            num_samples, num_thin, refresh, save_warmup,
+                            rngs[i], interrupt, logger, sample_writer[i],
+                            diagnostic_writer[i], init_chain_id + i,
+                            num_chains);
+                      }
+                    },
+                    tbb::simple_partitioner());
+  return error_codes::OK;
+}
+
+/**
+ * Runs multiple chains of HMC with NUTS with adaptation using diagonal
+ * Euclidean metric.
+ *
+ * @tparam Model Model class
+ * @tparam InitContextPtr A pointer with underlying type derived from
+ * `stan::io::var_context`
+ * @tparam SamplerWriter A type derived from `stan::callbacks::writer`
+ * @tparam DiagnosticWriter A type derived from `stan::callbacks::writer`
+ * @tparam InitWriter A type derived from `stan::callbacks::writer`
+ * @param[in] model Input model to test (with data already instantiated)
+ * @param[in] num_chains The number of chains to run in parallel. `init`,
+ * `init_writer`, `sample_writer`, and `diagnostic_writer` must be the same
+ * length as this value.
+ * @param[in] init An std vector of init var contexts for initialization of each
+ * chain.
+ * @param[in] random_seed random seed for the random number generator
+ * @param[in] init_chain_id first chain id. The pseudo random number generator
+ * will advance by for each chain by an integer sequence from `init_chain_id` to
+ * `init_chain_id+num_chains-1`
+ * @param[in] init_radius radius to initialize
+ * @param[in] num_warmup Number of warmup samples
+ * @param[in] num_samples Number of samples
+ * @param[in] num_thin Number to thin the samples
+ * @param[in] save_warmup Indicates whether to save the warmup iterations
+ * @param[in] refresh Controls the output
+ * @param[in] stepsize initial stepsize for discrete evolution
+ * @param[in] stepsize_jitter uniform random jitter of stepsize
+ * @param[in] max_depth Maximum tree depth
+ * @param[in] delta adaptation target acceptance statistic
+ * @param[in] gamma adaptation regularization scale
+ * @param[in] kappa adaptation relaxation exponent
+ * @param[in] t0 adaptation iteration offset
+ * @param[in] init_buffer width of initial fast adaptation interval
+ * @param[in] term_buffer width of final fast adaptation interval
+ * @param[in] window initial width of slow adaptation interval
+ * @param[in,out] interrupt Callback for interrupts
+ * @param[in,out] logger Logger for messages
+ * @param[in,out] init_writer std vector of Writer callbacks for unconstrained
+ * inits of each chain.
+ * @param[in,out] sample_writer std vector of Writers for draws of each chain.
+ * @param[in,out] diagnostic_writer std vector of Writers for diagnostic
+ * information of each chain.
+ * @return error_codes::OK if successful
+ */
+template <class Model, typename InitContextPtr, typename InitWriter,
+          typename SampleWriter, typename DiagnosticWriter>
+int hmc_nuts_diag_e_adapt(
+    Model& model, size_t num_chains, const std::vector<InitContextPtr>& init,
+    unsigned int random_seed, unsigned int init_chain_id, double init_radius,
+    int num_warmup, int num_samples, int num_thin, bool save_warmup,
+    int refresh, double stepsize, double stepsize_jitter, int max_depth,
+    double delta, double gamma, double kappa, double t0,
+    unsigned int init_buffer, unsigned int term_buffer, unsigned int window,
+    callbacks::interrupt& interrupt, callbacks::logger& logger,
+    std::vector<InitWriter>& init_writer,
+    std::vector<SampleWriter>& sample_writer,
+    std::vector<DiagnosticWriter>& diagnostic_writer) {
+  if (num_chains == 1) {
+    return hmc_nuts_diag_e_adapt(
+        model, *init[0], random_seed, init_chain_id, init_radius, num_warmup,
+        num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+        max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+        interrupt, logger, init_writer[0], sample_writer[0],
+        diagnostic_writer[0]);
+  }
+  std::vector<std::unique_ptr<stan::io::dump>> unit_e_metrics;
+  unit_e_metrics.reserve(num_chains);
+  for (size_t i = 0; i < num_chains; ++i) {
+    unit_e_metrics.emplace_back(std::make_unique<stan::io::dump>(
+        util::create_unit_e_diag_inv_metric(model.num_params_r())));
+  }
+  return hmc_nuts_diag_e_adapt(
+      model, num_chains, init, unit_e_metrics, random_seed, init_chain_id,
+      init_radius, num_warmup, num_samples, num_thin, save_warmup, refresh,
+      stepsize, stepsize_jitter, max_depth, delta, gamma, kappa, t0,
+      init_buffer, term_buffer, window, interrupt, logger, init_writer,
+      sample_writer, diagnostic_writer);
+}
+
 }  // namespace sample
 }  // namespace services
 }  // namespace stan
diff --git a/src/stan/services/util/create_rng.hpp b/src/stan/services/util/create_rng.hpp
index 9b32d4509c0..d63d6c79f1b 100644
--- a/src/stan/services/util/create_rng.hpp
+++ b/src/stan/services/util/create_rng.hpp
@@ -24,7 +24,7 @@ namespace util {
  */
 inline boost::ecuyer1988 create_rng(unsigned int seed, unsigned int chain) {
   using boost::uintmax_t;
-  static uintmax_t DISCARD_STRIDE = static_cast<uintmax_t>(1) << 50;
+  static constexpr uintmax_t DISCARD_STRIDE = static_cast<uintmax_t>(1) << 50;
   boost::ecuyer1988 rng(seed);
   rng.discard(DISCARD_STRIDE * chain);
   return rng;
diff --git a/src/stan/services/util/create_unit_e_dense_inv_metric.hpp b/src/stan/services/util/create_unit_e_dense_inv_metric.hpp
index 60de41e0875..ee9afc670a6 100644
--- a/src/stan/services/util/create_unit_e_dense_inv_metric.hpp
+++ b/src/stan/services/util/create_unit_e_dense_inv_metric.hpp
@@ -17,17 +17,12 @@ namespace util {
  * @return var_context
  */
 inline stan::io::dump create_unit_e_dense_inv_metric(size_t num_params) {
-  Eigen::MatrixXd inv_metric(num_params, num_params);
-  inv_metric.setIdentity();
-  size_t num_elements = num_params * num_params;
+  auto num_params_str = std::to_string(num_params);
+  std::string dims("),.Dim=c(" + num_params_str + ", " + num_params_str + "))");
+  Eigen::IOFormat RFmt(Eigen::StreamPrecision, Eigen::DontAlignCols, ", ", ",",
+                       "", "", "inv_metric <- structure(c(", dims);
   std::stringstream txt;
-  txt << "inv_metric <- structure(c(";
-  for (size_t i = 0; i < num_elements; i++) {
-    txt << inv_metric(i);
-    if (i < num_elements - 1)
-      txt << ", ";
-  }
-  txt << "),.Dim=c(" << num_params << ", " << num_params << "))";
+  txt << Eigen::MatrixXd::Identity(num_params, num_params).format(RFmt);
   return stan::io::dump(txt);
 }
 }  // namespace util
diff --git a/src/stan/services/util/create_unit_e_diag_inv_metric.hpp b/src/stan/services/util/create_unit_e_diag_inv_metric.hpp
index 18375e68f94..e14a80e3283 100644
--- a/src/stan/services/util/create_unit_e_diag_inv_metric.hpp
+++ b/src/stan/services/util/create_unit_e_diag_inv_metric.hpp
@@ -16,14 +16,11 @@ namespace util {
  * @return var_context
  */
 inline stan::io::dump create_unit_e_diag_inv_metric(size_t num_params) {
+  std::string dims("),.Dim=c(" + std::to_string(num_params) + "))");
+  Eigen::IOFormat RFmt(Eigen::StreamPrecision, Eigen::DontAlignCols, ", ", ",",
+                       "", "", "inv_metric <- structure(c(", dims);
   std::stringstream txt;
-  txt << "inv_metric <- structure(c(";
-  for (size_t i = 0; i < num_params; ++i) {
-    txt << "1.0";
-    if (i < num_params - 1)
-      txt << ", ";
-  }
-  txt << "),.Dim=c(" << num_params << "))";
+  txt << Eigen::VectorXd::Ones(num_params).format(RFmt);
   return stan::io::dump(txt);
 }
 }  // namespace util
diff --git a/src/stan/services/util/generate_transitions.hpp b/src/stan/services/util/generate_transitions.hpp
index 1516c6f0f99..5b1aa1a27dc 100644
--- a/src/stan/services/util/generate_transitions.hpp
+++ b/src/stan/services/util/generate_transitions.hpp
@@ -36,6 +36,8 @@ namespace util {
  * @param[in,out] base_rng random number generator
  * @param[in,out] callback interrupt callback called once an iteration
  * @param[in,out] logger logger for messages
+ * @param[in] num_chains The number of chains used in the program. This
+ *  is used in generate transitions to print out the chain number.
  */
 template <class Model, class RNG>
 void generate_transitions(stan::mcmc::base_mcmc& sampler, int num_iterations,
@@ -44,7 +46,8 @@ void generate_transitions(stan::mcmc::base_mcmc& sampler, int num_iterations,
                           util::mcmc_writer& mcmc_writer,
                           stan::mcmc::sample& init_s, Model& model,
                           RNG& base_rng, callbacks::interrupt& callback,
-                          callbacks::logger& logger) {
+                          callbacks::logger& logger, size_t chain_id = 1,
+                          size_t num_chains = 1) {
   for (int m = 0; m < num_iterations; ++m) {
     callback();
 
@@ -52,6 +55,9 @@ void generate_transitions(stan::mcmc::base_mcmc& sampler, int num_iterations,
         && (start + m + 1 == finish || m == 0 || (m + 1) % refresh == 0)) {
       int it_print_width = std::ceil(std::log10(static_cast<double>(finish)));
       std::stringstream message;
+      if (num_chains != 1) {
+        message << "Chain [" << chain_id << "] ";
+      }
       message << "Iteration: ";
       message << std::setw(it_print_width) << m + 1 + start << " / " << finish;
       message << " [" << std::setw(3)
diff --git a/src/stan/services/util/initialize.hpp b/src/stan/services/util/initialize.hpp
index fd180819b99..696013531c9 100644
--- a/src/stan/services/util/initialize.hpp
+++ b/src/stan/services/util/initialize.hpp
@@ -65,9 +65,10 @@ namespace util {
  *   std::domain_error)
  * @return valid unconstrained parameters for the model
  */
-template <bool Jacobian = true, class Model, class RNG>
-std::vector<double> initialize(Model& model, const stan::io::var_context& init,
-                               RNG& rng, double init_radius, bool print_timing,
+template <bool Jacobian = true, typename Model, typename InitContext,
+          typename RNG>
+std::vector<double> initialize(Model& model, const InitContext& init, RNG& rng,
+                               double init_radius, bool print_timing,
                                stan::callbacks::logger& logger,
                                stan::callbacks::writer& init_writer) {
   std::vector<double> unconstrained;
diff --git a/src/stan/services/util/run_adaptive_sampler.hpp b/src/stan/services/util/run_adaptive_sampler.hpp
index ba0135158dd..7115b138fc3 100644
--- a/src/stan/services/util/run_adaptive_sampler.hpp
+++ b/src/stan/services/util/run_adaptive_sampler.hpp
@@ -5,6 +5,7 @@
 #include <stan/callbacks/writer.hpp>
 #include <stan/services/util/generate_transitions.hpp>
 #include <stan/services/util/mcmc_writer.hpp>
+#include <tbb/parallel_for.h>
 #include <chrono>
 #include <vector>
 
@@ -33,8 +34,11 @@ namespace util {
  * @param[in,out] logger logger for messages
  * @param[in,out] sample_writer writer for draws
  * @param[in,out] diagnostic_writer writer for diagnostic information
+ * @param[in] chain_id The id for a given chain.
+ * @param[in] num_chains The number of chains used in the program. This
+ *  is used in generate transitions to print out the chain number.
  */
-template <class Sampler, class Model, class RNG>
+template <typename Sampler, typename Model, typename RNG>
 void run_adaptive_sampler(Sampler& sampler, Model& model,
                           std::vector<double>& cont_vector, int num_warmup,
                           int num_samples, int num_thin, int refresh,
@@ -42,7 +46,8 @@ void run_adaptive_sampler(Sampler& sampler, Model& model,
                           callbacks::interrupt& interrupt,
                           callbacks::logger& logger,
                           callbacks::writer& sample_writer,
-                          callbacks::writer& diagnostic_writer) {
+                          callbacks::writer& diagnostic_writer,
+                          size_t chain_id = 1, size_t num_chains = 1) {
   Eigen::Map<Eigen::VectorXd> cont_params(cont_vector.data(),
                                           cont_vector.size());
 
@@ -66,7 +71,8 @@ void run_adaptive_sampler(Sampler& sampler, Model& model,
   auto start_warm = std::chrono::steady_clock::now();
   util::generate_transitions(sampler, num_warmup, 0, num_warmup + num_samples,
                              num_thin, refresh, save_warmup, true, writer, s,
-                             model, rng, interrupt, logger);
+                             model, rng, interrupt, logger, chain_id,
+                             num_chains);
   auto end_warm = std::chrono::steady_clock::now();
   double warm_delta_t = std::chrono::duration_cast<std::chrono::milliseconds>(
                             end_warm - start_warm)
@@ -79,7 +85,8 @@ void run_adaptive_sampler(Sampler& sampler, Model& model,
   auto start_sample = std::chrono::steady_clock::now();
   util::generate_transitions(sampler, num_samples, num_warmup,
                              num_warmup + num_samples, num_thin, refresh, true,
-                             false, writer, s, model, rng, interrupt, logger);
+                             false, writer, s, model, rng, interrupt, logger,
+                             chain_id, num_chains);
   auto end_sample = std::chrono::steady_clock::now();
   double sample_delta_t = std::chrono::duration_cast<std::chrono::milliseconds>(
                               end_sample - start_sample)
@@ -87,6 +94,7 @@ void run_adaptive_sampler(Sampler& sampler, Model& model,
                           / 1000.0;
   writer.write_timing(warm_delta_t, sample_delta_t);
 }
+
 }  // namespace util
 }  // namespace services
 }  // namespace stan
diff --git a/src/test/test-models/good/optimization/rosenbrock.stan b/src/test/test-models/good/optimization/rosenbrock.stan
index db310524279..6cb7b8f8dec 100644
--- a/src/test/test-models/good/optimization/rosenbrock.stan
+++ b/src/test/test-models/good/optimization/rosenbrock.stan
@@ -5,4 +5,3 @@ parameters {
 model {
   target += -(pow(1 - x, 2) + 100 * pow(y - pow(x, 2), 2));
 }
-
diff --git a/src/test/unit/callbacks/unique_stream_writer_test.cpp b/src/test/unit/callbacks/unique_stream_writer_test.cpp
new file mode 100644
index 00000000000..6fca0f1808f
--- /dev/null
+++ b/src/test/unit/callbacks/unique_stream_writer_test.cpp
@@ -0,0 +1,50 @@
+#include <gtest/gtest.h>
+#include <boost/lexical_cast.hpp>
+#include <stan/callbacks/unique_stream_writer.hpp>
+
+class StanInterfaceCallbacksStreamWriter : public ::testing::Test {
+ public:
+  StanInterfaceCallbacksStreamWriter()
+      : writer(std::make_unique<std::stringstream>(std::stringstream{})) {}
+
+  void SetUp() {
+    static_cast<std::stringstream&>(writer.get_stream()).str(std::string());
+    static_cast<std::stringstream&>(writer.get_stream()).clear();
+  }
+  void TearDown() {}
+
+  stan::callbacks::unique_stream_writer<std::ostream> writer;
+};
+
+TEST_F(StanInterfaceCallbacksStreamWriter, double_vector) {
+  const int N = 5;
+  std::vector<double> x;
+  for (int n = 0; n < N; ++n)
+    x.push_back(n);
+
+  EXPECT_NO_THROW(writer(x));
+  EXPECT_EQ("0,1,2,3,4\n",
+            static_cast<std::stringstream&>(writer.get_stream()).str());
+}
+
+TEST_F(StanInterfaceCallbacksStreamWriter, string_vector) {
+  const int N = 5;
+  std::vector<std::string> x;
+  for (int n = 0; n < N; ++n)
+    x.push_back(boost::lexical_cast<std::string>(n));
+
+  EXPECT_NO_THROW(writer(x));
+  EXPECT_EQ("0,1,2,3,4\n",
+            static_cast<std::stringstream&>(writer.get_stream()).str());
+}
+
+TEST_F(StanInterfaceCallbacksStreamWriter, null) {
+  EXPECT_NO_THROW(writer());
+  EXPECT_EQ("\n", static_cast<std::stringstream&>(writer.get_stream()).str());
+}
+
+TEST_F(StanInterfaceCallbacksStreamWriter, string) {
+  EXPECT_NO_THROW(writer("message"));
+  EXPECT_EQ("message\n",
+            static_cast<std::stringstream&>(writer.get_stream()).str());
+}
diff --git a/src/test/unit/services/check_adaptation.hpp b/src/test/unit/services/check_adaptation.hpp
index e5d3f52dd4c..af863ea8139 100644
--- a/src/test/unit/services/check_adaptation.hpp
+++ b/src/test/unit/services/check_adaptation.hpp
@@ -13,11 +13,7 @@ namespace stan {
 namespace test {
 namespace unit {
 
-double stod(const std::string& val) {
-  char tmp[val.length()];
-  strcpy(tmp, val.c_str());
-  return atof(tmp);
-}
+double stod(const std::string& val) { return atof(val.c_str()); }
 
 void check_adaptation(const size_t& num_params,
                       const std::vector<double>& param_vals,
diff --git a/src/test/unit/services/instrumented_callbacks.hpp b/src/test/unit/services/instrumented_callbacks.hpp
index 5df61c17f74..797959ef0ee 100644
--- a/src/test/unit/services/instrumented_callbacks.hpp
+++ b/src/test/unit/services/instrumented_callbacks.hpp
@@ -9,7 +9,8 @@
 #include <string>
 #include <iostream>
 #include <exception>
-
+#include <atomic>
+#include <mutex>
 namespace stan {
 namespace test {
 namespace unit {
@@ -27,7 +28,7 @@ class instrumented_interrupt : public stan::callbacks::interrupt {
   unsigned int call_count() { return counter_; }
 
  private:
-  unsigned int counter_;
+  std::atomic<unsigned int> counter_;
 };
 
 /**
@@ -96,55 +97,54 @@ class instrumented_writer : public stan::callbacks::writer {
 
   unsigned int call_count() {
     unsigned int n = 0;
-    for (std::map<std::string, int>::iterator it = counter_.begin();
-         it != counter_.end(); ++it)
-      n += it->second;
+    for (auto& it : counter_)
+      n += it.second;
     return n;
   }
 
   unsigned int call_count(std::string s) { return counter_[s]; }
 
-  std::vector<std::pair<std::string, double> > string_double_values() {
+  std::vector<std::pair<std::string, double>> string_double_values() {
     return string_double;
   };
 
-  std::vector<std::pair<std::string, int> > string_int_values() {
+  std::vector<std::pair<std::string, int>> string_int_values() {
     return string_int;
   };
 
-  std::vector<std::pair<std::string, std::string> > string_string_values() {
+  std::vector<std::pair<std::string, std::string>> string_string_values() {
     return string_string;
   };
 
-  std::vector<std::pair<std::string, std::vector<double> > >
+  std::vector<std::pair<std::string, std::vector<double>>>
   string_pdouble_int_values() {
     return string_pdouble_int;
   };
 
-  std::vector<std::pair<std::string, Eigen::MatrixXd> >
+  std::vector<std::pair<std::string, Eigen::MatrixXd>>
   string_pdouble_int_int_values() {
     return string_pdouble_int_int;
   };
 
-  std::vector<std::vector<std::string> > vector_string_values() {
+  std::vector<std::vector<std::string>> vector_string_values() {
     return vector_string;
   };
 
-  std::vector<std::vector<double> > vector_double_values() {
+  std::vector<std::vector<double>> vector_double_values() {
     return vector_double;
   };
 
   std::vector<std::string> string_values() { return string; };
 
  private:
-  std::map<std::string, int> counter_;
-  std::vector<std::pair<std::string, double> > string_double;
-  std::vector<std::pair<std::string, int> > string_int;
-  std::vector<std::pair<std::string, std::string> > string_string;
-  std::vector<std::pair<std::string, std::vector<double> > > string_pdouble_int;
-  std::vector<std::pair<std::string, Eigen::MatrixXd> > string_pdouble_int_int;
-  std::vector<std::vector<std::string> > vector_string;
-  std::vector<std::vector<double> > vector_double;
+  std::map<std::string, std::atomic<int>> counter_;
+  std::vector<std::pair<std::string, double>> string_double;
+  std::vector<std::pair<std::string, int>> string_int;
+  std::vector<std::pair<std::string, std::string>> string_string;
+  std::vector<std::pair<std::string, std::vector<double>>> string_pdouble_int;
+  std::vector<std::pair<std::string, Eigen::MatrixXd>> string_pdouble_int_int;
+  std::vector<std::vector<std::string>> vector_string;
+  std::vector<std::vector<double>> vector_double;
   std::vector<std::string> string;
 };
 
@@ -156,35 +156,53 @@ class instrumented_writer : public stan::callbacks::writer {
  */
 class instrumented_logger : public stan::callbacks::logger {
  public:
+  std::mutex logger_guard;
   instrumented_logger() {}
 
-  void debug(const std::string& message) { debug_.push_back(message); }
+  void debug(const std::string& message) {
+    std::lock_guard<std::mutex> guard(logger_guard);
+    debug_.push_back(message);
+  }
 
   void debug(const std::stringstream& message) {
+    std::lock_guard<std::mutex> guard(logger_guard);
     debug_.push_back(message.str());
   }
 
-  void info(const std::string& message) { info_.push_back(message); }
+  void info(const std::string& message) {
+    std::lock_guard<std::mutex> guard(logger_guard);
+    info_.push_back(message);
+  }
 
   void info(const std::stringstream& message) {
+    std::lock_guard<std::mutex> guard(logger_guard);
     info_.push_back(message.str());
   }
 
-  void warn(const std::string& message) { warn_.push_back(message); }
+  void warn(const std::string& message) {
+    std::lock_guard<std::mutex> guard(logger_guard);
+    warn_.push_back(message);
+  }
 
   void warn(const std::stringstream& message) {
+    std::lock_guard<std::mutex> guard(logger_guard);
     warn_.push_back(message.str());
   }
 
-  void error(const std::string& message) { error_.push_back(message); }
+  void error(const std::string& message) {
+    std::lock_guard<std::mutex> guard(logger_guard);
+    error_.push_back(message);
+  }
 
   void error(const std::stringstream& message) {
+    std::lock_guard<std::mutex> guard(logger_guard);
     error_.push_back(message.str());
   }
 
   void fatal(const std::string& message) { fatal_.push_back(message); }
 
   void fatal(const std::stringstream& message) {
+    std::lock_guard<std::mutex> guard(logger_guard);
     fatal_.push_back(message.str());
   }
 
diff --git a/src/test/unit/services/sample/hmc_nuts_dense_e_adapt_parallel_match_test.cpp b/src/test/unit/services/sample/hmc_nuts_dense_e_adapt_parallel_match_test.cpp
new file mode 100644
index 00000000000..901afbc6ab7
--- /dev/null
+++ b/src/test/unit/services/sample/hmc_nuts_dense_e_adapt_parallel_match_test.cpp
@@ -0,0 +1,108 @@
+#include <stan/services/sample/hmc_nuts_dense_e_adapt.hpp>
+#include <stan/io/empty_var_context.hpp>
+#include <stan/callbacks/unique_stream_writer.hpp>
+#include <test/unit/util.hpp>
+#include <src/test/unit/services/util.hpp>
+#include <test/test-models/good/optimization/rosenbrock.hpp>
+#include <test/unit/services/instrumented_callbacks.hpp>
+#include <gtest/gtest.h>
+#include <iostream>
+
+auto&& blah = stan::math::init_threadpool_tbb();
+
+static constexpr size_t num_chains = 4;
+class ServicesSampleHmcNutsDenseEAdaptParMatch : public testing::Test {
+ public:
+  ServicesSampleHmcNutsDenseEAdaptParMatch()
+      : model(std::make_unique<rosenbrock_model_namespace::rosenbrock_model>(
+            data_context, 0, &model_log)) {
+    for (int i = 0; i < num_chains; ++i) {
+      init.push_back(stan::test::unit::instrumented_writer{});
+      par_parameters.emplace_back(std::make_unique<std::stringstream>(), "#");
+      seq_parameters.emplace_back(std::make_unique<std::stringstream>(), "#");
+      diagnostic.push_back(stan::test::unit::instrumented_writer{});
+      context.push_back(std::make_shared<stan::io::empty_var_context>());
+    }
+  }
+  stan::io::empty_var_context data_context;
+  std::stringstream model_log;
+  stan::test::unit::instrumented_logger logger;
+  std::vector<stan::test::unit::instrumented_writer> init;
+  using str_writer = stan::callbacks::unique_stream_writer<std::stringstream>;
+  std::vector<str_writer> par_parameters;
+  std::vector<str_writer> seq_parameters;
+  std::vector<stan::test::unit::instrumented_writer> diagnostic;
+  std::vector<std::shared_ptr<stan::io::empty_var_context>> context;
+  std::unique_ptr<rosenbrock_model_namespace::rosenbrock_model> model;
+};
+
+/**
+ * This test checks that running multiple chains in one call
+ * with the same initial id is the same as running multiple calls
+ * with incrementing chain ids.
+ */
+TEST_F(ServicesSampleHmcNutsDenseEAdaptParMatch, single_multi_match) {
+  constexpr unsigned int random_seed = 0;
+  constexpr unsigned int chain = 0;
+  constexpr double init_radius = 0;
+  constexpr int num_warmup = 200;
+  constexpr int num_samples = 400;
+  constexpr int num_thin = 5;
+  constexpr bool save_warmup = true;
+  constexpr int refresh = 0;
+  constexpr double stepsize = 0.1;
+  constexpr double stepsize_jitter = 0;
+  constexpr int max_depth = 8;
+  constexpr double delta = .1;
+  constexpr double gamma = .1;
+  constexpr double kappa = .1;
+  constexpr double t0 = .1;
+  constexpr unsigned int init_buffer = 50;
+  constexpr unsigned int term_buffer = 50;
+  constexpr unsigned int window = 100;
+  stan::test::unit::instrumented_interrupt interrupt;
+  EXPECT_EQ(interrupt.call_count(), 0);
+  int return_code = stan::services::sample::hmc_nuts_dense_e_adapt(
+      *model, num_chains, context, random_seed, chain, init_radius, num_warmup,
+      num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+      max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+      interrupt, logger, init, par_parameters, diagnostic);
+
+  EXPECT_EQ(0, return_code);
+
+  int num_output_lines = (num_warmup + num_samples) / num_thin;
+  EXPECT_EQ((num_warmup + num_samples) * num_chains, interrupt.call_count());
+  for (int i = 0; i < num_chains; ++i) {
+    stan::test::unit::instrumented_writer seq_init;
+    stan::test::unit::instrumented_writer seq_diagnostic;
+    return_code = stan::services::sample::hmc_nuts_dense_e_adapt(
+        *model, *(context[i]), random_seed, i, init_radius, num_warmup,
+        num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+        max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+        interrupt, logger, seq_init, seq_parameters[i], seq_diagnostic);
+    EXPECT_EQ(0, return_code);
+  }
+  std::vector<Eigen::MatrixXd> par_res;
+  for (int i = 0; i < num_chains; ++i) {
+    auto par_str = par_parameters[i].get_stream().str();
+    auto sub_par_str = par_str.substr(par_str.find("Elements") - 1);
+    std::istringstream sub_par_stream(sub_par_str);
+    Eigen::MatrixXd par_mat
+        = stan::test::read_stan_sample_csv(sub_par_stream, 80, 9);
+    par_res.push_back(par_mat);
+  }
+  std::vector<Eigen::MatrixXd> seq_res;
+  for (int i = 0; i < num_chains; ++i) {
+    auto seq_str = seq_parameters[i].get_stream().str();
+    auto sub_seq_str = seq_str.substr(seq_str.find("Elements") - 1);
+    std::istringstream sub_seq_stream(sub_seq_str);
+    Eigen::MatrixXd seq_mat
+        = stan::test::read_stan_sample_csv(sub_seq_stream, 80, 9);
+    seq_res.push_back(seq_mat);
+  }
+  for (int i = 0; i < num_chains; ++i) {
+    Eigen::MatrixXd diff_res
+        = (par_res[i].array() - seq_res[i].array()).matrix();
+    EXPECT_MATRIX_EQ(diff_res, Eigen::MatrixXd::Zero(80, 9));
+  }
+}
diff --git a/src/test/unit/services/sample/hmc_nuts_dense_e_adapt_parallel_test.cpp b/src/test/unit/services/sample/hmc_nuts_dense_e_adapt_parallel_test.cpp
new file mode 100644
index 00000000000..1f8d2a117b9
--- /dev/null
+++ b/src/test/unit/services/sample/hmc_nuts_dense_e_adapt_parallel_test.cpp
@@ -0,0 +1,175 @@
+#include <stan/services/sample/hmc_nuts_dense_e_adapt.hpp>
+#include <gtest/gtest.h>
+#include <stan/io/empty_var_context.hpp>
+#include <test/test-models/good/optimization/rosenbrock.hpp>
+#include <test/unit/services/instrumented_callbacks.hpp>
+#include <iostream>
+
+auto&& blah = stan::math::init_threadpool_tbb();
+
+static constexpr size_t num_chains = 4;
+class ServicesSampleHmcNutsDenseEAdaptPar : public testing::Test {
+ public:
+  ServicesSampleHmcNutsDenseEAdaptPar() : model(data_context, 0, &model_log) {
+    for (int i = 0; i < num_chains; ++i) {
+      init.push_back(stan::test::unit::instrumented_writer{});
+      parameter.push_back(stan::test::unit::instrumented_writer{});
+      diagnostic.push_back(stan::test::unit::instrumented_writer{});
+      context.push_back(std::make_shared<stan::io::empty_var_context>());
+    }
+  }
+  stan::io::empty_var_context data_context;
+  std::stringstream model_log;
+  stan::test::unit::instrumented_logger logger;
+  std::vector<stan::test::unit::instrumented_writer> init;
+  std::vector<stan::test::unit::instrumented_writer> parameter;
+  std::vector<stan::test::unit::instrumented_writer> diagnostic;
+  std::vector<std::shared_ptr<stan::io::empty_var_context>> context;
+  stan_model model;
+};
+
+TEST_F(ServicesSampleHmcNutsDenseEAdaptPar, call_count) {
+  unsigned int random_seed = 0;
+  unsigned int chain = 1;
+  double init_radius = 0;
+  int num_warmup = 200;
+  int num_samples = 400;
+  int num_thin = 5;
+  bool save_warmup = true;
+  int refresh = 0;
+  double stepsize = 0.1;
+  double stepsize_jitter = 0;
+  int max_depth = 8;
+  double delta = .1;
+  double gamma = .1;
+  double kappa = .1;
+  double t0 = .1;
+  unsigned int init_buffer = 50;
+  unsigned int term_buffer = 50;
+  unsigned int window = 100;
+  stan::test::unit::instrumented_interrupt interrupt;
+  EXPECT_EQ(interrupt.call_count(), 0);
+
+  int return_code = stan::services::sample::hmc_nuts_dense_e_adapt(
+      model, num_chains, context, random_seed, chain, init_radius, num_warmup,
+      num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+      max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+      interrupt, logger, init, parameter, diagnostic);
+
+  EXPECT_EQ(0, return_code);
+
+  int num_output_lines = (num_warmup + num_samples) / num_thin;
+  EXPECT_EQ((num_warmup + num_samples) * num_chains, interrupt.call_count());
+  for (int i = 0; i < num_chains; ++i) {
+    EXPECT_EQ(1, parameter[i].call_count("vector_string"));
+    EXPECT_EQ(num_output_lines, parameter[i].call_count("vector_double"));
+    EXPECT_EQ(1, diagnostic[i].call_count("vector_string"));
+    EXPECT_EQ(num_output_lines, diagnostic[i].call_count("vector_double"));
+  }
+}
+
+TEST_F(ServicesSampleHmcNutsDenseEAdaptPar, parameter_checks) {
+  unsigned int random_seed = 0;
+  unsigned int chain = 1;
+  double init_radius = 0;
+  int num_warmup = 200;
+  int num_samples = 400;
+  int num_thin = 5;
+  bool save_warmup = true;
+  int refresh = 0;
+  double stepsize = 0.1;
+  double stepsize_jitter = 0;
+  int max_depth = 8;
+  double delta = .1;
+  double gamma = .1;
+  double kappa = .1;
+  double t0 = .1;
+  unsigned int init_buffer = 50;
+  unsigned int term_buffer = 50;
+  unsigned int window = 100;
+  stan::test::unit::instrumented_interrupt interrupt;
+  EXPECT_EQ(interrupt.call_count(), 0);
+
+  int return_code = stan::services::sample::hmc_nuts_dense_e_adapt(
+      model, num_chains, context, random_seed, chain, init_radius, num_warmup,
+      num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+      max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+      interrupt, logger, init, parameter, diagnostic);
+
+  for (size_t i = 0; i < num_chains; ++i) {
+    std::vector<std::vector<std::string>> parameter_names;
+    parameter_names = parameter[i].vector_string_values();
+    std::vector<std::vector<double>> parameter_values;
+    parameter_values = parameter[i].vector_double_values();
+    std::vector<std::vector<std::string>> diagnostic_names;
+    diagnostic_names = diagnostic[i].vector_string_values();
+    std::vector<std::vector<double>> diagnostic_values;
+    diagnostic_values = diagnostic[i].vector_double_values();
+
+    // Expectations of parameter parameter names.
+    ASSERT_EQ(9, parameter_names[0].size());
+    EXPECT_EQ("lp__", parameter_names[0][0]);
+    EXPECT_EQ("accept_stat__", parameter_names[0][1]);
+    EXPECT_EQ("stepsize__", parameter_names[0][2]);
+    EXPECT_EQ("treedepth__", parameter_names[0][3]);
+    EXPECT_EQ("n_leapfrog__", parameter_names[0][4]);
+    EXPECT_EQ("divergent__", parameter_names[0][5]);
+    EXPECT_EQ("energy__", parameter_names[0][6]);
+    EXPECT_EQ("x", parameter_names[0][7]);
+    EXPECT_EQ("y", parameter_names[0][8]);
+
+    // Expect one name per parameter value.
+    EXPECT_EQ(parameter_names[0].size(), parameter_values[0].size());
+    EXPECT_EQ(diagnostic_names[0].size(), diagnostic_values[0].size());
+
+    EXPECT_EQ((num_warmup + num_samples) / num_thin, parameter_values.size());
+
+    // Expect one call to set parameter names, and one set of output per
+    // iteration.
+    EXPECT_EQ("lp__", diagnostic_names[0][0]);
+    EXPECT_EQ("accept_stat__", diagnostic_names[0][1]);
+  }
+  EXPECT_EQ(return_code, 0);
+}
+
+TEST_F(ServicesSampleHmcNutsDenseEAdaptPar, output_regression) {
+  unsigned int random_seed = 0;
+  unsigned int chain = 1;
+  double init_radius = 0;
+  int num_warmup = 200;
+  int num_samples = 400;
+  int num_thin = 5;
+  bool save_warmup = true;
+  int refresh = 0;
+  double stepsize = 0.1;
+  double stepsize_jitter = 0;
+  int max_depth = 8;
+  double delta = .1;
+  double gamma = .1;
+  double kappa = .1;
+  double t0 = .1;
+  unsigned int init_buffer = 50;
+  unsigned int term_buffer = 50;
+  unsigned int window = 100;
+  stan::test::unit::instrumented_interrupt interrupt;
+  EXPECT_EQ(interrupt.call_count(), 0);
+
+  stan::services::sample::hmc_nuts_dense_e_adapt(
+      model, num_chains, context, random_seed, chain, init_radius, num_warmup,
+      num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+      max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+      interrupt, logger, init, parameter, diagnostic);
+
+  for (auto&& init_it : init) {
+    std::vector<std::string> init_values;
+    init_values = init_it.string_values();
+
+    EXPECT_EQ(0, init_values.size());
+  }
+
+  EXPECT_EQ(num_chains, logger.find_info("Elapsed Time:"));
+  EXPECT_EQ(num_chains, logger.find_info("seconds (Warm-up)"));
+  EXPECT_EQ(num_chains, logger.find_info("seconds (Sampling)"));
+  EXPECT_EQ(num_chains, logger.find_info("seconds (Total)"));
+  EXPECT_EQ(0, logger.call_count_error());
+}
diff --git a/src/test/unit/services/sample/hmc_nuts_diag_e_adapt_parallel_match_test.cpp b/src/test/unit/services/sample/hmc_nuts_diag_e_adapt_parallel_match_test.cpp
new file mode 100644
index 00000000000..c4167419cbf
--- /dev/null
+++ b/src/test/unit/services/sample/hmc_nuts_diag_e_adapt_parallel_match_test.cpp
@@ -0,0 +1,108 @@
+#include <stan/services/sample/hmc_nuts_diag_e_adapt.hpp>
+#include <stan/io/empty_var_context.hpp>
+#include <stan/callbacks/unique_stream_writer.hpp>
+#include <test/unit/util.hpp>
+#include <src/test/unit/services/util.hpp>
+#include <test/test-models/good/optimization/rosenbrock.hpp>
+#include <test/unit/services/instrumented_callbacks.hpp>
+#include <gtest/gtest.h>
+#include <iostream>
+
+auto&& blah = stan::math::init_threadpool_tbb();
+
+static constexpr size_t num_chains = 4;
+class ServicesSampleHmcNutsDiagEAdaptParMatch : public testing::Test {
+ public:
+  ServicesSampleHmcNutsDiagEAdaptParMatch()
+      : model(std::make_unique<rosenbrock_model_namespace::rosenbrock_model>(
+            data_context, 0, &model_log)) {
+    for (int i = 0; i < num_chains; ++i) {
+      init.push_back(stan::test::unit::instrumented_writer{});
+      par_parameters.emplace_back(std::make_unique<std::stringstream>(), "#");
+      seq_parameters.emplace_back(std::make_unique<std::stringstream>(), "#");
+      diagnostic.push_back(stan::test::unit::instrumented_writer{});
+      context.push_back(std::make_shared<stan::io::empty_var_context>());
+    }
+  }
+  stan::io::empty_var_context data_context;
+  std::stringstream model_log;
+  stan::test::unit::instrumented_logger logger;
+  std::vector<stan::test::unit::instrumented_writer> init;
+  using str_writer = stan::callbacks::unique_stream_writer<std::stringstream>;
+  std::vector<str_writer> par_parameters;
+  std::vector<str_writer> seq_parameters;
+  std::vector<stan::test::unit::instrumented_writer> diagnostic;
+  std::vector<std::shared_ptr<stan::io::empty_var_context>> context;
+  std::unique_ptr<rosenbrock_model_namespace::rosenbrock_model> model;
+};
+
+/**
+ * This test checks that running multiple chains in one call
+ * with the same initial id is the same as running multiple calls
+ * with incrementing chain ids.
+ */
+TEST_F(ServicesSampleHmcNutsDiagEAdaptParMatch, single_multi_match) {
+  constexpr unsigned int random_seed = 0;
+  constexpr unsigned int chain = 0;
+  constexpr double init_radius = 0;
+  constexpr int num_warmup = 200;
+  constexpr int num_samples = 400;
+  constexpr int num_thin = 5;
+  constexpr bool save_warmup = true;
+  constexpr int refresh = 0;
+  constexpr double stepsize = 0.1;
+  constexpr double stepsize_jitter = 0;
+  constexpr int max_depth = 8;
+  constexpr double delta = .1;
+  constexpr double gamma = .1;
+  constexpr double kappa = .1;
+  constexpr double t0 = .1;
+  constexpr unsigned int init_buffer = 50;
+  constexpr unsigned int term_buffer = 50;
+  constexpr unsigned int window = 100;
+  stan::test::unit::instrumented_interrupt interrupt;
+  EXPECT_EQ(interrupt.call_count(), 0);
+  int return_code = stan::services::sample::hmc_nuts_diag_e_adapt(
+      *model, num_chains, context, random_seed, chain, init_radius, num_warmup,
+      num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+      max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+      interrupt, logger, init, par_parameters, diagnostic);
+
+  EXPECT_EQ(0, return_code);
+
+  int num_output_lines = (num_warmup + num_samples) / num_thin;
+  EXPECT_EQ((num_warmup + num_samples) * num_chains, interrupt.call_count());
+  for (int i = 0; i < num_chains; ++i) {
+    stan::test::unit::instrumented_writer seq_init;
+    stan::test::unit::instrumented_writer seq_diagnostic;
+    return_code = stan::services::sample::hmc_nuts_diag_e_adapt(
+        *model, *(context[i]), random_seed, i, init_radius, num_warmup,
+        num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+        max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+        interrupt, logger, seq_init, seq_parameters[i], seq_diagnostic);
+    EXPECT_EQ(0, return_code);
+  }
+  std::vector<Eigen::MatrixXd> par_res;
+  for (int i = 0; i < num_chains; ++i) {
+    auto par_str = par_parameters[i].get_stream().str();
+    auto sub_par_str = par_str.substr(par_str.find("Diagonal") - 1);
+    std::istringstream sub_par_stream(sub_par_str);
+    Eigen::MatrixXd par_mat
+        = stan::test::read_stan_sample_csv(sub_par_stream, 80, 9);
+    par_res.push_back(par_mat);
+  }
+  std::vector<Eigen::MatrixXd> seq_res;
+  for (int i = 0; i < num_chains; ++i) {
+    auto seq_str = seq_parameters[i].get_stream().str();
+    auto sub_seq_str = seq_str.substr(seq_str.find("Diagonal") - 1);
+    std::istringstream sub_seq_stream(sub_seq_str);
+    Eigen::MatrixXd seq_mat
+        = stan::test::read_stan_sample_csv(sub_seq_stream, 80, 9);
+    seq_res.push_back(seq_mat);
+  }
+  for (int i = 0; i < num_chains; ++i) {
+    Eigen::MatrixXd diff_res
+        = (par_res[i].array() - seq_res[i].array()).matrix();
+    EXPECT_MATRIX_EQ(diff_res, Eigen::MatrixXd::Zero(80, 9));
+  }
+}
diff --git a/src/test/unit/services/sample/hmc_nuts_diag_e_adapt_parallel_test.cpp b/src/test/unit/services/sample/hmc_nuts_diag_e_adapt_parallel_test.cpp
new file mode 100644
index 00000000000..453b03adcdd
--- /dev/null
+++ b/src/test/unit/services/sample/hmc_nuts_diag_e_adapt_parallel_test.cpp
@@ -0,0 +1,175 @@
+#include <stan/services/sample/hmc_nuts_diag_e_adapt.hpp>
+#include <gtest/gtest.h>
+#include <stan/io/empty_var_context.hpp>
+#include <test/test-models/good/optimization/rosenbrock.hpp>
+#include <test/unit/services/instrumented_callbacks.hpp>
+#include <iostream>
+
+auto&& blah = stan::math::init_threadpool_tbb();
+
+static constexpr size_t num_chains = 4;
+class ServicesSampleHmcNutsDiagEAdaptPar : public testing::Test {
+ public:
+  ServicesSampleHmcNutsDiagEAdaptPar() : model(data_context, 0, &model_log) {
+    for (int i = 0; i < num_chains; ++i) {
+      init.push_back(stan::test::unit::instrumented_writer{});
+      parameter.push_back(stan::test::unit::instrumented_writer{});
+      diagnostic.push_back(stan::test::unit::instrumented_writer{});
+      context.push_back(std::make_shared<stan::io::empty_var_context>());
+    }
+  }
+  stan::io::empty_var_context data_context;
+  std::stringstream model_log;
+  stan::test::unit::instrumented_logger logger;
+  std::vector<stan::test::unit::instrumented_writer> init;
+  std::vector<stan::test::unit::instrumented_writer> parameter;
+  std::vector<stan::test::unit::instrumented_writer> diagnostic;
+  std::vector<std::shared_ptr<stan::io::empty_var_context>> context;
+  stan_model model;
+};
+
+TEST_F(ServicesSampleHmcNutsDiagEAdaptPar, call_count) {
+  unsigned int random_seed = 0;
+  unsigned int chain = 1;
+  double init_radius = 0;
+  int num_warmup = 200;
+  int num_samples = 400;
+  int num_thin = 5;
+  bool save_warmup = true;
+  int refresh = 0;
+  double stepsize = 0.1;
+  double stepsize_jitter = 0;
+  int max_depth = 8;
+  double delta = .1;
+  double gamma = .1;
+  double kappa = .1;
+  double t0 = .1;
+  unsigned int init_buffer = 50;
+  unsigned int term_buffer = 50;
+  unsigned int window = 100;
+  stan::test::unit::instrumented_interrupt interrupt;
+  EXPECT_EQ(interrupt.call_count(), 0);
+
+  int return_code = stan::services::sample::hmc_nuts_diag_e_adapt(
+      model, num_chains, context, random_seed, chain, init_radius, num_warmup,
+      num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+      max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+      interrupt, logger, init, parameter, diagnostic);
+
+  EXPECT_EQ(0, return_code);
+
+  int num_output_lines = (num_warmup + num_samples) / num_thin;
+  EXPECT_EQ((num_warmup + num_samples) * num_chains, interrupt.call_count());
+  for (int i = 0; i < num_chains; ++i) {
+    EXPECT_EQ(1, parameter[i].call_count("vector_string"));
+    EXPECT_EQ(num_output_lines, parameter[i].call_count("vector_double"));
+    EXPECT_EQ(1, diagnostic[i].call_count("vector_string"));
+    EXPECT_EQ(num_output_lines, diagnostic[i].call_count("vector_double"));
+  }
+}
+
+TEST_F(ServicesSampleHmcNutsDiagEAdaptPar, parameter_checks) {
+  unsigned int random_seed = 0;
+  unsigned int chain = 1;
+  double init_radius = 0;
+  int num_warmup = 200;
+  int num_samples = 400;
+  int num_thin = 5;
+  bool save_warmup = true;
+  int refresh = 0;
+  double stepsize = 0.1;
+  double stepsize_jitter = 0;
+  int max_depth = 8;
+  double delta = .1;
+  double gamma = .1;
+  double kappa = .1;
+  double t0 = .1;
+  unsigned int init_buffer = 50;
+  unsigned int term_buffer = 50;
+  unsigned int window = 100;
+  stan::test::unit::instrumented_interrupt interrupt;
+  EXPECT_EQ(interrupt.call_count(), 0);
+
+  int return_code = stan::services::sample::hmc_nuts_diag_e_adapt(
+      model, num_chains, context, random_seed, chain, init_radius, num_warmup,
+      num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+      max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+      interrupt, logger, init, parameter, diagnostic);
+
+  for (size_t i = 0; i < num_chains; ++i) {
+    std::vector<std::vector<std::string>> parameter_names;
+    parameter_names = parameter[i].vector_string_values();
+    std::vector<std::vector<double>> parameter_values;
+    parameter_values = parameter[i].vector_double_values();
+    std::vector<std::vector<std::string>> diagnostic_names;
+    diagnostic_names = diagnostic[i].vector_string_values();
+    std::vector<std::vector<double>> diagnostic_values;
+    diagnostic_values = diagnostic[i].vector_double_values();
+
+    // Expectations of parameter parameter names.
+    ASSERT_EQ(9, parameter_names[0].size());
+    EXPECT_EQ("lp__", parameter_names[0][0]);
+    EXPECT_EQ("accept_stat__", parameter_names[0][1]);
+    EXPECT_EQ("stepsize__", parameter_names[0][2]);
+    EXPECT_EQ("treedepth__", parameter_names[0][3]);
+    EXPECT_EQ("n_leapfrog__", parameter_names[0][4]);
+    EXPECT_EQ("divergent__", parameter_names[0][5]);
+    EXPECT_EQ("energy__", parameter_names[0][6]);
+    EXPECT_EQ("x", parameter_names[0][7]);
+    EXPECT_EQ("y", parameter_names[0][8]);
+
+    // Expect one name per parameter value.
+    EXPECT_EQ(parameter_names[0].size(), parameter_values[0].size());
+    EXPECT_EQ(diagnostic_names[0].size(), diagnostic_values[0].size());
+
+    EXPECT_EQ((num_warmup + num_samples) / num_thin, parameter_values.size());
+
+    // Expect one call to set parameter names, and one set of output per
+    // iteration.
+    EXPECT_EQ("lp__", diagnostic_names[0][0]);
+    EXPECT_EQ("accept_stat__", diagnostic_names[0][1]);
+  }
+  EXPECT_EQ(return_code, 0);
+}
+
+TEST_F(ServicesSampleHmcNutsDiagEAdaptPar, output_regression) {
+  unsigned int random_seed = 0;
+  unsigned int chain = 1;
+  double init_radius = 0;
+  int num_warmup = 200;
+  int num_samples = 400;
+  int num_thin = 5;
+  bool save_warmup = true;
+  int refresh = 0;
+  double stepsize = 0.1;
+  double stepsize_jitter = 0;
+  int max_depth = 8;
+  double delta = .1;
+  double gamma = .1;
+  double kappa = .1;
+  double t0 = .1;
+  unsigned int init_buffer = 50;
+  unsigned int term_buffer = 50;
+  unsigned int window = 100;
+  stan::test::unit::instrumented_interrupt interrupt;
+  EXPECT_EQ(interrupt.call_count(), 0);
+
+  stan::services::sample::hmc_nuts_diag_e_adapt(
+      model, num_chains, context, random_seed, chain, init_radius, num_warmup,
+      num_samples, num_thin, save_warmup, refresh, stepsize, stepsize_jitter,
+      max_depth, delta, gamma, kappa, t0, init_buffer, term_buffer, window,
+      interrupt, logger, init, parameter, diagnostic);
+
+  for (auto&& init_it : init) {
+    std::vector<std::string> init_values;
+    init_values = init_it.string_values();
+
+    EXPECT_EQ(0, init_values.size());
+  }
+
+  EXPECT_EQ(num_chains, logger.find_info("Elapsed Time:"));
+  EXPECT_EQ(num_chains, logger.find_info("seconds (Warm-up)"));
+  EXPECT_EQ(num_chains, logger.find_info("seconds (Sampling)"));
+  EXPECT_EQ(num_chains, logger.find_info("seconds (Total)"));
+  EXPECT_EQ(0, logger.call_count_error());
+}
diff --git a/src/test/unit/services/util.hpp b/src/test/unit/services/util.hpp
new file mode 100644
index 00000000000..f283fa6f49f
--- /dev/null
+++ b/src/test/unit/services/util.hpp
@@ -0,0 +1,45 @@
+#ifndef STAN_SRC_TEST_UNIT_SERVICES_UTIL_HPP
+#define STAN_SRC_TEST_UNIT_SERVICES_UTIL_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <test/unit/util.hpp>
+#include <iostream>
+#include <string>
+
+namespace stan {
+namespace test {
+/**
+ * Read a CSV into an Eigen matrix.
+ * @param in An input string stream holding the CSV
+ * @param rows Number of rows
+ * @param cols Number of columns.
+ */
+Eigen::MatrixXd read_stan_sample_csv(std::istringstream& in, int rows,
+                                     int cols) {
+  std::string line;
+  int row = 0;
+  int col = 0;
+  Eigen::MatrixXd res = Eigen::MatrixXd(rows, cols);
+  while (std::getline(in, line)) {
+    if (line.find("#") != std::string::npos) {
+      continue;
+    }
+    const char* ptr = line.c_str();
+    int len = line.length();
+    col = 0;
+
+    const char* start = ptr;
+    for (int i = 0; i < len; i++) {
+      if (ptr[i] == ',') {
+        res(row, col++) = atof(start);
+        start = ptr + i + 1;
+      }
+    }
+    res(row, col) = atof(start);
+    row++;
+  }
+  return res;
+}
+}  // namespace test
+}  // namespace stan
+#endif