From 54588138a3bbb0e7a958c5649e14d41f40fccf57 Mon Sep 17 00:00:00 2001 From: Siu Wun Cheung Date: Thu, 6 Feb 2025 11:41:15 -0800 Subject: [PATCH 1/5] Minor change --- tests/PinnedH2O_3DOF/job.online_test1 | 2 +- tests/PinnedH2O_3DOF/job.online_test2 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/PinnedH2O_3DOF/job.online_test1 b/tests/PinnedH2O_3DOF/job.online_test1 index bf69a1ef..eee111cf 100644 --- a/tests/PinnedH2O_3DOF/job.online_test1 +++ b/tests/PinnedH2O_3DOF/job.online_test1 @@ -25,6 +25,6 @@ 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 +srun -p pdebug -n $ncpus $exe -c $cfg_online -i data/1.00_1.00_0.0/coords.in date diff --git a/tests/PinnedH2O_3DOF/job.online_test2 b/tests/PinnedH2O_3DOF/job.online_test2 index 88665f1e..2fc0c17b 100644 --- a/tests/PinnedH2O_3DOF/job.online_test2 +++ b/tests/PinnedH2O_3DOF/job.online_test2 @@ -25,6 +25,6 @@ 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 +srun -p pdebug -n $ncpus $exe -c $cfg_online -i data/1.00_1.00_0.0/coords.in date From f6fd1218d694f98d243c4d0ce29ae4311b84d067 Mon Sep 17 00:00:00 2001 From: Siu Wun Cheung Date: Mon, 10 Feb 2025 13:56:45 -0800 Subject: [PATCH 2/5] Fix calculations --- .../aux/coords_1.00_1.00_0.0.in | 3 + .../aux/coords_1.02_0.98_2.0.in | 3 + tests/PinnedH2O_3DOF/aux/rotation_test.cc | 125 +++++++++++------- tests/PinnedH2O_3DOF/aux/rotation_test.py | 83 +++++++----- 4 files changed, 132 insertions(+), 82 deletions(-) create mode 100644 tests/PinnedH2O_3DOF/aux/coords_1.00_1.00_0.0.in create mode 100644 tests/PinnedH2O_3DOF/aux/coords_1.02_0.98_2.0.in diff --git a/tests/PinnedH2O_3DOF/aux/coords_1.00_1.00_0.0.in b/tests/PinnedH2O_3DOF/aux/coords_1.00_1.00_0.0.in new file mode 100644 index 00000000..46f5681d --- /dev/null +++ b/tests/PinnedH2O_3DOF/aux/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/tests/PinnedH2O_3DOF/aux/coords_1.02_0.98_2.0.in b/tests/PinnedH2O_3DOF/aux/coords_1.02_0.98_2.0.in new file mode 100644 index 00000000..4c440443 --- /dev/null +++ b/tests/PinnedH2O_3DOF/aux/coords_1.02_0.98_2.0.in @@ -0,0 +1,3 @@ +O1 1 0.00 0.00 0.00 0 +H1 2 1.12 1.50 0.00 1 +H2 2 1.07 -1.44 0.00 1 diff --git a/tests/PinnedH2O_3DOF/aux/rotation_test.cc b/tests/PinnedH2O_3DOF/aux/rotation_test.cc index 0e5641ee..c3a5c4b6 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,42 +45,59 @@ 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 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_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 axis_to_align[3]; cross(plane_normal, target_plane_normal, axis_to_align); @@ -86,56 +106,61 @@ int main() { 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; + apply_rotation(rot_matrix_align_plane, H1, H1_temp); + apply_rotation(rot_matrix_align_plane, 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; + H2_rotated[1] *= -1; } - 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 (z=0 plane 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; + H2_rotated[1] *= -1; } - 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(rot_matrix_align_plane, H1_temp, H1_restored); + apply_transpose_rotation(rot_matrix_align_plane, 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..b05c0818 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) @@ -35,46 +48,52 @@ def rotation_matrix(axis, angle): 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) - -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}') - H1_rotated = np.dot(rot_matrix_align_plane, H1) H2_rotated = np.dot(rot_matrix_align_plane, H2) -fliped_bond = False + +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) + 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 + +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 about x=0 axis, with longer bondlength in H1)') +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}') + +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) -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 = np.dot(rot_matrix_align_plane.T, H1_restored) +H2_restored = np.dot(rot_matrix_align_plane.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}') From ee976ee2bff05522b620a32ecd3a097697a29240 Mon Sep 17 00:00:00 2001 From: Siu Wun Cheung Date: Mon, 10 Feb 2025 14:56:13 -0800 Subject: [PATCH 3/5] Rename --- tests/PinnedH2O_3DOF/aux/rotation_test.cc | 12 ++++++------ tests/PinnedH2O_3DOF/aux/rotation_test.py | 10 +++++----- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/tests/PinnedH2O_3DOF/aux/rotation_test.cc b/tests/PinnedH2O_3DOF/aux/rotation_test.cc index c3a5c4b6..8793cb82 100644 --- a/tests/PinnedH2O_3DOF/aux/rotation_test.cc +++ b/tests/PinnedH2O_3DOF/aux/rotation_test.cc @@ -106,10 +106,10 @@ int main() 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); - apply_rotation(rot_matrix_align_plane, H1, H1_temp); - apply_rotation(rot_matrix_align_plane, H2, H2_temp); + 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; @@ -148,8 +148,8 @@ int main() apply_transpose_rotation(planar_rotation_matrix, H1_rotated, H1_temp); apply_transpose_rotation(planar_rotation_matrix, H2_rotated, H2_temp); - apply_transpose_rotation(rot_matrix_align_plane, H1_temp, H1_restored); - apply_transpose_rotation(rot_matrix_align_plane, H2_temp, H2_restored); + 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); diff --git a/tests/PinnedH2O_3DOF/aux/rotation_test.py b/tests/PinnedH2O_3DOF/aux/rotation_test.py index b05c0818..41e20296 100644 --- a/tests/PinnedH2O_3DOF/aux/rotation_test.py +++ b/tests/PinnedH2O_3DOF/aux/rotation_test.py @@ -47,9 +47,9 @@ 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) +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) theta1 = np.arctan2(H1_rotated[1], H1_rotated[0]) theta2 = np.arctan2(H2_rotated[1], H2_rotated[0]) @@ -84,8 +84,8 @@ def rotation_matrix(axis, angle): H1_restored = np.dot(planar_rotation_matrix.T, H1_restored) H2_restored = np.dot(planar_rotation_matrix.T, H2_restored) -H1_restored = np.dot(rot_matrix_align_plane.T, H1_restored) -H2_restored = np.dot(rot_matrix_align_plane.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) From 674b58acc554b5462ddac1c540f09dbc775ae1e9 Mon Sep 17 00:00:00 2001 From: Siu Wun Cheung Date: Tue, 11 Feb 2025 09:20:34 -0800 Subject: [PATCH 4/5] Add rotation --- examples/PinnedH2O/mgmol.cfg | 2 +- examples/PinnedH2O/mgmol_offline.cfg | 2 +- tests/PinnedH2O_3DOF/aux/rotation_test.cc | 11 +- tests/PinnedH2O_3DOF/coords.in_rotate1 | 3 + tests/PinnedH2O_3DOF/coords.in_rotate2 | 3 + tests/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in | 3 - tests/PinnedH2O_3DOF/coords_1.02_0.98_2.0.in | 3 - tests/PinnedH2O_3DOF/job.basis | 2 +- tests/PinnedH2O_3DOF/job.basis_test1 | 34 --- tests/PinnedH2O_3DOF/job.basis_test2 | 39 --- tests/PinnedH2O_3DOF/job.offline | 2 +- tests/PinnedH2O_3DOF/job.online | 2 +- tests/PinnedH2O_3DOF/job.online_test1 | 30 -- tests/PinnedH2O_3DOF/job.online_test2 | 30 -- tests/PinnedH2O_3DOF/mgmol_online_test1.cfg | 32 --- tests/PinnedH2O_3DOF/mgmol_online_test2.cfg | 36 --- tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc | 286 +++++++++++++++++++ 17 files changed, 302 insertions(+), 218 deletions(-) create mode 100644 tests/PinnedH2O_3DOF/coords.in_rotate1 create mode 100644 tests/PinnedH2O_3DOF/coords.in_rotate2 delete mode 100644 tests/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in delete mode 100644 tests/PinnedH2O_3DOF/coords_1.02_0.98_2.0.in delete mode 100644 tests/PinnedH2O_3DOF/job.basis_test1 delete mode 100644 tests/PinnedH2O_3DOF/job.basis_test2 delete mode 100644 tests/PinnedH2O_3DOF/job.online_test1 delete mode 100644 tests/PinnedH2O_3DOF/job.online_test2 delete mode 100644 tests/PinnedH2O_3DOF/mgmol_online_test1.cfg delete mode 100644 tests/PinnedH2O_3DOF/mgmol_online_test2.cfg 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/aux/rotation_test.cc b/tests/PinnedH2O_3DOF/aux/rotation_test.cc index 8793cb82..21cfd93a 100644 --- a/tests/PinnedH2O_3DOF/aux/rotation_test.cc +++ b/tests/PinnedH2O_3DOF/aux/rotation_test.cc @@ -91,7 +91,6 @@ 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]; @@ -121,15 +120,15 @@ int main() if (bondlength1 < bondlength2) { swap(H1_rotated, H2_rotated); - H1_rotated[1] *= -1; - H2_rotated[1] *= -1; + H1_rotated[1] *= -1.0; + H2_rotated[1] *= -1.0; } 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 (z=0 plane about x=0 axis, with longer bondlength in Q1)" << endl; + 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_rotated << endl; @@ -141,8 +140,8 @@ int main() if (bondlength1 < bondlength2) { swap(H1_rotated, H2_rotated); - H1_rotated[1] *= -1; - H2_rotated[1] *= -1; + H1_rotated[1] *= -1.0; + H2_rotated[1] *= -1.0; } apply_transpose_rotation(planar_rotation_matrix, H1_rotated, H1_temp); 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/coords_1.00_1.00_0.0.in b/tests/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in deleted file mode 100644 index 46f5681d..00000000 --- a/tests/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in +++ /dev/null @@ -1,3 +0,0 @@ -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/coords_1.02_0.98_2.0.in b/tests/PinnedH2O_3DOF/coords_1.02_0.98_2.0.in deleted file mode 100644 index 4c440443..00000000 --- a/tests/PinnedH2O_3DOF/coords_1.02_0.98_2.0.in +++ /dev/null @@ -1,3 +0,0 @@ -O1 1 0.00 0.00 0.00 0 -H1 2 1.12 1.50 0.00 1 -H2 2 1.07 -1.44 0.00 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_test1 deleted file mode 100644 index eee111cf..00000000 --- a/tests/PinnedH2O_3DOF/job.online_test1 +++ /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_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 data/1.00_1.00_0.0/coords.in - -date diff --git a/tests/PinnedH2O_3DOF/job.online_test2 b/tests/PinnedH2O_3DOF/job.online_test2 deleted file mode 100644 index 2fc0c17b..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 data/1.00_1.00_0.0/coords.in - -date diff --git a/tests/PinnedH2O_3DOF/mgmol_online_test1.cfg b/tests/PinnedH2O_3DOF/mgmol_online_test1.cfg deleted file mode 100644 index e032350d..00000000 --- a/tests/PinnedH2O_3DOF/mgmol_online_test1.cfg +++ /dev/null @@ -1,32 +0,0 @@ -verbosity=2 -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=quench -[Quench] -max_steps=100 -atol=1.e-8 -[Orbitals] -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/mgmol_online_test2.cfg b/tests/PinnedH2O_3DOF/mgmol_online_test2.cfg deleted file mode 100644 index 8ce3e32d..00000000 --- a/tests/PinnedH2O_3DOF/mgmol_online_test2.cfg +++ /dev/null @@ -1,36 +0,0 @@ -verbosity=2 -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=quench -[Quench] -max_steps=1000 -atol=1.e-8 -[Orbitals] -initial_type=Random -initial_width=2. -nempty=4 -[Restart] -output_level=4 -[DensityMatrix] -solver=MVP -nb_inner_it=20 -[ROM.offline] -basis_file=data/PinnedH2O_3DOF_orbitals_basis_test2 -[ROM.basis] -compare_md=false -number_of_orbital_basis=8 diff --git a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc index e7813777..25b14916 100644 --- a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc +++ b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc @@ -22,10 +22,244 @@ #include #include #include +#include +#include +#include #include 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& positions, std::vector& anumbers, int atom_order[3], + 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_idx, H2_idx); + 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; + + atom_order[0] = H2_idx; + atom_order[1] = O1_idx; + atom_order[2] = H1_idx; +} + +void transpose_rotate_PinnedH2O(std::vector& positions, std::vector& forces, std::vector& anumbers, + const int atom_order[3], const double out_of_plane_rotation_matrix[3][3], + const double planar_rotation_angle, const bool flipped_bond) +{ + 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); + + int H2_idx = 0; + int O1_idx = 1; + int H1_idx = 2; + + //int H2_idx = atom_order[0]; + //int O1_idx = atom_order[1]; + //int H1_idx = atom_order[2]; + + 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; + std::swap(H1_rotated, H2_rotated); + f_O1_rotated[1] *= -1.0; + f_H1_rotated[1] *= -1.0; + f_H2_rotated[1] *= -1.0; + std::swap(f_H1_rotated, f_H2_rotated); + } + + 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]; + + 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]; + + anumbers[H2_idx] = 1; + anumbers[O1_idx] = 8; + anumbers[H1_idx] = 1; +} + int main(int argc, char** argv) { int mpirc = MPI_Init(&argc, &argv); @@ -127,6 +361,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, atom_order, 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 +413,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::iterator ita = anumbers.begin(); + for (std::vector::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::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, anumbers, atom_order, out_of_plane_rotation_matrix, planar_rotation_angle, flipped_bond); + + // print out results + if (MPIdata::onpe0) + { + std::cout << "Positions:" << std::endl; + std::vector::iterator ita = anumbers.begin(); + for (std::vector::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::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++; } } From 121e7fb42312b857c1f3ea53893d20b05e5d4e48 Mon Sep 17 00:00:00 2001 From: Siu Wun Cheung Date: Tue, 11 Feb 2025 10:02:07 -0800 Subject: [PATCH 5/5] Clean up --- tests/PinnedH2O_3DOF/job.online_rotate | 36 ++++++++++ tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg | 36 ++++++++++ tests/PinnedH2O_3DOF/mgmol_ref_rotate.cfg | 27 ++++++++ tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc | 72 +++++++------------- 4 files changed, 124 insertions(+), 47 deletions(-) create mode 100644 tests/PinnedH2O_3DOF/job.online_rotate create mode 100644 tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg create mode 100644 tests/PinnedH2O_3DOF/mgmol_ref_rotate.cfg diff --git a/tests/PinnedH2O_3DOF/job.online_rotate b/tests/PinnedH2O_3DOF/job.online_rotate new file mode 100644 index 00000000..11dd7419 --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.online_rotate @@ -0,0 +1,36 @@ +#!/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-20250206 + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +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 + +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/mgmol_online_rotate.cfg b/tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg new file mode 100644 index 00000000..72288ac3 --- /dev/null +++ b/tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg @@ -0,0 +1,36 @@ +verbosity=2 +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=quench +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +nempty=68 +[Restart] +output_level=4 +[DensityMatrix] +solver=MVP +nb_inner_it=100 +[ROM.offline] +basis_file=data/PinnedH2O_3DOF_orbitals_basis_2_2 +[ROM.basis] +compare_md=false +number_of_orbital_basis=72 diff --git a/tests/PinnedH2O_3DOF/mgmol_ref_rotate.cfg b/tests/PinnedH2O_3DOF/mgmol_ref_rotate.cfg new file mode 100644 index 00000000..96be6b8d --- /dev/null +++ b/tests/PinnedH2O_3DOF/mgmol_ref_rotate.cfg @@ -0,0 +1,27 @@ +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=quench +[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/testPinnedH2O_3DOF.cc b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc index 25b14916..a0cd51d3 100644 --- a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc +++ b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc @@ -95,7 +95,7 @@ void apply_transpose_rotation(const double matrix[3][3], const double vec[3], do result[2] = matrix[0][2] * vec[0] + matrix[1][2] * vec[1] + matrix[2][2] * vec[2]; } -void rotate_PinnedH2O(std::vector& positions, std::vector& anumbers, int atom_order[3], +void rotate_PinnedH2O(std::vector& positions, std::vector& anumbers, double out_of_plane_rotation_matrix[3][3], double& planar_rotation_angle, bool& flipped_bond) { int O1_idx = -1; @@ -159,7 +159,6 @@ void rotate_PinnedH2O(std::vector& positions, std::vector& anumbe { H1_rotated[1] *= -1.0; H2_rotated[1] *= -1.0; - std::swap(H1_idx, H2_idx); std::swap(H1_rotated, H2_rotated); } @@ -176,28 +175,12 @@ void rotate_PinnedH2O(std::vector& positions, std::vector& anumbe anumbers[0] = 1; anumbers[1] = 8; anumbers[2] = 1; - - atom_order[0] = H2_idx; - atom_order[1] = O1_idx; - atom_order[2] = H1_idx; } -void transpose_rotate_PinnedH2O(std::vector& positions, std::vector& forces, std::vector& anumbers, - const int atom_order[3], const double out_of_plane_rotation_matrix[3][3], +void transpose_rotate_PinnedH2O(std::vector& positions, std::vector& forces, + const double out_of_plane_rotation_matrix[3][3], const double planar_rotation_angle, const bool flipped_bond) { - 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); - - int H2_idx = 0; - int O1_idx = 1; - int H1_idx = 2; - - //int H2_idx = atom_order[0]; - //int O1_idx = atom_order[1]; - //int H1_idx = atom_order[2]; - 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]}; @@ -210,13 +193,15 @@ void transpose_rotate_PinnedH2O(std::vector& positions, std::vector& positions, std::vectorgrid(); @@ -439,7 +417,7 @@ int main(int argc, char** argv) } // rotate the forces to the original coordinate system - transpose_rotate_PinnedH2O(positions, forces, anumbers, atom_order, out_of_plane_rotation_matrix, planar_rotation_angle, flipped_bond); + transpose_rotate_PinnedH2O(positions, forces, out_of_plane_rotation_matrix, planar_rotation_angle, flipped_bond); // print out results if (MPIdata::onpe0)