diff --git a/examples/PinnedH2O/get_ROM_table.py b/examples/PinnedH2O/get_ROM_table.py
new file mode 100644
index 00000000..fa4289f3
--- /dev/null
+++ b/examples/PinnedH2O/get_ROM_table.py
@@ -0,0 +1,28 @@
+import subprocess
+import re
+
+pattern = r"For energy fraction: \d+\.\d+, take first (\d+) of \d+ basis vectors"
+
+print("\\begin{tabular}{|c||c|c|c|c|c|c|c|}")
+print("\\hline")
+print("$k$ & $\\varepsilon = 10^{-1}$ & $\\varepsilon = 10^{-2}$ & $\\varepsilon = 10^{-3}$ & $\\varepsilon = 10^{-4}$ & $\\varepsilon = 10^{-5}$ & Snapshots \\\\")
+print("\\hline")
+
+for t in range(10):
+    k = 50*(t+1)
+    snapshots = 4*k
+    grep_command = f"grep 'take first' basis_1_{k}_Pinned_H2O.out"
+    result = subprocess.run(grep_command, shell=True, capture_output=True, text=True)
+    matches = re.findall(pattern, result.stdout)
+    energy_fractions = {
+        "0.9": matches[0],
+        "0.99": matches[1],
+        "0.999": matches[2],
+        "0.9999": matches[3],
+        "0.99999": matches[4],
+    }
+    line = f"{k} & {energy_fractions['0.9']} & {energy_fractions['0.99']} & {energy_fractions['0.999']} & {energy_fractions['0.9999']} & {energy_fractions['0.99999']} & {snapshots} \\\\"
+    print(line)
+
+print("\\hline")
+print("\\end{tabular}")
diff --git a/examples/PinnedH2O_3DOF/PinnedH2O_3DOF.png b/examples/PinnedH2O_3DOF/PinnedH2O_3DOF.png
new file mode 100644
index 00000000..99135eb6
Binary files /dev/null and b/examples/PinnedH2O_3DOF/PinnedH2O_3DOF.png differ
diff --git a/examples/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in b/examples/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in
new file mode 100644
index 00000000..46f5681d
--- /dev/null
+++ b/examples/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in
@@ -0,0 +1,3 @@
+O1  1   0.00    0.00    0.00  0
+H1  2   1.12    1.45    0.00  1
+H2  2   1.12    -1.45    0.00  1
diff --git a/examples/PinnedH2O/generate_coord.py b/examples/PinnedH2O_3DOF/generate_coord.py
similarity index 100%
rename from examples/PinnedH2O/generate_coord.py
rename to examples/PinnedH2O_3DOF/generate_coord.py
diff --git a/examples/PinnedH2O/generate_coord_simple.py b/examples/PinnedH2O_3DOF/generate_coord_simple.py
similarity index 100%
rename from examples/PinnedH2O/generate_coord_simple.py
rename to examples/PinnedH2O_3DOF/generate_coord_simple.py
diff --git a/examples/PinnedH2O_3DOF/get_ROM_table.py b/examples/PinnedH2O_3DOF/get_ROM_table.py
new file mode 100644
index 00000000..5fea1e38
--- /dev/null
+++ b/examples/PinnedH2O_3DOF/get_ROM_table.py
@@ -0,0 +1,31 @@
+import subprocess
+import re
+
+bondlength_num_increments = (2, 5, 10)
+bondangle_num_increments = (2, 5, 10)
+
+pattern = r"For energy fraction: \d+\.\d+, take first (\d+) of \d+ basis vectors"
+
+print("\\begin{tabular}{|c|c||c|c|c|c|c|c|c|}")
+print("\\hline")
+print("$N_L$ & $N_\\theta$ & $\\varepsilon = 10^{-1}$ & $\\varepsilon = 10^{-2}$ & $\\varepsilon = 10^{-3}$ & $\\varepsilon = 10^{-4}$ & $\\varepsilon = 10^{-5}$ & Snapshots \\\\")
+print("\\hline")
+
+for _, N_L in enumerate(bondlength_num_increments): 
+    for _, N_theta in enumerate(bondangle_num_increments): 
+        snapshots = 2*(N_L+1)*(N_L+2)*(N_theta+1)
+        grep_command = f"grep 'take first' basis_{N_L}_{N_theta}_PinnedH2O_3DOF.out"
+        result = subprocess.run(grep_command, shell=True, capture_output=True, text=True)
+        matches = re.findall(pattern, result.stdout)
+        energy_fractions = {
+            "0.9": matches[0],
+            "0.99": matches[1],
+            "0.999": matches[2],
+            "0.9999": matches[3],
+            "0.99999": matches[4],
+        }
+        line = f"{N_L} & {N_theta} & {energy_fractions['0.9']} & {energy_fractions['0.99']} & {energy_fractions['0.999']} & {energy_fractions['0.9999']} & {energy_fractions['0.99999']} & {snapshots} \\\\"
+        print(line)
+
+print("\\hline")
+print("\\end{tabular}")
diff --git a/examples/PinnedH2O_3DOF/job.basis b/examples/PinnedH2O_3DOF/job.basis
new file mode 100644
index 00000000..30660259
--- /dev/null
+++ b/examples/PinnedH2O_3DOF/job.basis
@@ -0,0 +1,47 @@
+#!/bin/tcsh
+#SBATCH -N 2
+#SBATCH -t 0:10: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     = ${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 = 10
+set bondangle_min = -5.0 
+set bondangle_max = 5.0
+set bondangle_num_increments = 10
+
+foreach i (`seq 0 $bondlength_num_increments`)
+  set bondlength_one = `echo "scale=2; $bondlength_min + $i * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc`
+  set bondlength_one = `printf "%.2f" $bondlength_one` 
+  foreach j (`seq 0 $i`)
+    set bondlength_two = `echo "scale=2; $bondlength_min + $j * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc`
+    set bondlength_two = `printf "%.2f" $bondlength_two`
+    foreach k (`seq 0 $bondangle_num_increments`)
+      set bondangle = `echo "scale=1; $bondangle_min + $k * ($bondangle_max - $bondangle_min) / ($bondangle_num_increments)" | bc`
+      set bondangle = `printf "%.1f" $bondangle`
+      echo "bondlength1: $bondlength_one, bondlength2: $bondlength_two, bondangle: $bondangle"
+      set tag = ${bondlength_one}_${bondlength_two}_${bondangle}
+      set snapshot_files = "$snapshot_files results_${tag}/orbital_snapshot"
+    end
+  end
+end
+echo "Snapshot files: $snapshot_files"
+
+set basis_file = "PinnedH2O_3DOF_orbitals_basis_${bondlength_num_increments}_${bondangle_num_increments}"
+srun -n $ncpus $exe -f $basis_file $snapshot_files > basis_${bondlength_num_increments}_${bondangle_num_increments}_PinnedH2O_3DOF.out
+
+date
diff --git a/examples/PinnedH2O_3DOF/job.offline b/examples/PinnedH2O_3DOF/job.offline
new file mode 100644
index 00000000..e45fd5d5
--- /dev/null
+++ b/examples/PinnedH2O_3DOF/job.offline
@@ -0,0 +1,63 @@
+#!/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_3DOF
+
+set cfg_offline = mgmol_offline.cfg
+cp $datadir/$cfg_offline .
+
+# coords.in files will be generated by Python script
+#cp $datadir/generate_coords_simple.py .
+#python3 generate_coords_simple.py
+
+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
+
+set bondlength_min = 0.95
+set bondlength_max = 1.05
+set bondlength_num_increments = 10
+set bondangle_min = -5.0 
+set bondangle_max = 5.0
+set bondangle_num_increments = 10
+
+foreach i (`seq 0 $bondlength_num_increments`)
+  set bondlength_one = `echo "scale=2; $bondlength_min + $i * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc`
+  set bondlength_one = `printf "%.2f" $bondlength_one` 
+  foreach j (`seq 0 $i`)
+    set bondlength_two = `echo "scale=2; $bondlength_min + $j * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc`
+    set bondlength_two = `printf "%.2f" $bondlength_two`
+    foreach k (`seq 0 $bondangle_num_increments`)
+      set bondangle = `echo "scale=1; $bondangle_min + $k * ($bondangle_max - $bondangle_min) / ($bondangle_num_increments)" | bc`
+      set bondangle = `printf "%.1f" $bondangle`
+      echo "bondlength1: $bondlength_one, bondlength2: $bondlength_two, bondangle: $bondangle"
+      set tag = ${bondlength_one}_${bondlength_two}_${bondangle}
+      cp $cfg_offline ${tag}_${cfg_offline}
+      sed -i "s/CUSTOM_TAG/results_${tag}/g" ${tag}_${cfg_offline}
+      srun -n $ncpus $exe -c ${tag}_${cfg_offline} -i coords_${tag}.in > offline_PinnedH2O_${tag}.out
+      mv ${tag}_${cfg_offline} results_${tag}/${cfg_offline}
+      mv coords_${tag}.in results_${tag}/coords.in
+      mv offline_PinnedH2O_${tag}.out results_${tag}/offline_PinnedH2O.out
+    end
+  end
+end
+
+date
diff --git a/examples/PinnedH2O_3DOF/mgmol_offline.cfg b/examples/PinnedH2O_3DOF/mgmol_offline.cfg
new file mode 100644
index 00000000..42fabd17
--- /dev/null
+++ b/examples/PinnedH2O_3DOF/mgmol_offline.cfg
@@ -0,0 +1,34 @@
+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=quench
+[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
+output_filename=CUSTOM_TAG
+[ROM.offline]
+save_librom_snapshot=true
diff --git a/examples/PinnedH2O_3DOF/plot_PinnedH2O_3DOF.py b/examples/PinnedH2O_3DOF/plot_PinnedH2O_3DOF.py
new file mode 100644
index 00000000..9c859f0e
--- /dev/null
+++ b/examples/PinnedH2O_3DOF/plot_PinnedH2O_3DOF.py
@@ -0,0 +1,69 @@
+import matplotlib
+import matplotlib.pyplot as plt
+import numpy as np
+
+ref_bondlength = 1.83
+ref_bondangle = 104.5
+
+# factors and increments for bond lengths and bond angle
+bondlength_factor = np.linspace(0.95, 1.05, 3)
+bondangle_increment = np.linspace(-5, 5, 3)
+
+fig, ax = plt.subplots()
+
+for d_bondangle in bondangle_increment:
+    bondangle = ref_bondangle + d_bondangle
+    x = ref_bondlength * np.cos(np.radians(bondangle / 2))
+    y = ref_bondlength * np.sin(np.radians(bondangle / 2))
+    for i, f_bondlength1 in enumerate(bondlength_factor):
+        H1_x = f_bondlength1*x
+        H1_y = f_bondlength1*y
+        ax.plot(H1_x, H1_y, 'ro', markersize=8, alpha=0.1)
+        for f_bondlength2 in bondlength_factor[:(i+1)]:
+            H2_x = f_bondlength2*x
+            H2_y = -f_bondlength2*y
+            ax.plot(H2_x, H2_y, 'bo', markersize=8, alpha=0.1)
+
+ref_bondangle = np.radians(ref_bondangle)
+x = ref_bondlength * np.cos(ref_bondangle / 2)
+y = ref_bondlength * np.sin(ref_bondangle / 2)
+
+ax.plot([0, x], [0, y], 'r-', lw=2)
+ax.plot([0, x], [0, -y], 'b-', lw=2)
+ax.plot(0, 0, 'ko', markersize=8)
+ax.plot(x, y, 'ro', markersize=8)
+ax.plot(x, -y, 'bo', markersize=8)
+
+arc1 = matplotlib.patches.Arc(xy=(0, 0), 
+                             width=0.5,
+                             height=0.5, 
+                             angle=0, 
+                             theta1=0, 
+                             theta2=ref_bondangle/2,
+                             color='red',
+                             linewidth=2)
+ax.add_patch(arc1)
+ax.text(0.4, 0.25, r"$\frac{\theta}{2}$", ha='center', va='center', fontsize=12, color='red') 
+
+arc2 = matplotlib.patches.Arc(xy=(0, 0), 
+                             width=0.5,
+                             height=0.5, 
+                             angle=0, 
+                             theta1=-ref_bondangle/2,
+                             theta2=0, 
+                             color='blue',
+                             linewidth=2)
+ax.add_patch(arc2)
+ax.text(0.4, -0.25, r"$\frac{\theta}{2}$", ha='center', va='center', fontsize=12, color='blue') 
+
+ax.text(x / 2 - 0.1, y / 2, r"$L_1$", ha='right', color='red')
+ax.text(x / 2 - 0.1, -y / 2, r"$L_2$", ha='right', color='blue')
+
+ax.plot([-2, 2], [0, 0], 'k--', lw=1)
+ax.set_xlim(-2.0, 2.0)
+ax.set_ylim(-2.0, 2.0)
+ax.set_aspect('equal')
+ax.set_xlabel("x")
+ax.set_ylabel("y")
+ax.grid(True)
+plt.savefig("PinnedH2O_3DOF.png")
diff --git a/examples/PinnedH2O_3DOF/rotation_test.cc b/examples/PinnedH2O_3DOF/rotation_test.cc
new file mode 100644
index 00000000..0e5641ee
--- /dev/null
+++ b/examples/PinnedH2O_3DOF/rotation_test.cc
@@ -0,0 +1,141 @@
+#include <iostream>
+#include <cmath>
+#include <algorithm>
+#include <cstring>
+
+using namespace std;
+
+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 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);
+
+    if (!radian) {
+        angle = angle * 180.0 / M_PI;
+    }
+    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]);
+    if (norm > 0) {
+        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];
+}
+
+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);
+
+    cout << "Original system" << endl;
+    cout << "H1 = (" << H1[0] << ", " << H1[1] << ", " << H1[2] << ")" << endl;
+    cout << "H2 = (" << H2[0] << ", " << H2[1] << ", " << H2[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;
+
+    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;
+        swap(H1_rotated, H2_rotated);
+    }
+    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;
+    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;
+
+    if (flipped_bond) {
+        swap(H1_rotated, H2_rotated);
+    }
+    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);
+
+    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;
+
+    return 0;
+}
diff --git a/examples/PinnedH2O/rotation_test.py b/examples/PinnedH2O_3DOF/rotation_test.py
similarity index 73%
rename from examples/PinnedH2O/rotation_test.py
rename to examples/PinnedH2O_3DOF/rotation_test.py
index 79f1aaef..fbecb532 100644
--- a/examples/PinnedH2O/rotation_test.py
+++ b/examples/PinnedH2O_3DOF/rotation_test.py
@@ -1,8 +1,8 @@
 import numpy as np
 
 O1 = np.array([0.00, 0.00, 0.00])
