Orbital projection tests

examples/PinnedH2O/get_result.sh

@@ -0,0 +1,27 @@
+filename="offline_PinnedH2O.out" # FOM
+#filename="rom39_PinnedH2O.out" # compare MD
+#filename="39_force_PinnedH2O.out" # compare force
+
+# Extracting H1, H2, F1, F2 from MGmgol output log
+# if FOM, these files contain the FOM results
+# if compare MD, these files contain the results with projected orbitals
+awk '/H1 / {print $3, $4, $5}' $filename > H1_$filename
+awk '/H2 / {print $3, $4, $5}' $filename > H2_$filename
+awk '/F1 / {print $6, $7, $8}' $filename > F1_$filename
+awk '/F2 / {print $6, $7, $8}' $filename > F2_$filename
+
+# if compare force, files with "_fom" contain the FOM results
+# files with "_rom" contain the results with projected orbitals
+if [[ "$filename" == *"force_"* ]]; then
+    sed -n '1~2p' H1_$filename > H1_rom$filename
+    sed -n '1~2p' H2_$filename > H2_rom$filename
+    sed -n '1~2p' F1_$filename > F1_rom$filename
+    sed -n '1~2p' F2_$filename > F2_rom$filename
+
+    sed -n '2~2p' H1_$filename > H1_fom$filename
+    sed -n '2~2p' H2_$filename > H2_fom$filename
+    sed -n '2~2p' F1_$filename > F1_fom$filename
+    sed -n '2~2p' F2_$filename > F2_fom$filename
+fi
+
+rm -rf snapshot0_*

examples/PinnedH2O/job.basis → examples/PinnedH2O/job.basis_1_50

@@ -1,7 +1,7 @@
 #!/bin/tcsh
 #SBATCH -N 2
 #SBATCH -t 0:10:00
-#SBATCH -p pdebug
+#SBATCH -p pbatch
 
 date
 
@@ -10,23 +10,23 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 64
 
-set maindir = /p/lustre2/cheung26/mgmol-20240815
+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 increment_md_steps = 1
-set num_md_steps = 500
+set num_md_steps = 50
 
 foreach k (`seq $increment_md_steps $increment_md_steps $num_md_steps`)
     set snapshot_files = "$snapshot_files MD_mdstep${k}_snapshot"
 end
 echo "Snapshot files: $snapshot_files"
 
-set basis_file = "PinnedH2O_orbitals_basis"
+set basis_file = "PinnedH2O_orbitals_basis_${increment_md_steps}_${num_md_steps}"
 
-srun -n $ncpus $exe -f $basis_file $snapshot_files > basis_${increment_md_steps}_${num_md_steps}_Pinned_H2O.out
+srun -n $ncpus $exe -f $basis_file $snapshot_files > basis_${increment_md_steps}_${num_md_steps}_PinnedH2O.out
 
 date

examples/PinnedH2O/job.rom → examples/PinnedH2O/job.offline

@@ -10,7 +10,7 @@ setenv OMP_NUM_THREADS 1
 
 set ncpus = 64
 
-set maindir = /p/lustre2/cheung26/mgmol-20240815
+set maindir = /p/lustre2/cheung26/mgmol-20241012
 
 setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH
 
@@ -20,8 +20,8 @@ cp $maindir/install_quartz/bin/$exe .
 
 set datadir = $maindir/examples/PinnedH2O
 
-set cfg_rom = mgmol_rom.cfg
-#cp $datadir/$cfg_rom .
+set cfg_offline = mgmol_offline.cfg
+cp $datadir/$cfg_offline .
 
 cp $datadir/coords.in .
 
@@ -30,6 +30,6 @@ ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 .
 
 source $maindir/scripts/modules.quartz
 
-srun -n $ncpus $exe -c $cfg_rom -i coords.in > rom_PinnedH2O.out
+srun -n $ncpus $exe -c $cfg_offline -i coords.in > offline_PinnedH2O.out
 
 date

examples/PinnedH2O/job.rom_1_50

@@ -0,0 +1,39 @@
+#!/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
+
+set increment_md_steps = 1
+set num_md_steps = 50
+set basis_file = PinnedH2O_orbitals_basis_${increment_md_steps}_${num_md_steps}
+
+set cfg_rom = mgmol_rom_${increment_md_steps}_${num_md_steps}.cfg
+cp $datadir/$cfg_rom .
+
+cp $datadir/coords.in .
+
+ln -s -f $maindir/potentials/pseudo.O_ONCV_PBE_SG15 .
+ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 .
+
+source $maindir/scripts/modules.quartz
+
+srun -n $ncpus $exe -c $cfg_rom -i coords.in > rom_${increment_md_steps}_${num_md_steps}_PinnedH2O.out
+
+date

