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

examples/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

examples/PinnedH2O_3DOF/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}")

examples/PinnedH2O_3DOF/job.basis

@@ -0,0 +1,47 @@
+#!/bin/tcsh
+#SBATCH -N 2
+#SBATCH -t 0:10:00
+#SBATCH -p pbatch
+
+date
+
+setenv OMP_NUM_THREADS 1
+#setenv KMP_DETERMINISTIC_REDUCTION 1
+
+set ncpus = 64
+
+set maindir = /p/lustre2/cheung26/mgmol-20241012
+
+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 = 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}
+      set snapshot_files = "$snapshot_files results_${tag}/orbital_snapshot"
+    end
+  end
+end
+echo "Snapshot files: $snapshot_files"
+
+set basis_file = "PinnedH2O_3DOF_orbitals_basis_${bondlength_num_increments}_${bondangle_num_increments}"
+srun -n $ncpus $exe -f $basis_file $snapshot_files > basis_${bondlength_num_increments}_${bondangle_num_increments}_PinnedH2O_3DOF.out
+
+date

examples/PinnedH2O_3DOF/job.offline

@@ -0,0 +1,63 @@
+#!/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 = 64
+
+set maindir = /p/lustre2/cheung26/mgmol-20241012
+
+setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
+
+set exe     = mgmol-opt
+
+cp $maindir/install_quartz/bin/$exe .
+
+set datadir = $maindir/examples/PinnedH2O_3DOF
+
+set cfg_offline = mgmol_offline.cfg
+cp $datadir/$cfg_offline .
+
+# coords.in files will be generated by Python script
+#cp $datadir/generate_coords_simple.py .
+#python3 generate_coords_simple.py
+
+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}
+      cp $cfg_offline ${tag}_${cfg_offline}
+      sed -i "s/CUSTOM_TAG/results_${tag}/g" ${tag}_${cfg_offline}
+      srun -n $ncpus $exe -c ${tag}_${cfg_offline} -i coords_${tag}.in > offline_PinnedH2O_${tag}.out
+      mv ${tag}_${cfg_offline} results_${tag}/${cfg_offline}
+      mv coords_${tag}.in results_${tag}/coords.in
+      mv offline_PinnedH2O_${tag}.out results_${tag}/offline_PinnedH2O.out
+    end
+  end
+end
+
+date

examples/PinnedH2O_3DOF/mgmol_offline.cfg

@@ -0,0 +1,34 @@
+verbosity=1
+xcFunctional=PBE
+FDtype=Mehrstellen
+[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
+[Thermostat]
+type=Berendsen
+temperature=1000.
+relax_time=800.
+[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

examples/PinnedH2O_3DOF/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")
+ax.set_ylabel("y")
+ax.grid(True)
+plt.savefig("PinnedH2O_3DOF.png")