From fa36d2f75db123b5935b6cae8c13b1ed6e2ad7d0 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 14 Nov 2024 11:55:33 -0800
Subject: [PATCH 01/14] set up offline procedure

---
 .../PinnedH2O_3DOF/coords_1.00_1.00_0.0.in    |  3 +
 .../generate_coord.py                         |  0
 .../generate_coord_simple.py                  |  0
 examples/PinnedH2O_3DOF/job.offline           | 57 +++++++++++++++++++
 examples/PinnedH2O_3DOF/mgmol_offline.cfg     | 33 +++++++++++
 .../rotation_test.py                          |  0
 6 files changed, 93 insertions(+)
 create mode 100644 examples/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in
 rename examples/{PinnedH2O => PinnedH2O_3DOF}/generate_coord.py (100%)
 rename examples/{PinnedH2O => PinnedH2O_3DOF}/generate_coord_simple.py (100%)
 create mode 100644 examples/PinnedH2O_3DOF/job.offline
 create mode 100644 examples/PinnedH2O_3DOF/mgmol_offline.cfg
 rename examples/{PinnedH2O => PinnedH2O_3DOF}/rotation_test.py (100%)

diff --git a/examples/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in b/examples/PinnedH2O_3DOF/coords_1.00_1.00_0.0.in
new file mode 100644
index 00000000..46f5681d
--- /dev/null
+++ b/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
diff --git a/examples/PinnedH2O/generate_coord.py b/examples/PinnedH2O_3DOF/generate_coord.py
similarity index 100%
rename from examples/PinnedH2O/generate_coord.py
rename to examples/PinnedH2O_3DOF/generate_coord.py
diff --git a/examples/PinnedH2O/generate_coord_simple.py b/examples/PinnedH2O_3DOF/generate_coord_simple.py
similarity index 100%
rename from examples/PinnedH2O/generate_coord_simple.py
rename to examples/PinnedH2O_3DOF/generate_coord_simple.py
diff --git a/examples/PinnedH2O_3DOF/job.offline b/examples/PinnedH2O_3DOF/job.offline
new file mode 100644
index 00000000..6391fe07
--- /dev/null
+++ b/examples/PinnedH2O_3DOF/job.offline
@@ -0,0 +1,57 @@
+#!/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-20241031
+
+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"
+      srun -n $ncpus $exe -c $cfg_offline -i coords_${bondlength_one}_${bondlength_two}_${bondangle}.in > offline_PinnedH2O_${bondlength_one}_${bondlength_two}_${bondangle}.out
+    end
+  end
+end
+
+date
diff --git a/examples/PinnedH2O_3DOF/mgmol_offline.cfg b/examples/PinnedH2O_3DOF/mgmol_offline.cfg
new file mode 100644
index 00000000..64b5132b
--- /dev/null
+++ b/examples/PinnedH2O_3DOF/mgmol_offline.cfg
@@ -0,0 +1,33 @@
+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
+[ROM.offline]
+save_librom_snapshot=false
diff --git a/examples/PinnedH2O/rotation_test.py b/examples/PinnedH2O_3DOF/rotation_test.py
similarity index 100%
rename from examples/PinnedH2O/rotation_test.py
rename to examples/PinnedH2O_3DOF/rotation_test.py

From 152697416694f094e9e7faee6734a58a8838a8ea Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 14 Nov 2024 13:35:32 -0800
Subject: [PATCH 02/14] Parametric file

---
 examples/PinnedH2O_3DOF/job.offline       | 10 ++++++++--
 examples/PinnedH2O_3DOF/mgmol_offline.cfg |  3 ++-
 2 files changed, 10 insertions(+), 3 deletions(-)

diff --git a/examples/PinnedH2O_3DOF/job.offline b/examples/PinnedH2O_3DOF/job.offline
index 6391fe07..e45fd5d5 100644
--- a/examples/PinnedH2O_3DOF/job.offline
+++ b/examples/PinnedH2O_3DOF/job.offline
@@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 64
 
