diff --git a/examples/PinnedH2O/mgmol.cfg b/examples/PinnedH2O/mgmol.cfg
index 32325e47..ef82f07d 100644
--- a/examples/PinnedH2O/mgmol.cfg
+++ b/examples/PinnedH2O/mgmol.cfg
@@ -1,6 +1,6 @@
 verbosity=1
 xcFunctional=PBE
-FDtype=Mehrstellen
+FDtype=4th
 [Mesh]
 nx=64
 ny=64
diff --git a/examples/PinnedH2O/mgmol_offline.cfg b/examples/PinnedH2O/mgmol_offline.cfg
index 332dd77d..f535cff2 100644
--- a/examples/PinnedH2O/mgmol_offline.cfg
+++ b/examples/PinnedH2O/mgmol_offline.cfg
@@ -1,6 +1,6 @@
 verbosity=1
 xcFunctional=PBE
-FDtype=Mehrstellen
+FDtype=4th
 [Mesh]
 nx=64
 ny=64
diff --git a/tests/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in b/tests/PinnedH2O_3DOF/aux/coords_1.00_1.00_0.0.in
similarity index 100%
rename from tests/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in
rename to tests/PinnedH2O_3DOF/aux/coords_1.00_1.00_0.0.in
diff --git a/tests/PinnedH2O_3DOF/coords_1.02_0.98_2.0.in b/tests/PinnedH2O_3DOF/aux/coords_1.02_0.98_2.0.in
similarity index 100%
rename from tests/PinnedH2O_3DOF/coords_1.02_0.98_2.0.in
rename to tests/PinnedH2O_3DOF/aux/coords_1.02_0.98_2.0.in
diff --git a/tests/PinnedH2O_3DOF/aux/rotation_test.cc b/tests/PinnedH2O_3DOF/aux/rotation_test.cc
index 0e5641ee..21cfd93a 100644
--- a/tests/PinnedH2O_3DOF/aux/rotation_test.cc
+++ b/tests/PinnedH2O_3DOF/aux/rotation_test.cc
@@ -5,11 +5,13 @@
 
 using namespace std;
 
-double calculate_bondlength(const double atom1[3], const double atom2[3]) {
+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], bool radian) {
+double calculate_bondangle(const double atom1[3], const double atom2[3], const double atom3[3], bool radian)
+{
     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]};
 
@@ -24,7 +26,8 @@ double calculate_bondangle(const double atom1[3], const double atom2[3], const d
     return angle;
 }
 
-void rotation_matrix(const double axis[3], double angle, double matrix[3][3]) {
+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];
@@ -42,54 +45,41 @@ void rotation_matrix(const double axis[3], double angle, double matrix[3][3]) {
     matrix[2][2] = cos_theta + uz * uz * (1 - cos_theta);
 }
 
-void normalize(double vec[3]) {
+void normalize(double vec[3])
+{
     double norm = sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]);
-    if (norm > 0) {
-        vec[0] /= norm;
-        vec[1] /= norm;
-        vec[2] /= norm;
-    }
+    vec[0] /= norm;
+    vec[1] /= norm;
+    vec[2] /= norm;
 }
 
-void cross(const double a[3], const double b[3], double result[3]) {
+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]) {
+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]) {
+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];
 }
 
