From a733973bc917d0fd02a6d0d56dee9179e039374f Mon Sep 17 00:00:00 2001 From: Siu Wun Cheung Date: Tue, 5 Nov 2024 10:14:51 -0800 Subject: [PATCH] Add new files --- examples/PinnedH2O/generate_coord_simple.py | 32 +++++++++++ examples/PinnedH2O/rotation_test.py | 62 +++++++++++++++++++++ 2 files changed, 94 insertions(+) create mode 100644 examples/PinnedH2O/generate_coord_simple.py create mode 100644 examples/PinnedH2O/rotation_test.py diff --git a/examples/PinnedH2O/generate_coord_simple.py b/examples/PinnedH2O/generate_coord_simple.py new file mode 100644 index 00000000..181691a5 --- /dev/null +++ b/examples/PinnedH2O/generate_coord_simple.py @@ -0,0 +1,32 @@ +import numpy as np +import os + +O1 = np.array([0.00, 0.00, 0.00]) + +ref_bondlength = 1.83 +ref_bondangle = 104.5 + +# factors and increments for bond lengths and bond angle +bondlength1_factor = np.linspace(0.95, 1.05, 11) +bondlength2_factor = np.linspace(0.95, 1.05, 11) +bondangle_increment = np.linspace(-5, 5, 11) + +# output directory +output_dir = "PinnedH2O_3dof_coords" + +# generation +os.makedirs(output_dir, exist_ok=True) + +for d_bondangle in bondangle_increment: + bondangle = ref_bondangle + d_bondangle + x = ref_bondlength * np.cos(np.radians(bondangle / 2)) + y = ref_bondlength * np.cos(np.radians(bondangle / 2)) + for f_bondlength1 in bondlength1_factor: + for f_bondlength2 in bondlength2_factor: + H1 = np.array([f_bondlength1*x, f_bondlength1*y, 0.0]) + H2 = np.array([f_bondlength2*x, -f_bondlength2*y, 0.0]) + filename = f"{output_dir}/coords_{f_bondlength1:.2f}_{f_bondlength2:.2f}_{d_bondangle}.in" + with open(filename, "w") as file: + file.write(f"O1 1 {O1[0]:.2f} {O1[1]:.2f} {O1[2]:.2f} 0\n") + file.write(f"H1 2 {H1[0]:.2f} {H1[1]:.2f} {H1[2]:.2f} 1\n") + file.write(f"H2 2 {H2[0]:.2f} {H2[1]:.2f} {H2[2]:.2f} 1\n") diff --git a/examples/PinnedH2O/rotation_test.py b/examples/PinnedH2O/rotation_test.py new file mode 100644 index 00000000..7f74c55b --- /dev/null +++ b/examples/PinnedH2O/rotation_test.py @@ -0,0 +1,62 @@ +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]) + +def calculate_bondlength(atom1, atom2): + return np.linalg.norm(atom1 - atom2) + +def calculate_bondangle(atom1, atom2, atom3, radian): + vector1 = atom1 - atom2 + vector2 = atom3 - atom2 + dot_product = np.dot(vector1, vector2) + magnitude_product = np.linalg.norm(vector1) * np.linalg.norm(vector2) + angle = np.arccos(dot_product / magnitude_product) + if not radian: + angle = np.degrees(angle) + return angle + +def rotation_matrix(axis, angle): + cos_theta = np.cos(angle) + sin_theta = np.sin(angle) + ux, uy, uz = axis + return np.array([ + [cos_theta + ux**2 * (1 - cos_theta), ux * uy * (1 - cos_theta) - uz * sin_theta, ux * uz * (1 - cos_theta) + uy * sin_theta], + [uy * ux * (1 - cos_theta) + uz * sin_theta, cos_theta + uy**2 * (1 - cos_theta), uy * uz * (1 - cos_theta) - ux * sin_theta], + [uz * ux * (1 - cos_theta) - uy * sin_theta, uz * uy * (1 - cos_theta) + ux * sin_theta, cos_theta + uz**2 * (1 - cos_theta)] + ]) + +plane_normal = np.cross(H1, H2) +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) +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) +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}') + +bondlength1 = calculate_bondlength(H1_rotated, O1) +bondlength2 = calculate_bondlength(H2_rotated, O1) +bondangle = calculate_bondangle(H1_rotated, O1, H2_rotated, False) + +print('Aligned system') +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}')