diff --git a/examples/PinnedH2O/coords.in b/examples/PinnedH2O/coords_test1.in
similarity index 100%
rename from examples/PinnedH2O/coords.in
rename to examples/PinnedH2O/coords_test1.in
diff --git a/examples/PinnedH2O/coords_test2.in b/examples/PinnedH2O/coords_test2.in
new file mode 100644
index 00000000..e1e31f8b
--- /dev/null
+++ b/examples/PinnedH2O/coords_test2.in
@@ -0,0 +1,3 @@
+O1  1   0.00    0.00    0.00  0
+H1  2  -0.45    1.57   -1.07  1
+H2  2  -0.45   -1.48   -0.97  1
diff --git a/examples/PinnedH2O/get_result.sh b/examples/PinnedH2O/get_result.sh
index 443e4cdd..b737960d 100644
--- a/examples/PinnedH2O/get_result.sh
+++ b/examples/PinnedH2O/get_result.sh
@@ -1,27 +1,37 @@
-filename="offline_PinnedH2O.out" # FOM
-#filename="rom39_PinnedH2O.out" # compare MD
-#filename="39_force_PinnedH2O.out" # compare force
+#filename="offline_PinnedH2O" # FOM
+#filename="rom39_PinnedH2O" # ROM compare MD
+#filename="39_force_PinnedH2O" # ROM compare force
+
+#filename="PinnedH2O_test2_ref" # FOM
+filename="PinnedH2O_rom_3DOF_test2_2_2_34" # ROM PinnedH2O 3DOF MD
+
+# Extracting kinetic energy, total energy, temperature from MGmgol output log
+awk '/Kinetic/ {print $3}' $filename.out > ke_$filename.txt
+awk '/Kinetic/ {print $5}' $filename.out >temp_$filename.txt
+awk '/Total/ {print $3}' $filename.out > te_$filename.txt
 
 # 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
+awk '/O1 / {print $4, $5, $6}' $filename.out > O1_$filename.txt
+awk '/H1 / {print $3, $4, $5}' $filename.out > H1_$filename.txt
+awk '/H2 / {print $3, $4, $5}' $filename.out > H2_$filename.txt
+awk '/O1 / {print $7, $8, $9}' $filename.out > f_O1_$filename.txt
+awk '/H1 / {print $6, $7, $8}' $filename.out > f_H1_$filename.txt
+awk '/H2 / {print $6, $7, $8}' $filename.out > f_H2_$filename.txt
 
 # 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 '1~2p' H1_$filename.out > H1_rom$filename.txt
+    sed -n '1~2p' H2_$filename.out > H2_rom$filename.txt
+    sed -n '1~2p' f_H1_$filename.out > f_H1_rom$filename.txt
+    sed -n '1~2p' f_H2_$filename.out > f_H2_rom$filename.txt
 
-    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
+    sed -n '2~2p' H1_$filename.out > H1_fom$filename.txt
+    sed -n '2~2p' H2_$filename.out > H2_fom$filename.txt
+    sed -n '2~2p' f_H1_$filename.out > f_H1_fom$filename.txt
+    sed -n '2~2p' f_H2_$filename.out > f_H2_fom$filename.txt
 fi
 
-rm -rf snapshot0_*
+rm -rf snapshot_*
diff --git a/examples/PinnedH2O/job.basis_1_50 b/examples/PinnedH2O/job.basis_1_50
index a5069570..bf0a4c49 100644
--- a/examples/PinnedH2O/job.basis_1_50
+++ b/examples/PinnedH2O/job.basis_1_50
@@ -1,16 +1,16 @@
 #!/bin/tcsh
-#SBATCH -N 2
+#SBATCH -N 1
 #SBATCH -t 0:10:00
-#SBATCH -p pbatch
+#SBATCH -p pdebug
 
 date
 
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 64
+set ncpus = 8
 
-set maindir = /p/lustre2/cheung26/mgmol-20241012
+set maindir = /p/lustre2/cheung26/mgmol-20250219
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
diff --git a/examples/PinnedH2O/job.offline b/examples/PinnedH2O/job.offline
index 0f8ef90e..64997e09 100644
--- a/examples/PinnedH2O/job.offline
+++ b/examples/PinnedH2O/job.offline
@@ -1,16 +1,16 @@
 #!/bin/tcsh
-#SBATCH -N 2
+#SBATCH -N 1
 #SBATCH -t 1:00:00
-#SBATCH -p pbatch
+#SBATCH -p pdebug
 
 date
 
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 64
+set ncpus = 8
 
-set maindir = /p/lustre2/cheung26/mgmol-20241012
+set maindir = /p/lustre2/cheung26/mgmol-20250219
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
@@ -18,12 +18,7 @@ set exe     = mgmol-opt
 
 cp $maindir/install_quartz/bin/$exe .
 
-set datadir = $maindir/examples/PinnedH2O
-
 set cfg_offline = mgmol_offline.cfg
-cp $datadir/$cfg_offline .
-
-cp $datadir/coords.in .
 
 ln -s -f $maindir/potentials/pseudo.O_ONCV_PBE_SG15 .
 ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 .
diff --git a/examples/PinnedH2O/job.ref b/examples/PinnedH2O/job.ref
new file mode 100644
index 00000000..2c4660cf
--- /dev/null
+++ b/examples/PinnedH2O/job.ref
@@ -0,0 +1,36 @@
+#!/bin/tcsh
+#SBATCH -N 1
+#SBATCH -t 1:00:00
+#SBATCH -p pdebug
+
+date
+
+setenv OMP_NUM_THREADS 1
+#setenv KMP_DETERMINISTIC_REDUCTION 1
+
+set ncpus = 8
+set case = 2
+
+set maindir = /p/lustre2/cheung26/mgmol-20250219
+
+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 cfg = mgmol_ref_test${case}.cfg
+cp $datadir/$cfg .
+
+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 -i coords_test${case}.in > PinnedH2O_test${case}_ref.out
+
+date
diff --git a/examples/PinnedH2O/job.rom_1_50 b/examples/PinnedH2O/job.rom
similarity index 79%
rename from examples/PinnedH2O/job.rom_1_50
rename to examples/PinnedH2O/job.rom
index aa6df335..fc4ae1b0 100644
--- a/examples/PinnedH2O/job.rom_1_50
+++ b/examples/PinnedH2O/job.rom
@@ -1,16 +1,16 @@
 #!/bin/tcsh
-#SBATCH -N 2
+#SBATCH -N 1
 #SBATCH -t 1:00:00
-#SBATCH -p pbatch
+#SBATCH -p pdebug
 
 date
 
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 64
+set ncpus = 8
 
