ROM calculations of Pinned H2O with 3 DOF

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}")

src/DensityMatrix.cc

@@ -43,10 +43,6 @@ DensityMatrix<MatrixType>::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.;
 

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

src/MGmol.cc

@@ -1511,6 +1511,8 @@ double MGmol<OrbitalsType>::evaluateDMandEnergyAndForces(Orbitals* orbitals,
 
     force(*dorbitals, *ions_);
 
+    ions_->getForces(forces);
+
     return eks;
 }
 

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

src/md.cc

@@ -492,7 +492,6 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
         force(**orbitals, ions);
 
 #ifdef MGMOL_HAS_LIBROM
-        // TODO: cleanup
         if (ct.getROMOptions().num_orbbasis > 0)
         {
             if (onpe0)

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<OrbitalsType>::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<ORBDTYPE> 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<LocGridOrbitals>;
 template class MGmol<ExtendedGridOrbitals>;
 

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

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

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."

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}")

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")