examples/PinnedH2O/mgmol_rom_1_50.cfg

@@ -0,0 +1,40 @@
+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=MD
+[MD]
+num_steps=500
+dt=40.
+thermostat=ON
+[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]
+basis_file=PinnedH2O_orbitals_basis_1_50
+[ROM.basis]
+compare_md=false
+number_of_orbital_basis=39

examples/PinnedH2O/plot_PinnedH2O_force.m

@@ -0,0 +1,84 @@
+clc; clear all; close all;
+
+%%
+plot_fom = 0;
+plot_rom = 0;
+rdim = 77;
+
+%%
+load F_fom.mat
+fprintf(1, 'Force statistics using FOM orbitals\n');
+fprintf(1, 'Mean of force on H1: %6.4e, %6.4e, %6.4e\n', mean(F1_fom));
+fprintf(1, 'Variance of force on H1: %6.4e, %6.4e, %6.4e\n', var(F1_fom));
+fprintf(1, 'Mean of force on H2: %6.4e, %6.4e, %6.4e\n', mean(F2_fom));
+fprintf(1, 'Variance of force on H2: %6.4e, %6.4e, %6.4e\n', var(F2_fom));
+
+if plot_fom
+    plotForce(F1_fom, 'F_H1_fom');
+    plotForce(F2_fom, 'F_H2_fom');
+    plotForceHistograms(F1_fom, 'H1_fom');
+    plotForceHistograms(F2_fom, 'H2_fom');
+end
+
+%%
+load(['F_rom' int2str(rdim) '.mat'])
+fprintf(1, 'Force statistics using projected orbitals\n');
+fprintf(1, 'Mean of force on H1: %6.4e, %6.4e, %6.4e\n', mean(F1_rom));
+fprintf(1, 'Variance of force on H1: %6.4e, %6.4e, %6.4e\n', var(F1_rom));
+%H1_correlation = sum(F1_fom(:) .* F1_rom(:)) / (norm(F1_fom(:)) * norm(F1_rom(:)))
+fprintf(1, 'Mean of force on H2: %6.4e, %6.4e, %6.4e\n', mean(F2_rom));
+fprintf(1, 'Variance of force on H2: %6.4e, %6.4e, %6.4e\n', var(F2_rom));
+%H2_correlation = sum(F2_fom(:) .* F2_rom(:)) / (norm(F2_fom(:)) * norm(F2_rom(:)))
+
+if plot_rom
+    plotForce(F1_rom, ['F_H1_rom' int2str(rdim)]);
+    plotForce(F2_rom, ['F_H2_rom' int2str(rdim)]);
+    plotForceHistograms(F1_rom, ['H1_rom' int2str(rdim)]);
+    plotForceHistograms(F2_rom, ['H2_rom' int2str(rdim)]);
+    plotForceHistogram(abs(F1_fom - F1_rom), ['H1_rom' int2str(rdim)], 'Fdiff');
+    plotForceHistogram(abs(F2_fom - F2_rom), ['H2_rom' int2str(rdim)], 'Fdiff');
+end
+
+%% 
+function plotForce(F, suffix)
+    figure;
+    imagesc(F');
+    axis tight;
+    axis equal;
+    colorbar;
+    saveas(gcf, suffix, 'jpg');
+end
+
+function plotForceHistogram(F, suffix, var)
+    figure;
+    if strcmp(var,'Fx')
+        X = F(:,1);
+        var_name = 'x-directional Force';
+    elseif strcmp(var,'Fy')
+        X = F(:,2);
+        var_name = 'y-directional Force';
+    elseif strcmp(var,'Fz')
+        X = F(:,3);
+        var_name = 'z-directional Force';
+    elseif strcmp(var,'Fmag')
+        X = sqrt(sum(F.^2, 2));
+        var_name = 'Force Magitude';
+    elseif strcmp(var,'Fdiff')
+        X = sqrt(sum(F.^2, 2));
+        var_name = 'Magitude of Difference in Force';
+    else
+        error('Invalid type');
+    end
+    histogram(X, 20);
+    xlabel(var_name);
+    ylabel('Frequency');
+    title(['Histogram of ' var_name]);
+    saveas(gcf, [var '_' suffix], 'jpg');
+end
+
+function plotForceHistograms(F, suffix)
+    plotForceHistogram(F, suffix, 'Fx');
+    plotForceHistogram(F, suffix, 'Fy');
+    plotForceHistogram(F, suffix, 'Fz');
+    plotForceHistogram(F, suffix, 'Fmag');
+end

examples/PinnedH2O/plot_PinnedH2O_md.m

@@ -0,0 +1,59 @@
+clc; clear all; close all;
+
+%%
+plot_fom = 0;
+rdims = [77, 39];
+
+%% 
+load md_fom.mat
+if plot_fom
+    plotAngleHistogram(H1_fom, H2_fom, 'fom');
+end
+
+%%
+
+all_H1_rom = zeros(length(rdims), size(H1_fom, 1), 3);
+all_H2_rom = zeros(length(rdims), size(H2_fom, 1), 3);
+k = 0;
+
+for rdim = rdims
+    k = k + 1;
+    load(['md_rom' int2str(rdim) '.mat'])
+    plotAngleHistogram(H1_rom, H2_rom, ['rom' int2str(rdim)]);
+    all_H1_rom(k, :, :) = H1_rom;
+    all_H2_rom(k, :, :) = H2_rom;
+end
+
+plotAtomTrajectory(H1_fom(:,1), all_H1_rom(:,:,1), rdims, 'x', 1)
+plotAtomTrajectory(H1_fom(:,2), all_H1_rom(:,:,2), rdims, 'y', 1)
+plotAtomTrajectory(H1_fom(:,3), all_H1_rom(:,:,3), rdims, 'z', 1)
+
+plotAtomTrajectory(H2_fom(:,1), all_H2_rom(:,:,1), rdims, 'x', 2)
+plotAtomTrajectory(H2_fom(:,2), all_H2_rom(:,:,2), rdims, 'y', 2)
+plotAtomTrajectory(H2_fom(:,3), all_H2_rom(:,:,3), rdims, 'z', 2)
+
+%% 
+function plotAtomTrajectory(X_fom, all_X_rom, rdims, var, idx)
+    figure;
+    hold on;
+    plot(X_fom, 'Linewidth', 2, 'DisplayName', 'FOM');
+    k = 0;
+    for rdim = rdims
+        k = k + 1;
+        X_rom = all_X_rom(k, :);
+        plot(X_rom, 'Linewidth', 2, 'DisplayName', ['ROM dim = ' int2str(rdim)]);
+    end
+    title([var '-coordinate of H' int2str(idx)])
+    legend;
+    saveas(gcf, [var '_H' int2str(idx)], 'jpg');    
+end
+
+function plotAngleHistogram(X1, X2, suffix)
+    figure;
+    A = acosd(sum(X1.*X2,2) ./ sqrt(sum(X1.^2,2)) ./ sqrt(sum(X2.^2,2)));
+    histogram(A, 20);
+    xlabel('Angle');
+    ylabel('Frequency');
+    title('Histogram of angle');
+    saveas(gcf, ['angle_' suffix], 'jpg');
+end

src/Control.cc

@@ -2057,6 +2057,8 @@ void Control::setROMOptions(const boost::program_options::variables_map& vm)
         rom_pri_option.save_librom_snapshot = vm["ROM.offline.save_librom_snapshot"].as<bool>();
         rom_pri_option.librom_snapshot_freq = vm["ROM.offline.librom_snapshot_freq"].as<int>();
 
+        rom_pri_option.compare_md = vm["ROM.basis.compare_md"].as<bool>();
+        rom_pri_option.num_orbbasis = vm["ROM.basis.number_of_orbital_basis"].as<int>();
         rom_pri_option.num_potbasis = vm["ROM.basis.number_of_potential_basis"].as<int>();
     }  // onpe0
 
@@ -2106,6 +2108,12 @@ void Control::syncROMOptions()
     rom_pri_option.rom_stage = static_cast<ROMStage>(rom_stage);
     rom_pri_option.variable = static_cast<ROMVariable>(rom_var);
 
+    mpirc = MPI_Bcast(&rom_pri_option.compare_md, 1, MPI_C_BOOL, 0, comm_global_);
+    bcast_check(mpirc);
+
+    mpirc = MPI_Bcast(&rom_pri_option.num_orbbasis, 1, MPI_INT, 0, comm_global_);
+    bcast_check(mpirc);
+
     mpirc = MPI_Bcast(&rom_pri_option.num_potbasis, 1, MPI_INT, 0, comm_global_);
     bcast_check(mpirc);
 }

src/MGmol.h

@@ -359,6 +359,7 @@ 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);
 #endif
 };
 // Instantiate static variables here to avoid clang warnings