diff --git a/docker/Dockerfile b/docker/Dockerfile
index a097df9b..8b1ffa94 100644
--- a/docker/Dockerfile
+++ b/docker/Dockerfile
@@ -14,7 +14,7 @@ RUN sudo apt-get install -yq libopenblas-dev libmpich-dev libblas-dev liblapack-
 RUN sudo apt-get install -yq libboost-all-dev
 RUN sudo apt-get install -yq vim
 RUN sudo apt-get install -yq git-lfs
-RUN sudo apt-get install -yq valgrind
+RUN sudo apt-get install -yq valgrind hdf5-tools
 RUN sudo apt-get install -yq wget
 ### clang-format seems to be updated to 14.0. Not using it for now.
 # RUN sudo apt-get install -yq clang-format
diff --git a/examples/Carbyne/carbyne.rom.cfg b/examples/Carbyne/carbyne.rom.cfg
new file mode 100644
index 00000000..cb0cd295
--- /dev/null
+++ b/examples/Carbyne/carbyne.rom.cfg
@@ -0,0 +1,34 @@
+verbosity=2
+xcFunctional=PBE
+FDtype=4th
+[Mesh]
+nx= 96
+ny= 96
+nz= 192
+[Domain]
+ox= -10.
+oy= -10.
+oz= -20.
+lx= 20.
+ly= 20.
+lz= 40.
+[Potentials]
+pseudopotential=pseudo.H_ONCV_PBE_SG15
+pseudopotential=pseudo.C_ONCV_PBE_SG15
+[Run]
+type=QUENCH
+[Quench]
+max_steps=5
+atol=1.e-8
+[Orbitals]
+initial_type=Fourier
+[Restart]
+output_level=4
+input_level=4
+input_filename=snapshot0_000
+
+[ROM.offline]
+restart_filefmt=snapshot0_%03d
+restart_min_idx=0
+restart_max_idx=1
+basis_file=carom
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 1bad52e7..e33d03d2 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -162,6 +162,11 @@ set(SOURCES
  rom.cc
 )
 
+if(USE_LIBROM)
+  list(APPEND SOURCES
+    rom_workflows.cc)
+endif(USE_LIBROM)
+
 add_library(mgmol_src ${SOURCES})
 
 target_include_directories(mgmol_src PRIVATE ${HDF5_INCLUDE_DIRS})
diff --git a/src/rom_main.cc b/src/rom_main.cc
index c7ec278e..dc877a58 100644
--- a/src/rom_main.cc
+++ b/src/rom_main.cc
@@ -23,76 +23,115 @@
 //   Potentials, eigenvalues and operators in Rydberg
 //   Energies in Hartree
 //
-#include <cassert>
-#include <iostream>
-#include <iterator>
-#include <vector>
-using namespace std;
-
-#ifdef _OPENMP
-#include <omp.h>
-#endif
-
-#ifdef USE_CNR
-#include <mkl.h>
-#endif
-
-#include <mpi.h>
-
-#include "Control.h"
-#include "DistMatrix.h"
-#include "ExtendedGridOrbitals.h"
-#include "LocGridOrbitals.h"
-#include "MGmol.h"
-#include "MGmol_MPI.h"
-#include "MPIdata.h"
-#include "MatricesBlacsContext.h"
-#include "Mesh.h"
-#include "PackedCommunicationBuffer.h"
-#include "ReplicatedWorkSpace.h"
-#include "SparseDistMatrix.h"
-#include "magma_singleton.h"
-#include "tools.h"
-
-#include <fenv.h>
-#include <sys/cdefs.h>
-#include <time.h>
-
-#include <boost/program_options.hpp>
-namespace po = boost::program_options;
-
-#include "librom.h"
-
-//#include "MemTrack.h"
+#include "rom_workflows.h"
 
-int main(int argc, char** argv)
+// A helper function
+template <class T>
+std::ostream& operator<<(std::ostream& os, const std::vector<T>& v)
 {
-    // change handling of memory allocation errors
-    set_new_handler(noMoreMemory);
-
-    cout.sync_with_stdio();
+    copy(v.begin(), v.end(), std::ostream_iterator<T>(std::cout, " "));
+    return os;
+}
 