-set maindir = /p/lustre2/cheung26/mgmol-20241012
+set maindir = /p/lustre2/cheung26/mgmol-20250219
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
@@ -18,16 +18,11 @@ 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 .
diff --git a/examples/PinnedH2O/job.rom_3DOF b/examples/PinnedH2O/job.rom_3DOF
new file mode 100644
index 00000000..8b328199
--- /dev/null
+++ b/examples/PinnedH2O/job.rom_3DOF
@@ -0,0 +1,35 @@
+#!/bin/tcsh
+#SBATCH -N 1
+#SBATCH -t 1:00:00
+#SBATCH -p pdebug
+
+date
+
+setenv OMP_NUM_THREADS 1
+#setenv KMP_DETERMINISTIC_REDUCTION 1
+
+set ncpus = 8
+set case = 2
+
+set maindir = /p/lustre2/cheung26/mgmol-20250219
+
+setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
+
+set exe     = mgmol-opt
+
+cp $maindir/install_quartz/bin/$exe .
+
+set bondlength_num_increments = 2
+set bondangle_num_increments = 2
+set num_basis = 34
+
+set cfg_rom = mgmol_rom_3DOF_test${case}.cfg
+
+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_test${case}.in > PinnedH2O_rom_3DOF_test${case}_${bondlength_num_increments}_${bondangle_num_increments}_${num_basis}.out
+
+date
diff --git a/examples/PinnedH2O/mgmol.cfg b/examples/PinnedH2O/mgmol_ref_test1.cfg
similarity index 89%
rename from examples/PinnedH2O/mgmol.cfg
rename to examples/PinnedH2O/mgmol_ref_test1.cfg
index ef82f07d..fd456f92 100644
--- a/examples/PinnedH2O/mgmol.cfg
+++ b/examples/PinnedH2O/mgmol_ref_test1.cfg
@@ -18,12 +18,12 @@ pseudopotential=pseudo.H_ONCV_PBE_SG15
 [Run]
 type=MD
 [MD]
-num_steps=50
+num_steps=500
 dt=40.
 thermostat=ON
 [Thermostat]
 type=Berendsen
-temperature=1000.
+temperature=300.
 relax_time=800.
 [Quench]
 max_steps=100
@@ -32,4 +32,4 @@ atol=1.e-8
 initial_type=Random
 initial_width=2.
 [Restart]
