diff --git a/examples/PinnedH2O/generate_coord.py b/examples/PinnedH2O/generate_coord.py new file mode 100644 index 00000000..c239eb87 --- /dev/null +++ b/examples/PinnedH2O/generate_coord.py @@ -0,0 +1,63 @@ +import numpy as np +import os + +# coords.in +O1 = np.array([0.00, 0.00, 0.00]) +ref_H1 = np.array([-0.45, 1.42, -1.07]) +ref_H2 = np.array([-0.45, -1.48, -0.97]) + +# 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" + +# utilities +def calculate_bondlength(atom1, atom2): + return np.linalg.norm(atom1 - atom2) + +def calculate_bondangle(atom1, atom2, atom3): + vector1 = atom1 - atom2 + vector2 = atom3 - atom2 + dot_product = np.dot(vector1, vector2) + magnitude_product = np.linalg.norm(vector1) * np.linalg.norm(vector2) + angle_rad = np.arccos(dot_product / magnitude_product) + angle_deg = np.degrees(angle_rad) + return angle_deg + +# Rodrigues' rotation formula +def rotation_matrix(axis, angle_degrees): + angle = np.radians(angle_degrees) + 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)] + ]) + +# generation +os.makedirs(output_dir, exist_ok=True) + +ref_bondlength1 = calculate_bondlength(ref_H1, O1) +ref_bondlength2 = calculate_bondlength(ref_H2, O1) +ref_bondangle = calculate_bondangle(ref_H1, O1, ref_H2) + +normal_vector = np.cross(ref_H1, ref_H2) +normal_unit_vector = normal_vector / np.linalg.norm(normal_vector) + +for d_bondangle in bondangle_increment: + Q = rotation_matrix(normal_unit_vector, d_bondangle) + Q_ref_H2 = np.dot(Q, ref_H2) + for f_bondlength1 in bondlength1_factor: + for f_bondlength2 in bondlength2_factor: + H1 = f_bondlength1 * ref_H1 + H2 = f_bondlength2 * Q_ref_H2 + 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/generate_coord_simple.py b/examples/PinnedH2O/generate_coord_simple.py new file mode 100644 index 00000000..7e799bdd --- /dev/null +++ b/examples/PinnedH2O/generate_coord_simple.py @@ -0,0 +1,31 @@ +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 +bondlength_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.sin(np.radians(bondangle / 2)) + for i, f_bondlength1 in enumerate(bondlength_factor): + for f_bondlength2 in bondlength_factor[:(i+1)]: + 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..79f1aaef --- /dev/null +++ b/examples/PinnedH2O/rotation_test.py @@ -0,0 +1,64 @@ +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(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) +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) +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(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}')