diff --git a/examples/PinnedH2O/get_result.sh b/examples/PinnedH2O/get_result.sh
new file mode 100644
index 00000000..443e4cdd
--- /dev/null
+++ b/examples/PinnedH2O/get_result.sh
@@ -0,0 +1,27 @@
+filename="offline_PinnedH2O.out" # FOM
+#filename="rom39_PinnedH2O.out" # compare MD
+#filename="39_force_PinnedH2O.out" # compare force
+
+# Extracting H1, H2, F1, F2 from MGmgol output log
+# if FOM, these files contain the FOM results
+# if compare MD, these files contain the results with projected orbitals
+awk '/H1 / {print $3, $4, $5}' $filename > H1_$filename
+awk '/H2 / {print $3, $4, $5}' $filename > H2_$filename
+awk '/F1 / {print $6, $7, $8}' $filename > F1_$filename
+awk '/F2 / {print $6, $7, $8}' $filename > F2_$filename
+
+# if compare force, files with "_fom" contain the FOM results
+# files with "_rom" contain the results with projected orbitals
+if [[ "$filename" == *"force_"* ]]; then
+    sed -n '1~2p' H1_$filename > H1_rom$filename
+    sed -n '1~2p' H2_$filename > H2_rom$filename
+    sed -n '1~2p' F1_$filename > F1_rom$filename
+    sed -n '1~2p' F2_$filename > F2_rom$filename
+
+    sed -n '2~2p' H1_$filename > H1_fom$filename
+    sed -n '2~2p' H2_$filename > H2_fom$filename
+    sed -n '2~2p' F1_$filename > F1_fom$filename
+    sed -n '2~2p' F2_$filename > F2_fom$filename
+fi
+
+rm -rf snapshot0_*
diff --git a/examples/PinnedH2O/job.basis b/examples/PinnedH2O/job.basis_1_50
similarity index 73%
rename from examples/PinnedH2O/job.basis
rename to examples/PinnedH2O/job.basis_1_50
index 8f8c79f6..a5069570 100644
--- a/examples/PinnedH2O/job.basis
+++ b/examples/PinnedH2O/job.basis_1_50
@@ -1,7 +1,7 @@
 #!/bin/tcsh
 #SBATCH -N 2
 #SBATCH -t 0:10:00
-#SBATCH -p pdebug
+#SBATCH -p pbatch
 
 date
 
@@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 64
 
-set maindir = /p/lustre2/cheung26/mgmol-20240815
+set maindir = /p/lustre2/cheung26/mgmol-20241012
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
@@ -18,15 +18,15 @@ set exe     = ${maindir}/build_quartz/libROM/build/examples/misc/combine_samples
 
 set snapshot_files = ""
 set increment_md_steps = 1
-set num_md_steps = 500
+set num_md_steps = 50
 
 foreach k (`seq $increment_md_steps $increment_md_steps $num_md_steps`)
     set snapshot_files = "$snapshot_files MD_mdstep${k}_snapshot"
 end
 echo "Snapshot files: $snapshot_files"
 
-set basis_file = "PinnedH2O_orbitals_basis"
+set basis_file = "PinnedH2O_orbitals_basis_${increment_md_steps}_${num_md_steps}"
 
-srun -n $ncpus $exe -f $basis_file $snapshot_files > basis_${increment_md_steps}_${num_md_steps}_Pinned_H2O.out
+srun -n $ncpus $exe -f $basis_file $snapshot_files > basis_${increment_md_steps}_${num_md_steps}_PinnedH2O.out
 
 date
diff --git a/examples/PinnedH2O/job.rom b/examples/PinnedH2O/job.offline
similarity index 73%
rename from examples/PinnedH2O/job.rom
rename to examples/PinnedH2O/job.offline
index 854c56d8..0f8ef90e 100644
--- a/examples/PinnedH2O/job.rom
+++ b/examples/PinnedH2O/job.offline
@@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 64
 
