Add restore check. Make H1 shorter than H2 in original system for ver…

…ification

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,13 +47,32 @@ 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)
 bondlength1 = calculate_bondlength(H1_rotated, O1)
 bondlength2 = calculate_bondlength(H2_rotated, O1)
 bondangle = calculate_bondangle(H1_rotated, O1, H2_rotated, False)
+fliped_bond = False
 if bondlength1 < bondlength2:
+    fliped_bond = True
     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_rotated = np.dot(rot_matrix_align_plane.T, H1_rotated)
+H2_rotated = np.dot(rot_matrix_align_plane.T, H2_rotated)
+bondlength1 = calculate_bondlength(H1_rotated, O1)
+bondlength2 = calculate_bondlength(H2_rotated, O1)
+bondangle = calculate_bondangle(H1_rotated, O1, H2_rotated, False)
+
+print('Restored system')
 print(f'H1 = {H1_rotated}')
 print(f'H2 = {H2_rotated}')
 print(f'Bondlength of O1-H1 = {bondlength1}')