-H1 = np.array([-0.45, 1.42, -1.07])
-H2 = np.array([-0.45, -1.48, -0.97])
+H1 = np.array([-0.45, -1.48, -0.97])
+H2 = np.array([-0.45, 1.42, -1.07])
 
 def calculate_bondlength(atom1, atom2):
     return np.linalg.norm(atom1 - atom2)
@@ -34,10 +34,7 @@ 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)
-H1_rotated = np.dot(rot_matrix_align_plane, H1)
-H2_rotated = np.dot(rot_matrix_align_plane, H2)
 
 bondlength1 = calculate_bondlength(H1, O1)
 bondlength2 = calculate_bondlength(H2, O1)
@@ -50,15 +47,34 @@ def rotation_matrix(axis, angle):
 print(f'Bondlength of O1-H2 = {bondlength2}')
 print(f'Angle between O1-H1 and O1-H2 = {bondangle}')
 
+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)
-if bondlength1 < bondlength2:
-    H1_rotated, H2_rotated, bondlength1, bondlength2 = H2_rotated, H1_rotated, bondlength2, bondlength1
 
-print('Rotated system in z=0 plane about x=0 axis, with longer bondlength in H1')
+print('Reference system (z=0 plane about x=0 axis, with longer bondlength in H1)')
 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}')