-set maindir = /p/lustre2/cheung26/mgmol-20240815
+set maindir = /p/lustre2/cheung26/mgmol-20241012
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
@@ -20,8 +20,8 @@ cp $maindir/install_quartz/bin/$exe .
 
 set datadir = $maindir/examples/PinnedH2O
 
-set cfg_rom = mgmol_rom.cfg
-#cp $datadir/$cfg_rom .
+set cfg_offline = mgmol_offline.cfg
+cp $datadir/$cfg_offline .
 
 cp $datadir/coords.in .
 
@@ -30,6 +30,6 @@ ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 .
 
 source $maindir/scripts/modules.quartz
 
-srun -n $ncpus $exe -c $cfg_rom -i coords.in > rom_PinnedH2O.out
+srun -n $ncpus $exe -c $cfg_offline -i coords.in > offline_PinnedH2O.out
 
 date
diff --git a/examples/PinnedH2O/job.rom_1_50 b/examples/PinnedH2O/job.rom_1_50
new file mode 100644
index 00000000..aa6df335
--- /dev/null
+++ b/examples/PinnedH2O/job.rom_1_50
@@ -0,0 +1,39 @@
+#!/bin/tcsh
+#SBATCH -N 2
+#SBATCH -t 1:00:00
+#SBATCH -p pbatch
+
+date
+
+setenv OMP_NUM_THREADS 1
+#setenv KMP_DETERMINISTIC_REDUCTION 1
+
+set ncpus = 64
+
+set maindir = /p/lustre2/cheung26/mgmol-20241012
+
+setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
+
+set exe     = mgmol-opt
+
+cp $maindir/install_quartz/bin/$exe .
+
+set datadir = $maindir/examples/PinnedH2O
+
+set increment_md_steps = 1
+set num_md_steps = 50
+set basis_file = PinnedH2O_orbitals_basis_${increment_md_steps}_${num_md_steps}
+
+set cfg_rom = mgmol_rom_${increment_md_steps}_${num_md_steps}.cfg
+cp $datadir/$cfg_rom .
+
+cp $datadir/coords.in .
+
+ln -s -f $maindir/potentials/pseudo.O_ONCV_PBE_SG15 .
+ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 .
+
+source $maindir/scripts/modules.quartz
+
+srun -n $ncpus $exe -c $cfg_rom -i coords.in > rom_${increment_md_steps}_${num_md_steps}_PinnedH2O.out
+
+date
diff --git a/examples/PinnedH2O/mgmol_rom.cfg b/examples/PinnedH2O/mgmol_offline.cfg
similarity index 100%
rename from examples/PinnedH2O/mgmol_rom.cfg
rename to examples/PinnedH2O/mgmol_offline.cfg
diff --git a/examples/PinnedH2O/mgmol_rom_1_50.cfg b/examples/PinnedH2O/mgmol_rom_1_50.cfg
new file mode 100644
index 00000000..337a86eb
--- /dev/null
+++ b/examples/PinnedH2O/mgmol_rom_1_50.cfg
@@ -0,0 +1,40 @@
+verbosity=1
+xcFunctional=PBE
+FDtype=Mehrstellen
+[Mesh]
+nx=64
+ny=64
+nz=64
+[Domain]
+ox=-6.
+oy=-6.
+oz=-6.
+lx=12.
+ly=12.
+lz=12.
+[Potentials]
+pseudopotential=pseudo.O_ONCV_PBE_SG15
+pseudopotential=pseudo.H_ONCV_PBE_SG15
+[Run]
+type=MD
+[MD]
+num_steps=500
+dt=40.
+thermostat=ON
+[Thermostat]
+type=Berendsen
+temperature=1000.
+relax_time=800.
+[Quench]
+max_steps=100
+atol=1.e-8
+[Orbitals]
+initial_type=Random
+initial_width=2.
+[Restart]
+output_level=4
+[ROM.offline]
+basis_file=PinnedH2O_orbitals_basis_1_50
+[ROM.basis]
+compare_md=false
+number_of_orbital_basis=39
diff --git a/examples/PinnedH2O/plot_PinnedH2O_force.m b/examples/PinnedH2O/plot_PinnedH2O_force.m
new file mode 100644
index 00000000..a41457a5
--- /dev/null
+++ b/examples/PinnedH2O/plot_PinnedH2O_force.m
@@ -0,0 +1,84 @@
+clc; clear all; close all;
+
+%%
+plot_fom = 0;
+plot_rom = 0;
+rdim = 77;
+
+%%
+load F_fom.mat
+fprintf(1, 'Force statistics using FOM orbitals\n');
+fprintf(1, 'Mean of force on H1: %6.4e, %6.4e, %6.4e\n', mean(F1_fom));
+fprintf(1, 'Variance of force on H1: %6.4e, %6.4e, %6.4e\n', var(F1_fom));
+fprintf(1, 'Mean of force on H2: %6.4e, %6.4e, %6.4e\n', mean(F2_fom));
+fprintf(1, 'Variance of force on H2: %6.4e, %6.4e, %6.4e\n', var(F2_fom));
+
+if plot_fom
+    plotForce(F1_fom, 'F_H1_fom');
+    plotForce(F2_fom, 'F_H2_fom');
+    plotForceHistograms(F1_fom, 'H1_fom');
+    plotForceHistograms(F2_fom, 'H2_fom');
+end
+
+%%
+load(['F_rom' int2str(rdim) '.mat'])
+fprintf(1, 'Force statistics using projected orbitals\n');
+fprintf(1, 'Mean of force on H1: %6.4e, %6.4e, %6.4e\n', mean(F1_rom));
+fprintf(1, 'Variance of force on H1: %6.4e, %6.4e, %6.4e\n', var(F1_rom));
+%H1_correlation = sum(F1_fom(:) .* F1_rom(:)) / (norm(F1_fom(:)) * norm(F1_rom(:)))
+fprintf(1, 'Mean of force on H2: %6.4e, %6.4e, %6.4e\n', mean(F2_rom));
+fprintf(1, 'Variance of force on H2: %6.4e, %6.4e, %6.4e\n', var(F2_rom));
+%H2_correlation = sum(F2_fom(:) .* F2_rom(:)) / (norm(F2_fom(:)) * norm(F2_rom(:)))
+
+if plot_rom
+    plotForce(F1_rom, ['F_H1_rom' int2str(rdim)]);
+    plotForce(F2_rom, ['F_H2_rom' int2str(rdim)]);
+    plotForceHistograms(F1_rom, ['H1_rom' int2str(rdim)]);
+    plotForceHistograms(F2_rom, ['H2_rom' int2str(rdim)]);
+    plotForceHistogram(abs(F1_fom - F1_rom), ['H1_rom' int2str(rdim)], 'Fdiff');
+    plotForceHistogram(abs(F2_fom - F2_rom), ['H2_rom' int2str(rdim)], 'Fdiff');
+end
+
+%% 
+function plotForce(F, suffix)
+    figure;
+    imagesc(F');
+    axis tight;
+    axis equal;
+    colorbar;
+    saveas(gcf, suffix, 'jpg');
+end
+
+function plotForceHistogram(F, suffix, var)
+    figure;
+    if strcmp(var,'Fx')
+        X = F(:,1);
+        var_name = 'x-directional Force';
+    elseif strcmp(var,'Fy')
+        X = F(:,2);
+        var_name = 'y-directional Force';
+    elseif strcmp(var,'Fz')
+        X = F(:,3);
+        var_name = 'z-directional Force';
+    elseif strcmp(var,'Fmag')
+        X = sqrt(sum(F.^2, 2));
+        var_name = 'Force Magitude';
+    elseif strcmp(var,'Fdiff')
+        X = sqrt(sum(F.^2, 2));
+        var_name = 'Magitude of Difference in Force';
+    else
+        error('Invalid type');
+    end
+    histogram(X, 20);
+    xlabel(var_name);
+    ylabel('Frequency');
+    title(['Histogram of ' var_name]);
+    saveas(gcf, [var '_' suffix], 'jpg');
+end
+
+function plotForceHistograms(F, suffix)
+    plotForceHistogram(F, suffix, 'Fx');
+    plotForceHistogram(F, suffix, 'Fy');
+    plotForceHistogram(F, suffix, 'Fz');
+    plotForceHistogram(F, suffix, 'Fmag');
+end
diff --git a/examples/PinnedH2O/plot_PinnedH2O_md.m b/examples/PinnedH2O/plot_PinnedH2O_md.m
new file mode 100644
index 00000000..47c309b2
--- /dev/null
+++ b/examples/PinnedH2O/plot_PinnedH2O_md.m
@@ -0,0 +1,59 @@
+clc; clear all; close all;
+
+%%
+plot_fom = 0;
+rdims = [77, 39];
+
+%% 
+load md_fom.mat
+if plot_fom
+    plotAngleHistogram(H1_fom, H2_fom, 'fom');
+end
+
+%%
+
+all_H1_rom = zeros(length(rdims), size(H1_fom, 1), 3);
+all_H2_rom = zeros(length(rdims), size(H2_fom, 1), 3);
+k = 0;
+
+for rdim = rdims
+    k = k + 1;
+    load(['md_rom' int2str(rdim) '.mat'])
+    plotAngleHistogram(H1_rom, H2_rom, ['rom' int2str(rdim)]);
+    all_H1_rom(k, :, :) = H1_rom;
+    all_H2_rom(k, :, :) = H2_rom;
+end
+
+plotAtomTrajectory(H1_fom(:,1), all_H1_rom(:,:,1), rdims, 'x', 1)
+plotAtomTrajectory(H1_fom(:,2), all_H1_rom(:,:,2), rdims, 'y', 1)
+plotAtomTrajectory(H1_fom(:,3), all_H1_rom(:,:,3), rdims, 'z', 1)
+
+plotAtomTrajectory(H2_fom(:,1), all_H2_rom(:,:,1), rdims, 'x', 2)
+plotAtomTrajectory(H2_fom(:,2), all_H2_rom(:,:,2), rdims, 'y', 2)
+plotAtomTrajectory(H2_fom(:,3), all_H2_rom(:,:,3), rdims, 'z', 2)
+
+%% 
+function plotAtomTrajectory(X_fom, all_X_rom, rdims, var, idx)
+    figure;
+    hold on;
+    plot(X_fom, 'Linewidth', 2, 'DisplayName', 'FOM');
+    k = 0;
+    for rdim = rdims
+        k = k + 1;
+        X_rom = all_X_rom(k, :);
+        plot(X_rom, 'Linewidth', 2, 'DisplayName', ['ROM dim = ' int2str(rdim)]);
+    end
+    title([var '-coordinate of H' int2str(idx)])
+    legend;
+    saveas(gcf, [var '_H' int2str(idx)], 'jpg');    
+end
+
+function plotAngleHistogram(X1, X2, suffix)
+    figure;
+    A = acosd(sum(X1.*X2,2) ./ sqrt(sum(X1.^2,2)) ./ sqrt(sum(X2.^2,2)));
+    histogram(A, 20);
+    xlabel('Angle');
+    ylabel('Frequency');
+    title('Histogram of angle');
+    saveas(gcf, ['angle_' suffix], 'jpg');
+end
\ No newline at end of file
diff --git a/src/Control.cc b/src/Control.cc
index bd8f5ced..3f8a11cc 100644
--- a/src/Control.cc
+++ b/src/Control.cc
@@ -2061,6 +2061,8 @@ void Control::setROMOptions(const boost::program_options::variables_map& vm)
         rom_pri_option.save_librom_snapshot = vm["ROM.offline.save_librom_snapshot"].as<bool>();
         rom_pri_option.librom_snapshot_freq = vm["ROM.offline.librom_snapshot_freq"].as<int>();
 
+        rom_pri_option.compare_md = vm["ROM.basis.compare_md"].as<bool>();
+        rom_pri_option.num_orbbasis = vm["ROM.basis.number_of_orbital_basis"].as<int>();
         rom_pri_option.num_potbasis = vm["ROM.basis.number_of_potential_basis"].as<int>();
         rom_pri_option.pot_rom_file = vm["ROM.potential_rom_file"].as<std::string>();
     }  // onpe0