+int main(int argc, char** argv)
+{
     int mpirc = MPI_Init(&argc, &argv);
     if (mpirc != MPI_SUCCESS)
     {
-        cerr << "MPI Initialization failed!!!" << endl;
+        std::cerr << "MPI Initialization failed!!!" << std::endl;
         MPI_Abort(MPI_COMM_WORLD, 0);
     }
-    MPI_Comm_rank(MPI_COMM_WORLD, &mype);
-    assert(mype > -1);
-    onpe0 = (mype == 0);
 
-    CAROM::Vector librom_vector(10, false);
+    MPI_Comm comm = MPI_COMM_WORLD;
+
+    mgmol_init(comm);
+
+    // read runtime parameters
+    std::string input_filename("");
+    std::string lrs_filename;
+    std::string constraints_filename("");
+    bool tcheck = false;
+
+    float total_spin = 0.;
+    bool with_spin   = false;
+
+    po::variables_map vm;
+
+    // use configure file if it can be found
+    // std::string config_filename("mgmol.cfg");
+
+    // read options from PE0 only
+    if (MPIdata::onpe0)
+    {
+        read_config(argc, argv, vm, input_filename, lrs_filename,
+            constraints_filename, total_spin, with_spin, tcheck);
+    }
+
+    MGmol_MPI::setup(comm, std::cout, with_spin);
+    MGmol_MPI& mmpi      = *(MGmol_MPI::instance());
+    MPI_Comm global_comm = mmpi.commGlobal();
+
+    Control::setup(global_comm, with_spin, total_spin);
+    Control& ct = *(Control::instance());
+
+    ct.setOptions(vm);
+
+    int ret = ct.checkOptions();
+    if (ret < 0) return ret;
+
+    mmpi.bcastGlobal(input_filename);
+    mmpi.bcastGlobal(lrs_filename);
+
+    // Enter main scope
+    {
+        MGmolInterface* mgmol;
+        if (ct.isLocMode())
+            mgmol = new MGmol<LocGridOrbitals>(global_comm, *MPIdata::sout);
+        else
+            mgmol
+                = new MGmol<ExtendedGridOrbitals>(global_comm, *MPIdata::sout);
+
+        unsigned ngpts[3]    = { ct.ngpts_[0], ct.ngpts_[1], ct.ngpts_[2] };
+        double origin[3]     = { ct.ox_, ct.oy_, ct.oz_ };
+        const double cell[3] = { ct.lx_, ct.ly_, ct.lz_ };
+        Mesh::setup(mmpi.commSpin(), ngpts, origin, cell, ct.lap_type);
+
+        mgmol->setupFromInput(input_filename);
+
+        if (ct.isLocMode() || ct.init_loc == 1) mgmol->setupLRs(lrs_filename);
+
+        mgmol->setupConstraintsFromInput(constraints_filename);
+
+        mgmol_setup();
+
+        if (!tcheck)
+        {
+            mgmol->setup();
+
+            if (ct.isLocMode())
+                readRestartFiles<LocGridOrbitals>(mgmol);
+            else
+                readRestartFiles<ExtendedGridOrbitals>(mgmol);
+        }
+        else
+        {
+            *MPIdata::sout << " Input parameters OK\n";
+        }
+        delete mgmol;
+
+    } // close main scope
+
+    mgmol_finalize();
 
     mpirc = MPI_Finalize();
     if (mpirc != MPI_SUCCESS)
     {
-        cerr << "MPI Finalize failed!!!" << endl;
+        std::cerr << "MPI Finalize failed!!!" << std::endl;
     }
 
     time_t tt;
     time(&tt);
-    if (onpe0) cout << " Run ended at " << ctime(&tt) << endl;
+    if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl;
 
     // MemTrack::TrackDumpBlocks();
 
diff --git a/src/rom_workflows.cc b/src/rom_workflows.cc
new file mode 100644
index 00000000..4594dd3a
--- /dev/null
+++ b/src/rom_workflows.cc
@@ -0,0 +1,71 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge
+// National Laboratory.
+// LLNL-CODE-743438
+// All rights reserved.
+// This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
+// Please also read this link https://github.com/llnl/mgmol/LICENSE
+
+#include "rom_workflows.h"
+#include <memory>
+#include <string>
+#include <stdexcept>
+
+template<typename ... Args>
+std::string string_format( const std::string& format, Args ... args )
+{
+    int size_s = std::snprintf( nullptr, 0, format.c_str(), args ... ) + 1; // Extra space for '\0'
+    if( size_s <= 0 ){ throw std::runtime_error( "Error during formatting." ); }
+    auto size = static_cast<size_t>( size_s );
+    std::unique_ptr<char[]> buf( new char[ size ] );
+    std::snprintf( buf.get(), size, format.c_str(), args ... );
+    return std::string( buf.get(), buf.get() + size - 1 ); // We don't want the '\0' inside
+}
+
+template <class OrbitalsType>
+void readRestartFiles(MGmolInterface *mgmol_)
+{
+    Control& ct              = *(Control::instance());
+    ROMPrivateOptions rom_options = ct.getROMOptions();
+    assert(rom_options.restart_file_minidx >= 0);
+    assert(rom_options.restart_file_maxidx >= 0);
+    const int minidx = rom_options.restart_file_minidx;
+    const int maxidx = rom_options.restart_file_maxidx;
+    const int num_restart = maxidx - minidx + 1;
+
+    MGmol<OrbitalsType> *mgmol = static_cast<MGmol<OrbitalsType> *>(mgmol_);
+    OrbitalsType *orbitals = nullptr;
+    std::string filename;
+
+    /* Read the first snapshot to determin dimension and number of snapshots */
+    filename = string_format(rom_options.restart_file_fmt, minidx);
+    orbitals = mgmol->loadOrbitalFromRestartFile(filename);
+    const int dim = orbitals->getLocNumpt();
+    const int chrom_num = orbitals->chromatic_number();
+    const int totalSamples = orbitals->chromatic_number() * num_restart;
+    delete orbitals;
+
+    /* Initialize libROM classes */
+    CAROM::Options svd_options(dim, totalSamples, 1);
+    CAROM::BasisGenerator basis_generator(svd_options, false, rom_options.basis_file);
+    
+    /* Collect the restart files */
+    for (int k = minidx; k <= maxidx; k++)
+    {
+        filename = string_format(rom_options.restart_file_fmt, k);
+        orbitals = mgmol->loadOrbitalFromRestartFile(filename);
+        assert(dim == orbitals->getLocNumpt());
+        assert(chrom_num == orbitals->chromatic_number());
+
+        for (int i = 0; i < chrom_num; ++i)
+            basis_generator.takeSample(orbitals->getPsi(i));
+
+        delete orbitals;
+    }
+    basis_generator.writeSnapshot();
+    basis_generator.endSamples();
+}
+
+template void readRestartFiles<LocGridOrbitals>(MGmolInterface *mgmol_);
+template void readRestartFiles<ExtendedGridOrbitals>(MGmolInterface *mgmol_);
\ No newline at end of file
diff --git a/src/rom_workflows.h b/src/rom_workflows.h
new file mode 100644
index 00000000..dfb22b79
--- /dev/null
+++ b/src/rom_workflows.h
@@ -0,0 +1,41 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge
+// National Laboratory.
+// LLNL-CODE-743438
+// All rights reserved.
+// This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
+// Please also read this link https://github.com/llnl/mgmol/LICENSE
+
+#ifndef ROM_WORKFLOWS_H
+#define ROM_WORKFLOWS_H
+
+#include "Control.h"
+#include "ExtendedGridOrbitals.h"
+#include "LocGridOrbitals.h"
+#include "MGmol.h"
+#include "MGmol_MPI.h"
+#include "MPIdata.h"
+#include "Mesh.h"
+#include "mgmol_run.h"
+#include "tools.h"
+
+#include <cassert>
+#include <fenv.h>
+#include <iostream>
+#include <iterator>
+#include <sys/cdefs.h>
+#include <time.h>
+#include <vector>
+
+#include <mpi.h>
+
+#include <boost/program_options.hpp>
+namespace po = boost::program_options;
+
+#include "librom.h"
+
+template <class OrbitalsType>
+void readRestartFiles(MGmolInterface *mgmol_);
+
+#endif  // ROM_WORKFLOWS_H