diff --git a/examples/PinnedH2O/generate_coord.py b/examples/PinnedH2O/generate_coord.py deleted file mode 100644 index c239eb87..00000000 --- a/examples/PinnedH2O/generate_coord.py +++ /dev/null @@ -1,63 +0,0 @@ -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/get_ROM_table.py b/examples/PinnedH2O/get_ROM_table.py new file mode 100644 index 00000000..fa4289f3 --- /dev/null +++ b/examples/PinnedH2O/get_ROM_table.py @@ -0,0 +1,28 @@ +import subprocess +import re + +pattern = r"For energy fraction: \d+\.\d+, take first (\d+) of \d+ basis vectors" + +print("\\begin{tabular}{|c||c|c|c|c|c|c|c|}") +print("\\hline") +print("$k$ & $\\varepsilon = 10^{-1}$ & $\\varepsilon = 10^{-2}$ & $\\varepsilon = 10^{-3}$ & $\\varepsilon = 10^{-4}$ & $\\varepsilon = 10^{-5}$ & Snapshots \\\\") +print("\\hline") + +for t in range(10): + k = 50*(t+1) + snapshots = 4*k + grep_command = f"grep 'take first' basis_1_{k}_Pinned_H2O.out" + result = subprocess.run(grep_command, shell=True, capture_output=True, text=True) + matches = re.findall(pattern, result.stdout) + energy_fractions = { + "0.9": matches[0], + "0.99": matches[1], + "0.999": matches[2], + "0.9999": matches[3], + "0.99999": matches[4], + } + line = f"{k} & {energy_fractions['0.9']} & {energy_fractions['0.99']} & {energy_fractions['0.999']} & {energy_fractions['0.9999']} & {energy_fractions['0.99999']} & {snapshots} \\\\" + print(line) + +print("\\hline") +print("\\end{tabular}") diff --git a/src/DensityMatrix.cc b/src/DensityMatrix.cc index 2f5e7c1c..78d7008e 100644 --- a/src/DensityMatrix.cc +++ b/src/DensityMatrix.cc @@ -43,10 +43,6 @@ DensityMatrix::DensityMatrix(const int ndim) { assert(ndim > 0); - occ_uptodate_ = false; - stripped_ = false; - uniform_occ_ = false; - MGmol_MPI& mmpi = *(MGmol_MPI::instance()); orbital_occupation_ = mmpi.nspin() > 1 ? 1. : 2.; diff --git a/src/ExtendedGridOrbitals.h b/src/ExtendedGridOrbitals.h index c1a718d0..74af1982 100644 --- a/src/ExtendedGridOrbitals.h +++ b/src/ExtendedGridOrbitals.h @@ -400,6 +400,9 @@ class ExtendedGridOrbitals : public Orbitals const pb::Grid& mygrid = mymesh->grid(); return mygrid.maxDomainSize(); } +#ifdef MGMOL_HAS_LIBROM + void set(std::string file_path, int rdim); +#endif }; #endif diff --git a/src/MGmol.cc b/src/MGmol.cc index acc2b944..0a5c4ad4 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -1511,6 +1511,8 @@ double MGmol::evaluateDMandEnergyAndForces(Orbitals* orbitals, force(*dorbitals, *ions_); + ions_->getForces(forces); + return eks; } diff --git a/src/MGmol.h b/src/MGmol.h index 1eb8b2d3..315a0ff9 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -359,8 +359,8 @@ class MGmol : public MGmolInterface } #ifdef MGMOL_HAS_LIBROM - int save_orbital_snapshot(std::string snapshot_dir, OrbitalsType& orbitals); - void project_orbital(std::string snapshot_dir, int rdim, OrbitalsType& orbitals); + int save_orbital_snapshot(std::string file_path, OrbitalsType& orbitals); + void project_orbital(std::string file_path, int rdim, OrbitalsType& orbitals); #endif }; // Instantiate static variables here to avoid clang warnings diff --git a/src/md.cc b/src/md.cc index 36236f85..01e34a87 100644 --- a/src/md.cc +++ b/src/md.cc @@ -492,7 +492,6 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) force(**orbitals, ions); #ifdef MGMOL_HAS_LIBROM - // TODO: cleanup if (ct.getROMOptions().num_orbbasis > 0) { if (onpe0) diff --git a/src/rom.cc b/src/rom.cc index 19ff66a3..986e136f 100644 --- a/src/rom.cc +++ b/src/rom.cc @@ -11,6 +11,7 @@ #ifdef MGMOL_HAS_LIBROM #include "LocGridOrbitals.h" +#include "ExtendedGridOrbitals.h" #include "MGmol.h" #include "librom.h" @@ -91,6 +92,25 @@ void MGmol::project_orbital(std::string file_path, int rdim, Orbit } } +void ExtendedGridOrbitals::set(std::string file_path, int rdim) +{ + const int dim = getLocNumpt(); + + CAROM::BasisReader reader(file_path); + CAROM::Matrix* orbital_basis = reader.getSpatialBasis(rdim); + + Control& ct = *(Control::instance()); + Mesh* mymesh = Mesh::instance(); + pb::GridFunc gf_psi(mymesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]); + CAROM::Vector psi; + for (int i = 0; i < rdim; ++i) + { + orbital_basis->getColumn(i, psi); + gf_psi.assign(psi.getData()); + setPsi(gf_psi, i); + } +} + template class MGmol; template class MGmol; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b29c1aad..1855eb1e 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -258,6 +258,8 @@ add_executable(testDMandEnergyAndForces ${CMAKE_SOURCE_DIR}/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc) add_executable(testRestartEnergyAndForces ${CMAKE_SOURCE_DIR}/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc) +add_executable(testPinnedH2O_3DOF + ${CMAKE_SOURCE_DIR}/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc) if(${MAGMA_FOUND}) add_executable(testOpenmpOffload @@ -382,6 +384,14 @@ add_test(NAME testRestartEnergyAndForces ${CMAKE_CURRENT_SOURCE_DIR}/RestartEnergyAndForces/restart.cfg ${CMAKE_CURRENT_SOURCE_DIR}/RestartEnergyAndForces/h2o.xyz ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +add_test(NAME testPinnedH2O_3DOF + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/PinnedH2O_3DOF/test.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/testPinnedH2O_3DOF + ${CMAKE_CURRENT_SOURCE_DIR}/PinnedH2O_3DOF/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/PinnedH2O_3DOF/coords.in + ${CMAKE_CURRENT_SOURCE_DIR}/PinnedH2O_3DOF/lrs.in + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) if(${MAGMA_FOUND}) add_test(NAME testOpenmpOffload @@ -579,6 +589,7 @@ target_link_libraries(testEnergyAndForces PRIVATE mgmol_src) target_link_libraries(testWFEnergyAndForces PRIVATE mgmol_src) target_link_libraries(testDMandEnergyAndForces PRIVATE mgmol_src) target_link_libraries(testRestartEnergyAndForces PRIVATE mgmol_src) +target_link_libraries(testPinnedH2O_3DOF PRIVATE mgmol_src) target_link_libraries(testIons PRIVATE mgmol_src) target_link_libraries(testDensityMatrix PRIVATE ${HDF5_LIBRARIES}) diff --git a/tests/PinnedH2O_3DOF/README.md b/tests/PinnedH2O_3DOF/README.md new file mode 100644 index 00000000..90551b47 --- /dev/null +++ b/tests/PinnedH2O_3DOF/README.md @@ -0,0 +1 @@ +python3 test.py srun -p debug -n 64 ../../build_quartz/tests/testPinnedH2O_3DOF mgmol_online.cfg coords_1.00_1.00_0.0.in ../../potentials diff --git a/tests/PinnedH2O_3DOF/assemble_results.sh b/tests/PinnedH2O_3DOF/assemble_results.sh new file mode 100644 index 00000000..6a9a4d0d --- /dev/null +++ b/tests/PinnedH2O_3DOF/assemble_results.sh @@ -0,0 +1,31 @@ +#!/bin/bash + +rm -rf PinnedH2O_3DOF_assembled_results +mkdir PinnedH2O_3DOF_assembled_results + +for d in results_*; do + echo "Processing $d" + rm $d/energy_and_forces_*.out + rm $d/*.png + + suffix=$(echo "$d" | cut -d'_' -f2,3,4) + N_l=$(echo "$suffix" | cut -d'_' -f1) + N_theta=$(echo "$suffix" | cut -d'_' -f2) + rdim=$(echo "$suffix" | cut -d'_' -f3) + python3 compare_energy_and_forces.py --N_l "$N_l" --N_theta "$N_theta" --rdim "$rdim" + + cp "$d/Eks_difference_histogram_reproductive.png" "PinnedH2O_3DOF_assembled_results/Eks_difference_histogram_${suffix}_reproductive.png" + cp "$d/f_H1_difference_histogram_reproductive.png" "PinnedH2O_3DOF_assembled_results/f_H1_difference_histogram_${suffix}_reproductive.png" + cp "$d/f_H2_difference_histogram_reproductive.png" "PinnedH2O_3DOF_assembled_results/f_H2_difference_histogram_${suffix}_reproductive.png" + cp "$d/f_O1_difference_histogram_reproductive.png" "PinnedH2O_3DOF_assembled_results/f_O1_difference_histogram_${suffix}_reproductive.png" + + cp "$d/Eks_difference_histogram_predictive.png" "PinnedH2O_3DOF_assembled_results/Eks_difference_histogram_${suffix}_predictive.png" + cp "$d/f_H1_difference_histogram_predictive.png" "PinnedH2O_3DOF_assembled_results/f_H1_difference_histogram_${suffix}_predictive.png" + cp "$d/f_H2_difference_histogram_predictive.png" "PinnedH2O_3DOF_assembled_results/f_H2_difference_histogram_${suffix}_predictive.png" + cp "$d/f_O1_difference_histogram_predictive.png" "PinnedH2O_3DOF_assembled_results/f_O1_difference_histogram_${suffix}_predictive.png" +done + +tar -cvf PinnedH2O_3DOF_assembled_results.tar PinnedH2O_3DOF_assembled_results +mv PinnedH2O_3DOF_assembled_results.tar /p/lustre2/cheung26/scp_local + +echo "Finished assembling results." diff --git a/tests/PinnedH2O_3DOF/aux/PinnedH2O_3DOF.png b/tests/PinnedH2O_3DOF/aux/PinnedH2O_3DOF.png new file mode 100644 index 00000000..b484f4c9 Binary files /dev/null and b/tests/PinnedH2O_3DOF/aux/PinnedH2O_3DOF.png differ diff --git a/tests/PinnedH2O_3DOF/aux/get_ROM_table.py b/tests/PinnedH2O_3DOF/aux/get_ROM_table.py new file mode 100644 index 00000000..5fea1e38 --- /dev/null +++ b/tests/PinnedH2O_3DOF/aux/get_ROM_table.py @@ -0,0 +1,31 @@ +import subprocess +import re + +bondlength_num_increments = (2, 5, 10) +bondangle_num_increments = (2, 5, 10) + +pattern = r"For energy fraction: \d+\.\d+, take first (\d+) of \d+ basis vectors" + +print("\\begin{tabular}{|c|c||c|c|c|c|c|c|c|}") +print("\\hline") +print("$N_L$ & $N_\\theta$ & $\\varepsilon = 10^{-1}$ & $\\varepsilon = 10^{-2}$ & $\\varepsilon = 10^{-3}$ & $\\varepsilon = 10^{-4}$ & $\\varepsilon = 10^{-5}$ & Snapshots \\\\") +print("\\hline") + +for _, N_L in enumerate(bondlength_num_increments): + for _, N_theta in enumerate(bondangle_num_increments): + snapshots = 2*(N_L+1)*(N_L+2)*(N_theta+1) + grep_command = f"grep 'take first' basis_{N_L}_{N_theta}_PinnedH2O_3DOF.out" + result = subprocess.run(grep_command, shell=True, capture_output=True, text=True) + matches = re.findall(pattern, result.stdout) + energy_fractions = { + "0.9": matches[0], + "0.99": matches[1], + "0.999": matches[2], + "0.9999": matches[3], + "0.99999": matches[4], + } + line = f"{N_L} & {N_theta} & {energy_fractions['0.9']} & {energy_fractions['0.99']} & {energy_fractions['0.999']} & {energy_fractions['0.9999']} & {energy_fractions['0.99999']} & {snapshots} \\\\" + print(line) + +print("\\hline") +print("\\end{tabular}") diff --git a/tests/PinnedH2O_3DOF/aux/plot_PinnedH2O_3DOF.py b/tests/PinnedH2O_3DOF/aux/plot_PinnedH2O_3DOF.py new file mode 100644 index 00000000..e679c691 --- /dev/null +++ b/tests/PinnedH2O_3DOF/aux/plot_PinnedH2O_3DOF.py @@ -0,0 +1,69 @@ +import matplotlib +import matplotlib.pyplot as plt +import numpy as np + +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, 3) +bondangle_increment = np.linspace(-5, 5, 3) + +fig, ax = plt.subplots() + +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): + H1_x = f_bondlength1*x + H1_y = f_bondlength1*y + ax.plot(H1_x, H1_y, 'ro', markersize=8, alpha=0.1) + for f_bondlength2 in bondlength_factor[:(i+1)]: + H2_x = f_bondlength2*x + H2_y = -f_bondlength2*y + ax.plot(H2_x, H2_y, 'bo', markersize=8, alpha=0.1) + +ref_bondangle = np.radians(ref_bondangle) +x = ref_bondlength * np.cos(ref_bondangle / 2) +y = ref_bondlength * np.sin(ref_bondangle / 2) + +ax.plot([0, x], [0, y], 'r-', lw=2) +ax.plot([0, x], [0, -y], 'b-', lw=2) +ax.plot(0, 0, 'ko', markersize=8) +ax.plot(x, y, 'ro', markersize=8) +ax.plot(x, -y, 'bo', markersize=8) + +arc1 = matplotlib.patches.Arc(xy=(0, 0), + width=0.5, + height=0.5, + angle=0, + theta1=0, + theta2=ref_bondangle/2, + color='red', + linewidth=2) +ax.add_patch(arc1) +ax.text(0.4, 0.25, r"$\frac{\theta}{2}$", ha='center', va='center', fontsize=12, color='red') + +arc2 = matplotlib.patches.Arc(xy=(0, 0), + width=0.5, + height=0.5, + angle=0, + theta1=-ref_bondangle/2, + theta2=0, + color='blue', + linewidth=2) +ax.add_patch(arc2) +ax.text(0.4, -0.25, r"$\frac{\theta}{2}$", ha='center', va='center', fontsize=12, color='blue') + +ax.text(x / 2 - 0.1, y / 2, r"$L_1$", ha='right', color='red') +ax.text(x / 2 - 0.1, -y / 2, r"$L_2$", ha='right', color='blue') + +ax.plot([-2, 2], [0, 0], 'k--', lw=1) +ax.set_xlim(-2.0, 2.0) +ax.set_ylim(-2.0, 2.0) +ax.set_aspect('equal') +ax.set_xlabel("x (Bohr)") +ax.set_ylabel("y (Bohr)") +ax.grid(True) +plt.savefig("PinnedH2O_3DOF.png") diff --git a/tests/PinnedH2O_3DOF/aux/rotation_test.cc b/tests/PinnedH2O_3DOF/aux/rotation_test.cc new file mode 100644 index 00000000..0e5641ee --- /dev/null +++ b/tests/PinnedH2O_3DOF/aux/rotation_test.cc @@ -0,0 +1,141 @@ +#include +#include +#include +#include + +using namespace std; + +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 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); + + if (!radian) { + angle = angle * 180.0 / M_PI; + } + 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]); + if (norm > 0) { + 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]; +} + +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); + + 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; + swap(H1_rotated, H2_rotated); + } + 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; + 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; + + if (flipped_bond) { + swap(H1_rotated, H2_rotated); + } + 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); + + 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; + + return 0; +} diff --git a/examples/PinnedH2O/rotation_test.py b/tests/PinnedH2O_3DOF/aux/rotation_test.py similarity index 73% rename from examples/PinnedH2O/rotation_test.py rename to tests/PinnedH2O_3DOF/aux/rotation_test.py index 79f1aaef..fbecb532 100644 --- a/examples/PinnedH2O/rotation_test.py +++ b/tests/PinnedH2O_3DOF/aux/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,15 +47,34 @@ 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) +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) -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('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_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) + +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}') diff --git a/tests/PinnedH2O_3DOF/compare_energy_and_forces.py b/tests/PinnedH2O_3DOF/compare_energy_and_forces.py new file mode 100644 index 00000000..b11dbe91 --- /dev/null +++ b/tests/PinnedH2O_3DOF/compare_energy_and_forces.py @@ -0,0 +1,233 @@ +import numpy as np +import os +import matplotlib.pyplot as plt +import matplotlib.ticker as ticker +import argparse + +parser = argparse.ArgumentParser(description="Calculate and plot differences in energies and forces.") + +parser.add_argument("--N_l", type=int, default=2, help="Sampling frequency of bond length points") +parser.add_argument("--N_theta", type=int, default=2, help="Sampling frequency of bond angle points") +parser.add_argument("--rdim", type=int, default=18, help="Dimension of the ROM basis") + +args = parser.parse_args() + +N_l = args.N_l +N_theta = args.N_theta +rdim = args.rdim + +bondlength_min = 0.95 +bondlength_max = 1.05 +bondlength_num_increments = 10 + +bondangle_min = -5.0 +bondangle_max = 5.0 +bondangle_num_increments = 10 + +output_dir = f'results_{N_l}_{N_theta}_{rdim}' + +def process_data(s_l1, s_l2, s_theta, N_l, N_theta, rdim): + try: + data_dir = f'data/{s_l1}_{s_l2}_{s_theta}' + offline_file = os.path.join(data_dir, 'offline_PinnedH2O.out') + log_file = f'{output_dir}/{s_l1}_{s_l2}_{s_theta}.log' + + f_O1_fom = None + f_H1_fom = None + f_H2_fom = None + + if not os.path.exists(offline_file): + print(f"Error: File not found: {offline_file}") + return None + + with open(offline_file, 'r') as f: + eks_lines = [line for line in f if " SC ENERGY [Ha]" in line] + if eks_lines: + Eks_fom = float(eks_lines[-1].split()[-1]) + else: + print(f"Error: Could not find 'SC ENERGY [Ha]' line in {offline_file}") + return None + + with open(offline_file, 'r') as f: + for line in f: + if "O1" in line: + f_O1_fom = np.array(line.split()[6:9], dtype=float) + elif "H1" in line: + f_H1_fom = np.array(line.split()[5:8], dtype=float) + elif "H2" in line: + f_H2_fom = np.array(line.split()[5:8], dtype=float) + + if f_O1_fom is not None and f_H1_fom is not None and f_H2_fom is not None: + break # Exit the loop once we found all lines + + if f_O1_fom is None or f_H1_fom is None or f_H2_fom is None: + print(f"Error: Could not find O1, H1, or H2 lines in {offline_file}") + return None + + f_H2_rom = None + f_O1_rom = None + f_H1_rom = None + + if not os.path.exists(log_file): + print(f"Error: File not found: {log_file}") + return None + + with open(log_file, 'r') as f: + lines = f.readlines() + + eks_index = -1 + for i in range(len(lines) - 1, -1, -1): # Search from the end + if "Eks:" in lines[i]: + eks_index = i + break + + if eks_index == -1: + print(f"Error: Could not find 'Eks:' line in {log_file}") + return None + + Eks_rom = float(lines[eks_index].split()[-1]) # Get the last element + + force_index = -1 + for i in range(len(lines) - 1, -1, -1): # Search from the end + if "Forces:" in lines[i]: + force_index = i + break + + if force_index == -1: + print(f"Error: Could not find 'Forces:' line in {log_file}") + return None + + if force_index + 3 >= len(lines): + print(f"Error: Not enough lines after 'Forces:' in {log_file}") + return None + + f_H2_rom = np.array(lines[force_index + 1].split(), dtype=float) + f_O1_rom = np.array(lines[force_index + 2].split(), dtype=float) + f_H1_rom = np.array(lines[force_index + 3].split(), dtype=float) + + return Eks_fom, f_O1_fom, f_H1_fom, f_H2_fom, Eks_rom, f_O1_rom, f_H1_rom, f_H2_rom + + except Exception as e: + print(f"An error occurred: {e}") + return None + +Eks_diff_reproductive = [] +f_O1_diff_reproductive = [] +f_H1_diff_reproductive = [] +f_H2_diff_reproductive = [] + +Eks_diff_predictive = [] +f_O1_diff_predictive = [] +f_H1_diff_predictive = [] +f_H2_diff_predictive = [] + +for i in range(bondlength_num_increments + 1): + bondlength_one = round(bondlength_min + i * (bondlength_max - bondlength_min) / bondlength_num_increments, 2) + for j in range(i + 1): + bondlength_two = round(bondlength_min + j * (bondlength_max - bondlength_min) / bondlength_num_increments, 2) + for k in range(bondangle_num_increments + 1): + bondangle = round(bondangle_min + k * (bondangle_max - bondangle_min) / bondangle_num_increments, 1) + + s_l1 = f"{bondlength_one:.2f}" + s_l2 = f"{bondlength_two:.2f}" + s_theta = f"{bondangle:.1f}" + tag = f'{s_l1}_{s_l2}_{s_theta}' + + output_filename = f'{output_dir}/energy_and_forces_{tag}.out' + + with open(output_filename, 'w') as outfile: + print("s_l1:", s_l1, file=outfile) + print("s_l2:", s_l2, file=outfile) + print("s_theta:", s_theta, file=outfile) + print("N_l:", N_l, file=outfile) + print("N_theta:", N_theta, file=outfile) + print("rdim:", rdim, file=outfile) + + results = process_data(s_l1, s_l2, s_theta, N_l, N_theta, rdim) + + if results: + Eks_fom, f_O1_fom, f_H1_fom, f_H2_fom, Eks_rom, f_O1_rom, f_H1_rom, f_H2_rom = results + + print("Eks_fom:", Eks_fom, file=outfile) + print("Eks_rom:", Eks_rom, file=outfile) + + print("f_O1_fom:", f_O1_fom, file=outfile) + print("f_O1_rom:", f_O1_rom, file=outfile) + + print("f_H1_fom:", f_H1_fom, file=outfile) + print("f_H1_rom:", f_H1_rom, file=outfile) + + print("f_H2_fom:", f_H2_fom, file=outfile) + print("f_H2_rom:", f_H2_rom, file=outfile) + + def calculate_differences(fom, rom, name): + abs_diff = np.linalg.norm(fom - rom) + rel_diff = abs_diff / np.linalg.norm(fom) if np.linalg.norm(fom)!= 0 else float('inf') + print(f"Absolute difference in {name}:", abs_diff, file=outfile) + print(f"Relative difference in {name}:", rel_diff, file=outfile) + return abs_diff + + Eks_diff = calculate_differences(Eks_fom, Eks_rom, "Eks") + f_O1_diff = calculate_differences(f_O1_fom, f_O1_rom, "f_O1") + f_H1_diff = calculate_differences(f_H1_fom, f_H1_rom, "f_H1") + f_H2_diff = calculate_differences(f_H2_fom, f_H2_rom, "f_H2") + + if i * N_l % bondlength_num_increments == 0 and j * N_l % bondlength_num_increments == 0 and k * N_theta % bondangle_num_increments == 0: + Eks_diff_reproductive.append(Eks_diff) + f_O1_diff_reproductive.append(f_O1_diff) + f_H1_diff_reproductive.append(f_H1_diff) + f_H2_diff_reproductive.append(f_H2_diff) + else: + Eks_diff_predictive.append(Eks_diff) + f_O1_diff_predictive.append(f_O1_diff) + f_H1_diff_predictive.append(f_H1_diff) + f_H2_diff_predictive.append(f_H2_diff) + else: + print("Error occurred during data processing. Differences cannot be calculated.", file=outfile) + +def plot_histogram(data, quantity, test_case): + plt.figure(figsize=(8, 6)) + plt.hist(data, bins=20, color='skyblue', edgecolor='black') + + if quantity == "Eks": + quantity_name = "absolute difference in total energy" + elif quantity.startswith("f_"): + quantity_name = f"magnitude of difference in force on {quantity[2:]}" + else: + raise ValueError("Invalid input quantity") + + plt.title(f'Histogram of {quantity_name}') + plt.xlabel('Difference') + plt.ylabel('Frequency') + + min_val, max_val = np.min(data), np.max(data) + plt.xlim(min_val, max_val) + num_ticks = 8 + xticks = np.linspace(min_val, max_val, num_ticks) + plt.xticks(xticks) + formatter = ticker.ScalarFormatter(useMathText=True) + formatter.set_scientific(True) + formatter.set_powerlimits((0, 0)) + plt.gca().xaxis.set_major_formatter(formatter) + + total_count = len(data) + mean_val = np.mean(data) + max_val = np.max(data) + stats_text = (f"Total {test_case} cases: {total_count}\n" + f"Mean: {mean_val:.3e}") + plt.text(0.95, 0.95, stats_text, transform=plt.gca().transAxes, + fontsize=12, verticalalignment='top', horizontalalignment='right', + bbox=dict(facecolor='white', alpha=0.7, edgecolor='black')) + + plt.tight_layout() + plt.savefig(f"{output_dir}/{quantity}_difference_histogram_{test_case}.png") + +plot_histogram(Eks_diff_reproductive, "Eks", "reproductive") +plot_histogram(f_O1_diff_reproductive, "f_O1", "reproductive") +plot_histogram(f_H1_diff_reproductive, "f_H1", "reproductive") +plot_histogram(f_H2_diff_reproductive, "f_H2", "reproductive") + +plot_histogram(Eks_diff_predictive, "Eks", "predictive") +plot_histogram(f_O1_diff_predictive, "f_O1", "predictive") +plot_histogram(f_H1_diff_predictive, "f_H1", "predictive") +plot_histogram(f_H2_diff_predictive, "f_H2", "predictive") 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 new file mode 100644 index 00000000..46f5681d --- /dev/null +++ b/tests/PinnedH2O_3DOF/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/coords_1.02_0.98_2.0.in b/tests/PinnedH2O_3DOF/coords_1.02_0.98_2.0.in new file mode 100644 index 00000000..4c440443 --- /dev/null +++ b/tests/PinnedH2O_3DOF/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/examples/PinnedH2O/generate_coord_simple.py b/tests/PinnedH2O_3DOF/generate_coord.py similarity index 96% rename from examples/PinnedH2O/generate_coord_simple.py rename to tests/PinnedH2O_3DOF/generate_coord.py index 7e799bdd..194c2df7 100644 --- a/examples/PinnedH2O/generate_coord_simple.py +++ b/tests/PinnedH2O_3DOF/generate_coord.py @@ -11,7 +11,7 @@ bondangle_increment = np.linspace(-5, 5, 11) # output directory -output_dir = "PinnedH2O_3dof_coords" +output_dir = "data" # generation os.makedirs(output_dir, exist_ok=True) diff --git a/tests/PinnedH2O_3DOF/job.basis b/tests/PinnedH2O_3DOF/job.basis new file mode 100644 index 00000000..ccc03c9d --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.basis @@ -0,0 +1,47 @@ +#!/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 snapshot_files = "" +set bondlength_min = 0.95 +set bondlength_max = 1.05 +set bondlength_num_increments = 5 +set bondangle_min = -5.0 +set bondangle_max = 5.0 +set bondangle_num_increments = 5 + +foreach i (`seq 0 $bondlength_num_increments`) + set bondlength_one = `echo "scale=2; $bondlength_min + $i * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_one = `printf "%.2f" $bondlength_one` + foreach j (`seq 0 $i`) + set bondlength_two = `echo "scale=2; $bondlength_min + $j * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_two = `printf "%.2f" $bondlength_two` + foreach k (`seq 0 $bondangle_num_increments`) + set bondangle = `echo "scale=1; $bondangle_min + $k * ($bondangle_max - $bondangle_min) / ($bondangle_num_increments)" | bc` + set bondangle = `printf "%.1f" $bondangle` + 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" + end + end +end +echo "Snapshot files: $snapshot_files" + +set basis_file = "data/PinnedH2O_3DOF_orbitals_basis_${bondlength_num_increments}_${bondangle_num_increments}" +srun -n $ncpus $exe -f $basis_file $snapshot_files > data/basis_${bondlength_num_increments}_${bondangle_num_increments}_PinnedH2O_3DOF.out + +date diff --git a/tests/PinnedH2O_3DOF/job.basis_test1 b/tests/PinnedH2O_3DOF/job.basis_test1 new file mode 100644 index 00000000..3955cdc6 --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.basis_test1 @@ -0,0 +1,34 @@ +#!/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 new file mode 100644 index 00000000..d0ac6301 --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.basis_test2 @@ -0,0 +1,39 @@ +#!/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 new file mode 100644 index 00000000..7fc51dd7 --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.offline @@ -0,0 +1,61 @@ +#!/bin/tcsh +#SBATCH -N 2 +#SBATCH -t 4:00:00 +#SBATCH -p pbatch + +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 = mgmol-opt + +rm $exe +cp $maindir/install_quartz/bin/$exe . + +set cfg_offline = mgmol_offline.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 + +set bondlength_min = 0.95 +set bondlength_max = 1.05 +set bondlength_num_increments = 10 +set bondangle_min = -5.0 +set bondangle_max = 5.0 +set bondangle_num_increments = 10 + +foreach i (`seq 0 $bondlength_num_increments`) + set bondlength_one = `echo "scale=2; $bondlength_min + $i * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_one = `printf "%.2f" $bondlength_one` + foreach j (`seq 0 $i`) + set bondlength_two = `echo "scale=2; $bondlength_min + $j * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_two = `printf "%.2f" $bondlength_two` + foreach k (`seq 0 $bondangle_num_increments`) + set bondangle = `echo "scale=1; $bondangle_min + $k * ($bondangle_max - $bondangle_min) / ($bondangle_num_increments)" | bc` + set bondangle = `printf "%.1f" $bondangle` + echo "bondlength1: $bondlength_one, bondlength2: $bondlength_two, bondangle: $bondangle" + set tag = ${bondlength_one}_${bondlength_two}_${bondangle} + if (-e data/${tag}) then + continue + endif + cp $cfg_offline ${tag}_${cfg_offline} + sed -i "s/CUSTOM_TAG/${tag}/g" ${tag}_${cfg_offline} + srun -n $ncpus $exe -c ${tag}_${cfg_offline} -i data/coords_${tag}.in > ${tag}_offline_PinnedH2O.out + mv data/coords_${tag}.in ${tag}/coords.in + mv ${tag}_${cfg_offline} ${tag}/${cfg_offline} + mv ${tag}_offline_PinnedH2O.out ${tag}/offline_PinnedH2O.out + mv -f ${tag} data + end + end +end + +date diff --git a/tests/PinnedH2O_3DOF/job.online b/tests/PinnedH2O_3DOF/job.online new file mode 100644 index 00000000..e703fc29 --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.online @@ -0,0 +1,70 @@ +#!/bin/tcsh +#SBATCH -N 2 +#SBATCH -t 1:00: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 + +rm $exe +cp $maindir/build_quartz/tests/$exe . + +set cfg_online = mgmol_online.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 + +set sample_bondlength_num_increments = 5 +set sample_bondangle_num_increments = 5 +set num_orbital_basis = 36 + +set bondlength_min = 0.95 +set bondlength_max = 1.05 +set bondlength_num_increments = 10 +set bondangle_min = -5.0 +set bondangle_max = 5.0 +set bondangle_num_increments = 10 + +set case = ${sample_bondlength_num_increments}_${sample_bondangle_num_increments}_${num_orbital_basis} +mkdir results_${case} +cp ${cfg_online} ${case}_${cfg_online} + +@ num_empty_orbital = $num_orbital_basis - 4 +sed -i "s/CUSTOM_NE/${num_empty_orbital}/g" ${case}_${cfg_online} +sed -i "s/CUSTOM_NL/${sample_bondlength_num_increments}/g" ${case}_${cfg_online} +sed -i "s/CUSTOM_NT/${sample_bondangle_num_increments}/g" ${case}_${cfg_online} +sed -i "s/CUSTOM_NB/${num_orbital_basis}/g" ${case}_${cfg_online} + +foreach i (`seq 0 $bondlength_num_increments`) + set bondlength_one = `echo "scale=2; $bondlength_min + $i * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_one = `printf "%.2f" $bondlength_one` + foreach j (`seq 0 $i`) + set bondlength_two = `echo "scale=2; $bondlength_min + $j * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_two = `printf "%.2f" $bondlength_two` + foreach k (`seq 0 $bondangle_num_increments`) + set bondangle = `echo "scale=1; $bondangle_min + $k * ($bondangle_max - $bondangle_min) / ($bondangle_num_increments)" | bc` + set bondangle = `printf "%.1f" $bondangle` + echo "bondlength1: $bondlength_one, bondlength2: $bondlength_two, bondangle: $bondangle" + set tag = ${bondlength_one}_${bondlength_two}_${bondangle} + + echo $ncpus $exe $case $cfg_online $tag + srun -p pdebug -n $ncpus $exe -c ${case}_${cfg_online} -i data/${tag}/coords.in >& results_${case}/${tag}.log + end + end +end + +mv ${case}_${cfg_online} results_${case} + +date diff --git a/tests/PinnedH2O_3DOF/job.online_test1 b/tests/PinnedH2O_3DOF/job.online_test1 new file mode 100644 index 00000000..bf69a1ef --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.online_test1 @@ -0,0 +1,30 @@ +#!/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 coords_1.00_1.00_0.0.in + +date diff --git a/tests/PinnedH2O_3DOF/job.online_test2 b/tests/PinnedH2O_3DOF/job.online_test2 new file mode 100644 index 00000000..88665f1e --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.online_test2 @@ -0,0 +1,30 @@ +#!/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 coords_1.00_1.00_0.0.in + +date diff --git a/tests/PinnedH2O_3DOF/mgmol_offline.cfg b/tests/PinnedH2O_3DOF/mgmol_offline.cfg new file mode 100644 index 00000000..7b2dbcfa --- /dev/null +++ b/tests/PinnedH2O_3DOF/mgmol_offline.cfg @@ -0,0 +1,30 @@ +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 +output_filename=CUSTOM_TAG +[ROM.offline] +save_librom_snapshot=true diff --git a/tests/PinnedH2O_3DOF/mgmol_online.cfg b/tests/PinnedH2O_3DOF/mgmol_online.cfg new file mode 100644 index 00000000..5bc5b8b8 --- /dev/null +++ b/tests/PinnedH2O_3DOF/mgmol_online.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=CUSTOM_NE +[Restart] +output_level=4 +[DensityMatrix] +solver=MVP +nb_inner_it=100 +[ROM.offline] +basis_file=data/PinnedH2O_3DOF_orbitals_basis_CUSTOM_NL_CUSTOM_NT +[ROM.basis] +compare_md=false +number_of_orbital_basis=CUSTOM_NB diff --git a/tests/PinnedH2O_3DOF/mgmol_online_test1.cfg b/tests/PinnedH2O_3DOF/mgmol_online_test1.cfg new file mode 100644 index 00000000..e032350d --- /dev/null +++ b/tests/PinnedH2O_3DOF/mgmol_online_test1.cfg @@ -0,0 +1,32 @@ +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 new file mode 100644 index 00000000..8ce3e32d --- /dev/null +++ b/tests/PinnedH2O_3DOF/mgmol_online_test2.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=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 new file mode 100644 index 00000000..e7813777 --- /dev/null +++ b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc @@ -0,0 +1,203 @@ +// Copyright (c) 2017, Lawrence Livermore National Security, LLC and +// UT-Battelle, LLC. +// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge +// National Laboratory. +// LLNL-CODE-743438 +// All rights reserved. +// This file is part of MGmol. For details, see https://github.com/llnl/mgmol. +// Please also read this link https://github.com/llnl/mgmol/LICENSE + +#include "Control.h" +#include "ExtendedGridOrbitals.h" +#include "LocGridOrbitals.h" +#include "MGmol.h" +#include "MGmol_MPI.h" +#include "MPIdata.h" +#include "mgmol_run.h" + +#ifdef MGMOL_HAS_LIBROM +#include "librom.h" + +#include +#include +#include +#include + +#include +namespace po = boost::program_options; + +int main(int argc, char** argv) +{ + int mpirc = MPI_Init(&argc, &argv); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Initialization failed!!!" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + } + + MPI_Comm comm = MPI_COMM_WORLD; + + /* + * Initialize general things, like magma, openmp, IO, ... + */ + mgmol_init(comm); + + /* + * read runtime parameters + */ + std::string input_filename(""); + std::string lrs_filename; + std::string constraints_filename(""); + + float total_spin = 0.; + bool with_spin = false; + + po::variables_map vm; + + // read from PE0 only + if (MPIdata::onpe0) + { + read_config(argc, argv, vm, input_filename, lrs_filename, + constraints_filename, total_spin, with_spin); + } + + MGmol_MPI::setup(comm, std::cout, with_spin); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + MPI_Comm global_comm = mmpi.commGlobal(); + + /* + * Setup control struct with run time parameters + */ + Control::setup(global_comm, with_spin, total_spin); + Control& ct = *(Control::instance()); + + ct.setOptions(vm); + + int ret = ct.checkOptions(); + if (ret < 0) return ret; + + mmpi.bcastGlobal(input_filename); + mmpi.bcastGlobal(lrs_filename); + + // Enter main scope + { + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "Construct MGmol object..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + + MGmolInterface* mgmol = new MGmol(global_comm, + *MPIdata::sout, input_filename, lrs_filename, constraints_filename); + + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "MGmol setup..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + mgmol->setup(); + + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "Setup done..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + + // here we just use the atomic positions read in and used + // to initialize MGmol + std::vector positions; + mgmol->getAtomicPositions(positions); + std::vector anumbers; + mgmol->getAtomicNumbers(anumbers); + 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++; + } + } + + Mesh* mymesh = Mesh::instance(); + const pb::Grid& mygrid = mymesh->grid(); + const pb::PEenv& myPEenv = mymesh->peenv(); + + // compute energy and forces again with projected problem onto ROM subspace + const int rdim = ct.getROMOptions().num_orbbasis; + if (rdim != ct.numst) + { + std::cerr << "The number of functions in the ROM basis file, " + << rdim << " is not equal to ct.numst, " << ct.numst + << std::endl; + MPI_Abort(mmpi.commSameSpin(), 0); + } + + std::shared_ptr projmatrices + = mgmol->getProjectedMatrices(); + + ExtendedGridOrbitals orbitals("new_orbitals", mygrid, mymesh->subdivx(), + ct.numst, ct.bcWF, projmatrices.get(), nullptr, nullptr, nullptr, + nullptr); + + orbitals.set(ct.getROMOptions().basis_file, ct.numst); + orbitals.orthonormalizeLoewdin(); + orbitals.setDataWithGhosts(true); + + // set the iterative index to 1 to differentiate it from first instance + // in MGmol initial() function. This is not very clean and could be + // better designed, but works for now + orbitals.setIterativeIndex(10); + + // set initial DM with uniform occupations + projmatrices->setDMuniform(ct.getNelSpin(), 0); + projmatrices->printDM(std::cout); + + // + // evaluate energy and forces with ROM bases just read + // + std::vector forces; + double eks = mgmol->evaluateDMandEnergyAndForces( + &orbitals, positions, anumbers, forces); + + // print out results + if (MPIdata::onpe0) + { + std::cout << "Eks: " << eks << std::endl; + std::cout << "Forces:" << std::endl; + for (std::vector::iterator it = forces.begin(); + it != forces.end(); it += 3) + { + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + } + } + + delete mgmol; + + } // close main scope + + mgmol_finalize(); + + mpirc = MPI_Finalize(); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Finalize failed!!!" << std::endl; + } + + time_t tt; + time(&tt); + if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl; + + return 0; +} +#endif // MGMOL_HAS_LIBROM diff --git a/tests/RestartEnergyAndForces/job.test b/tests/RestartEnergyAndForces/job.test new file mode 100644 index 00000000..89b4865f --- /dev/null +++ b/tests/RestartEnergyAndForces/job.test @@ -0,0 +1,28 @@ +#!/bin/tcsh +#SBATCH -N 2 +#SBATCH -t 1:00:00 +#SBATCH -p pbatch + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 2 + +set maindir = /p/lustre2/cheung26/mgmol-20241220 + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = testRestartEnergyAndForces + +cp $maindir/build_quartz/tests/$exe . + +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 + +python3 test.py srun -n $ncpus ${maindir}/install_quartz/bin/mgmol-opt ${exe} mgmol.cfg restart.cfg h2o.xyz . + +date diff --git a/tests/WFEnergyAndForces/job.test b/tests/WFEnergyAndForces/job.test new file mode 100644 index 00000000..7a9e49ee --- /dev/null +++ b/tests/WFEnergyAndForces/job.test @@ -0,0 +1,27 @@ +#!/bin/tcsh +#SBATCH -N 2 +#SBATCH -t 1:00:00 +#SBATCH -p pbatch + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 2 + +set maindir = /p/lustre2/cheung26/mgmol-20241220 + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = testWFEnergyAndForces + +cp $maindir/build_quartz/tests/$exe . + +set cfg_online = mgmol.cfg + +source $maindir/scripts/modules.quartz + +python3 test.py srun -p pdebug -n $ncpus $exe $cfg_online sih4.xyz $maindir/potentials + +date