@@ -2112,6 +2114,12 @@ void Control::syncROMOptions()
     rom_pri_option.rom_stage = static_cast<ROMStage>(rom_stage);
     rom_pri_option.variable = static_cast<ROMVariable>(rom_var);
 
+    mpirc = MPI_Bcast(&rom_pri_option.compare_md, 1, MPI_C_BOOL, 0, comm_global_);
+    bcast_check(mpirc);
+
+    mpirc = MPI_Bcast(&rom_pri_option.num_orbbasis, 1, MPI_INT, 0, comm_global_);
+    bcast_check(mpirc);
+
     mpirc = MPI_Bcast(&rom_pri_option.num_potbasis, 1, MPI_INT, 0, comm_global_);
     bcast_check(mpirc);
 }
diff --git a/src/MGmol.h b/src/MGmol.h
index a6c383f3..f7525600 100644
--- a/src/MGmol.h
+++ b/src/MGmol.h
@@ -360,6 +360,7 @@ class MGmol : public MGmolInterface
 
 #ifdef MGMOL_HAS_LIBROM
     int save_orbital_snapshot(std::string snapshot_dir, OrbitalsType& orbitals);
+    void project_orbital(std::string snapshot_dir, int rdim, OrbitalsType& orbitals);
 #endif
 };
 // Instantiate static variables here to avoid clang warnings