-output_level=0
+output_level=4
diff --git a/examples/PinnedH2O/mgmol_ref_test2.cfg b/examples/PinnedH2O/mgmol_ref_test2.cfg
new file mode 100644
index 00000000..67aa443c
--- /dev/null
+++ b/examples/PinnedH2O/mgmol_ref_test2.cfg
@@ -0,0 +1,31 @@
+verbosity=1
+xcFunctional=PBE
+FDtype=4th
+[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=OFF
+[Quench]
+max_steps=100
+atol=1.e-8
+[Orbitals]
+initial_type=Random
+initial_width=2.
+[Restart]
+output_level=4
diff --git a/examples/PinnedH2O/mgmol_rom_1_50.cfg b/examples/PinnedH2O/mgmol_rom_1_50.cfg
index 337a86eb..13cc3921 100644
--- a/examples/PinnedH2O/mgmol_rom_1_50.cfg
+++ b/examples/PinnedH2O/mgmol_rom_1_50.cfg
@@ -1,6 +1,6 @@
 verbosity=1
 xcFunctional=PBE
-FDtype=Mehrstellen
+FDtype=4th
 [Mesh]
 nx=64
 ny=64
@@ -23,7 +23,7 @@ dt=40.
 thermostat=ON
 [Thermostat]
 type=Berendsen
-temperature=1000.
+temperature=300.
 relax_time=800.
 [Quench]
 max_steps=100
@@ -33,6 +33,8 @@ initial_type=Random
 initial_width=2.
 [Restart]
 output_level=4
+[ROM]
+stage=test_orbital
 [ROM.offline]
 basis_file=PinnedH2O_orbitals_basis_1_50
 [ROM.basis]
diff --git a/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg b/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg
new file mode 100644
index 00000000..adfdfc43
--- /dev/null
+++ b/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg
@@ -0,0 +1,46 @@
+verbosity=1
+xcFunctional=PBE
+FDtype=4th
+[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=300.
+relax_time=800.
+[Quench]
+max_steps=100
+atol=1.e-8
+[Orbitals]
+initial_type=Random
+initial_width=2.
+nempty=30
+[DensityMatrix]
+solver=MVP
+nb_inner_it=100
+[Restart]
+output_level=4
+[ROM]
+stage=online_pinned_H2O_3dof
+[ROM.offline]
+basis_file=/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_8/PinnedH2O_3DOF_orbitals_basis_2_2
+[ROM.basis]
+compare_md=false
+number_of_orbital_basis=34
diff --git a/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg b/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg
new file mode 100644
index 00000000..f0368790
--- /dev/null
+++ b/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg
@@ -0,0 +1,42 @@
+verbosity=1
+xcFunctional=PBE
+FDtype=4th
+[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=OFF
+[Quench]
+max_steps=100
+atol=1.e-8
+[Orbitals]
+initial_type=Random
+initial_width=2.
+nempty=30
+[DensityMatrix]
+solver=MVP
+nb_inner_it=100
+[Restart]
+output_level=4
+[ROM]
+stage=online_pinned_H2O_3dof
+[ROM.offline]
+basis_file=/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_8/PinnedH2O_3DOF_orbitals_basis_2_2
+[ROM.basis]
+compare_md=false
+number_of_orbital_basis=34
diff --git a/src/Control.cc b/src/Control.cc
index d0b1a224..151fe920 100644
--- a/src/Control.cc
+++ b/src/Control.cc
@@ -1953,6 +1953,10 @@ void Control::setROMOptions(const boost::program_options::variables_map& vm)
             rom_pri_option.rom_stage = ROMStage::ONLINE;
         else if (str.compare("build") == 0)
             rom_pri_option.rom_stage = ROMStage::BUILD;
+        else if (str.compare("online_pinned_H2O_3dof") == 0)
+            rom_pri_option.rom_stage = ROMStage::ONLINE_PINNED_H2O_3DOF;
+        else if (str.compare("test_orbital") == 0)
+            rom_pri_option.rom_stage = ROMStage::TEST_ORBITAL;
         else if (str.compare("online_poisson") == 0)
             rom_pri_option.rom_stage = ROMStage::ONLINE_POISSON;
         else if (str.compare("test_poisson") == 0)
diff --git a/src/Ions.cc b/src/Ions.cc
index cb6dddf5..072f0950 100644
--- a/src/Ions.cc
+++ b/src/Ions.cc
@@ -1315,6 +1315,22 @@ void Ions::setLocalForces(
         {
             if (ion->compareName(*s))
             {
+                //std::cout << "Ion found: " << ion->name() << std::endl;
+                //std::cout << "Ion force: (" << *it << ", " << *(it + 1) << ", " << *(it + 2) << ")" << std::endl;
+                //std::cout << "names: ";
+                //for (int i = 0; i < names.size(); i++)
+                //{
+                //    std::cout << names[i];
+                //    if (i == forces.size() - 1) std::cout << std::endl;
+                //    else std::cout << ", ";
+                //}
+                //std::cout << "forces: ";
+                //for (int i = 0; i < forces.size(); i++)
+                //{
+                //    std::cout << forces[i];
+                //    if (i == forces.size() - 1) std::cout << ")" << std::endl;
+                //    else std::cout << ", ";
+                //}
                 ion->set_force(0, *it);
                 ion->set_force(1, *(it + 1));
                 ion->set_force(2, *(it + 2));
@@ -1343,11 +1359,14 @@ void Ions::setLocalForces(
             double p[3];
             ion->getPosition(&p[0]);
             double d2 = (p[0] - (*cit)) * (p[0] - (*cit))
-                        + (p[1] - (*(cit + 1))) * (p[0] - (*(cit + 1)))
-                        + (p[2] - (*(cit + 2))) * (p[0] - (*(cit + 2)));
+                        + (p[1] - (*(cit + 1))) * (p[1] - (*(cit + 1)))
+                        + (p[2] - (*(cit + 2))) * (p[2] - (*(cit + 2)));
             double d = std::sqrt(d2);
             if (d < tol)
             {
+                std::cout << "Ion found: " << ion->name() << std::endl;
+                std::cout << "Ion position:( " << *cit << ", " << *(cit + 1) << ", " << *(cit + 2) << ")" << std::endl;
+                std::cout << "Ion force: (" << *fit << ", " << *(fit + 1) << ", " << *(fit + 2) << ")" << std::endl;
                 ion->set_force(0, *fit);
                 ion->set_force(1, *(fit + 1));
                 ion->set_force(2, *(fit + 2));
diff --git a/src/MGmol.cc b/src/MGmol.cc
index 4ccb8e90..facd362e 100644
--- a/src/MGmol.cc
+++ b/src/MGmol.cc
@@ -1422,6 +1422,33 @@ void MGmol<OrbitalsType>::getAtomicNumbers(std::vector<short>& an)
     ions_->getAtomicNumbers(an);
 }
 
+template <class OrbitalsType>
+void MGmol<OrbitalsType>::updateDMandEnergy(OrbitalsType& orbitals, Ions& ions, double& eks)
+{
+    // initialize electronic density
+    rho_->update(orbitals);
+
+    // initialize potential
+    update_pot(ions);
+
+    // initialize projected matrices
+    updateHmatrix(orbitals, ions);
+    proj_matrices_->updateThetaAndHB();
+
+    // compute DM
+    std::shared_ptr<DMStrategy<OrbitalsType>> dm_strategy(
+        DMStrategyFactory<OrbitalsType,
+            dist_matrix::DistMatrix<double>>::create(comm_, os_, ions,
+            rho_.get(), energy_.get(), electrostat_.get(), this,
+            proj_matrices_.get(), &orbitals));
+
+    dm_strategy->update(orbitals);
+
+    // evaluate energy and forces
+    double ts = 0.;
+    eks = energy_->evaluateTotal(ts, proj_matrices_.get(), ions, orbitals, 2, os_);
+}
+
 template <class OrbitalsType>
 double MGmol<OrbitalsType>::evaluateEnergyAndForces(
     const std::vector<double>& tau, const std::vector<short>& atnumbers,
diff --git a/src/MGmol.h b/src/MGmol.h
index 799a7edf..c7777a62 100644
--- a/src/MGmol.h
+++ b/src/MGmol.h
@@ -368,6 +368,7 @@ class MGmol : public MGmolInterface
     int save_orbital_snapshot(std::string file_path, OrbitalsType& orbitals);
     void project_orbital(std::string file_path, int rdim, OrbitalsType& orbitals);
 #endif
+    void updateDMandEnergy(OrbitalsType& orbitals, Ions& ions, double& eks);
 };
 // Instantiate static variables here to avoid clang warnings
 template <class OrbitalsType>
diff --git a/src/md.cc b/src/md.cc
index 8ad3d7b8..bc986955 100644
--- a/src/md.cc
+++ b/src/md.cc
@@ -35,6 +35,12 @@
 #include "mgmol_Signal.h"
 #include "tools.h"
 
+#ifdef MGMOL_HAS_LIBROM
+#include "KBPsiMatrixSparse.h"
+#include "PinnedH2O.h"
+#include "rom_workflows.h"
+#endif
+
 #include <sstream>
 #include <string>
 #include <vector>
@@ -399,8 +405,27 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
         h5f_file_.reset();
     }
 
+    bool ROM_MVP = (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF);
+#ifdef MGMOL_HAS_LIBROM
+    // ROM - initialize orbitals and density matrix
+    if (ROM_MVP)
+    {
+        if (onpe0) os_ << "Setup ROM MVP solver..." << std::endl;
+        ExtendedGridOrbitals** extended_orbitals = reinterpret_cast<ExtendedGridOrbitals**>(orbitals);
+        (*extended_orbitals)->set(ct.getROMOptions().basis_file, ct.numst); 
+        (*extended_orbitals)->orthonormalizeLoewdin();
+        (*extended_orbitals)->setDataWithGhosts(true);
+        (*extended_orbitals)->setIterativeIndex(10);
+
+        std::shared_ptr<ProjectedMatricesInterface> projmatrices
+            = getProjectedMatrices();
+        projmatrices->setDMuniform(ct.getNelSpin(), 0);
+        projmatrices->printDM(os_);
+    }
+#endif
+
     // additional SC steps to compensate random start
-    if (ct.restart_info < 3)
+    if (ct.restart_info < 3 && !ROM_MVP)
     {
         double eks = 0.;
         quench(**orbitals, ions, ct.max_electronic_steps, 20, eks);
@@ -432,10 +457,43 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
         int retval              = 0;
         bool small_move         = true;
         bool last_move_is_small = true;
-        do
+
+        bool force_on_ions = true;
+#ifdef MGMOL_HAS_LIBROM
+        // Ions for ROM
+        Mesh* mymesh            = Mesh::instance();
+        const pb::Grid& mygrid  = mymesh->grid();
+        const double lattice[3] = { mygrid.ll(0), mygrid.ll(1), mygrid.ll(2) };
+        std::vector<double> positions;
+        std::vector<short> anumbers;
+        getAtomicPositions(positions);
+        getAtomicNumbers(anumbers);
+        Ions* ROM_ions;
+
+        // Pinned H2O 3 DOF
+        PinnedH2O H2O_molecule;
+        if (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF)
+        {
+            if (onpe0) os_ << "Rotate Pinned H2O molecule in timestep " << mdstep << std::endl;
+            H2O_molecule.rotate(positions, anumbers);
+            if (onpe0) H2O_molecule.print(os_);
+            ROM_ions = new Ions(positions, anumbers, lattice, ions_->getSpecies());
+            setupPotentials(*ROM_ions);
+            force_on_ions = false;
+        }
+#endif
+
+        if (ROM_MVP)
+        {
+            updateDMandEnergy(**orbitals, *ROM_ions, eks);
+        }
+        else
         {
             retval = quench(**orbitals, ions, ct.max_electronic_steps, 0, eks);
+        }
 
+        do
+        {
             // update localization regions
             if (ct.adaptiveLRs())
             {
@@ -495,10 +553,23 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
                 << std::endl;
 
         // Compute forces
-        force(**orbitals, ions);
+        if (force_on_ions)
+            force(**orbitals, ions);
 
 #ifdef MGMOL_HAS_LIBROM
-        if (ct.getROMOptions().num_orbbasis > 0)
+        if (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF)
+        {
+            force(**orbitals, *ROM_ions);
+            // Pinned H2O 3 DOF
+            if (onpe0) os_ << "Transpose rotate the PinnedH2O molecule" << std::endl;
+            std::vector<double> forces;
+            ROM_ions->getForces(forces);
+            H2O_molecule.transpose_rotate(positions, anumbers, forces);
+            ions.setLocalForces(forces, positions);
+            delete ROM_ions;
+        }
+
+        if (ct.getROMOptions().rom_stage == ROMStage::TEST_ORBITAL)
         {
             if (onpe0)
             {
@@ -507,26 +578,30 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
                 os_ << "Loading ROM basis " << ct.getROMOptions().basis_file << std::endl;
                 os_ << "ROM basis dimension = " << ct.getROMOptions().num_orbbasis << std::endl;
             }
+
+            // Project orbitals to ROM subspace
             project_orbital(ct.getROMOptions().basis_file, ct.getROMOptions().num_orbbasis, **orbitals);
             if (ct.getROMOptions().compare_md)
             {
+                // overwrite ions force for MD time stepping
                 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);
+                // write ROM_ions force for one-step comparison
+                ROM_ions = new Ions(positions, anumbers, lattice, ions_->getSpecies());
+                force(**orbitals, *ROM_ions);
                 std::string zero = "0";
                 if (ions_->getNumIons() < 256 || ct.verbose > 2)
                 {
-                    if (ct.verbose > 0) ROM_ions.printForcesGlobal(os_);
+                    if (ct.verbose > 0) ROM_ions->printForcesGlobal(os_);
+
                 }
                 else if (zero.compare(ct.md_print_filename) == 0)
                 {
-                    ROM_ions.printForcesLocal(os_);
+                    ROM_ions->printForcesLocal(os_);
                 }
+                delete ROM_ions;
             }
         }
 #endif
diff --git a/src/rom_Control.h b/src/rom_Control.h
index 200061cb..2d143d84 100644
--- a/src/rom_Control.h
+++ b/src/rom_Control.h
@@ -22,6 +22,8 @@ enum class ROMStage
     ONLINE,
     RESTORE,    // TODO(kevin): what stage is this?
     BUILD,
+    ONLINE_PINNED_H2O_3DOF,
+    TEST_ORBITAL,
     ONLINE_POISSON,
     TEST_POISSON,
     TEST_RHO,
diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt
index baae6cda..ca187a0e 100644
--- a/src/tools/CMakeLists.txt
+++ b/src/tools/CMakeLists.txt
@@ -8,6 +8,7 @@ set(SOURCES mgmol_mpi_tools.cc
        random.cc 
        coloring.cc 
        SymmetricMatrix.cc
+       PinnedH2O.cc
 )
 
 add_library(mgmol_tools ${SOURCES})
diff --git a/src/tools/PinnedH2O.cc b/src/tools/PinnedH2O.cc
new file mode 100644
index 00000000..22f4677a
--- /dev/null
+++ b/src/tools/PinnedH2O.cc
@@ -0,0 +1,237 @@
+// 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 "PinnedH2O.h"
+#include <iostream>
+
+PinnedH2O::PinnedH2O()
+    : flipped_bond(false),
+      O1_idx(-1),
+      H1_idx(-1),
+      H2_idx(-1),
+      planar_rotation_angle(0.0)
+{
+    for (int i = 0; i < 3; ++i) {
+        for (int j = 0; j < 3; ++j) {
+            out_of_plane_rotation_matrix[i][j] = 0.0;
+        }
+    }
+}
+
+double PinnedH2O::calculate_bondlength(const double atom1[3], const double atom2[3]) const
+{
+    return sqrt(pow(atom1[0] - atom2[0], 2) + pow(atom1[1] - atom2[1], 2) + pow(atom1[2] - atom2[2], 2));
+}
+
+double PinnedH2O::calculate_bondangle(const double atom1[3], const double atom2[3], const double atom3[3]) const
+{
+    double vector1[3] = {atom1[0] - atom2[0], atom1[1] - atom2[1], atom1[2] - atom2[2]};
+    double vector2[3] = {atom3[0] - atom2[0], atom3[1] - atom2[1], atom3[2] - atom2[2]};
+
+    double dot_product = vector1[0] * vector2[0] + vector1[1] * vector2[1] + vector1[2] * vector2[2];
+    double magnitude_product = sqrt(pow(vector1[0], 2) + pow(vector1[1], 2) + pow(vector1[2], 2)) *
+                             sqrt(pow(vector2[0], 2) + pow(vector2[1], 2) + pow(vector2[2], 2));
+    double angle = acos(dot_product / magnitude_product);
+
+    return angle;
+}
+
+void PinnedH2O::rotation_matrix(const double axis[3], double angle, double matrix[3][3]) const
+{
+    double cos_theta = cos(angle);
+    double sin_theta = sin(angle);
+    double ux = axis[0], uy = axis[1], uz = axis[2];
+
+    matrix[0][0] = cos_theta + ux * ux * (1 - cos_theta);
+    matrix[0][1] = ux * uy * (1 - cos_theta) - uz * sin_theta;
+    matrix[0][2] = ux * uz * (1 - cos_theta) + uy * sin_theta;
+
+    matrix[1][0] = uy * ux * (1 - cos_theta) + uz * sin_theta;
+    matrix[1][1] = cos_theta + uy * uy * (1 - cos_theta);
+    matrix[1][2] = uy * uz * (1 - cos_theta) - ux * sin_theta;
+
+    matrix[2][0] = uz * ux * (1 - cos_theta) - uy * sin_theta;
+    matrix[2][1] = uz * uy * (1 - cos_theta) + ux * sin_theta;
+    matrix[2][2] = cos_theta + uz * uz * (1 - cos_theta);
+}
+
+void PinnedH2O::normalize(double vec[3]) const
+{
+    double norm = sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
+    vec[0] /= norm;
+    vec[1] /= norm;
+    vec[2] /= norm;
+}
+
+void PinnedH2O::cross(const double a[3], const double b[3], double result[3]) const
+{
+    result[0] = a[1] * b[2] - a[2] * b[1];
+    result[1] = a[2] * b[0] - a[0] * b[2];
+    result[2] = a[0] * b[1] - a[1] * b[0];
+}
+
+void PinnedH2O::apply_rotation(const double matrix[3][3], const double vec[3], double result[3]) const
+{
+    result[0] = matrix[0][0] * vec[0] + matrix[0][1] * vec[1] + matrix[0][2] * vec[2];
+    result[1] = matrix[1][0] * vec[0] + matrix[1][1] * vec[1] + matrix[1][2] * vec[2];
+    result[2] = matrix[2][0] * vec[0] + matrix[2][1] * vec[1] + matrix[2][2] * vec[2];
+}
+
+void PinnedH2O::apply_transpose_rotation(const double matrix[3][3], const double vec[3], double result[3]) const
+{
+    result[0] = matrix[0][0] * vec[0] + matrix[1][0] * vec[1] + matrix[2][0] * vec[2];
+    result[1] = matrix[0][1] * vec[0] + matrix[1][1] * vec[1] + matrix[2][1] * vec[2];
+    result[2] = matrix[0][2] * vec[0] + matrix[1][2] * vec[1] + matrix[2][2] * vec[2];
+}
+
+void PinnedH2O::rotate(std::vector<double>& positions, std::vector<short>& anumbers)
+{
+    for (int i = 0; i < 3; i++) {
+        if (positions[3 * i] == 0.0 && positions[3 * i + 1] == 0.0 && positions[3 * i + 2] == 0.0) {
+            O1_idx = i;
+            break;
+        }
+    }
+    if (O1_idx == -1) return;
+
+    H1_idx = (O1_idx + 1) % 3;
+    H2_idx = (O1_idx + 2) % 3;
+
+    double O1[3] = {positions[3 * O1_idx], positions[3 * O1_idx + 1], positions[3 * O1_idx + 2]};
+    double H1[3] = {positions[3 * H1_idx], positions[3 * H1_idx + 1], positions[3 * H1_idx + 2]};
+    double H2[3] = {positions[3 * H2_idx], positions[3 * H2_idx + 1], positions[3 * H2_idx + 2]};
+
+    bondlength1 = calculate_bondlength(H1, O1);
+    bondlength2 = calculate_bondlength(H2, O1);
+    bondangle = calculate_bondangle(H1, O1, H2);
+
+    double H1_temp[3], H2_temp[3];
+    double H1_rotated[3], H2_rotated[3];
+
+    double plane_normal[3];
+    cross(H2, H1, plane_normal);
+    normalize(plane_normal);
+    double target_plane_normal[3] = {0, 0, 1};
+    double dot_product = plane_normal[0] * target_plane_normal[0] +
+                       plane_normal[1] * target_plane_normal[1] +
+                       plane_normal[2] * target_plane_normal[2];
+    double angle_to_align = std::acos(std::min(std::max(dot_product, -1.0), 1.0));
+    double axis_to_align[3];
+    if (abs(dot_product) > 1.0 - 1e-8)
+    {
+        axis_to_align[0] = 1.0;
+        axis_to_align[1] = 0.0;
+        axis_to_align[2] = 0.0;
+    }
+    else
+    {
+        cross(plane_normal, target_plane_normal, axis_to_align);
+        normalize(axis_to_align);
+    }
+
+    rotation_matrix(axis_to_align, angle_to_align, out_of_plane_rotation_matrix);
+    apply_rotation(out_of_plane_rotation_matrix, H1, H1_temp);
+    apply_rotation(out_of_plane_rotation_matrix, H2, H2_temp);
+
+    double theta1 = std::atan2(H1_temp[1], H1_temp[0]);
+    planar_rotation_angle = -theta1 + bondangle / 2.0;
+    double planar_rotation_matrix[3][3];
+    rotation_matrix(target_plane_normal, planar_rotation_angle, planar_rotation_matrix);
+    apply_rotation(planar_rotation_matrix, H1_temp, H1_rotated);
+    apply_rotation(planar_rotation_matrix, H2_temp, H2_rotated);
+
+    flipped_bond = (bondlength1 < bondlength2);
+    if (flipped_bond)
+    {
+        H1_rotated[1] *= -1.0;
+        H2_rotated[1] *= -1.0;
+        std::swap(H1_rotated, H2_rotated);
+        std::swap(bondlength1, bondlength2);
+    }
+
+    positions[0] = H2_rotated[0];
+    positions[1] = H2_rotated[1];
+    positions[2] = H2_rotated[2];
+    positions[3] = 0.0;
+    positions[4] = 0.0;
+    positions[5] = 0.0;
+    positions[6] = H1_rotated[0];
+    positions[7] = H1_rotated[1];
+    positions[8] = H1_rotated[2];
+
+    anumbers[0] = 1;
+    anumbers[1] = 8;
+    anumbers[2] = 1;
+}
+
+void PinnedH2O::transpose_rotate(std::vector<double>& positions, std::vector<short>& anumbers, std::vector<double>& forces)
+{
+    double H2_rotated[3] = {positions[0], positions[1], positions[2]};
+    double O1_rotated[3] = {positions[3], positions[4], positions[5]};
+    double H1_rotated[3] = {positions[6], positions[7], positions[8]};
+
+    double f_H2_rotated[3] = {forces[0], forces[1], forces[2]};
+    double f_O1_rotated[3] = {forces[3], forces[4], forces[5]};
+    double f_H1_rotated[3] = {forces[6], forces[7], forces[8]};
+
+    if (flipped_bond)
+    {
+        H1_rotated[1] *= -1.0;
+        H2_rotated[1] *= -1.0;
+        f_O1_rotated[1] *= -1.0;
+        f_H1_rotated[1] *= -1.0;
+        f_H2_rotated[1] *= -1.0;
+    }
+
+    double planar_rotation_matrix[3][3];
+    double target_plane_normal[3] = {0, 0, 1};
+    rotation_matrix(target_plane_normal, planar_rotation_angle, planar_rotation_matrix);
+
+    double H1_temp[3], H2_temp[3];
+    apply_transpose_rotation(planar_rotation_matrix, H1_rotated, H1_temp);
+    apply_transpose_rotation(planar_rotation_matrix, H2_rotated, H2_temp);
+
+    double f_O1_temp[3], f_H1_temp[3], f_H2_temp[3];
+    apply_transpose_rotation(planar_rotation_matrix, f_O1_rotated, f_O1_temp);
+    apply_transpose_rotation(planar_rotation_matrix, f_H1_rotated, f_H1_temp);
+    apply_transpose_rotation(planar_rotation_matrix, f_H2_rotated, f_H2_temp);
+
+    double H1_restored[3], H2_restored[3];
+    apply_transpose_rotation(out_of_plane_rotation_matrix, H1_temp, H1_restored);
+    apply_transpose_rotation(out_of_plane_rotation_matrix, H2_temp, H2_restored);
+
+    double f_O1_restored[3], f_H1_restored[3], f_H2_restored[3];
+    apply_transpose_rotation(out_of_plane_rotation_matrix, f_O1_temp, f_O1_restored);
+    apply_transpose_rotation(out_of_plane_rotation_matrix, f_H1_temp, f_H1_restored);
+    apply_transpose_rotation(out_of_plane_rotation_matrix, f_H2_temp, f_H2_restored);
+
+    positions[3*H2_idx] = H2_restored[0];
+    positions[3*H2_idx+1] = H2_restored[1];
+    positions[3*H2_idx+2] = H2_restored[2];
+    positions[3*O1_idx] = 0.0;
+    positions[3*O1_idx+1] = 0.0; 
+    positions[3*O1_idx+2] = 0.0; 
+    positions[3*H1_idx] = H1_restored[0];
+    positions[3*H1_idx+1] = H1_restored[1];
+    positions[3*H1_idx+2] = H1_restored[2];
+
+    anumbers[H2_idx] = 1;
+    anumbers[O1_idx] = 8;
+    anumbers[H1_idx] = 1;
+
+    forces[3*H2_idx] = f_H2_restored[0];
+    forces[3*H2_idx+1] = f_H2_restored[1];
+    forces[3*H2_idx+2] = f_H2_restored[2];
+    forces[3*O1_idx] = f_O1_restored[0];
+    forces[3*O1_idx+1] = f_O1_restored[1];
+    forces[3*O1_idx+2] = f_O1_restored[2];
+    forces[3*H1_idx] = f_H1_restored[0];
+    forces[3*H1_idx+1] = f_H1_restored[1];
+    forces[3*H1_idx+2] = f_H1_restored[2];
+}
diff --git a/src/tools/PinnedH2O.h b/src/tools/PinnedH2O.h
new file mode 100644
index 00000000..bee880d8
--- /dev/null
+++ b/src/tools/PinnedH2O.h
@@ -0,0 +1,54 @@
+// 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 PINNED_H2O_H
+#define PINNED_H2O_H
+
+#include <vector>
+#include <cmath>
+#include <algorithm>
+#include <iostream>
+
+class PinnedH2O
+{
+
+public:
+    PinnedH2O();
+    ~PinnedH2O() = default;
+
+    void rotate(std::vector<double>& positions, std::vector<short>& anumbers);
+    void transpose_rotate(std::vector<double>& positions, std::vector<short>& anumbers, std::vector<double>& forces);
+    void print(std::ostream& os)
+    {
+        os << "Bondlengths = " << bondlength1 << ", " << bondlength2 << " Bohrs; "
+           << "Bondangle = " << bondangle * 180.0 / M_PI << " degrees." << std::endl;
+    }
+
+private:
+    double calculate_bondlength(const double atom1[3], const double atom2[3]) const;
+    double calculate_bondangle(const double atom1[3], const double atom2[3], const double atom3[3]) const;
+    void rotation_matrix(const double axis[3], double angle, double matrix[3][3]) const;
+    void normalize(double vec[3]) const;
+    void cross(const double a[3], const double b[3], double result[3]) const;
+    void apply_rotation(const double matrix[3][3], const double vec[3], double result[3]) const;
+    void apply_transpose_rotation(const double matrix[3][3], const double vec[3], double result[3]) const;
+
+    double bondlength1;
+    double bondlength2;
+    double bondangle;
+
+    bool flipped_bond;
+    int O1_idx;
+    int H1_idx;
+    int H2_idx;
+    double planar_rotation_angle;
+    double out_of_plane_rotation_matrix[3][3];
+};
+
+#endif // PINNED_H2O_H
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 02edbe8b..b4844557 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -291,7 +291,8 @@ add_executable(testDMandEnergyAndForces
 add_executable(testRestartEnergyAndForces
                ${CMAKE_SOURCE_DIR}/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc)
 add_executable(testPinnedH2O_3DOF
-               ${CMAKE_SOURCE_DIR}/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc)
+               ${CMAKE_SOURCE_DIR}/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
+               ${CMAKE_SOURCE_DIR}/src/tools/PinnedH2O.cc)
 
 if(${MAGMA_FOUND})
   add_executable(testOpenmpOffload
diff --git a/tests/PinnedH2O_3DOF/coords.in_rotate1 b/tests/PinnedH2O_3DOF/coords_rotate1.in
similarity index 100%
rename from tests/PinnedH2O_3DOF/coords.in_rotate1
rename to tests/PinnedH2O_3DOF/coords_rotate1.in
diff --git a/tests/PinnedH2O_3DOF/coords.in_rotate2 b/tests/PinnedH2O_3DOF/coords_rotate2.in
similarity index 100%
rename from tests/PinnedH2O_3DOF/coords.in_rotate2
rename to tests/PinnedH2O_3DOF/coords_rotate2.in
diff --git a/tests/PinnedH2O_3DOF/job.basis b/tests/PinnedH2O_3DOF/job.basis
index eaf8bfde..cab929d1 100644
--- a/tests/PinnedH2O_3DOF/job.basis
+++ b/tests/PinnedH2O_3DOF/job.basis
@@ -1,5 +1,5 @@
 #!/bin/tcsh
-#SBATCH -N 2
+#SBATCH -N 1
 #SBATCH -t 0:10:00
 #SBATCH -p pdebug
 
@@ -8,9 +8,9 @@ date
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 64
+set ncpus = 8
 
-set maindir = /p/lustre2/cheung26/mgmol-20250206
+set maindir = /p/lustre2/cheung26/mgmol-20250219
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
@@ -19,10 +19,10 @@ set exe     = ${maindir}/build_quartz/libROM/build/examples/misc/combine_samples
 set snapshot_files = ""
 set bondlength_min = 0.95
 set bondlength_max = 1.05
-set bondlength_num_increments = 5
+set bondlength_num_increments = 2
 set bondangle_min = -5.0 
 set bondangle_max = 5.0
-set bondangle_num_increments = 5
+set bondangle_num_increments = 2
 
 foreach i (`seq 0 $bondlength_num_increments`)
   set bondlength_one = `echo "scale=2; $bondlength_min + $i * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc`
diff --git a/tests/PinnedH2O_3DOF/job.offline b/tests/PinnedH2O_3DOF/job.offline
index 3870a3a3..243f563a 100644
--- a/tests/PinnedH2O_3DOF/job.offline
+++ b/tests/PinnedH2O_3DOF/job.offline
@@ -1,5 +1,5 @@
 #!/bin/tcsh
-#SBATCH -N 2
+#SBATCH -N 1
 #SBATCH -t 4:00:00
 #SBATCH -p pbatch
 
@@ -8,15 +8,14 @@ date
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 64
+set ncpus = 8
 
-set maindir = /p/lustre2/cheung26/mgmol-20250206
+set maindir = /p/lustre2/cheung26/mgmol-20250219
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
 set exe     = mgmol-opt
 
-rm $exe
 cp $maindir/install_quartz/bin/$exe .
 
 set cfg_offline = mgmol_offline.cfg
@@ -46,14 +45,13 @@ foreach i (`seq 0 $bondlength_num_increments`)
       set tag = ${bondlength_one}_${bondlength_two}_${bondangle}
       if (-e data/${tag}) then
         continue
+      else
+        mkdir data/${tag}
       endif
-      cp $cfg_offline ${tag}_${cfg_offline}
-      sed -i "s/CUSTOM_TAG/${tag}/g" ${tag}_${cfg_offline}
-      srun -n $ncpus $exe -c ${tag}_${cfg_offline} -i data/coords_${tag}.in > ${tag}_offline_PinnedH2O.out
-      mv data/coords_${tag}.in ${tag}/coords.in
-      mv ${tag}_${cfg_offline} ${tag}/${cfg_offline}
-      mv ${tag}_offline_PinnedH2O.out ${tag}/offline_PinnedH2O.out
-      mv -f ${tag} data
+      cp $cfg_offline data/${tag}/${cfg_offline}
+      mv data/coords_${tag}.in data/${tag}/coords.in
+      sed -i "s/CUSTOM_TAG/${tag}/g" data/${tag}/${cfg_offline}
+      srun -n $ncpus $exe -c data/${tag}/${cfg_offline} -i data/${tag}/coords.in > data/${tag}/offline_PinnedH2O.out
     end
   end
 end
diff --git a/tests/PinnedH2O_3DOF/job.online b/tests/PinnedH2O_3DOF/job.online
index 77b8fde0..1373e011 100644
--- a/tests/PinnedH2O_3DOF/job.online
+++ b/tests/PinnedH2O_3DOF/job.online
@@ -1,5 +1,5 @@
 #!/bin/tcsh
-#SBATCH -N 2
+#SBATCH -N 1
 #SBATCH -t 1:00:00
 #SBATCH -p pdebug
 
@@ -8,15 +8,14 @@ date
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 64
+set ncpus = 8
 
-set maindir = /p/lustre2/cheung26/mgmol-20250206
+set maindir = /p/lustre2/cheung26/mgmol-20250219
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
 set exe     = testPinnedH2O_3DOF
 
-rm $exe
 cp $maindir/build_quartz/tests/$exe .
 
 set cfg_online = mgmol_online.cfg
@@ -26,9 +25,9 @@ ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 .
 
 source $maindir/scripts/modules.quartz
 
-set sample_bondlength_num_increments = 5
-set sample_bondangle_num_increments = 5
-set num_orbital_basis = 36
+set sample_bondlength_num_increments = 2
+set sample_bondangle_num_increments = 2
+set num_orbital_basis = 34
 
 set bondlength_min = 0.95
 set bondlength_max = 1.05
@@ -58,9 +57,7 @@ foreach i (`seq 0 $bondlength_num_increments`)
       set bondangle = `printf "%.1f" $bondangle`
       echo "bondlength1: $bondlength_one, bondlength2: $bondlength_two, bondangle: $bondangle"
       set tag = ${bondlength_one}_${bondlength_two}_${bondangle}
- 
-      echo $ncpus $exe $case $cfg_online $tag
-      srun -p pdebug -n $ncpus $exe -c ${case}_${cfg_online} -i data/${tag}/coords.in >& results_${case}/${tag}.log
+      srun -p pdebug -n $ncpus $exe -c ${case}_${cfg_online} -i /usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_${ncpus}/${tag}/coords.in >& results_${case}/${tag}.log
     end
   end
 end
diff --git a/tests/PinnedH2O_3DOF/job.online_rotate b/tests/PinnedH2O_3DOF/job.online_rotate
index 11dd7419..3f8a65ad 100644
--- a/tests/PinnedH2O_3DOF/job.online_rotate
+++ b/tests/PinnedH2O_3DOF/job.online_rotate
@@ -1,5 +1,5 @@
 #!/bin/tcsh
-#SBATCH -N 2
+#SBATCH -N 1
 #SBATCH -t 0:10:00
 #SBATCH -p pdebug
 
@@ -8,9 +8,9 @@ date
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 64
+set ncpus = 8
 
-set maindir = /p/lustre2/cheung26/mgmol-20250206
+set maindir = /p/lustre2/cheung26/mgmol-20250219
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
@@ -21,7 +21,7 @@ source $maindir/scripts/modules.quartz
 
 rm -rf snapshot0_*
 
-set coord_file = coords.in_rotate2
+set coord_file = coords_rotate2.in
 
 set exe     = mgmol-opt
 set cfg_online = mgmol_ref_rotate.cfg
diff --git a/tests/PinnedH2O_3DOF/mgmol_online.cfg b/tests/PinnedH2O_3DOF/mgmol_online.cfg
index 5bc5b8b8..6f66067d 100644
--- a/tests/PinnedH2O_3DOF/mgmol_online.cfg
+++ b/tests/PinnedH2O_3DOF/mgmol_online.cfg
@@ -30,7 +30,7 @@ output_level=4
 solver=MVP
 nb_inner_it=100
 [ROM.offline]
-basis_file=data/PinnedH2O_3DOF_orbitals_basis_CUSTOM_NL_CUSTOM_NT
+basis_file=/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_8/PinnedH2O_3DOF_orbitals_basis_CUSTOM_NL_CUSTOM_NT
 [ROM.basis]
 compare_md=false
 number_of_orbital_basis=CUSTOM_NB
diff --git a/tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg b/tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg
index 72288ac3..8c598248 100644
--- a/tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg
+++ b/tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg
@@ -23,14 +23,14 @@ atol=1.e-8
 [Orbitals]
 initial_type=Random
 initial_width=2.
-nempty=68
+nempty=30
 [Restart]
 output_level=4
 [DensityMatrix]
 solver=MVP
 nb_inner_it=100
 [ROM.offline]
-basis_file=data/PinnedH2O_3DOF_orbitals_basis_2_2
+basis_file=/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_8/PinnedH2O_3DOF_orbitals_basis_2_2
 [ROM.basis]
 compare_md=false
-number_of_orbital_basis=72
+number_of_orbital_basis=34
diff --git a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
index a0cd51d3..83bfe8a6 100644
--- a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
+++ b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
@@ -14,6 +14,7 @@
 #include "MGmol_MPI.h"
 #include "MPIdata.h"
 #include "mgmol_run.h"
+#include "PinnedH2O.h"
 
 #ifdef MGMOL_HAS_LIBROM
 #include "librom.h"
@@ -29,215 +30,6 @@
 #include <boost/program_options.hpp>
 namespace po = boost::program_options;
 
-double calculate_bondlength(const double atom1[3], const double atom2[3])
-{
-    return sqrt(pow(atom1[0] - atom2[0], 2) + pow(atom1[1] - atom2[1], 2) + pow(atom1[2] - atom2[2], 2));
-}
-
-double calculate_bondangle(const double atom1[3], const double atom2[3], const double atom3[3])
-{
-    double vector1[3] = {atom1[0] - atom2[0], atom1[1] - atom2[1], atom1[2] - atom2[2]};
-    double vector2[3] = {atom3[0] - atom2[0], atom3[1] - atom2[1], atom3[2] - atom2[2]};
-
-    double dot_product = vector1[0] * vector2[0] + vector1[1] * vector2[1] + vector1[2] * vector2[2];
-    double magnitude_product = sqrt(pow(vector1[0], 2) + pow(vector1[1], 2) + pow(vector1[2], 2)) *
-                               sqrt(pow(vector2[0], 2) + pow(vector2[1], 2) + pow(vector2[2], 2));
-    double angle = acos(dot_product / magnitude_product);
-
-    return angle;
-}
-
-void rotation_matrix(const double axis[3], double angle, double matrix[3][3])
-{
-    double cos_theta = cos(angle);
-    double sin_theta = sin(angle);
-    double ux = axis[0], uy = axis[1], uz = axis[2];
-
-    matrix[0][0] = cos_theta + ux * ux * (1 - cos_theta);
-    matrix[0][1] = ux * uy * (1 - cos_theta) - uz * sin_theta;
-    matrix[0][2] = ux * uz * (1 - cos_theta) + uy * sin_theta;
-
-    matrix[1][0] = uy * ux * (1 - cos_theta) + uz * sin_theta;
-    matrix[1][1] = cos_theta + uy * uy * (1 - cos_theta);
-    matrix[1][2] = uy * uz * (1 - cos_theta) - ux * sin_theta;
-
-    matrix[2][0] = uz * ux * (1 - cos_theta) - uy * sin_theta;
-    matrix[2][1] = uz * uy * (1 - cos_theta) + ux * sin_theta;
-    matrix[2][2] = cos_theta + uz * uz * (1 - cos_theta);
-}
-
-void normalize(double vec[3])
-{
-    double norm = sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
-    vec[0] /= norm;
-    vec[1] /= norm;
-    vec[2] /= norm;
-}
-
-void cross(const double a[3], const double b[3], double result[3])
-{
-    result[0] = a[1] * b[2] - a[2] * b[1];
-    result[1] = a[2] * b[0] - a[0] * b[2];
-    result[2] = a[0] * b[1] - a[1] * b[0];
-}
-
-void apply_rotation(const double matrix[3][3], const double vec[3], double result[3])
-{
-    result[0] = matrix[0][0] * vec[0] + matrix[0][1] * vec[1] + matrix[0][2] * vec[2];
-    result[1] = matrix[1][0] * vec[0] + matrix[1][1] * vec[1] + matrix[1][2] * vec[2];
-    result[2] = matrix[2][0] * vec[0] + matrix[2][1] * vec[1] + matrix[2][2] * vec[2];
-}
-
-void apply_transpose_rotation(const double matrix[3][3], const double vec[3], double result[3])
-{
-    result[0] = matrix[0][0] * vec[0] + matrix[1][0] * vec[1] + matrix[2][0] * vec[2];
-    result[1] = matrix[0][1] * vec[0] + matrix[1][1] * vec[1] + matrix[2][1] * vec[2];
-    result[2] = matrix[0][2] * vec[0] + matrix[1][2] * vec[1] + matrix[2][2] * vec[2];
-}
-
-void rotate_PinnedH2O(std::vector<double>& positions, std::vector<short>& anumbers,
-                      double out_of_plane_rotation_matrix[3][3], double& planar_rotation_angle, bool& flipped_bond)
-{
-    int O1_idx = -1;
-    for (int i = 0; i < 3; i++)
-    {
-        if (positions[3*i] == 0.0 && positions[3*i+1] == 0.0 && positions[3*i+2] == 0.0)
-        {
-            O1_idx = i;
-            break;
-        }
-    }
-    if (O1_idx == -1) return;
-
-    int H1_idx = (O1_idx + 1) % 3;
-    int H2_idx = (O1_idx + 2) % 3;
-
-    double O1[3] = {positions[3*O1_idx], positions[3*O1_idx+1], positions[3*O1_idx+2]};
-    double H1[3] = {positions[3*H1_idx], positions[3*H1_idx+1], positions[3*H1_idx+2]};
-    double H2[3] = {positions[3*H2_idx], positions[3*H2_idx+1], positions[3*H2_idx+2]};
-
-    double bondlength1 = calculate_bondlength(H1, O1);
-    double bondlength2 = calculate_bondlength(H2, O1);
-    double bondangle = calculate_bondangle(H1, O1, H2);
-
-    double H1_temp[3], H2_temp[3];
-    double H1_rotated[3], H2_rotated[3];
-
-    double plane_normal[3];
-    cross(H2, H1, plane_normal);
-    normalize(plane_normal);
-    double target_plane_normal[3] = {0, 0, 1};
-    double dot_product = plane_normal[0] * target_plane_normal[0] +
-                         plane_normal[1] * target_plane_normal[1] +
-                         plane_normal[2] * target_plane_normal[2];
-    double angle_to_align = std::acos(std::min(std::max(dot_product, -1.0), 1.0));
-    double axis_to_align[3];
-    if (abs(dot_product) > 1.0 - 1e-8)
-    {
-        axis_to_align[0] = 1.0; 
-        axis_to_align[1] = 0.0;
-        axis_to_align[2] = 0.0;
-    }
-    else
-    {
-        cross(plane_normal, target_plane_normal, axis_to_align);
-        normalize(axis_to_align);
-    }
-    rotation_matrix(axis_to_align, angle_to_align, out_of_plane_rotation_matrix);
-    apply_rotation(out_of_plane_rotation_matrix, H1, H1_temp);
-    apply_rotation(out_of_plane_rotation_matrix, H2, H2_temp);
-
-    double theta1 = std::atan2(H1_temp[1], H1_temp[0]);
-    planar_rotation_angle = -theta1 + bondangle / 2.0;
-    double planar_rotation_matrix[3][3];
-    rotation_matrix(target_plane_normal, planar_rotation_angle, planar_rotation_matrix);
-    apply_rotation(planar_rotation_matrix, H1_temp, H1_rotated);
-    apply_rotation(planar_rotation_matrix, H2_temp, H2_rotated);
-
-    flipped_bond = (bondlength1 < bondlength2);
-    if (flipped_bond)
-    {
-        H1_rotated[1] *= -1.0;
-        H2_rotated[1] *= -1.0;
-        std::swap(H1_rotated, H2_rotated);
-    }
-
-    positions[0] = H2_rotated[0];
-    positions[1] = H2_rotated[1];
-    positions[2] = H2_rotated[2];
-    positions[3] = 0.0;
-    positions[4] = 0.0;
-    positions[5] = 0.0;
-    positions[6] = H1_rotated[0];
-    positions[7] = H1_rotated[1];
-    positions[8] = H1_rotated[2];
-
-    anumbers[0] = 1;
-    anumbers[1] = 8;
-    anumbers[2] = 1;
-}
-
-void transpose_rotate_PinnedH2O(std::vector<double>& positions, std::vector<double>& forces,
-                                const double out_of_plane_rotation_matrix[3][3], 
-                                const double planar_rotation_angle, const bool flipped_bond)
-{
-    double H2_rotated[3] = {positions[0], positions[1], positions[2]};
-    double O1_rotated[3] = {positions[3], positions[4], positions[5]};
-    double H1_rotated[3] = {positions[6], positions[7], positions[8]};
-
-    double f_H2_rotated[3] = {forces[0], forces[1], forces[2]};
-    double f_O1_rotated[3] = {forces[3], forces[4], forces[5]};
-    double f_H1_rotated[3] = {forces[6], forces[7], forces[8]};
-
-    if (flipped_bond)
-    {
-        H1_rotated[1] *= -1.0;
-        H2_rotated[1] *= -1.0;
-        f_O1_rotated[1] *= -1.0;
-        f_H1_rotated[1] *= -1.0;
-        f_H2_rotated[1] *= -1.0;
-    }
-
-    double planar_rotation_matrix[3][3];
-    double target_plane_normal[3] = {0, 0, 1};
-    rotation_matrix(target_plane_normal, planar_rotation_angle, planar_rotation_matrix);
-
-    double H1_temp[3], H2_temp[3];
-    apply_transpose_rotation(planar_rotation_matrix, H1_rotated, H1_temp);
-    apply_transpose_rotation(planar_rotation_matrix, H2_rotated, H2_temp);
-
-    double f_O1_temp[3], f_H1_temp[3], f_H2_temp[3];
-    apply_transpose_rotation(planar_rotation_matrix, f_O1_rotated, f_O1_temp);
-    apply_transpose_rotation(planar_rotation_matrix, f_H1_rotated, f_H1_temp);
-    apply_transpose_rotation(planar_rotation_matrix, f_H2_rotated, f_H2_temp);
-
-    double H1_restored[3], H2_restored[3];
-    apply_transpose_rotation(out_of_plane_rotation_matrix, H1_temp, H1_restored);
-    apply_transpose_rotation(out_of_plane_rotation_matrix, H2_temp, H2_restored);
-
-    double f_O1_restored[3], f_H1_restored[3], f_H2_restored[3];
-    apply_transpose_rotation(out_of_plane_rotation_matrix, f_O1_temp, f_O1_restored);
-    apply_transpose_rotation(out_of_plane_rotation_matrix, f_H1_temp, f_H1_restored);
-    apply_transpose_rotation(out_of_plane_rotation_matrix, f_H2_temp, f_H2_restored);
-
-    positions[0] = H2_restored[0];
-    positions[1] = H2_restored[1];
-    positions[2] = H2_restored[2];
-    positions[6] = H1_restored[0];
-    positions[7] = H1_restored[1];
-    positions[8] = H1_restored[2];
-
-    forces[0] = f_H2_restored[0];
-    forces[1] = f_H2_restored[1];
-    forces[2] = f_H2_restored[2];
-    forces[3] = f_O1_restored[0];
-    forces[4] = f_O1_restored[1];
-    forces[5] = f_O1_restored[2];
-    forces[6] = f_H1_restored[0];
-    forces[7] = f_H1_restored[1];
-    forces[8] = f_H1_restored[2];
-}
-
 int main(int argc, char** argv)
 {
     int mpirc = MPI_Init(&argc, &argv);
@@ -340,11 +132,8 @@ int main(int argc, char** argv)
         }
 
         // rotate the molecule to the reference coordinate system
-        int atom_order[3];
-        double out_of_plane_rotation_matrix[3][3];
-        double planar_rotation_angle;
-        bool flipped_bond;
-        rotate_PinnedH2O(positions, anumbers, out_of_plane_rotation_matrix, planar_rotation_angle, flipped_bond);
+        PinnedH2O H2O_molecule;
+        H2O_molecule.rotate(positions, anumbers);
 
         Mesh* mymesh             = Mesh::instance();
         const pb::Grid& mygrid   = mymesh->grid();
@@ -417,7 +206,7 @@ int main(int argc, char** argv)
         }
 
         // rotate the forces to the original coordinate system
-        transpose_rotate_PinnedH2O(positions, forces, out_of_plane_rotation_matrix, planar_rotation_angle, flipped_bond);
+        H2O_molecule.transpose_rotate(positions, anumbers, forces);
 
         // print out results
         if (MPIdata::onpe0)
diff --git a/tests/PinnedH2O_3DOF/transfer_nlrom.sh b/tests/PinnedH2O_3DOF/transfer_nlrom.sh
new file mode 100644
index 00000000..62824de2
--- /dev/null
+++ b/tests/PinnedH2O_3DOF/transfer_nlrom.sh
@@ -0,0 +1,20 @@
+#!/bin/bash
+
+ncpus=8
+target_dir="/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_${ncpus}"
+
+if [ ! -d $target_dir ]; then
+  mkdir -p $target_dir
+fi
+
+for dir in data/*/; do
+  dirname=$(basename "$dir")
+  new_dir="$target_dir/${dirname}"
+  if [ ! -d $new_dir ]; then
+    mkdir -p $new_dir
+  fi
+  mv -f $dir/* $new_dir/
+  rm -rf $dir
+done
+
+rm -rf data