-int main() {
+int main()
+{
     double O1[3] = {0.00, 0.00, 0.00};
     double H1[3] = {-0.45, -1.48, -0.97};
     double H2[3] = {-0.45, 1.42, -1.07};
 
-    double plane_normal[3];
-    cross(H2, H1, plane_normal);
-    normalize(plane_normal);
-
-    double target_plane_normal[3] = {0, 0, 1};
-    double axis_to_align[3];
-    cross(plane_normal, target_plane_normal, axis_to_align);
-    normalize(axis_to_align);
-    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 = acos(min(max(dot_product, -1.0), 1.0));
-
-    double rot_matrix_align_plane[3][3];
-    rotation_matrix(axis_to_align, angle_to_align, rot_matrix_align_plane);
-
     double bondlength1 = calculate_bondlength(H1, O1);
     double bondlength2 = calculate_bondlength(H2, O1);
     double bondangle = calculate_bondangle(H1, O1, H2, false);
@@ -101,41 +91,75 @@ int main() {
     cout << "Bondlength of O1-H2 = " << bondlength2 << endl;
     cout << "Angle between O1-H1 and O1-H2 = " << bondangle << endl;
 
+    double H1_temp[3], H2_temp[3];
     double H1_rotated[3], H2_rotated[3];
-    apply_rotation(rot_matrix_align_plane, H1, H1_rotated);
-    apply_rotation(rot_matrix_align_plane, H2, H2_rotated);
-    bool flipped_bond = false;
-    if (bondlength1 < bondlength2) {
-        flipped_bond = true;
+
+    double plane_normal[3];
+    cross(H2, H1, plane_normal);
+    normalize(plane_normal);
+    double target_plane_normal[3] = {0, 0, 1};
+    double axis_to_align[3];
+    cross(plane_normal, target_plane_normal, axis_to_align);
+    normalize(axis_to_align);
+    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 = acos(min(max(dot_product, -1.0), 1.0));
+    double out_of_plane_rotation_matrix[3][3];
+    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 = atan2(H1_temp[1], H1_temp[0]);
+    double planar_rotation_angle = -theta1 + bondangle / 360.0 * M_PI;
+    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);
+
+    if (bondlength1 < bondlength2)
+    {
         swap(H1_rotated, H2_rotated);
+        H1_rotated[1] *= -1.0;
+        H2_rotated[1] *= -1.0;
     }
-    bondlength1 = calculate_bondlength(H1_rotated, O1);
-    bondlength2 = calculate_bondlength(H2_rotated, O1);
-    bondangle = calculate_bondangle(H1_rotated, O1, H2_rotated, false);
 
-    cout << "Reference system (z=0 plane about x=0 axis, with longer bondlength in H1)" << endl;
+    double bondlength1_rotated = calculate_bondlength(H1_rotated, O1);
+    double bondlength2_rotated = calculate_bondlength(H2_rotated, O1);
+    double bondangle_rotated = calculate_bondangle(H1_rotated, O1, H2_rotated, false);
+
+    cout << "Reference system (in z=0 plane, symmetric about x=0 axis, with longer bondlength in Q1)" << endl;
     cout << "H1 = (" << H1_rotated[0] << ", " << H1_rotated[1] << ", " << H1_rotated[2] << ")" << endl;
     cout << "H2 = (" << H2_rotated[0] << ", " << H2_rotated[1] << ", " << H2_rotated[2] << ")" << endl;
-    cout << "Bondlength of O1-H1 = " << bondlength1 << endl;
-    cout << "Bondlength of O1-H2 = " << bondlength2 << endl;
-    cout << "Angle between O1-H1 and O1-H2 = " << bondangle << endl;
+    cout << "Bondlength of O1-H1 = " << bondlength1_rotated << endl;
+    cout << "Bondlength of O1-H2 = " << bondlength2_rotated << endl;
+    cout << "Angle between O1-H1 and O1-H2 = " << bondangle_rotated << endl;
+
+    double H1_restored[3], H2_restored[3];
 
-    if (flipped_bond) {
+    if (bondlength1 < bondlength2)
+    {
         swap(H1_rotated, H2_rotated);
+        H1_rotated[1] *= -1.0;
+        H2_rotated[1] *= -1.0;
     }
-    double H1_restored[3], H2_restored[3];
-    apply_transpose_rotation(rot_matrix_align_plane, H1_rotated, H1_restored);
-    apply_transpose_rotation(rot_matrix_align_plane, H2_rotated, H2_restored);
-    bondlength1 = calculate_bondlength(H1_restored, O1);
-    bondlength2 = calculate_bondlength(H2_restored, O1);
-    bondangle = calculate_bondangle(H1_restored, O1, H2_restored, false);
+
+    apply_transpose_rotation(planar_rotation_matrix, H1_rotated, H1_temp);
+    apply_transpose_rotation(planar_rotation_matrix, H2_rotated, H2_temp);
+
+    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 bondlength1_restored = calculate_bondlength(H1_restored, O1);
+    double bondlength2_restored = calculate_bondlength(H2_restored, O1);
+    double bondangle_restored = calculate_bondangle(H1_restored, O1, H2_restored, false);
 
     cout << "Restored system" << endl;
     cout << "H1 = (" << H1_restored[0] << ", " << H1_restored[1] << ", " << H1_restored[2] << ")" << endl;
     cout << "H2 = (" << H2_restored[0] << ", " << H2_restored[1] << ", " << H2_restored[2] << ")" << endl;
-    cout << "Bondlength of O1-H1 = " << bondlength1 << endl;
-    cout << "Bondlength of O1-H2 = " << bondlength2 << endl;
-    cout << "Angle between O1-H1 and O1-H2 = " << bondangle << endl;
+    cout << "Bondlength of O1-H1 = " << bondlength1_restored << endl;
+    cout << "Bondlength of O1-H2 = " << bondlength2_restored << endl;
+    cout << "Angle between O1-H1 and O1-H2 = " << bondangle_restored << endl;
 
     return 0;
 }
diff --git a/tests/PinnedH2O_3DOF/aux/rotation_test.py b/tests/PinnedH2O_3DOF/aux/rotation_test.py
index fbecb532..41e20296 100644
--- a/tests/PinnedH2O_3DOF/aux/rotation_test.py
+++ b/tests/PinnedH2O_3DOF/aux/rotation_test.py
@@ -17,6 +17,17 @@ def calculate_bondangle(atom1, atom2, atom3, radian):
         angle = np.degrees(angle)
     return angle
 
+bondlength1 = calculate_bondlength(H1, O1)
+bondlength2 = calculate_bondlength(H2, O1)
+bondangle = calculate_bondangle(H1, O1, H2, False)
+
+print('Original system')
+print(f'H1 = {H1}')
+print(f'H2 = {H2}')
+print(f'Bondlength of O1-H1 = {bondlength1}')
+print(f'Bondlength of O1-H2 = {bondlength2}')
+print(f'Angle between O1-H1 and O1-H2 = {bondangle}')
+
 def rotation_matrix(axis, angle):
     cos_theta = np.cos(angle)
     sin_theta = np.sin(angle)
@@ -27,6 +38,8 @@ def rotation_matrix(axis, angle):
         [uz * ux * (1 - cos_theta) - uy * sin_theta, uz * uy * (1 - cos_theta) + ux * sin_theta, cos_theta + uz**2 * (1 - cos_theta)]
     ])
 
+H1_rotated, H2_rotated = H1, H2
+
 plane_normal = np.cross(H2, H1)
 plane_normal = plane_normal / np.linalg.norm(plane_normal)
 
@@ -34,47 +47,53 @@ def rotation_matrix(axis, angle):
 axis_to_align = np.cross(plane_normal, target_plane_normal)
 axis_to_align /= np.linalg.norm(axis_to_align)
 angle_to_align = np.arccos(np.clip(np.dot(plane_normal, target_plane_normal), -1.0, 1.0))
-rot_matrix_align_plane = rotation_matrix(axis_to_align, angle_to_align)
-
-bondlength1 = calculate_bondlength(H1, O1)
-bondlength2 = calculate_bondlength(H2, O1)
-bondangle = calculate_bondangle(H1, O1, H2, False)
+out_of_plane_rotation_matrix = rotation_matrix(axis_to_align, angle_to_align)
+H1_rotated = np.dot(out_of_plane_rotation_matrix, H1)
+H2_rotated = np.dot(out_of_plane_rotation_matrix, H2)
 
-print('Original system')
-print(f'H1 = {H1}')
-print(f'H2 = {H2}')
-print(f'Bondlength of O1-H1 = {bondlength1}')
-print(f'Bondlength of O1-H2 = {bondlength2}')
-print(f'Angle between O1-H1 and O1-H2 = {bondangle}')
+theta1 = np.arctan2(H1_rotated[1], H1_rotated[0]) 
+theta2 = np.arctan2(H2_rotated[1], H2_rotated[0])
+planar_rotation_angle = -theta1 + np.radians(bondangle) / 2 
+planar_rotation_matrix = rotation_matrix(target_plane_normal, planar_rotation_angle)
+H1_rotated = np.dot(planar_rotation_matrix, H1_rotated)
+H2_rotated = np.dot(planar_rotation_matrix, H2_rotated)
 
-H1_rotated = np.dot(rot_matrix_align_plane, H1)
-H2_rotated = np.dot(rot_matrix_align_plane, H2)
-fliped_bond = False
 if bondlength1 < bondlength2:
-    fliped_bond = True
-    H1_rotated, H2_rotated = H2_rotated, H1_rotated
-bondlength1 = calculate_bondlength(H1_rotated, O1)
-bondlength2 = calculate_bondlength(H2_rotated, O1)
-bondangle = calculate_bondangle(H1_rotated, O1, H2_rotated, False)
+    H1_rotated, H2_rotated = H2_rotated, H1_rotated 
+    H1_rotated[1] *= -1
+    H2_rotated[1] *= -1
 
-print('Reference system (z=0 plane about x=0 axis, with longer bondlength in H1)')
+bondlength1_rotated = calculate_bondlength(H1_rotated, O1)
+bondlength2_rotated = calculate_bondlength(H2_rotated, O1)
+bondangle_rotated = calculate_bondangle(H1_rotated, O1, H2_rotated, False)
+
+print('Reference system (z=0 plane, symmetric about x=0 axis, with longer bondlength in Q1)')
 print(f'H1 = {H1_rotated}')
 print(f'H2 = {H2_rotated}')
-print(f'Bondlength of O1-H1 = {bondlength1}')
-print(f'Bondlength of O1-H2 = {bondlength2}')
-print(f'Angle between O1-H1 and O1-H2 = {bondangle}')
+print(f'Bondlength of O1-H1 = {bondlength1_rotated}')
+print(f'Bondlength of O1-H2 = {bondlength2_rotated}')
+print(f'Angle between O1-H1 and O1-H2 = {bondangle_rotated}')
 
-if fliped_bond:
-    H1_rotated, H2_rotated = H2_rotated, H1_rotated
-H1_restored = np.dot(rot_matrix_align_plane.T, H1_rotated)
-H2_restored = np.dot(rot_matrix_align_plane.T, H2_rotated)
-bondlength1 = calculate_bondlength(H1_restored, O1)
-bondlength2 = calculate_bondlength(H2_restored, O1)
-bondangle = calculate_bondangle(H1_restored, O1, H2_restored, False)
+H1_restored, H2_restored = H1_rotated, H2_rotated
+
+if bondlength1 < bondlength2:
+    H1_rotated, H2_rotated = H2_rotated, H1_rotated 
+    H1_rotated[1] *= -1
+    H2_rotated[1] *= -1
+
+H1_restored = np.dot(planar_rotation_matrix.T, H1_restored)
+H2_restored = np.dot(planar_rotation_matrix.T, H2_restored)
+
+H1_restored = np.dot(out_of_plane_rotation_matrix.T, H1_restored)
+H2_restored = np.dot(out_of_plane_rotation_matrix.T, H2_restored)
+
+bondlength1_restored = calculate_bondlength(H1_restored, O1)
+bondlength2_restored = calculate_bondlength(H2_restored, O1)
+bondangle_restored = calculate_bondangle(H1_restored, O1, H2_restored, False)
 
 print('Restored system')
 print(f'H1 = {H1_restored}')
 print(f'H2 = {H2_restored}')
-print(f'Bondlength of O1-H1 = {bondlength1}')
-print(f'Bondlength of O1-H2 = {bondlength2}')
-print(f'Angle between O1-H1 and O1-H2 = {bondangle}')
+print(f'Bondlength of O1-H1 = {bondlength1_restored}')
+print(f'Bondlength of O1-H2 = {bondlength2_restored}')
+print(f'Angle between O1-H1 and O1-H2 = {bondangle_restored}')
diff --git a/tests/PinnedH2O_3DOF/coords.in_rotate1 b/tests/PinnedH2O_3DOF/coords.in_rotate1
new file mode 100644
index 00000000..7498ff55
--- /dev/null
+++ b/tests/PinnedH2O_3DOF/coords.in_rotate1
@@ -0,0 +1,3 @@
+O1  1   0.00000  0.00000 0.00  0
+H1  2  -1.08248 -1.47039 0.00  1
+H2  2  -1.14768  1.43060 0.00  1
diff --git a/tests/PinnedH2O_3DOF/coords.in_rotate2 b/tests/PinnedH2O_3DOF/coords.in_rotate2
new file mode 100644
index 00000000..307e170c
--- /dev/null
+++ b/tests/PinnedH2O_3DOF/coords.in_rotate2
@@ -0,0 +1,3 @@
+O1  1   0.00    0.00    0.00  0
+H1  2  -0.45    1.42   -1.07  1
+H2  2  -0.45   -1.48   -0.97  1
diff --git a/tests/PinnedH2O_3DOF/job.basis b/tests/PinnedH2O_3DOF/job.basis
index ccc03c9d..eaf8bfde 100644
--- a/tests/PinnedH2O_3DOF/job.basis
+++ b/tests/PinnedH2O_3DOF/job.basis
@@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 64
 
-set maindir = /p/lustre2/cheung26/mgmol-20241220
+set maindir = /p/lustre2/cheung26/mgmol-20250206
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
diff --git a/tests/PinnedH2O_3DOF/job.basis_test1 b/tests/PinnedH2O_3DOF/job.basis_test1
deleted file mode 100644
index 3955cdc6..00000000
--- a/tests/PinnedH2O_3DOF/job.basis_test1
+++ /dev/null
@@ -1,34 +0,0 @@
-#!/bin/tcsh
-#SBATCH -N 2
-#SBATCH -t 0:10:00
-#SBATCH -p pdebug
-
-date
-
-setenv OMP_NUM_THREADS 1
-#setenv KMP_DETERMINISTIC_REDUCTION 1
-
-set ncpus = 64
-
-set maindir = /p/lustre2/cheung26/mgmol-20241220
-
-setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
-
-set exe     = ${maindir}/build_quartz/libROM/build/examples/misc/combine_samples
-
-set basis_tag = "test1"
-set snapshot_files = ""
-
-set bondlength_one = 1.00
-set bondlength_two = 1.00
-set bondangle = 0.0
-echo "bondlength1: $bondlength_one, bondlength2: $bondlength_two, bondangle: $bondangle"
-set tag = ${bondlength_one}_${bondlength_two}_${bondangle}
-set snapshot_files = "$snapshot_files data/${tag}/orbital_snapshot"
-
-echo "Snapshot files: $snapshot_files"
-
-set basis_file = "data/PinnedH2O_3DOF_orbitals_basis_${basis_tag}"
-srun -n $ncpus $exe -f $basis_file $snapshot_files > data/basis_${basis_tag}_PinnedH2O_3DOF.out
-
-date
diff --git a/tests/PinnedH2O_3DOF/job.basis_test2 b/tests/PinnedH2O_3DOF/job.basis_test2
deleted file mode 100644
index d0ac6301..00000000
--- a/tests/PinnedH2O_3DOF/job.basis_test2
+++ /dev/null
@@ -1,39 +0,0 @@
-#!/bin/tcsh
-#SBATCH -N 2
-#SBATCH -t 0:10:00
-#SBATCH -p pdebug
-
-date
-
-setenv OMP_NUM_THREADS 1
-#setenv KMP_DETERMINISTIC_REDUCTION 1
-
-set ncpus = 64
-
-set maindir = /p/lustre2/cheung26/mgmol-20241220
-
-setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
-
-set exe     = ${maindir}/build_quartz/libROM/build/examples/misc/combine_samples
-
-set basis_tag = "test2"
-set snapshot_files = ""
-
-set bondlength_one = 1.00
-set bondlength_two = 1.00
-set bondangle = 0.0
-echo "bondlength1: $bondlength_one, bondlength2: $bondlength_two, bondangle: $bondangle"
-set tag = ${bondlength_one}_${bondlength_two}_${bondangle}
-set snapshot_files = "$snapshot_files data/${tag}/orbital_snapshot"
-
-set bondangle = 1.0
-echo "bondlength1: $bondlength_one, bondlength2: $bondlength_two, bondangle: $bondangle"
-set tag = ${bondlength_one}_${bondlength_two}_${bondangle}
-set snapshot_files = "$snapshot_files data/${tag}/orbital_snapshot"
-
-echo "Snapshot files: $snapshot_files"
-
-set basis_file = "data/PinnedH2O_3DOF_orbitals_basis_${basis_tag}"
-srun -n $ncpus $exe -f $basis_file $snapshot_files > data/basis_${basis_tag}_PinnedH2O_3DOF.out
-
-date
diff --git a/tests/PinnedH2O_3DOF/job.offline b/tests/PinnedH2O_3DOF/job.offline
index 7fc51dd7..3870a3a3 100644
--- a/tests/PinnedH2O_3DOF/job.offline
+++ b/tests/PinnedH2O_3DOF/job.offline
@@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 64
 
-set maindir = /p/lustre2/cheung26/mgmol-20241220
+set maindir = /p/lustre2/cheung26/mgmol-20250206
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
diff --git a/tests/PinnedH2O_3DOF/job.online b/tests/PinnedH2O_3DOF/job.online
index e703fc29..77b8fde0 100644
--- a/tests/PinnedH2O_3DOF/job.online
+++ b/tests/PinnedH2O_3DOF/job.online
@@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 64
 
-set maindir = /p/lustre2/cheung26/mgmol-20241220
+set maindir = /p/lustre2/cheung26/mgmol-20250206
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
diff --git a/tests/PinnedH2O_3DOF/job.online_test1 b/tests/PinnedH2O_3DOF/job.online_rotate
similarity index 55%
rename from tests/PinnedH2O_3DOF/job.online_test1
rename to tests/PinnedH2O_3DOF/job.online_rotate
index bf69a1ef..11dd7419 100644
--- a/tests/PinnedH2O_3DOF/job.online_test1
+++ b/tests/PinnedH2O_3DOF/job.online_rotate
@@ -10,21 +10,27 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 64
 
-set maindir = /p/lustre2/cheung26/mgmol-20241220
+set maindir = /p/lustre2/cheung26/mgmol-20250206
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
-set exe     = testPinnedH2O_3DOF
-
-cp $maindir/build_quartz/tests/$exe .
-
-set cfg_online = mgmol_online_test1.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 -p pdebug -n $ncpus $exe -c $cfg_online -i coords_1.00_1.00_0.0.in
+rm -rf snapshot0_*
+
+set coord_file = coords.in_rotate2
+
+set exe     = mgmol-opt
+set cfg_online = mgmol_ref_rotate.cfg
+cp $maindir/install_quartz/bin/$exe .
+srun -p pdebug -n $ncpus $exe -c $cfg_online -i $coord_file
+
+set exe     = testPinnedH2O_3DOF
+set cfg_online = mgmol_online_rotate.cfg
+cp $maindir/build_quartz/tests/$exe .
+srun -p pdebug -n $ncpus $exe -c $cfg_online -i $coord_file
 
 date
diff --git a/tests/PinnedH2O_3DOF/job.online_test2 b/tests/PinnedH2O_3DOF/job.online_test2
deleted file mode 100644
index 88665f1e..00000000
--- a/tests/PinnedH2O_3DOF/job.online_test2
+++ /dev/null
@@ -1,30 +0,0 @@
-#!/bin/tcsh
-#SBATCH -N 2
-#SBATCH -t 0:10:00
-#SBATCH -p pdebug
-
-date
-
-setenv OMP_NUM_THREADS 1
-#setenv KMP_DETERMINISTIC_REDUCTION 1
-
-set ncpus = 64
-
-set maindir = /p/lustre2/cheung26/mgmol-20241220
-
-setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
-
-set exe     = testPinnedH2O_3DOF
-
-cp $maindir/build_quartz/tests/$exe .
-
-set cfg_online = mgmol_online_test2.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 -p pdebug -n $ncpus $exe -c $cfg_online -i coords_1.00_1.00_0.0.in
-
-date
diff --git a/tests/PinnedH2O_3DOF/mgmol_online_test2.cfg b/tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg
similarity index 76%
rename from tests/PinnedH2O_3DOF/mgmol_online_test2.cfg
rename to tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg
index 8ce3e32d..72288ac3 100644
--- a/tests/PinnedH2O_3DOF/mgmol_online_test2.cfg
+++ b/tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg
@@ -18,19 +18,19 @@ pseudopotential=pseudo.H_ONCV_PBE_SG15
 [Run]
 type=quench
 [Quench]
-max_steps=1000
+max_steps=100
 atol=1.e-8
 [Orbitals]
 initial_type=Random
 initial_width=2.
-nempty=4
+nempty=68
 [Restart]
 output_level=4
 [DensityMatrix]
 solver=MVP
-nb_inner_it=20
+nb_inner_it=100
 [ROM.offline]
-basis_file=data/PinnedH2O_3DOF_orbitals_basis_test2
+basis_file=data/PinnedH2O_3DOF_orbitals_basis_2_2
 [ROM.basis]
 compare_md=false
-number_of_orbital_basis=8
+number_of_orbital_basis=72
diff --git a/tests/PinnedH2O_3DOF/mgmol_online_test1.cfg b/tests/PinnedH2O_3DOF/mgmol_ref_rotate.cfg
similarity index 70%
rename from tests/PinnedH2O_3DOF/mgmol_online_test1.cfg
rename to tests/PinnedH2O_3DOF/mgmol_ref_rotate.cfg
index e032350d..96be6b8d 100644
--- a/tests/PinnedH2O_3DOF/mgmol_online_test1.cfg
+++ b/tests/PinnedH2O_3DOF/mgmol_ref_rotate.cfg
@@ -1,4 +1,4 @@
-verbosity=2
+verbosity=1
 xcFunctional=PBE
 FDtype=4th
 [Mesh]
@@ -25,8 +25,3 @@ initial_type=Random
 initial_width=2.
 [Restart]
 output_level=4
-[ROM.offline]
-basis_file=data/PinnedH2O_3DOF_orbitals_basis_test1
-[ROM.basis]
-compare_md=false
-number_of_orbital_basis=4
diff --git a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
index e7813777..a0cd51d3 100644
--- a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
+++ b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
@@ -22,10 +22,222 @@
 #include <iostream>
 #include <time.h>
 #include <vector>
+#include <cmath>
+#include <algorithm>
+#include <cstring>
 
 #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);
@@ -127,6 +339,13 @@ 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);
+
         Mesh* mymesh             = Mesh::instance();
         const pb::Grid& mygrid   = mymesh->grid();
         const pb::PEenv& myPEenv = mymesh->peenv();