diff --git a/src/md.cc b/src/md.cc
index ed1e5969..314fd0dc 100644
--- a/src/md.cc
+++ b/src/md.cc
@@ -492,6 +492,41 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
         // Compute forces
         force(**orbitals, ions);
 
+#ifdef MGMOL_HAS_LIBROM
+        // TODO: cleanup
+        if (ct.getROMOptions().num_orbbasis > 0)
+        {
+            if (onpe0)
+            {
+                os_ << "Projecting orbitals onto ROM subspaces to compare " 
+                    << ((ct.getROMOptions().compare_md) ? "MD dynamics" : "force") << std::endl;
+                os_ << "Loading ROM basis " << ct.getROMOptions().basis_file << std::endl;
+                os_ << "ROM basis dimension = " << ct.getROMOptions().num_orbbasis << std::endl;
+            }
+            project_orbital(ct.getROMOptions().basis_file, ct.getROMOptions().num_orbbasis, **orbitals);
+            if (ct.getROMOptions().compare_md)
+            {
+                force(**orbitals, ions);
+            }
+            else
+            {
+                double shift[3];
+                for (short i = 0; i < 3; i++) shift[i] = 0.;
+                Ions ROM_ions(ions, shift);
+                force(**orbitals, ROM_ions);
+                std::string zero = "0";
+                if (ions_->getNumIons() < 256 || ct.verbose > 2)
+                {
+                    if (ct.verbose > 0) ROM_ions.printForcesGlobal(os_);
+                }
+                else if (zero.compare(ct.md_print_filename) == 0)
+                {
+                    ROM_ions.printForcesLocal(os_);
+                }
+            }
+        }
+#endif
+
         // set fion
         ions.getLocalForces(fion);
 
