ROM on PinnedH2O with geometric transformation

examples/PinnedH2O/mgmol.cfg

@@ -1,6 +1,6 @@
 verbosity=1
 xcFunctional=PBE
-FDtype=Mehrstellen
+FDtype=4th
 [Mesh]
 nx=64
 ny=64

examples/PinnedH2O/mgmol_offline.cfg

@@ -1,6 +1,6 @@
 verbosity=1
 xcFunctional=PBE
-FDtype=Mehrstellen
+FDtype=4th
 [Mesh]
 nx=64
 ny=64

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;
 }

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,54 +38,62 @@ 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)
 
 target_plane_normal = np.array([0, 0, 1])
 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}')

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

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

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