-set maindir = /p/lustre2/cheung26/mgmol-20241031
+set maindir = /p/lustre2/cheung26/mgmol-20241012
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
@@ -49,7 +49,13 @@ foreach i (`seq 0 $bondlength_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"
-      srun -n $ncpus $exe -c $cfg_offline -i coords_${bondlength_one}_${bondlength_two}_${bondangle}.in > offline_PinnedH2O_${bondlength_one}_${bondlength_two}_${bondangle}.out
+      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
diff --git a/examples/PinnedH2O_3DOF/mgmol_offline.cfg b/examples/PinnedH2O_3DOF/mgmol_offline.cfg
index 64b5132b..42fabd17 100644
--- a/examples/PinnedH2O_3DOF/mgmol_offline.cfg
+++ b/examples/PinnedH2O_3DOF/mgmol_offline.cfg
@@ -29,5 +29,6 @@ initial_type=Random
 initial_width=2.
 [Restart]
 output_level=4
+output_filename=CUSTOM_TAG
 [ROM.offline]
-save_librom_snapshot=false
+save_librom_snapshot=true

From 65939c549466fe407c150b404c5a45cf2c3feaae Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 14 Nov 2024 13:45:47 -0800
Subject: [PATCH 03/14] Add basis script

---
 examples/PinnedH2O_3DOF/job.basis | 47 +++++++++++++++++++++++++++++++
 1 file changed, 47 insertions(+)
 create mode 100644 examples/PinnedH2O_3DOF/job.basis

diff --git a/examples/PinnedH2O_3DOF/job.basis b/examples/PinnedH2O_3DOF/job.basis
new file mode 100644
index 00000000..30660259
--- /dev/null
+++ b/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

From 03e59eb5184a914c10c4a48b0b16404abe7c0962 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Thu, 14 Nov 2024 14:32:23 -0800
Subject: [PATCH 04/14] Add postprocessing scripts

---
 examples/PinnedH2O/get_ROM_table.py      | 28 +++++++++++++++++++++
 examples/PinnedH2O_3DOF/get_ROM_table.py | 31 ++++++++++++++++++++++++
 2 files changed, 59 insertions(+)
 create mode 100644 examples/PinnedH2O/get_ROM_table.py
 create mode 100644 examples/PinnedH2O_3DOF/get_ROM_table.py

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/examples/PinnedH2O_3DOF/get_ROM_table.py b/examples/PinnedH2O_3DOF/get_ROM_table.py
new file mode 100644
index 00000000..5fea1e38
--- /dev/null
+++ b/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}")

From 318325cfbd74956392942745d0975e175ef5040e Mon Sep 17 00:00:00 2001
From: "Siu Wun \"Tony\" Cheung" <cheung26@llnl.gov>
Date: Fri, 15 Nov 2024 09:06:38 -0800
Subject: [PATCH 05/14] Add files via upload

---
 examples/PinnedH2O_3DOF/PinnedH2O_3DOF.png    | Bin 0 -> 22126 bytes
 .../PinnedH2O_3DOF/plot_PinnedH2O_3DOF.py     |  69 ++++++++++++++++++
 2 files changed, 69 insertions(+)
 create mode 100644 examples/PinnedH2O_3DOF/PinnedH2O_3DOF.png
 create mode 100644 examples/PinnedH2O_3DOF/plot_PinnedH2O_3DOF.py

diff --git a/examples/PinnedH2O_3DOF/PinnedH2O_3DOF.png b/examples/PinnedH2O_3DOF/PinnedH2O_3DOF.png
new file mode 100644
index 0000000000000000000000000000000000000000..99135eb68d56a47ab5505347bf23c53a91b6c819
GIT binary patch
literal 22126
zcmeFZbyQVt*FL%^3l$X+K?RYLmX^k#Q&M7sN;d+VE)xq-x<MtSl<tj*v~Eh+q)6;d
zN;hXNp7-~@@A$@d{yP7hGsbZ|&ro1<ueI)&^P1On&Gl4UQ|ZVd#zP2#98poeu7e=t
z!3aV&d|)4ZNAzd^ANWhcL*b@}u8XyY*KIc|MB}!HtD}pDqurhJo>p$|b}r7s{FnJf
zc+cB<c(}St3J5s;`vd$gZZ-nU?5tQA<e;mvkvoFW+(!Q=%a+TsLlAv}%JnOH-bph9
zex~~RJG*~ZM~+@Nap3-{iJeamzJvpd|B=x2^K9{RyG)EnV0>7Hm{E6B_n$Enqw#W;
z{6_vb!AQ|%1=V~uzG$vi<<^?}T!#)EID2b%=O)<&dNPR#*UN<yqj#C6OI;>DGf39i
z%+z+58O2x~gxg{I+}uS4f7e$Y2u8nrj2u9a;oJZJ7yr-Z;1|-8DratJdiN2dxO3k3
z@5-U41+O+XHcrU!i-{Q|OL-Z{%Li9sr!vFkc77e~B1aHIRzwUrBktT8#ihHn++)ro
zD{JOijp@O7B^MPJ^T=2)j<+&OdFF9wCS{+K^0fW(<T#_?orj*Y{lqU#WJoa2FW%-N
zQw$}k`|oH29Rov7=oxWRf<r-$YpS5dw~JCzCR5${JVst)(T4AA+Ea|n9OH{sVdeSd
zX>L(Oy?Vu>a)_Eb*CX}Im+MCugwk88;^I!Tb8twOh{1mFl`Y^$u({dU${m>+9r}it
zi~>xxL6y0^{ksHF2Lo$Nv2|-yLW1$gj*pK|Q*$%VqAaIYik`Q3g`=3sz)2d!m9e?V
zV`Irq)Me1r%xr1D4<q?43Sq0JlQaF5y5@oFR>>(|J6r1s;x1;%DYS6GkfY~xnq&Ad
zn0-kvkYB;-$&&hT|8VuhwoPVtJ-w)~b5aJ@10O6JPSVlokL(oNw0AENJJk9Z7#Qpa
zs!9Vlh*{hMqCMg;hp9n3mHTLzc+Jht4P$j#XN{$>QuJ7<)avfQg+D47w?=mR{)!K|
zYWJv*A!m^0qNNGBT_a}y`FmoHztTKs`jS`UcD6TpSXi#W^?8iwnVAdWYr_g0{Mp*x
zZqBEnrq=cL_Re?k60064NJ~5aC?ur1i@x+vU8!!l`;GF&O81$FiqM8@$B{F^T)LUI
z%cQQ29^6jGx6m`O{S7_oDzsGl_Nhj|mszLyjBQ~iBGUIW`)gC|Kic`}5jFPn5%I}n
zufd<?LxH(t&9T+ft8iA^K4Maj)bkY9ym|5B@DWB4<)1%)vKM*{1(aK+KG;vg&CZ^Y
zE_16M$v3MtZjKc&3}=!lv2JDj`t_^6gexV2h#X#mE!P#_bocY)BV*m6OftrAmBKyc
z6%Zs%_U-;-5EMQiL&y+B`Y_r5>%Um5s6wZtxhTEgl-R<MldZh45J66n5sWliyJNAz
z1QvuGA(g`q*Y0xhsB!a(j1<CWR$}fy$ar#qmXb{G3_PTE(Ac=DQkb%j56AR+(sH}O
zm3tl=B~hsYrM6cNP}&Fk!E`J~TE(Z}YNWYDN*vcB)ulu_9pwW1r)koRA!r6~aEbSr
zt{gd9;ENl+Rc!8D;7%fEMEJ=l5K@+~O^Zy>#L8B3Gr^>(Ehda9Izi-i@B!^3l=6uD
zww`ofId^5Hj&Xv=blkxgK0b$!k-fVI;W#^pxA|eEZkBplaj*JIscJ(3aYI!N(FgKd
z+ZpdKxof3;L5>egSL}aJ_N&A$(z)iZB6`Zi5;gALy{kWqyZw1zuVa~8VP~h_#L2^m
zeq(F&QRbMyZGFkRpO$9NNm+NFld&%|i6O2xx5wVU&nAZ}VhJa^ARj$R=7>P(>3XX~
zGZ3&kz$oEnSsJdUscE3BtgIp>B!oXcMkc>JHl8pX8*Agcbm&mab^r^@cwqgsW23Y9
z2tTKiedqT)htBNZW_88;pHS|T-0mSG@xZ~q<iF;6gju%cTZ4k4B6rD}+d-tFud<^~
zS;xELpnrUPW|jck7v8M&IC)#6%vAyDnyp}(W}$`@X**xa&z~92l5ag$t7poCF!n4g
ztoz67Ys+1R=9N`v5vDdCxuK=;WSJQ?`QVduR=#U(Nw4r78uMx@f~#wy_LgO-gw5DQ
z2Yz~HSx0`>-Hlwk-G%2y4N-#gu1kn-+T^4ypGlc|ta)JN!otE*pF^SDkF+a&MOM#!
zS7$LJg8zzcb~e`5su9zk)$pG#set-xX_ap4*ZxvZk|)1c0YOU6e5K#CquI-IdvMTN
z#KkNnDuth*#N0R)sil7N#WS57>M!yVUHFgO^YO`eaP9v67M{yD?fZY3>d0c7Ol6%1
zYc%(TR4m#YZMQ1soqge?>?oEfyEaS#IpmFE!;!O>wewBNyAxgNj7C;CQ@pQ2X5c9g
zgqU*Z(4lu4If(XS7sbb0!ovDLvTr;hu;SJgF@^<EbMeAOj$*yMxqq@`E@?|~s{beq
z*!aHtN5yAi{N1IC%4MIf`ixERl-udjt*=KXHFNT)QQf=O<kZ(`A+Ueq&!5RMg0rro
zqGJ5!#}S4D!AI1cD27|&g_p--&A9~yKYW^@TSm(4yHvTjxDL_K;9f=~H#wK1S?pTK
z(dfd$!rbPZ<6im6RmshKQ)RrpTB+TJA3jvQyS9H_Q86J!B4{_l)p~Qy!bdCqP7KeQ
zP4k8m?i0J3C6;(+c6sry;3~D4J)443>=9oos!aO#8i~|QXMe|}R;HJfh$tv2(IPS@
zV}n^#R8>`VJUl!y2;1A+V{>!n&b4r9lYkLA_ns>D)xFC*k=j~CDdg}SW0BG9L(`7y
zI~%zSzmA`@+~<5VX6ZYt^{(m$uVV~*uj7Kx(hXwmYRlMQ(pmQ~>(-X_w(%*4SuWm-
zmzuHlPO^(b6}tJ^ew6k5&hzl7`K`|K+06~%%vDs-T&fHCO4M(~@tlO)b(fmIkJ;=R
zBiQYG3X-clh;RD)O>0;B95$wl+n0iNSAz;R=n(HW5S_jscW*b-e!1}Zllslge2mao
zbEbxMeZ%_AnD_;&8DqO8TjHDN>cM%_vT_0f!i^%s_h2dN(Y)~S!T#8iEUi?2m8z+N
z>Jqmx)r|D?$>EO=*a?e^i=r0a4x@K+<;oS(U~P4Gp}G0^S!7syLUa92d7@g$?c3J(
zC#D@~$=%)3!)1jC@7x<AS~&tX^A8435O9*=2^5SM(UWWY<GenWXQC;Zr@4eVaP79t
z*8JzA1&5G#FB2Vnr`l2$6u6UD*!fytbjt+vNmHd)mFJqoROXrZAG|DD`rGafM2X2e
z93hWGe(!D|NYY#jgK%o97W?+GulL#sgu}nHXO%ID`Af%LSy-&O<4Q^v$cireFTA+r
zeUbb;Qoqh!y=5#NVy+oh8B&OKu>1ZpVHw#Op|Sp<WSz;Nk`$>c`Pco*y?FJ$p`jMO
z$=AB|h$;C#<g&Q}#U<Y6rS7qb#KG!2e8zfnNn(k_uIjN<#Yr*NpN4XDt<znUpAc_2
z_Wyb)=es-n5e2{#)Q3fG=cT3n5qi%jH$+)h{CVgpE&e{$gq-<Zt-()^GFF#GU7US)
z`~zBTho2lsp*oDlCG__WT&l$y8L{Wq%y_z+r<Io*-MD#CW+R({@`B1a$+$tO9Fcmu
z-NQ#?%mVbQMjqo=-}kByQy&x@5#VhpIT(b<z{Tm|dl#GA1!_6xrZ>r{52eSNi|JM=
z>vUw|&$z^CS2!n^1lBUJoIS;H<gD%NX6~ep;M;(EhmJ#rDmZ@Z)>VRw7WrQpcg3df
zmssWrE*is;>VAzUicN3Y&i<{Qm2_`U3lC>tPfEu8jS?_2%ItZ$xg9K2L3tQ)O%6tK
ztkiNe&Wv4GOfdUcBTbCAQB~S$OJK*isSt->{vPhBsbrWjk+xgv)vdeKN1uHHz2TmX
z8ME4&_m^n9nY`<lC{^+AZ{r^x9Dj8;lrbJFAz@vqLQ8x8xzc9e%y{6H#8yX}5&>;4
zSk$z@_q~6XS8;09p7EZ9GW_qanM&{9|5%!~vQ99=&AZYrk&~+`ZH|wta%v^4QzVL`
zzet7tqM2><%ky(|JnF3H>mTColXsl)-tZvmUitL7d4)t5zQV^eLP$UOKsI>sJS2Zp
z0WK98vw+f)ZBnj|`U@Y+BplAjP`du@JG0uRluP7khtVBY3Sq&0D&P8otxDCLtX==*
zV$=R&L+x7}9tGT6f@9R$PFq84NbuB2I+f9HjVp`onY~$B#3+HU+;JjRSB;FU<@N<2
z6AB0=Qqkt#uEHF6QKiYH|K_8)8eQMU@|%Ctf9~v9+e4l$BE&CHWecLnz;v{@z2)s7
zn~_*zU9v33u*6RHZ=H=6V{){WB9)&{Jkxa-mrQ>$`zHVfT<boEWSF>!l4<Jvv<N49
z7n<oa>q^^sI2SOS)sY@%z4VgHJecwKiAzw?Ei7ipZ8t8BtKjEHxD$l4d>c+_t)Dn$
zokDZ=l!4?nb+A1v%uABPxH^~MPu_An*4lLr-Eu_pThLwD#6A<}h>z%7ic(WctV82}
zF$C+Pd#mm5XBxcBYE+rIuS(UHTTd5-_a$`R{%WSiD27E96>Ioytrguq_Rpuotc>tq
zp6a+2*{x2p@L`(jkJLXxPgcW5>h7yM7OifP9OHGkwDUdPl-`%oVeZy#;a_a<t$aJ9
z(-2mWyFZu>{q7AN-B$mgui?&#GA30-4*0bE2|#PZ2dJpVRhYTBGPgH})U;%YeF;bu
zy40k>x8bU1#1qEGTQkasV*I@~E^3>jP!>D&fF^w+MYvg3)*4Z^Om}b~*WQQXKChK&
zLH8{>lL~FdOZ8S`NcL&;yW&N~t(NV8M*^nu@u|DhFuoqT2o<g|>;qn@CGJ^v=Ds;g
zxZP&Uc6k|e(M>B1&7*hLwTjQ4<;}-91d?aRqNtF%7@NB^{2*hd$2?Zj%M&Amv^__W
zT`etVx{33vg+-0c*CFXY-(pN?gBLHr6WC|H*RI<3B>rIOs3UO`du`C4*s}m$PZV>)
zhl&qokGCc^TNiWj^Uu~aIR1KnrPgmX2I>(>&*|)F9z)3*P2@;Rq?Rh=7HbJ19p5>b
z#-5&1QL^k)5XE06C+9QC_}UWj8C`uryF1$ugwBX~&);HaXSeGsv8$R$@mnHhrAk)l
zAV<D6FvrUFwKl8yu3{Fy&>+tDCKg1H&mlBD^CO>gYj=yRnmY?HQhcTr+7TYaj95~$
z0Iu`**XMM!w4>QVWFjXwSJk}wZ}g?enKO!fN9R_8JZ-<Kx_akJmTOIU&!0a1Sz_10
zEyE!opauAfeWqoI!?%<$Q$(0Kn$vOonFBT};d!OQOy`X0(LyP4L)o{Ajh>_^ww%e#
zM7_&S2_E<6CJIu^q1LK4sm(T_l&678^qqBULbOp&o{_MXwKWM~%Av!DdA(Y;))xUA
z*SJl|?QXeKQBg_c87Mss5BI_zEFCXEYc+{?vrKBS?`FczMWp9V+B%&TR;1cp*Gf(v
z!pONaZrisvHr~Y^bY&qw|ITV<x;H^Kz*j{|SXejTxD*emmX3*O_EVGnxBmh@i354F
z{|sH^On9fiu+I1A=Dm9hyYtsDJ+*$xNd$D@00%8KHTS)N)Am1m_6!eoZBw3<=XBiF
zt5><Fw<Fntf`ZmJH)BoarlzJ0%)J3@I>T>eIZwP>tcr|ilJ4nQA>}U2bY#u8?Vj=4
zu3=xz9Z0rU7_d0`em`jlhdXiN1k81KxFM`IM(=pA?@T9`UFSupR!Zh@=2eTDaHxjw
zL;V|^m&Xg+08pGtMEp4Xt+JnnDaw+6DC~@Qf#vSzY_%3lUNIBo{_VBVaNqT@*yTaL
zIpXUOCO3t3gjD?I0<oYvJvB&dhNs}=DFII>4TS$R*?Y3J6@p7wuAu<nRJ90^fDI4y
zAMO)v(QJTIq8?sRL9vnS#w6?4hK5|Fa3<pl*Vkr#GoqQnWZHwp?W)XywlSYSXNWqa
zA9CJLHf)9#Jf}r$U-CA;^YBy9T1Wrwxg{1sgK$T^JionMXx(caIq0iqKhTHnd)g&9
zR!vg=HyNe6P03@@Ia6L%s3Ib#(b~&seH}mf6aRNKQ^w!hPV0Q|KD6$-I3*DDs&1zb
zjf~rispCF4gz@ms^L<E^@1_&+D}5`3#-*K2ottrJ7@0uhnne#kF~v`nAW+LDajVD%
zDNur#pd24>T{Q7~;<IL)w%B2ETW$E{Fp*a+;qnbtO|cm!ne+^`)=mM05sh`;RqSdl
zl<g6n=}#-?3L_@&Q`nw{&ucSdu+btuMWaI#&jr%sh1CcyH7wZ4;k8w&VI6pcY?o)x
z(sy?ZNdtrrbCC($p$x2arxIX?;x)u2R%FOJsrNn<5t%XV*s!nuJ~^@5BqQo%N8O4s
zF(>QnY9K<QgR6In6Ux<OWp!3(!~DV}+|xtPw7NY%g(r1Il$VbV(#Y3dMJR(etT#Q0
zMa}8K!3%_)H6_5K5W{k|*Lw=0cXqsrr$<@8ep7w*s#J~{DZlp5=U18pc4gPn7_x7e
z$ZfNQGMo3Al_<luT$!L4jzJ%9h>cx;hcfKSgY<8q4qk)Pj>5f*#EwiE&pZ*V{Y)EG
zeLf_gG}+#AEmf_;O|tX<iNBulCH@kTG89Zr4JeV>RuoXZPMpCocC$6MYgPS;=u7qL
z6(XkYiZp-x2+J@$+?PrdYyyc-K0xOzQdm<oy^&ie7Ae;Atvi2w@Z+6A96xr}wA!7_
z_85wWmK(R$$K7Ts3;Q19XgOP++^!xdVfIyKl&P8}Q4PmTl2IaFLY7Z_RFyr;I*T&m
zKV=l%sCo1V?^Lc5ax{LW%TTRdcB;mhtP=tg;>tPIV;D5#vzWf-UEH4+zcgDNQCMg}
zA|EzDuXW_*%xJh8^UU~+(<k%A@dCa>^0ilB0{3xD{u1=jiqo#Tv*T|!EZ;*%nKlS~
zhWxdN#9W?~0E03h2mJmh69t0q!f7L&gOwCHKiHAUwKS2cG5Y&?Ym|Wc%HiR&a0|8<
z&4Msy@Q^A+<P?zHG!iN*E+WgvA>A?FKPQ#WrK_6oNC3x+of`jiAL+C~@7gO%Grlp!
z-JM-!Gvr7lc5{ZG|B!ruB5Z9h13k9TGo0I=L*_r}3~eXAJl*>dy6e5lZrv+!PBuD9
zV>m)^$$iq8LYA~2uF~JNv5dC}kGD$HFs?QJ=Wx2ByVF%EghtJxAwjJ^!e^R*?Wru|
zK%5`IIh^Hy%pqe8bVW%a&dhh|2862md`Rc@11)hZ0o$cU&IAMBMUBSih3-gKgZMtA
zBC;ztBi2l)d1nWoAo|Y%9rlDm=?gwXt*gw*go=_f+Z-3zSA})R*eo!u)DcU}Z{azA
zzW!+!d3H3WI$b4l?+`B&N2`313O>l0^Rv}+_8Ut?->tt5r2(6s%acE^2?`2spj534
z2OHaUh?wk!K-gXWa0d~MiRTJE%Qf~VZsp;@6N}qq+3YaHd)uEEG!jHqXgSn>SjA#8
zt(s%FE?!hF++9)s@$8qQn1X^rbPe#S<6oYUbw06D@fzaWUW&9^o#x`h*k=Zk1sl)b
z#gbKAT3=u9D|1S;Zco7!wKX@ZU3+}^#+^HNG8(L`ta3CG#kje-b>;?Zms^C}v(;kx
zqi0iNW3y^%Y9c(Xtuxsgoy(h@%X>b${6701-lERcOfBa7=JqPPih*tRp+O&|g}<-3
zv-4%UFtoGA|4y!A%|*P1^!NIte3!C~OKe}aCW_|*eqf(MOF^gpiW??nj$B^qckaA_
z(DQb7pqhu3^;)~!ZjG2zUpl7RGY1}-uCAR{OZv6X;&`+P&Ue`!(s0*+0A<>X_n2yr
zE{pKSsq*!mx)fohhK&{MU47dFZf<TK1gCN?7>b33rGS6GrXu3uUJ8XNK!mC#yqc)z
z^i7%NPF=cL-5WHwat-pMA|l?*&(BluDkOdR?^G<q-;8m{nx%f3z@1tbV5O!J%|<7?
z@Z=fyx=YE;+U)96|72Cj#kw9IB4j_V-97lLRn)0ZZ*XuBfYr_T<%Y0xx_WwgrEMcq
zQ@LifzIHIr(X*p-a{@s@yNANDSpqe8{+TabuOq!P_;X&$xJPmrzx<liKw}uCOERgQ
zu->2V>Gh5d8-C@yb6b16R=~#64aoD2Rd3#0U}s~CA|B6Ey?y(2-R_PbFdI_u)Dclp
zn=s(^jGTLQAot`XjL;x$Ax!oAZo!5`t9=j1!nrO(b;d=OPl2)%z5R(&1z&M@wDXMu
zbt?!F5lep}wMb1og?jJ+(J3hfxZUkJN!Q`v*~$qf!yQKCSw4nGX1B&@WjXqe1rVZe
z#L}*yaMR-3i=QAUh<eSMW@Kif%E1N<rJCs};XR*L)8?z*nWLK}VuxpgpOW4@=hsP<
z4Jh$hw8nV$TA`%I$xR}s-PRg@<&6S9+OR6`EOLaF?ESqsN$k5!&1B<|vLk22qU;3B
zYHmKK7x28QkJz(GnJ|i~N(S_mG=xQGs-qE&k!_X?@=kqGy5C&1=d8#KKPG*5+nmF1
z-|#Vr9N9ExJ_^25d?@Iy-L31)egz5n$1TVh^<8C-uYj|6D!j>pCH*OKLqz=1P?Is*
zkYVD?6_7_*B4!2@2-LmwT^6dnFFwGQEHT<INhWd;{Z#a)tHc#T#^(@?O;T<K6?M~G
zxlK63v4v06)^{e~DrJZdc{t}Ij3*)7Ry5ISDgXUzU1=nY5=J#a6A7{y-#^ERUbsj`
z@fgJyFYC6u<D<!KZ=uQm5OV<A_Kb^lFM;q7_bla``$9St^FF+$AKqNxI$tk=CAWPU
z2LybQ)9y~uJ?D?)!~_TKL2;Ik{mj1#FvKchoAj(yPsdP%QN$M>GTzM{A~}(mDdUzd
zVxc&;$$a_N?5>)YK!AdhnxJ4<A)WjKv^cRzkks${@R`@#|3myj#yd2mXQZX!J+IzB
z$;b@>{|BITNA$sCg0?=_u+fsb3!<-HeJ7i##Z0?_$EOcex4s>T;m7*JLyN$QAU3}~
zyq$Irh!u#f@|4-VD=28~8`qtmaf%OWtjL0@rn^}3`2*;HCAwKH0YN@KncMXg);l=T
z{6d65frtQ-&5F`9O5EHaeF(D#vFEs$+x@(N#`jTl^<Is)$9CrnH;;@!7>~f=U|EO%
zqM0)61hat5N}N$yZ#vX>#-;CiCEhWzn6?zzo1H<z(afZD*+aQ%rVQ`8@+T7i<F%+9
z&B)PU3fF_D5t4CWorrIV$M1)BL#_`*Y{FpF;9^K8X$2ao{7$`hVhCl?gZ>>a#`8mg
z`<Sh66}iO10-IO$=~Ju$H`0C`QgYQ2;3rzq%}upAoB+n-#!NeK1XusWnkO{WYb0_E
zc)yK`*5Ch?d|nS8KbuF&OEr=+whUsi?`qw9cVIGGlO24K+M09|m*}mn^gn;7r3ChF
zJh<v*yd8&MTchAi^;`^xX5Qv%Z?Tf3j}GbAhuQWc^%iiW!&i-B4*JV&TlUrZstDBT
zt<JfRzwl@{hoA1%iR8cr1tLsn-lVtxdBMFfmF>~vgvJyKim4TuX;0yz3JWE7mD316
z0=by~OQr^%H9V?p$Y0lfsLpxlZMHI$0>WHN@ASJ8&oUzK>aYJmv`x6U@Nt4Vw_n@I
zZDcQn$aQXS_WKq36l%#u`>*xwa~%N`6x{cops}W=Av(iyE!}64v=eaj>`OC0jlL><
z`#!xxN5u*S>asUi$%h5ka1n_1nPQvR$M0kKGP-ivUFNESQpKFS&z=pIRH6Y8+L<7V
z|MbYZx5w~VX!L~(rY5vVCwjZ~7bJSs*L58K-dXe(LtI(WsV(?FP%(hugr}gX0ML+)
zl0$kSR`CrrGksQlm@4JwSYdW*q#Rv1lJ~lm@Q-T|jN)Rna^m)XprlWWIZbsKv1x2=
z)pvBv>n_0LLM?4S))dXyTU)ItxA=nmnV0d0tEA4xXe;u1Dga*fE<4+EK@MAMe<iV`
zD9Egk64D_BUIE3$R~yf*BHk8vdzOL{5&y<MyzdT_+4h98d~$N~&O&oJCr~Gwwk<3y
zx+>jmP#|Jb?z{n&RI(nIwj2-{yWLNYvqq()@T))`=Ux=X;9DbhWNrnkktl+zt|N}f
zV8ALyBUKgy(s4zZ%TNJS@)tQc6mdHn4v#jWx;3(B2<!Utl(y8O;V3FQvnsy$(UvTg
zY$~BG_h)?E9&S5=EjL0XmQSUqS7ax}FDfx{s-zF8KYkz>J6TXYiwBi8X{YY_DLy;^
zArVL>Xf#I^<NNpTFKtW}8dmJ?Znf{){(65qLBc(gzjjFn&A+69>T)w*B0Hzy(N3BV
zA3nh3=w*x)##tGm)6|`3G@9Oo)C)r4cGSw!()RmH7W<COP!9FDwCdR^H9q4K*Ga$e
zi3tZnnNxy{pJyX@C<H)@zW(&-Q#z)|xVUU$M^>y6JuU6qii!%=A<lmj1}*9;_%bKL
zNk&HFB3Tg;y_vpJU1#S4&$&U&taPpKas(SYCnpalXI4uY0FjHlyy^|%%snHYLz*3_
zj~?y*N^hR`>C;uvxV*;N2&E2&>l+(z90J%&=^O@YDglvY>Sk%$wWrARmAj;tIS+7y
zk_;<1iB}VdW}AT^*zY6dHJevoFMsZ`t2`**yv2M%LOSq7Kj4;2z2*jm2>$;5o%tr>
zD)`l*phB2OAp(^5*`WANeSgI|=4=1?)2Gf%jWfv-?uv3dYZi^Zpd^EEZ3OytL<vX>
zZZbEV2xVOg^q6#5U*hD?w^n739zM(lbp>D$-O393w16=BifzJ+ii(;X9eRs0E?l_a
z0TMUdRa#3K=+f8?C;96NskZ?Iu;+YQR8&yll9V)#kB^53Gmm!)3=AaUGoG^%tLN%;
zZ{NNhKSK|>;{JmN*tPlF=i~xwTFU0<T}zz@j3p!_R5dmEi(9RWa{v=W?CpMFZ+y3x
zd#s9k+s+{Di$jKRocJsMM7R-rmnK<YF}IgGw3^!5_>1%Pa}DhKN^&%lF26?Qq2BK7
z8?2}h4?7+ow7WUvHZcP?vbD9fwAd~OAo6WX8EA6=khx@JFq^A$u>;Nm>t|?S@c>A3
zpsEk|>L!-dsRH_jBp?$!9^d5n{K*q-@DM~tM`y_Rtwi;X&&=clCRK|7O+N$1h9&Tr
zOndqA7!VK_OV^AVaXHgcZz^Vg<6H~_l(rU{1+2>SC!{_#G$^5mt8rCT&?ZRfB+pKY
z%#ha?8o4d)*_+<SsIjCT^PtjG{HC_pxxs;jzkfeGK61wT>$h*l0Mw8A;8>NyZbsw*
zRE{;0b01enaT*=d=+u-02=i4fRT*lr#;;zzimq!P4VQCnG^NU25Op1v=i%pX^98d8
zcRv$7y$y)=M#+*`?s6xf=L}3ulSF+M{=mGhOxvM9$~7!3n6ClkTlD+)?>wG@iZ}gb
z2WdHOLOkjQ@~hd&0k*S>hODgYSQ&2q6ZPH*mfv5x5fnoit<Vv)WR~?Rokm7R*e8zA
z&<K=o=VYbq4R+(LVl+X_^82y98ya~;-{iFT_b-=(%jU{Vvr{S@e>%87jkrs%`^(_M
zpPoH`UTlzWtUqE?0pK3}B~ov(O{G&Lr&bwyb9Nn>(S&XQjS&gll5^l*xPx9E>L%WD
zr<OQDxHHp7LFPAZ+<+T7c#x|+;NyKN^mefw8scyxpd6!P0y%WM(&^s4d+ZZVUGMZu
zcM5nml6j1o&zx~Uhr?8PaFshXMRDyts8v&x8v#_ey(>XnqDxsV?J+Ax?}m(GUlhT!
z;hfjrR3%BcH#s@Lc#U(6r%&H>-URPOenrJ)hqd{Uv9bY=4h?j?(pw+zyAc4etG1m1
z?}_s!x^C!I0V{Krht56OSUm>|dUs^^^v9WByw9^<NHwj@5K8ZSJ%nn0hhWp5^Xh-S
zpUGPuj@RAf+J>KRmAU|_J3;K;-dq(h+MWMQ<1qRy6ohBBN(<PjCH7t1q35K^EC+yG
z&4bd3`0IlOU7ePKOfbqkW~By_JY|BR__~1$-0m;h?+V&u)#{+%7|IsP1e^Q+eZ*h4
z_S=_0vG$!c+!luJnkC!%%)FDP#2ir4U$%j#xVa@h>%2F{PY8sZ5aG(F$o!}B%RPFn
zVk<E;x75{N2nh+DJaIy2R@!#hG*RwM0oX;1sZa5~hwyGG3SvKaoszIKr9H;5(YoI#
zJSR;DfQ-)T;^G3lyX?i{!2XK6-5TNpdLx2x+eN1Y%uKBZ`YYYhQ)##vAS@)rtJ3|H
zmebz2NNsFl!UVSO(WzHhHn04FbELr0P;vWCinZ}Xm#V5is+FH7FrX(Biu%h}XgM{1
zP)vS2O{o(Ciu_>EZou{kO^`NF<tNXa$%P$<hD+m8`)H6o3t?Kcv}D&UlWVh5Q&XXU
zdi(I<!w4+I7FA{Ciyf@qWnk`!?KO}1Jvz!Wz-&9I7|xXcoWDxBY*D%!;vMi^{CK|d
z{@;7a487)?w{P_U003ue5A)6*8{a18JwGfIt3LG!EJ;|fa1kMv7{U=OwQM|LWn<&;
z>EXeM^t+?qaulAQN{|bxV<Udt&)ivmkD}Cd<m!Oy=Oexo$zDq{<wN%C3u9Pg`eA+m
z59lTs@Lz1zKfLclV-b*shFs^*r`?}HhGijf+y4Fil~KmmBYOAp!^e;9f#$_-)cx<(
zTPM)sDtOR~$Hw&gdZljtV-X74o-!A;X*-MY|N9jVv2ITRW~zk0KDpm}G!(w(S!H|q
zJ5?t;UEg<e=8UNIi$7nNX%Xia$jy-1nnlgn+Eu6oMIi(qDtG$z{&jy9Y-)$4$)6w&
zy9$U8ctnOV$-EPF=$_nO9qJjZ_3Qdzp#)Y2RV^9pbIp{)KSG#QsZS}pk5^O~BiS;5
z;F8$SpO^RQzZ>%H1;bmwIyWK8*7^qoV0*1VT8yw#376qmp5Zyd@S9{$mgZUUgbe9&
z9l+vie5@<8l^LZBByBiOAg*Unl$fc+99@7L<W0TI63WDvEKJ3ORP3R|$k^DW@kAX$
zStcYd{0jztzJrMUBLL5$$4+tN8?ml$Ca$OnsN@<-?F8V~t^#9i->D1|6j+QJaGb%9
zQ58qDGb{XH;8zX}$SZFBn=`!@vl9LzR9vF^lD?fw*U9CbQ10>+hh|h@q60=k0{{IL
zULrg$U$`SnEnND-_N<rA=r{b6<M{F+$?;mmdmlV}=lxa=J`6XV-26mVs`|yed6o<G
z!+%ulwtb0VQmQ?4bSly>DC-}A*d@<LwNDj_5IVVmyRD3Ne}7j6gg>NgbdK#S6`rGI
z<d)ek>8!#I)X^aAD0xf$qpo)8+nJ1+5N0t7F;de#it!!X(h}a#qvRy5uGbt{{UTan
zkPmoU5Z{g*qDB?O!`<1`U`*21wv%G{xJ2$bWgtGWe>fiHZTF}6>r|~67bpEgIiECL
z^VwVyd@aDI(N}AdAoR8059*NnFiF^Y0-<I9uQ!kET>98rQvyaNLZ7}WsO#%r^jd5M
zO&t%W!ZjZ<Z!z@YUwu~=KK7-jrzADV%(!~xQyX!NJ>Qhe$|mvEH&#Vc+h{}{9z_y6
z7L9+SkTLUjRKt}bc;g{@8ONrR!fu8Us2Uk6B;$?j>3(Tpy^i>WF}zi5{7PFgK47)W
zIl?w^W(x{XRDr%M@A|%SSK#um5|Da`(T4VcOwE$J)QJ5v81X1RC$i*@f}LNVh+Xrk
z*k}rb9WCMI%lj{Aq&mmUcy`DJ)@(FPi~(KelzH#qLBR*=Cuq;<@*wTN-6H!<-q*K`
zajF2fSw7^Pa$4sc6%}kZ{J$zi>lEpT=&`grWZsnMjl>leE)DwW6kBJWV-{#lG_18a
zBN{Q@CMpwbW`@r%R}BmVEp?Pc^<y3_j}N|it8rjyMQkQ3B|>&1dwk;Gm7?<$#jp{2
zn`}L%7X16^t4EL8?e376gT`M-Sg$zds6<v3(gaVy_6eqKjbjNtYZuee(aQF-vT{PM
z%{f*vTudg&eL}M2I>HnTbD(KitOArV#{Yxe+%J~x^>)yl_1Fe{vLrmaxH*EL0z|K6
zsgV$;^5(&P_07^A*s?9PtKbFFam@QOI-11zFJ|Fh0j72I_cNLKVB!(!o_iD==K)QN
zI@zc2o;b@J@b}T8Jv<vq?J79q1}Nnk)z_r$qdZ7@y+8*EhWkmme{0z`+}aqJOEMP-
zevDSR{3W)pI?D@<265($;%)tl`;joH^AP%lRq@r8M<IzTRamR6XJm*H`moPSo%w1(
zR-GJdV2bQ6^g_9%y&Xq!EkspI)-LAxGYI1u*FvD?{O52le@T+SFZU#l0r~%>eSN(k
z|DgP5w2<`u1d8aQMvjw7$6(P{X$eY%#7Op|HCft62i!Q-01u!Xc?Zdin};X2$ztCT
zMz+LOgK_Loz^X*f9Yp^ib>v{s(BY%6(qqjc27Q}Bqq6lr5bQh)G8u&4>XYCy#ZGlY
z^3q6<_F3==PZY6>(9zLRHEWkyO`?5ArG2i}Pfaa0ri@c7#&gyKG-XC6^0YIATEBAH
z^*^*I)FAZ~rWVv~*w?Q8evW1o4~WOy;YSV~(t(1tk%x`#i;6x)cGy$ejOsbDX7Ff4
zOtz?_)ZfT?r0N4Wfge47OoG%UU?dIz5a7KXz-fJ^&2AHy9Jn!3UuiCo(D_M8XO&{9
z(;l1os*rj+cz2f~4?c_;!$%+@$D|GVLxCCO77{9Ikz{6Oe)Q-Om#Ande$=^oMy%AX
zLlIzc4_H>Sl$d<3gFo81+S4hc;m_C@mw-TFdh)u7;llxE0Q<Q5p)WrFk%BRPva)!o
zIb}mL;@c&U#K9nQ-<|cV8i@w$Z?BTB4U~=)VAW<#s;j%6S1BNQAs_~4*rW!OH3gtW
z($r#VK=`e4=Yf;w&gDB4)*1kk)djyXx3u-YEXbA1sAU+6G<(eR+%F9c)m0CT-Fy0S
zs#h}A?J~1c>vSSS9AY#)H&0=S5e8#|%dcnqTv2~nZZeNtN~^VzsdlzY7_(eKNy$50
zKcU>EFlc8nH9^qgv1an}(%Sv|<WN0c1S5FR)+lrCa$s+7FNz8|Q~lq>#l;<3&B@;X
z9x5bN6_p#_e}(h<g9n@q{rzhm!#x5n`#&D$p3qA7+LBxM%GtQMoA5DwRL{crapnQJ
z)m8sZW0wGEtDxzp_(^e-M&bJPix)2ZuGw6k8!{&niKvPI0Ekru^`W60M&R~_{^oW>
zL<F14Pr7Y~g)^w^Oy90{_Y_>WeB%?(_};`yVa?-k+Zr3=XSi6dsA`l#PVxsUj=q6$
zKu&cgdi8C4uPoymtEVY{e?!>;AbHgzIW+Ly`6j4U(~i`cQ?=Y}(pBT*u7ZEXs?-To
z2=F5M5hX4|vM7~+-dW)K82ug}0z|e)t7uO)RKW#XSVZq5fgmkVVWrlw^P76{&oLzX
z2zgqfA3#-f7>wUcY1Pig6c#Sz0O*OGFymPr5gDmAyF9<R*s-xZg&n~5<bd#ut_PGq
z8Boam(8rMn+dq5q@{Wa?-*!sSpgQ9bMD#2BEphr&r!qmcHms_J>RWPCtQ~3_JHXFv
z7mox}nEGUtcz+PVzO&PcJCB{=28y|_3aOw+N!eEE`CyPOtk;+9gjGB|%TVsp3=Nak
zj$(G}ZaHf79fPIK()?i(=MX)5R2b8Ks|dOHTd0QLvOQ!Wl=fOhQ;2xie`lb?W)c|+
z_z-Oxkc2Dq{J;;PB5WOJRgD2<+k}B3FpTeLFaZ@fWgtjzui`)p>LLW~_yZR$+#Etf
zD6zf4Om|AMKo*$|hXA^)?pDOmWio(FPVv-k7Ifyt+wHF6UimJl)D4p9g}H!D#iq<9
z*55xpK{|<3(ipXpfK!huMQ#XM3M8`S5$C;OT2^Kt!K#pvS@keWvK?WxEW^$@ryO!C
zb0~O~;**)39x#!NZ<mv*Gec%W(Cfl4(3$Tn|9ut1`wF<T4cxBq+rai@#J8{H0>)?T
zhKi1Fl}9rlK3+Ji;B;^`DtO97MzJjubpr!DN;-?PjKZEpB#1?8rtnsWw{mRbcFV7w
zlS)wgDCMafv{QVonQHhkI#1!#1RLMkHgUv~AeQr9+rI1F7H}ctwX0_L?nAO--n}G%
zTm;TI*w=E?bFqV9(pOo$wA`VUv6&>z6?(>oX}@870jAkGRVZVr9~mAcLAOT6vul0P
zhB26+<t{wEIrzm<{FY~g4<^%3odV6^%XK@Vj(dOhjm(Klye%vl6t1k`KaiIZ3aYa-
z8Z$Tdza1FHEU>+^YW>g!%WHaZdCDdxX#MNb*;91-J_isQz&D8a-!4PzPZJfFg-@s3
zXc0H(a%H~_4XM-9GCB7#2UuR^K-#~t-$Gnz4b3stgVhH4C>MU`j=I2ai*S{IP1S25
z-2tIUKLwcchRGfujZX1x;x&m8@yzp>PeXU3Bp~DYW5qUV;{CT`I`CjWY8^g?q@ZiV
z^k2!XN<TsHa|Sz^>v2FR-Y=i?NBonwT22bO)Q$mvdYKt^ueAEyAih~37P668qz1)2
zsx{BcgMz7%WE%sX(pRh*d#R<~2R#F?&HQhWx|(_CdRaB%ZYdt1uwC2xHW>E7-|<7c
z!>roCxa2!})f4Dd->=NHejmOAXxXXn4=A#iQ6*WrpLu7`=TZ4yn^9O%RRxub&Z3C?
z2Z%A(9><M!>Jhbqz^5G>Ls<_#Sp@O%Q_D3}<1<NE0sO@og4o8P-&N7M{*dA1?<t2q
zYWf7R;K%3;T>M#>T4JeU2Y#1~Kb8Rd>aV~@gN0_nEVJ!d?7vYG&nC%(W^L4Rgvc|9
zB`$%lFqS`a`!mh>-Mi{1j)<UhATJThr~>d@)Y-Am1&3K#|A_SObsD_ZOiricQ6<pG
zuz`d<MIzJoLbE@(uDXNCC{CU{nf~MVghX<Pl~IC}R{@lmbf-^evs`;@3p8XH-^pN9
zB->~<Zz&q!hmz);+$^kw>17Cn&F$?yQ0`)Pw^t1}<G{(KDNr+iKZegF6RhpoK<7Mr
z`t$~{qZB+;!v!F)6jXlWWtQa)J?D^en!G;XI5N;*<FmN*;~Gs4Fw<RqLq)uNMn%`r
z@({b&DxM2b*cVE1z!Eor^pGgKi8wblDEUsDgDg}+dm_Do8HH8+hsl0?o(3~?PEHP^
z&Q=Sj3Mk(wI}sHZH~FavoCQCB;{E30IUg(IJG9iMx>32mbOH|7uqpgXMEahAYbejC
zIKqjpvAJ2-%*<@hX%%5vPO$QNyZ`a29+IXgJ(`vYm1m7sL7!Qi@FR|dXdHraRLwjH
zleRJ=u+=}i;q!kB1@eIwm&&<d_dnBql7SNDf7-&H3&d=!;;ue@+8UkW+%Z4WYy5tg
zdVQ0=t8sl`c`>HHa0OtHVJt#_0qMRbY~Au>?eA!Ae^vnBu7Fpki}>HOQl9E)gU(d=
zIjPJP89xKx{^{OgQP6((WO-mGGN2;ok(RD8D6P2vpM*d3dw)vR!GDbN4A@;L7I$}I
zR;E7yp4Z-<kQ29EU!Xoi_xm3QTvkR#Wff(`F|zke&w%{Ae*NLKwgE8bp~NW*%WJlI
zaKqW77%q%$Wpy<chYPA)uB)qaUVdcQf|-ixFGo8Vl)BF`@q>9W<xEl}YrZKlBe!q`
zT%7cQ&j0NzJN6_|3W!#~3Lf*rSLjZi;s&c2{ArS$`j`w2Zr#wobRGl-#{@B_gq@Xg
z9C2uOrySL^p)pMpg?%7%m@b(huHu+#)!j{8#tfHk`V0YMNE4j!djF5s)?_I|U^Kby
zR#C%PrML@pVnl#aW(nFs)bt`yWIWK+2DGAquWuEKasw6`Iic^-)qV!)2fu*I)2G1C
ztN8d><0qX6>F+NY8txnRpCt2e=YehrE2A;9uOnJQ?z_nWTCZEJJ0i7YbD%|L3<Q$x
zjb$(nsU?|9Bjx3;uPaaTs^GRRK)Iu*Cq8-wxqJ+QkG2EUVvc#8B{_CMTOXEMgrP8~
zUx%s;p-<Q6A<mSu7Y^+x#7h~v6*JsG!odfHn7#*TEGtE*ImRl5-hVF!YfC_MOVeyg
zU<bm=T7pl5QKs+kQ8MjaD1i_x8iYrW<Y1&a2CMWX9d^8t<pb!YZGepO-3`d`9$bqx
zj{&zS3fZOiMp`CW-8cWx7@0AL^MZ<?2{rhV`Ukw!R{Jet+LHt$ucwMK%OLdV_s{<R
zN}n-U+q#e<14eD@hnJ9^y>DMQFK=9|c<6}e42|q1t)pbaz7WTPKbJTddLH`@Iqq91
zYRpBE$>q)?;>tnWxkmiFQQuy?3S+`LX#{8O4Kn+bGyV0o)`bri9F%=CK&#B#?v1h+
zffFb8si75c^GC!Kt>IjsWol&nx(N=|5_&$N2iN~eDZ?V$8!BosQDvO!FC!u`21kcK
zpwT+3YBvx7XN#LbjEGzM`oX=7f@Br-rdu_q7hp6tLzvSy=T@yFZ0n!;n_>{8eo@4p
z4SWG$lE4!?@YQp^Tx{vd&K0joqrh);>Mb_4$ASi%eQ|MN#Zz-e{u-KmE55MDD)Sg7
zOi-%fpV48k6Uhotb#f2NiRy!-De3=9on>^CmF{HR6?Q%o)qu@}^9P2v_U!ZVkm_(-
zUwP?IdfQFkblKi;Sl`ys&3H%&W?n8yV^Y5n>DSM#7hz&kv?p2ukv0@Ox#p}dW(;i?
zC`Wkwd7^=uplZJ!qSk0P_pLe@aOH|sIL#rUDX7o`ZJ4Rk?w}K?w!Wt9Qmdc+oWCiW
zEnp?y$D{uhD~MotC1w@#0QU7w!*jC<;^MI(L@LU*aI|lYpPr6=drV-+cD6rP|E6un
zeX2~QLU$Wl%GK#tVB4;r>V&Qpo}jR_^X@LU?B8|1dDIaEzPle3008N$JS1`?+@GoT
zf$dP#?#`#%9wi3@8p728YC-wQWT?kxWvQ$buF9#Yyt$e<qo6=na}ZHM^Pdc~J}gh(
zNH?ohm*P~~APlv*-FQ_=*<VvhJgBZIwwWsH1a^a;l(?vq#CkHjMQdWBi@Lyt-~}|@
z`m)FjXo33tg2CFDUZt;G>)G?jnws<)>0*BA{-vdJl9pvQmdc@L)S&Oo=QE}%E+BBb
zYw6Sp@_Xp7BeBrtx4I5;sV+NbVsp7}Rsxfp%EoVM+7->DOHM&S!n@{fMNx?kSAhtw
zls?N~_bC;JtE9Qa(>&cUb}Nk+j#RxBc+W2mVW-jqHm^oS++jN6HTkkhn(j}BR%^x(
ztCG6bt-!jdcst*1=rlXZB+`tjjtRH`ZL%Dv@ze*RzFxcqYjsj<HYtT))PX}mA^BCS
zHUIATm`V&goK^nOz8QSb4(ottV#(dNB@VGf2~*sl^bO1TTRC95SZg!rH}NLvePs<H
z(?T5_mj@~%rRYy8fBd-K)GRzw=v}WgRF}DbZvOV-;#c}+v$odnmmbL0ZK#4v&hFRy
zg}JdAr#6~2=!_ymG-DI&uVNO*Gp6o8pjO@Iee(H>sDfI9LkBL+xX+Ch^DYW=HmEag
zr|w1##K4q*Rjjgnz9k-SfoP$Iipv=LcTjDtSPazEw2bwC3ssw2UR)K|I+W8kI>4*u
zJI+80Tt8@&>T|cOGkprdA7dUU_)%EA=QfUU_|)FfsM_+Uoo=H$#j!kXx9|R0i95b4
z5fSRaxnW6XiRQZz;21d@8Od5|$6D&_L~adzJW`kBhF+ZDpEaVCK7#g`WGyb{Zq-Sf
zo1>=FAKF}(Z>p;0?)W%~3BPBPBB+zmp16RWy!9#7@~Vw(Hjza-Fsb?DyeXx;_$nxm
z4)wug&>^Q|YinB?#sF$2^d`smyni;Ev=7<;^qJ5V0nV)HjXYlAO0b?<i@0WK(yg8}
z_*SHTi#LV&y@TQ6<w|#|{fVo;buLNM0o4Nv1xj#N?e1(kc=W%&{TXa$It###m!H#P
zL89pCtjGs{zRj3=_DprIUCPgG;^nc9V#9M{i3j~<eB%riMzpwV3t@(-cQJ_97nxvZ
zX#DeuNC)?R_vB?Lydh~zq+LKR%elA+e{F`yHCSt10?mZwLP9NLL&7$ypND?D9=Bua
zx{cNynO=q2-=*o&i#E4vYP_!=kfxUjX0nIQ-e)gf=w)dpdjKBf$>$L7MdZ&tc_yS=
zQ(|wWI6XTwbW<rT?jtj&W?bXz!iuX<mKVK!l#9=P<uMwnb1Mx_gJ|hObB%&R#_B%u
zl*df+G$M|_E`^?z&;gabqb^`m1T+I^<}?MBenEDVLDgmtt^JMs#SUQOpDxL+b18;2
zC;tnt_2<vdyMl8q{aVPhr1vFV78WY|uBJQ?9>J0|-7D6bJhe?o3w>6w-vr`1XDV+L
z4wrGCoR(2oSfk=a7{4DwWm^5xsf{(gpoe#@`e%d%#+{enIQQ~C7t~&E1<<$zgoZ(7
zJLpQ_vY3cifqpgc@rh*g{Js~9+SesP_+6grnSu__eAiTHqI~nmZL!y&&@4ZVv}5j$
zWM5TPeHpmUe)H8UpD}n>0%$R>fQddaHnFK(WO+Yro#B{v+TBslcjh<?Ngup<W#!<Z
zL)?0|Uu&$do7wk?7zPGr^!ODQ`JX&}Gc8lYJvSie_O-~!$hQRNY7J;nKXmZmU95zw
z5cw?`5Ab6wnngGb)zyO1r`kt<0%7O6aN%muCSk~9u6EU+qV^3_J<twZ?Cc64(3bjY
z8NMgXJPQx+fzDS{;Ul@(e*s^k^Mu<(RZqdT>Aib6!|M;V;x|_o7Fs%_-P~FnMimpq
zTvl|h+^h6*U7NoU%3#5D>AlBpr=D<+x%SJ-lMwqdgSakUn)$8N+VxdBcy=Qmuq2oR
z5|GK)6_lp`qK-b<&3S<q8X<7MX>D_}d1wgN@B2TPA*ooY?t$zZ_by!U_aH?T-ha?q
zRU^H4K)iZDnuSe74_ZJDAAb5mM7MEm-KVj=omXw#ra8u1GVb?}xM(YWPCviID;&p<
zW91H(R0oXD>hX}6%5hbdXg~Wutq6MHj0Dgb6&019mzNh2APL3`8re1Zh$xWx<ArTV
zQ-!#00o>M`hYu-HwG4byZ^0aEy6k`iL*%Ln0vq^wW#<AQmF?DeF22;!X^{_KnbC-Y
z{yUifJ3A{*?9l_4Wc;B><g&WJs!{D^QXIw^T;K7NVsd=EEjbs6x#$wVWfHr}#@1rH
zC^ST>=(xER0nNsy0{)8)@Cl<zT}D<`E(o<%vsF_|ow{23zdr9r+Q4TAQ9m<4!v7a#
zUA??;+-GnARvtaqa`s*bO=9I1@DT#gGp@={p=SLfq_EJZ#uU(EIax{=GC6KgVCpj{
zfwb4P3SYXUof^22^B|WzxR3gXXa`Qj$^H}Ot)g{_Gjvpqjam6_KAW3^C%Q7<k6esb
z!-byd7~P+=*K_1mQ=|UNq2n9nC-+`Q>x9$Q14~DZ4SLn;si-)%K=0bzmd2(FXfvRA
z(oAwFQj(m!#S(94H?DJ6QIRM(IG8!ne*A$<@MBaN$c`na(0Ivn@ltFs&xvEV_Ia<Z
z#ij1-a8T~sm*BLtMMj5$K@xP~3tw1UB0dwKpl@yygH~P_&0&gn&+Tp@$Z`ZE#F357
z#HVT0hnk!PE?%0dY^TYV01gdVHjobn9OAK6!O(f*=J$X=c?_`l2&ojt3x|2~7=Ljv
zsDEO=>1p(42cbQF*h0uA&Dj-PWJqub+VYQhqd)uq$QNm6HldqZ<nHh5r4Bus;BFxO
zc>RDWXzMPd8bcpbDB3gC<H6bNO0|DK^wr-0PjX!3yL|{!0FG!2unoP1-XWkN3Zaz-
zV_f6?0j>Ha+{g8N`mfQ*q@#`%us`%w+<l9RvB0}Z=GY_1us!;Ak5X|6>%i0Ir4Q5v
z_<>()D%#NOC%A9SM+6~YZ8XU3yapn4Y8AIz%VziW86E1pUIzF36xfPARA&HZW5Hps
zuBiC<)2C0|19eUWvECU^^j#vBn2EOJF7VpJ`$P;Fk;R3DsMuJ;@c=gBpT<b6UZwl5
z8PBmle=hRzX+SYrz@Me~Pb)piQ`RfLc5RbT#7l2TZg?IMI}hYT*K>ion|5=lsN=W^
z+|#dLzqaZYcj#6J+8~-q9EeURE(`lN)c7mYyO+_vPZ$@-hr3vYy?86)(9J8eIbEzj
zTkN;$jGgVTOb__;{WwZdgP)-r9B-tWOt7pt05i@$v$W)a$`sH8<&i58*3zP03N7yF
z%Sn`#pZvR{{drY#JI0~?`*&^7|8r~XAA&A?lcqIJ@Y-lYc%fS=xIbcfnPq(SAdI2x
zPKp2e0(bcHXU}v31MBb|9nDMROh6BS=_Q(J;vX-$+Z5<MztK)FP?GwK{?+XD#(4q=
zkJ_5kG^D2of|t<(A1q$)O6OoL(f++<VuuosLEh!LpzlB|yn@zu==&0}@4UDdqK-02
zXM?xF+!YopZHING`U6;#Y@w`1Uzm~-wBrL7FO+<-njMpJ`||{itni`%RUsKqFj#bg
zm0@{(B6ZK5cj@;1jn%pK18{z&DFM0gLYqIjG^nl*{eVJ*XG&OD{C~ObED6vKI5swx
z0pt_tI$F>+fwqk~EDi4L?4bP=JT<(gsg(``5u4qb18oWeZ4(8h)eZyMXT+ZW9YWx(
zM32<!(ZBxpmlI{jzzbv);1>0&J!eh=OLp{}<mBH5X2`*T(G6%p4@Mw(C4r%<>wA>%
z2m72vj{WF=m^QuH<KEH#0fe_dpEa43u-7gB->J6}axnOHf;Xc8oMRx6+5i6Xw5h4-
z-9p@z2m7KoTjr|gWJ9_uJ%~D`P6PuoGBQx%%o{yaqMwF_>cc^xJ9aG1d7z32t|3p`
zfan~R0pKHeS1b2XQeOUGYyo;-)!%PtP-&#sQ4AP1c>4_O&1^6cF08H^f^)AFTKhwK
zd<qh}Yj1yk?3P;IYMr}qBQfj}Y<<k#uXO0kY%Z{}ir)Eh02Uy1wPIsI{AdF!^2RuD
z>#q_W9UT!Qm@beQOT7QyL0@;_KHieEP!Yo9s|+62*LN%|O7D&;o`dZQ+rvYrOA$6H
z>Vesu!{Old6L3!UT3tZpBwV30$P>fLt@gY7f-|hdu~*2cw@4wBQB()qLT>Kv=Wg2h
zplT>;xLaHEhTWuz_OHURZ};nmlh))vj;fiT?ft=f30#<m%+A_qUtK^o7zSU%3y8MC
zY)UtagMI>VKN+CA5MKKu0uP4zTETMZaWfQqqueq@`|iqw?{GaCr-7rTjm=*7V3mhm
z>B<l;2<KcevhN6!^iANVP$~-TcnE-(105|W$xG<Q@JdKjw%9N+Fc8xY@sz=HYs3pF
zm<Mf_g8Rp+oPcUYTvAficnebV(B-5r<o3kMiIn0#Fo_aZAT~7j^z=Ln4dpW`hghN(
z5e?+#ixbT^KiqlL1&?Px;|Jy+Bk>{s;%iem1*C3_o}MMZdLz{0sMZY-!30h~H^Bur
zQ-7#=V<dzb{7>*{-J)5!-Su!VQs%)D==t`8GbkVUv>ylQ7y~v}N+ApA;Q4uZOW=J;
z&d}e1`?S5glqF{*-tT%X9t-qU@hqqtA>HU6!9?QxY+hLg%zxakKZ~+vCa5{C*ctWO
z!6H%c<Qij1JAh>PR8Uu5l8!D6FjW}Ch8obKXLZeYzKxt+Lf&)(z<KBYWDR__z`%-v
zc`yQtzQ+n8-ak(*@EvTAC*j3Yn#q#Hoy|EzW8;KewRd=kIMMKGsPRQ-eEL(6-zGZS
zr33|UK-jJVi=|pO>hJ|}7xlVDH`IYDh#EtyJZIiP433Sjui;1QEU*My$@2W;2TsSg
zqc0q3*@-v_EykVgWnx%09Rq_Ha7Mlbm`ntN7?`-LlWxpX<<O+AKAP+;&I0~i)MHXR
z07St8ekI6f5qe`@Wz|9B*$gR^_l<k2y$a3z|0>ArZY@U7MkXY1%gN!~Rtv$e)IwbS
z+}R1$orv2XWi%~=THFCfK@o5wLOyFE5C~mk`q^v3PGE43UvAr}(bd&$#>|2b^v2Db
zH>0uro&+b<=L~CAFju#!87{r>1oX&KI1`eEf5E5C&C8oNU%UYki$|J6RJ3a>)yau8
z+nEz-7PO58ODP6!7v${;(OdNh)Lgm!v=2#<O|o;4j8L~Qzy$txdG^39qYiC;+?E@J
zS|S)6JJPN^D6wijwLCjOU>AjRzt`GLdY9V?kD<DnOA|Ee%H)g4^h3SGU2}VVx~FpY
z{aJWN2o#MSpdg4sS5J9PHHUNUvH|LSBz5N}D!w@V8e9p$(bWBf<tAYV!v3#Tt~I8~
zGmKU-GyZS}BbqTCP6ZhZaCk#hz|AcJa#5u8V!2F6m8-W#Q<Z|-!ViLiKpRA@A!3oX
zK!KuGE;2SlS6zwFuYv=nb%_12Wy(cbOLtypmi=D#kL%a>doRy(-sd^zJR2QAKvdwg
z;ulfmv=5y$EhPJmkcwB#F1;cXIcEtq5Ny>O`7ZO*@e-&`_sHwp&E|*h*T{Z>f;b+^
zp$&M&HZ*VD*CpqH9XQr_kpLcSGc(+;j<6Tjh(_m39R62z-9F{1gI3r8?KPPmzohoo
zUrbe1Rb^6X{lvsXU|SvnCKOwJKa}xInt**+JOk_%UV?9EOgxpP)SJ3Bg`vGj3ApQM
zQXSy%^}$MvbNO|qLAl+>r}nqJ@6oOYtE2H7=bdTV7<1u8mIe2&B3c=NQ2octm%Qwq
z!gx_7lJ~}eVb}Dy>DwTw5ur0a#OUftRbZ)#i;IPi?FkjQE)#kDaqX#?%faf@mDG&z
zG~I>3QotB}f=s@|)pd{Z%=zV8s?ZCN$4Emsh$OA8wU#hTuk1XYTV7mz6<IVqaJA`u
zHs7kWXmF)VsjjXrU%*s9>DGY3E~}~WCQ7n-?bsv;$<_B(U-jKPH~{7}rsX(V*--=^
zpm5mFCn|#l0z*kz*_!ZWZz~dh`bUv=7OVZKcE4BTu3bGK_xWzSj)eqCvIlzo2nsUo
zr19aSQGM4f$58$^BRv=a2d#G_0``M_sEgY1bWDGzZjGnsmS@&mKH{_&dUR!pj}iER
zJFo|~TECcllCwPenR$=*;Y9_DJX(Nf50d^p1!}$$H@pkRe*1CgE=X+9YS4<)6QnZP
z&zG`?bL~n)he$-^htWz00wjghK&Y|o&d$zvU&AGjMf`C3Wnk}92-$wRMn$hIhAW!U
zB`Lse$-povL?k6(3zsn6%PT9nIWQ;EaT=2I-Q3;N(erXPw%aXRX#|xPJc5rCckhsZ
zD^5cmyl*zhC2E{78f)4Ng6+2%V~<pxkM#~N2`G~g%5lgz)gosu9%PdkfUIeF@*QCn
zSzBM9gFjJ@a{ujzB=9S=rE9r}iA-4PrLo1RH=p3ZRL`g1iDS|T@qpsFR5UaIa?d2R
zB3$5P3PJ*FgF9ViAumm?gKCf`NTj5cH;2-}FkvC@9!(!4b1?I<I2Ce)?g7<-0s#^w
zbetK8R}_vFyU_`SW+r5q?DhO#o57ID^5bdxUQJE4fJam+<V^#;y<<V6^GH54XZxGq
zzY|L&5_PUaD5K7xu*3zvR>3VVcElbGnMJD9X50MX^Ob-cOa&_?ymmZ36LT$IL*qUX
zz&<f@|06no{wEFon$cCh-rn6dlZTV#wIR5RcVaL;Wsq<DOiMOC=Apc_^Z<S+U;^o)
z)R6%!e(&OjLwmKM4U_6@U0_`2wy8%4n3$>0Q86Rs95rI=We%_#6ZwZILw$pS9--W#
zg3ct?>(&+f&7<g7kyTLPo<jahj?+(A@=Dm~6;Lxz3j;U+DlI6+C~#gm00RPT3*r_B
z<4CA*44<-Yr!PeMtH>Yt^E51vGZQ6W(em9<<|~`BINJ6H7dJN5b<7Gg5Nx9BZ=V#F
z9ZULIJp>*wg>9-^(c8q2L6cTbBqkLdlQzzDlRhLvFxu0J?o|T@r2vp*OHtJCrR&D*
zwlRG}h_k5c&rER|``lOxocqEYKccP-15wgeVvMhw7xtgNdGml+MY_3F^TE$Z9fYF@
zQmqVilV!{g7=l9n9GAQ75NJ^d*s9!|!$cJk!7+wZRUrN?kw?1Z1!Dwx$94lV(Z+fs
zF0Xag4rk)<>DNUtM&$99GWtjz;gSi{OW0yS=GPE>NMpIqe*we1WSTMe?iO5+WVm|O
z!d|lD?W;n$ZqDCAuuGnyp-I3y@z{8V8*yt%8+useYy-IE@!a1(I&I{B1x`xQER<6|
zK|w4i3HE_%)SNvSoFKt#5J#r*N}Se%u-_7*f$g?<@hSj`wwaMi8`gRP*VkK`ZH$`r
m#2^CY0#bPR|11yRisouJCC0X`--RV&gnu2qv+m2}vws2alm@o|

literal 0
HcmV?d00001

diff --git a/examples/PinnedH2O_3DOF/plot_PinnedH2O_3DOF.py b/examples/PinnedH2O_3DOF/plot_PinnedH2O_3DOF.py
new file mode 100644
index 00000000..9c859f0e
--- /dev/null
+++ b/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")

From b616e8e24f5af01c7a25436aedc375e4fa2f459d Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Mon, 18 Nov 2024 19:22:49 -0800
Subject: [PATCH 06/14] Add restore check. Make H1 shorter than H2 in original
 system for verification

---
 examples/PinnedH2O_3DOF/rotation_test.py | 28 +++++++++++++++++++-----
 1 file changed, 22 insertions(+), 6 deletions(-)

diff --git a/examples/PinnedH2O_3DOF/rotation_test.py b/examples/PinnedH2O_3DOF/rotation_test.py
index 79f1aaef..c92194a3 100644
--- a/examples/PinnedH2O_3DOF/rotation_test.py
+++ b/examples/PinnedH2O_3DOF/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,13 +47,32 @@ 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)
 bondlength1 = calculate_bondlength(H1_rotated, O1)
 bondlength2 = calculate_bondlength(H2_rotated, O1)
 bondangle = calculate_bondangle(H1_rotated, O1, H2_rotated, False)
+fliped_bond = False
 if bondlength1 < bondlength2:
+    fliped_bond = True
     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_rotated = np.dot(rot_matrix_align_plane.T, H1_rotated)
+H2_rotated = np.dot(rot_matrix_align_plane.T, H2_rotated)
+bondlength1 = calculate_bondlength(H1_rotated, O1)
+bondlength2 = calculate_bondlength(H2_rotated, O1)
+bondangle = calculate_bondangle(H1_rotated, O1, H2_rotated, False)
+
+print('Restored system')
 print(f'H1 = {H1_rotated}')
 print(f'H2 = {H2_rotated}')
 print(f'Bondlength of O1-H1 = {bondlength1}')

From 38a0511a8af15bc02d7efe80bb77847595c8fe11 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Mon, 18 Nov 2024 19:28:23 -0800
Subject: [PATCH 07/14] Minor change

---
 examples/PinnedH2O_3DOF/rotation_test.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/examples/PinnedH2O_3DOF/rotation_test.py b/examples/PinnedH2O_3DOF/rotation_test.py
index c92194a3..eaf22fe4 100644
--- a/examples/PinnedH2O_3DOF/rotation_test.py
+++ b/examples/PinnedH2O_3DOF/rotation_test.py
@@ -49,13 +49,13 @@ def rotation_matrix(axis, angle):
 
 H1_rotated = np.dot(rot_matrix_align_plane, H1)
 H2_rotated = np.dot(rot_matrix_align_plane, H2)
-bondlength1 = calculate_bondlength(H1_rotated, O1)
-bondlength2 = calculate_bondlength(H2_rotated, O1)
-bondangle = calculate_bondangle(H1_rotated, O1, H2_rotated, False)
 fliped_bond = False
 if bondlength1 < bondlength2:
     fliped_bond = True
-    H1_rotated, H2_rotated, bondlength1, bondlength2 = H2_rotated, H1_rotated, bondlength2, bondlength1
+    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)
 
 print('Reference system (z=0 plane about x=0 axis, with longer bondlength in H1)')
 print(f'H1 = {H1_rotated}')

From 2f79fdb31dcc54956aae32ad0c1c3efb1865c842 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Tue, 19 Nov 2024 10:20:45 -0800
Subject: [PATCH 08/14] Add C++ translation

---
 examples/PinnedH2O_3DOF/rotation_test.cc | 120 +++++++++++++++++++++++
 1 file changed, 120 insertions(+)
 create mode 100644 examples/PinnedH2O_3DOF/rotation_test.cc

diff --git a/examples/PinnedH2O_3DOF/rotation_test.cc b/examples/PinnedH2O_3DOF/rotation_test.cc
new file mode 100644
index 00000000..c342f632
--- /dev/null
+++ b/examples/PinnedH2O_3DOF/rotation_test.cc
@@ -0,0 +1,120 @@
+#include <iostream>
+#include <cmath>
+#include <algorithm>
+#include <cstring>
+
+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];
+}
+
+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;
+
+    return 0;
+}

From 8f615237973070d980b8df7ac2e514ec7ae6cc94 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Tue, 19 Nov 2024 10:30:32 -0800
Subject: [PATCH 09/14] Add restoring test

---
 examples/PinnedH2O_3DOF/rotation_test.cc | 27 +++++++++++++++++++++---
 examples/PinnedH2O_3DOF/rotation_test.py | 14 ++++++------
 2 files changed, 31 insertions(+), 10 deletions(-)

diff --git a/examples/PinnedH2O_3DOF/rotation_test.cc b/examples/PinnedH2O_3DOF/rotation_test.cc
index c342f632..0e5641ee 100644
--- a/examples/PinnedH2O_3DOF/rotation_test.cc
+++ b/examples/PinnedH2O_3DOF/rotation_test.cc
@@ -63,6 +63,12 @@ void apply_rotation(const double matrix[3][3], const double vec[3], double resul
     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};
@@ -99,22 +105,37 @@ int main() {
     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 << "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_3DOF/rotation_test.py b/examples/PinnedH2O_3DOF/rotation_test.py
index eaf22fe4..fbecb532 100644
--- a/examples/PinnedH2O_3DOF/rotation_test.py
+++ b/examples/PinnedH2O_3DOF/rotation_test.py
@@ -66,15 +66,15 @@ def rotation_matrix(axis, angle):
 
 if fliped_bond:
     H1_rotated, H2_rotated = H2_rotated, H1_rotated
-H1_rotated = np.dot(rot_matrix_align_plane.T, H1_rotated)
-H2_rotated = np.dot(rot_matrix_align_plane.T, H2_rotated)
-bondlength1 = calculate_bondlength(H1_rotated, O1)
-bondlength2 = calculate_bondlength(H2_rotated, O1)
-bondangle = calculate_bondangle(H1_rotated, O1, H2_rotated, False)
+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_rotated}')
-print(f'H2 = {H2_rotated}')
+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}')

From a09ca14943e0581c267c0fc462a9d1b847cdff8c Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Tue, 3 Dec 2024 09:20:43 -0800
Subject: [PATCH 10/14] Set up test

---
 tests/PinnedH2O_3DOF/coords.in |  3 ++
 tests/PinnedH2O_3DOF/mgmol.cfg | 31 ++++++++++++
 tests/PinnedH2O_3DOF/test.py   | 91 ++++++++++++++++++++++++++++++++++
 3 files changed, 125 insertions(+)
 create mode 100644 tests/PinnedH2O_3DOF/coords.in
 create mode 100644 tests/PinnedH2O_3DOF/mgmol.cfg
 create mode 100755 tests/PinnedH2O_3DOF/test.py

diff --git a/tests/PinnedH2O_3DOF/coords.in b/tests/PinnedH2O_3DOF/coords.in
new file mode 100644
index 00000000..46f5681d
--- /dev/null
+++ b/tests/PinnedH2O_3DOF/coords.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/mgmol.cfg b/tests/PinnedH2O_3DOF/mgmol.cfg
new file mode 100644
index 00000000..4a0b39ef
--- /dev/null
+++ b/tests/PinnedH2O_3DOF/mgmol.cfg
@@ -0,0 +1,31 @@
+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
diff --git a/tests/PinnedH2O_3DOF/test.py b/tests/PinnedH2O_3DOF/test.py
new file mode 100755
index 00000000..524027ba
--- /dev/null
+++ b/tests/PinnedH2O_3DOF/test.py
@@ -0,0 +1,91 @@
+#!/usr/bin/env python
+import sys
+import os
+import subprocess
+import string
+
+print("Test PinnedH2O_3DOF...")
+
+nargs=len(sys.argv)
+
+mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3]
+for i in range(4,nargs-4):
+  mpicmd = mpicmd + " "+sys.argv[i]
+print("MPI run command: {}".format(mpicmd))
+
+exe = sys.argv[nargs-4]
+inp = sys.argv[nargs-3]
+coords = sys.argv[nargs-2]
+print("coordinates file: %s"%coords)
+
+#create links to potentials files
+dst1 = 'pseudo.O_ONCV_PBE_SG15'
+dst2 = 'pseudo.H_ONCV_PBE_SG15'
+src1 = sys.argv[nargs-1] + '/' + dst1
+src2 = sys.argv[nargs-1] + '/' + dst2
+
+if not os.path.exists(dst1):
+  print("Create link to %s"%dst1)
+  os.symlink(src1, dst1)
+if not os.path.exists(dst2):
+  print("Create link to %s"%dst2)
+  os.symlink(src2, dst2)
+
+#run
+command = "{} {} -c {} -i {}".format(mpicmd,exe,inp,coords)
+print("Run command: {}".format(command))
+output = subprocess.check_output(command,shell=True)
+lines=output.split(b'\n')
+
+#analyse output
+energies=[]
+for line in lines:
+  if line.count(b'%%'):
+    print(line)
+    words=line.split()
+    words=words[5].split(b',')[0]
+    energy = words.decode()
+  if line.count(b'achieved'):
+    energies.append(energy)
+
+flag=0
+forces=[]
+for line in lines:
+  if flag>0:
+    print(line)
+    words=line.split(b'    ')
+    forces.append(words[1].decode())
+    forces.append(words[2].decode())
+    forces.append(words[3].decode())
+    flag=flag-1
+  if line.count(b'Forces:'):
+    flag=2
+
+
+print("Check energies...")
+print( energies )
+if len(energies)<2:
+  print("Expected two converged energies")
+  sys.exit(1)
+
+tol = 1.e-6
+diff=eval(energies[1])-eval(energies[0])
+print(diff)
+if abs(diff)>tol:
+  print("Energies differ: {} vs {} !!!".format(energies[0],energies[1]))
+  sys.exit(1)
+
+print("Check forces...")
+print(forces)
+flag=0
+for i in range(6):
+  diff=eval(forces[i+6])-eval(forces[i])
+  print(diff)
+  if abs(diff)>1.e-3:
+    print("Forces difference larger than tol")
+    flag=1
+if flag>0:
+  sys.exit(1)
+
+print("Test SUCCESSFUL!")
+sys.exit(0)

From 47e63350f01676f4992316b2b9c6d8bcf0786b52 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Tue, 3 Dec 2024 10:02:11 -0800
Subject: [PATCH 11/14] Update cfg

---
 tests/PinnedH2O_3DOF/mgmol.cfg | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/tests/PinnedH2O_3DOF/mgmol.cfg b/tests/PinnedH2O_3DOF/mgmol.cfg
index 4a0b39ef..fb7959db 100644
--- a/tests/PinnedH2O_3DOF/mgmol.cfg
+++ b/tests/PinnedH2O_3DOF/mgmol.cfg
@@ -29,3 +29,8 @@ initial_type=Random
 initial_width=2.
 [Restart]
 output_level=4
+[ROM.offline]
+basis_file=PinnedH2O_3DOF_orbitals_basis_5_5
+[ROM.basis]
+compare_md=false
+number_of_orbital_basis=36

From cdf38c8e6597231ee7538f84f4b87a544ff6f2d9 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Fri, 6 Dec 2024 07:54:23 -0800
Subject: [PATCH 12/14] Fix conversions

---
 src/MGmol.h |  6 ++++++
 src/rom.cc  | 59 ++++++++++++++++++++++++++++++-----------------------
 2 files changed, 40 insertions(+), 25 deletions(-)

diff --git a/src/MGmol.h b/src/MGmol.h
index f7525600..68ce2592 100644
--- a/src/MGmol.h
+++ b/src/MGmol.h
@@ -21,6 +21,10 @@
 #include "MGmol_prototypes.h"
 #include <iostream>
 
+#ifdef MGMOL_HAS_LIBROM
+#include "librom.h"
+#endif  // MGMOL_HAS_LIBROM
+
 // inline double one(const double r){ return 1.; }
 
 inline double linear(const double r) { return 1. - r; }
@@ -359,6 +363,8 @@ class MGmol : public MGmolInterface
     }
 
 #ifdef MGMOL_HAS_LIBROM
+    const CAROM::Matrix* orbitals_to_carom_matrix(const OrbitalsType& orbitals);
+    void carom_matrix_to_orbitals(const CAROM::Matrix* Psi, OrbitalsType& orbitals);
     int save_orbital_snapshot(std::string snapshot_dir, OrbitalsType& orbitals);
     void project_orbital(std::string snapshot_dir, int rdim, OrbitalsType& orbitals);
 #endif
diff --git a/src/rom.cc b/src/rom.cc
index 19ff66a3..878f72cf 100644
--- a/src/rom.cc
+++ b/src/rom.cc
@@ -20,7 +20,36 @@
 #include <fstream>
 #include <sys/stat.h>
 
-// Save the wavefunction snapshots
+template <class OrbitalsType>
+const CAROM::Matrix* MGmol<OrbitalsType>::orbitals_to_carom_matrix(const OrbitalsType& orbitals)
+{
+    const int dim = orbitals.getLocNumpt();
+    const int num_orbitals = orbitals.chromatic_number();
+
+    CAROM::Options svd_options(dim, num_orbitals, 1);
+    CAROM::BasisGenerator basis_generator(svd_options, false, "foo");
+
+    for (int i = 0; i < num_orbitals; ++i)
+        basis_generator.takeSample(orbitals.getPsi(i));
+
+    return basis_generator.getSnapshotMatrix();
+}
+
+template <class OrbitalsType>
+void MGmol<OrbitalsType>::carom_matrix_to_orbitals(const CAROM::Matrix* Psi, OrbitalsType& orbitals)
+{
+    Control& ct = *(Control::instance());
+    Mesh* mesh = Mesh::instance();
+    pb::GridFunc<ORBDTYPE> gf_psi(mesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]);
+    CAROM::Vector psi;
+    for (int i = 0; i < Psi->numColumns(); ++i)
+    {
+        Psi->getColumn(i, psi);
+        gf_psi.assign(psi.getData());
+        orbitals.setPsi(gf_psi, i);
+    }
+}
+
 template <class OrbitalsType>
 int MGmol<OrbitalsType>::save_orbital_snapshot(std::string file_path, OrbitalsType& orbitals)
 {
@@ -60,35 +89,15 @@ int MGmol<OrbitalsType>::save_orbital_snapshot(std::string file_path, OrbitalsTy
 template <class OrbitalsType>
 void MGmol<OrbitalsType>::project_orbital(std::string file_path, int rdim, OrbitalsType& orbitals)
 {
-    const int dim = orbitals.getLocNumpt();
-    const int totalSamples = orbitals.chromatic_number();
-
-    CAROM::Options svd_options(dim, totalSamples, 1);
-    CAROM::BasisGenerator basis_generator(svd_options, false, "foo");
-
-    for (int i = 0; i < totalSamples; ++i)
-        basis_generator.takeSample(orbitals.getPsi(i));
-    const CAROM::Matrix* orbital_snapshots = basis_generator.getSnapshotMatrix();
+    const CAROM::Matrix* Psi = orbitals_to_carom_matrix(orbitals);
 
     CAROM::BasisReader reader(file_path);
     CAROM::Matrix* orbital_basis = reader.getSpatialBasis(rdim);
 
-    CAROM::Matrix* proj_orbital_coeff = orbital_basis->transposeMult(orbital_snapshots);
-    CAROM::Matrix* proj_orbital_snapshots = orbital_basis->mult(proj_orbital_coeff);
+    CAROM::Matrix* Psi_reduced = orbital_basis->transposeMult(Psi);
+    CAROM::Matrix* Psi_projected = orbital_basis->mult(Psi_reduced);
     
-    Control& ct = *(Control::instance());
-    Mesh* mesh = Mesh::instance();
-    pb::GridFunc<ORBDTYPE> gf_psi(mesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]);
-    CAROM::Vector snapshot, proj_snapshot;
-    for (int i = 0; i < totalSamples; ++i)
-    {
-        orbital_snapshots->getColumn(i, snapshot);
-        proj_orbital_snapshots->getColumn(i, proj_snapshot);
-        gf_psi.assign(proj_snapshot.getData());
-        orbitals.setPsi(gf_psi, i);
-        snapshot -= proj_snapshot;
-        std::cout << "Error for orbital " << i << " = " << snapshot.norm() << std::endl;
-    }
+    carom_matrix_to_orbitals(Psi_projected, orbitals);
 }
 
 template class MGmol<LocGridOrbitals>;

From 15555f49a9530a79488efe6815a3cec41c1eb05b Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Fri, 6 Dec 2024 08:13:59 -0800
Subject: [PATCH 13/14] Add new script

---
 tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc | 215 +++++++++++++++++++++
 1 file changed, 215 insertions(+)
 create mode 100644 tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc

diff --git a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
new file mode 100644
index 00000000..f13c93d2
--- /dev/null
+++ b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
@@ -0,0 +1,215 @@
+// 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"
+#endif  // MGMOL_HAS_LIBROM
+
+#include <cassert>
+#include <iostream>
+#include <time.h>
+#include <vector>
+
+#include <boost/program_options.hpp>
+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;
+        if (ct.isLocMode())
+            mgmol = new MGmol<LocGridOrbitals>(global_comm, *MPIdata::sout,
+                input_filename, lrs_filename, constraints_filename);
+        else
+            mgmol = new MGmol<ExtendedGridOrbitals>(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<double> positions;
+        mgmol->getAtomicPositions(positions);
+        std::vector<short> anumbers;
+        mgmol->getAtomicNumbers(anumbers);
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Positions:" << std::endl;
+            std::vector<short>::iterator ita = anumbers.begin();
+            for (std::vector<double>::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++;
+            }
+        }
+
+        // compute energy and forces using all MPI tasks
+        // expect positions to be replicated on all MPI tasks
+        std::vector<double> forces;
+        double eks
+            = mgmol->evaluateEnergyAndForces(positions, anumbers, forces);
+        mgmol->dumpRestart();
+
+        // print out results
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Eks: " << eks << std::endl;
+            std::cout << "Forces:" << std::endl;
+            for (std::vector<double>::iterator it = forces.begin();
+                 it != forces.end(); it += 3)
+            {
+                for (int i = 0; i < 3; i++)
+                    std::cout << "    " << *(it + i);
+                std::cout << std::endl;
+            }
+        }
+
+        // compute energy and forces again using wavefunctions
+        // from previous call
+        Mesh* mymesh           = Mesh::instance();
+        const pb::Grid& mygrid = mymesh->grid();
+
+        std::shared_ptr<ProjectedMatricesInterface> projmatrices
+            = mgmol->getProjectedMatrices();
+
+        ExtendedGridOrbitals orbitals("new_orbitals", mygrid, mymesh->subdivx(),
+            ct.numst, ct.bcWF, projmatrices.get(), nullptr, nullptr, nullptr,
+            nullptr);
+
+        int rdim = 36;
+        CAROM::BasisReader reader(file_path);
+        CAROM::Matrix* orbital_basis = reader.getSpatialBasis(rdim);
+        carom_matrix_to_orbitals(Psi_projected, orbitals);
+
+        //
+        // evaluate energy and forces again
+        //
+
+        // convergence should be really quick since we start with an initial
+        // guess which is the solution
+        ct.max_electronic_steps = 300;
+        eks                     = mgmol->evaluateEnergyAndForces(
+            &orbitals, positions, anumbers, forces);
+
+        // print out results
+        if (MPIdata::onpe0)
+        {
+            std::cout << "Eks: " << eks << std::endl;
+            std::cout << "Forces:" << std::endl;
+            for (std::vector<double>::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;
+}

From f2c963c346ab4cb7bfdeab3e3b1ac633891a1d29 Mon Sep 17 00:00:00 2001
From: Siu Wun Cheung <cheung26@llnl.gov>
Date: Fri, 6 Dec 2024 08:48:57 -0800
Subject: [PATCH 14/14] Add CMake

---
 tests/CMakeLists.txt                       | 11 +++++++++++
 tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc | 20 +++++++++++++-------
 2 files changed, 24 insertions(+), 7 deletions(-)

diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 6de34971..6712c47d 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -254,6 +254,8 @@ add_executable(testWFEnergyAndForces
                ${CMAKE_SOURCE_DIR}/tests/WFEnergyAndForces/testWFEnergyAndForces.cc)
 add_executable(testDMandEnergyAndForces
                ${CMAKE_SOURCE_DIR}/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc)
+add_executable(testPinnedH2O_3DOF
+               ${CMAKE_SOURCE_DIR}/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc)
 
 if(${MAGMA_FOUND})
   add_executable(testOpenmpOffload
@@ -369,6 +371,14 @@ add_test(NAME testDMandEnergyAndForces
          ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/coords.in
          ${CMAKE_CURRENT_SOURCE_DIR}/DMandEnergyAndForces/lrs.in
          ${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
@@ -557,6 +567,7 @@ target_link_libraries(testDirectionalReduce PRIVATE MPI::MPI_CXX)
 target_link_libraries(testEnergyAndForces PRIVATE mgmol_src)
 target_link_libraries(testWFEnergyAndForces PRIVATE mgmol_src)
 target_link_libraries(testDMandEnergyAndForces PRIVATE mgmol_src)
+target_link_libraries(testPinnedH2O_3DOF PRIVATE mgmol_src)
 target_link_libraries(testIons PRIVATE mgmol_src)
 
 if(${MAGMA_FOUND})
diff --git a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
index f13c93d2..1cf93fdf 100644
--- a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
+++ b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc
@@ -166,18 +166,24 @@ int main(int argc, char** argv)
             ct.numst, ct.bcWF, projmatrices.get(), nullptr, nullptr, nullptr,
             nullptr);
 
-        int rdim = 36;
-        CAROM::BasisReader reader(file_path);
-        CAROM::Matrix* orbital_basis = reader.getSpatialBasis(rdim);
-        carom_matrix_to_orbitals(Psi_projected, orbitals);
+#ifdef MGMOL_HAS_LIBROM
+        CAROM::BasisReader reader(ct.getROMOptions().basis_file);
+        CAROM::Matrix* Psi = reader.getSpatialBasis(ct.getROMOptions().num_orbbasis);
+        //mgmol->carom_matrix_to_orbitals(Psi, orbitals);
+        pb::GridFunc<ORBDTYPE> gf_psi(mymesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]);
+        CAROM::Vector psi;
+        for (int i = 0; i < Psi->numColumns(); ++i)
+        {
+            Psi->getColumn(i, psi);
+            gf_psi.assign(psi.getData());
+            orbitals.setPsi(gf_psi, i);
+        }
+#endif  // MGMOL_HAS_LIBROM
 
         //
         // evaluate energy and forces again
         //
 
-        // convergence should be really quick since we start with an initial
-        // guess which is the solution
-        ct.max_electronic_steps = 300;
         eks                     = mgmol->evaluateEnergyAndForces(
             &orbitals, positions, anumbers, forces);