@@ -633,8 +668,8 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
 
 #ifdef MGMOL_HAS_LIBROM
         // Save orbital snapshots
-        if (md_iteration_ % librom_snapshot_freq == 0
-            && ct.getROMOptions().save_librom_snapshot > 0)
+        if (ct.getROMOptions().save_librom_snapshot > 0 && 
+            md_iteration_ % librom_snapshot_freq == 0)
         {
             int ierr = save_orbital_snapshot(
                 ct.md_print_filename + "_mdstep" + std::to_string(mdstep), **orbitals);
diff --git a/src/read_config.cc b/src/read_config.cc
index 4f3b34c1..884aa8d3 100644
--- a/src/read_config.cc
+++ b/src/read_config.cc
@@ -430,6 +430,10 @@ void setupROMConfigOption(po::options_description &rom_cfg)
             "Frequency of saving libROM snapshot file at FOM simulation.")
         ("ROM.offline.variable", po::value<std::string>()->default_value(""),
             "FOM variable to perform POD: either orbitals or potential.")
+        ("ROM.basis.compare_md", po::value<bool>()->default_value(false),
+            "Compare MD or single-step force.")
+        ("ROM.basis.number_of_orbital_basis", po::value<int>()->default_value(-1),
+            "Number of orbital POD basis.")
         ("ROM.basis.number_of_potential_basis", po::value<int>()->default_value(-1),
             "Number of potential POD basis to build Hartree potential ROM operator.")
         ("ROM.potential_rom_file", po::value<std::string>()->default_value(""),
diff --git a/src/rom.cc b/src/rom.cc
index a554d2cf..19ff66a3 100644
--- a/src/rom.cc
+++ b/src/rom.cc
@@ -57,6 +57,40 @@ int MGmol<OrbitalsType>::save_orbital_snapshot(std::string file_path, OrbitalsTy
     return 0;
 }
 
+template <class OrbitalsType>
+void MGmol<OrbitalsType>::project_orbital(std::string file_path, int rdim, OrbitalsType& orbitals)
+{
+    const int dim = orbitals.getLocNumpt();
+    const int totalSamples = orbitals.chromatic_number();
+
+    CAROM::Options svd_options(dim, totalSamples, 1);
+    CAROM::BasisGenerator basis_generator(svd_options, false, "foo");
+
+    for (int i = 0; i < totalSamples; ++i)
+        basis_generator.takeSample(orbitals.getPsi(i));
+    const CAROM::Matrix* orbital_snapshots = basis_generator.getSnapshotMatrix();
+
+    CAROM::BasisReader reader(file_path);
+    CAROM::Matrix* orbital_basis = reader.getSpatialBasis(rdim);
+
+    CAROM::Matrix* proj_orbital_coeff = orbital_basis->transposeMult(orbital_snapshots);
+    CAROM::Matrix* proj_orbital_snapshots = orbital_basis->mult(proj_orbital_coeff);
+    
+    Control& ct = *(Control::instance());
+    Mesh* mesh = Mesh::instance();
+    pb::GridFunc<ORBDTYPE> gf_psi(mesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]);
+    CAROM::Vector snapshot, proj_snapshot;
+    for (int i = 0; i < totalSamples; ++i)
+    {
+        orbital_snapshots->getColumn(i, snapshot);
+        proj_orbital_snapshots->getColumn(i, proj_snapshot);
+        gf_psi.assign(proj_snapshot.getData());
+        orbitals.setPsi(gf_psi, i);
+        snapshot -= proj_snapshot;
+        std::cout << "Error for orbital " << i << " = " << snapshot.norm() << std::endl;
+    }
+}
+
 template class MGmol<LocGridOrbitals>;
 template class MGmol<ExtendedGridOrbitals>;
 
diff --git a/src/rom_Control.h b/src/rom_Control.h
index 8218bbb4..200061cb 100644
--- a/src/rom_Control.h
+++ b/src/rom_Control.h
@@ -45,13 +45,15 @@ struct ROMPrivateOptions
     int restart_file_minidx = -1;
     int restart_file_maxidx = -1;
     std::string basis_file = "";
-    ROMVariable variable=ROMVariable::NONE;
+    ROMVariable variable = ROMVariable::NONE;
 
-    /* save librom snapshot matrix at FOM simulation. */
+    /* save librom orbital snapshot matrix at FOM simulation. */
     bool save_librom_snapshot = false;
     int librom_snapshot_freq = -1;
 
     /* options for ROM building */
+    bool compare_md = false;
+    int num_orbbasis = -1;
     int num_potbasis = -1;
     std::string pot_rom_file = "";
 };