@@ -172,13 +391,58 @@ int main(int argc, char** argv)
         if (MPIdata::onpe0)
         {
             std::cout << "Eks: " << eks << std::endl;
+            std::cout << "Positions in the reference domain:" << std::endl;
+            std::vector<short>::iterator ita = anumbers.begin();
+            for (std::vector<double>::iterator it = positions.begin();
+                 it != positions.end(); it += 3)
+            {
+                std::cout << *ita;
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << std::endl;
+                ita++;
+            }
+
+            std::cout << "Forces in the reference domain:" << std::endl;
+            ita = anumbers.begin();
+            for (std::vector<double>::iterator it = forces.begin();
+                 it != forces.end(); it += 3)
+            {
+                std::cout << *ita;
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << std::endl;
+                ita++;
+            }
+        }
+
+        // rotate the forces to the original coordinate system
+        transpose_rotate_PinnedH2O(positions, forces, out_of_plane_rotation_matrix, planar_rotation_angle, flipped_bond);
+
+        // print out results
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Positions:" << std::endl;
+            std::vector<short>::iterator ita = anumbers.begin();
+            for (std::vector<double>::iterator it = positions.begin();
+                 it != positions.end(); it += 3)
+            {
+                std::cout << *ita;
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << std::endl;
+                ita++;
+            }
             std::cout << "Forces:" << std::endl;
+            ita = anumbers.begin();
             for (std::vector<double>::iterator it = forces.begin();
                  it != forces.end(); it += 3)
             {
+                std::cout << *ita;
                 for (int i = 0; i < 3; i++)
                     std::cout << "    " << *(it + i);
                 std::cout << std::endl;
+                ita++;
             }
         }