+
+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)
+
+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}')
diff --git a/tests/PinnedH2O_3DOF/coords.in b/tests/PinnedH2O_3DOF/coords.in
new file mode 100644
index 00000000..46f5681d
--- /dev/null
+++ b/tests/PinnedH2O_3DOF/coords.in
@@ -0,0 +1,3 @@
+O1  1   0.00    0.00    0.00  0
+H1  2   1.12    1.45    0.00  1
+H2  2   1.12    -1.45    0.00  1
diff --git a/tests/PinnedH2O_3DOF/mgmol.cfg b/tests/PinnedH2O_3DOF/mgmol.cfg
new file mode 100644
index 00000000..4a0b39ef
--- /dev/null
+++ b/tests/PinnedH2O_3DOF/mgmol.cfg
@@ -0,0 +1,31 @@
+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=quench
+[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
diff --git a/tests/PinnedH2O_3DOF/test.py b/tests/PinnedH2O_3DOF/test.py
new file mode 100755
index 00000000..524027ba
--- /dev/null
+++ b/tests/PinnedH2O_3DOF/test.py
@@ -0,0 +1,91 @@
+#!/usr/bin/env python
+import sys
+import os
+import subprocess
+import string
+
+print("Test PinnedH2O_3DOF...")
+
+nargs=len(sys.argv)
+
+mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3]
+for i in range(4,nargs-4):
+  mpicmd = mpicmd + " "+sys.argv[i]
+print("MPI run command: {}".format(mpicmd))
+
+exe = sys.argv[nargs-4]
+inp = sys.argv[nargs-3]
+coords = sys.argv[nargs-2]
+print("coordinates file: %s"%coords)
+
+#create links to potentials files
+dst1 = 'pseudo.O_ONCV_PBE_SG15'
+dst2 = 'pseudo.H_ONCV_PBE_SG15'
+src1 = sys.argv[nargs-1] + '/' + dst1
+src2 = sys.argv[nargs-1] + '/' + dst2
+
+if not os.path.exists(dst1):
+  print("Create link to %s"%dst1)
+  os.symlink(src1, dst1)
+if not os.path.exists(dst2):
+  print("Create link to %s"%dst2)
+  os.symlink(src2, dst2)
+
+#run
+command = "{} {} -c {} -i {}".format(mpicmd,exe,inp,coords)
+print("Run command: {}".format(command))
+output = subprocess.check_output(command,shell=True)
+lines=output.split(b'\n')
+
+#analyse output
+energies=[]
+for line in lines:
+  if line.count(b'%%'):
+    print(line)
+    words=line.split()
+    words=words[5].split(b',')[0]
+    energy = words.decode()
+  if line.count(b'achieved'):
+    energies.append(energy)
+
+flag=0
+forces=[]
+for line in lines:
+  if flag>0:
+    print(line)
+    words=line.split(b'    ')
+    forces.append(words[1].decode())
+    forces.append(words[2].decode())
+    forces.append(words[3].decode())
+    flag=flag-1
+  if line.count(b'Forces:'):
+    flag=2
+
+
+print("Check energies...")
+print( energies )
+if len(energies)<2:
+  print("Expected two converged energies")
+  sys.exit(1)
+
+tol = 1.e-6
+diff=eval(energies[1])-eval(energies[0])
+print(diff)
+if abs(diff)>tol:
+  print("Energies differ: {} vs {} !!!".format(energies[0],energies[1]))
+  sys.exit(1)
+
+print("Check forces...")
+print(forces)
+flag=0
+for i in range(6):
+  diff=eval(forces[i+6])-eval(forces[i])
+  print(diff)
+  if abs(diff)>1.e-3:
+    print("Forces difference larger than tol")
+    flag=1
+if flag>0:
+  sys.exit(1)
+
+print("Test SUCCESSFUL!")
+sys.exit(0)