diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 00000000..43953585 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,82 @@ +name: CI +on: + workflow_dispatch: {} + pull_request: + types: [opened, labeled, synchronize] + branches: + - 'ROMFPMD' + # push: + # branches: + # - release + +jobs: + docker-image: + uses: ./.github/workflows/docker_image.yml + build: + runs-on: ubuntu-latest + needs: [docker-image] + container: + image: ghcr.io/llnl/mgmol/mgmol_env:latest + options: --user 1001 --privileged + volumes: + - /mnt:/mnt + steps: + - name: Cancel previous runs + uses: styfle/cancel-workflow-action@0.11.0 + with: + access_token: ${{ github.token }} + # - name: Set Swap Space + # uses: pierotofy/set-swap-space@master + # with: + # swap-size-gb: 10 + - name: Check out mgmol + uses: actions/checkout@v1 + with: + submodules: 'true' + - name: cmake + run: | + mkdir ${GITHUB_WORKSPACE}/build + cd ${GITHUB_WORKSPACE}/build + cmake .. -DCMAKE_CXX_COMPILER=mpic++ -DCMAKE_Fortran_COMPILER=mpif90 -DMPIEXEC_PREFLAGS="--oversubscribe" -DUSE_LIBROM=On -DLIBROM_PATH=/env/dependencies/libROM + - name: make + run: | + cd ${GITHUB_WORKSPACE}/build && make -j 4 + - name: test + run: | + cd ${GITHUB_WORKSPACE}/build && ctest --no-compress-output -V -T Test -I 1,20,1 + - name: test ROM Poisson operator + run: | + cd ${GITHUB_WORKSPACE}/tests/ROM/test_rom_poisson + ln -s ${GITHUB_WORKSPACE}/build/src/mgmol-rom . + ln -s ${GITHUB_WORKSPACE}/potentials/* . + mpirun -n 3 --oversubscribe ./mgmol-rom -c carbyne.poisson.cfg -i carbyne.in + - name: test ROM ion density evaluation + run: | + cd ${GITHUB_WORKSPACE}/tests/ROM/test_rom_poisson + mpirun -n 3 --oversubscribe ./mgmol-rom -c carbyne.ion.cfg -i carbyne.in + # code-style: + # runs-on: ubuntu-latest + # needs: [docker-image] + # container: + # image: ghcr.io/llnl/mgmol/mgmol_env:latest + # options: --user 1001 --privileged + # volumes: + # - /mnt:/mnt + # steps: + # - name: Cancel previous runs + # uses: styfle/cancel-workflow-action@0.11.0 + # with: + # access_token: ${{ github.token }} + # - name: Check out mgmol + # uses: actions/checkout@v1 + # with: + # submodules: 'true' + # - name: cmake + # run: | + # mkdir build + # cd build + # cmake .. -DCMAKE_CXX_COMPILER=mpic++ -DCMAKE_Fortran_COMPILER=mpif90 -DMGMOL_WITH_CLANG_FORMAT=ON + # - name: make + # run: | + # cd build && make format + diff --git a/.github/workflows/docker_image.yml b/.github/workflows/docker_image.yml new file mode 100644 index 00000000..90a2447b --- /dev/null +++ b/.github/workflows/docker_image.yml @@ -0,0 +1,59 @@ +name: docker-image +on: + workflow_call: + +env: + REGISTRY: ghcr.io + # github.repository as / + IMAGE_NAME: llnl/mgmol/mgmol_env + DOCKERPATH: docker + +jobs: + docker-ci: + runs-on: ubuntu-latest + name: "docker env" + steps: + - uses: actions/checkout@v3 + with: + fetch-depth: 0 + - uses: Ana06/get-changed-files@v2.2.0 + id: files + - name: DockerPATH configuration + run: echo "DOCKERPATH=$DOCKERPATH" + - name: DockerPATH - check if files in docker path changed + if: contains(steps.files.outputs.all,env.DOCKERPATH) || contains(steps.files.outputs.all,'docker_image.yml') + run: | + echo "CI container needs rebuilding..." + echo "CI_NEEDS_REBUILD=true" >> $GITHUB_ENV + - name: Log into registry ${{ env.REGISTRY }} + if: env.CI_NEEDS_REBUILD + uses: docker/login-action@v2 + with: + registry: ${{ env.REGISTRY }} + username: ${{ github.actor }} + password: ${{ secrets.GITHUB_TOKEN }} + - name: Extract metadata (tags, labels) for Docker + id: meta + if: env.CI_NEEDS_REBUILD + uses: docker/metadata-action@v4 + with: + images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }} + tags: type=sha + flavor: latest=true + - name: Build Container motd + if: env.CI_NEEDS_REBUILD + run: | + echo "#!/bin/bash" > ${{env.DOCKERPATH}}/motd.sh + echo "echo --------------------------" >> ${{env.DOCKERPATH}}/motd.sh + echo "echo mgmol_env/CI Development Container" >> ${{env.DOCKERPATH}}/motd.sh + echo "echo \"Revision: `echo ${GITHUB_SHA} | cut -c1-8`\"" >> ${{env.DOCKERPATH}}/motd.sh + echo "echo --------------------------" >> ${{env.DOCKERPATH}}/motd.sh + chmod 755 ${{env.DOCKERPATH}}/motd.sh + cat ${{env.DOCKERPATH}}/motd.sh + - name: Docker Image - Build and push + if: env.CI_NEEDS_REBUILD + uses: docker/build-push-action@v3 + with: + push: true + context: ${{ env.DOCKERPATH }} + tags: ${{ steps.meta.outputs.tags }} diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..81c6b154 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +src/mgmol_config.h diff --git a/CMakeLists.txt b/CMakeLists.txt index c1a853f0..23d058f0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -139,6 +139,22 @@ if (${MGMOL_WITH_SCALAPACK} OR DEFINED SCALAPACK_ROOT) endif(${SCALAPACK_FOUND}) endif(${MGMOL_WITH_SCALAPACK} OR DEFINED SCALAPACK_ROOT) +# libROM (optional) +set(USE_LIBROM False CACHE BOOL "Build with libROM") +set(LIBROM_PATH "" CACHE STRING "Path of libROM") +if(USE_LIBROM) + message(STATUS "LIBROM_PATH: ${LIBROM_PATH}") + if(NOT LIBROM_PATH) + message(FATAL_ERROR "Cmake is asked to use libROM, but LIBROM_PATH not specified.") + endif(NOT LIBROM_PATH) + + find_package(libROM REQUIRED) + + if(libROM_FOUND) + add_definitions(-DMGMOL_HAS_LIBROM) + endif(libROM_FOUND) +endif(USE_LIBROM) + # ARPACK (optional) set(MGMOL_WITH_ARPACK FALSE CACHE BOOL "Compile with ARPACK package") if(${MGMOL_WITH_ARPACK} OR DEFINED ARPACK_ROOT) @@ -263,6 +279,9 @@ include_directories("${PROJECT_SOURCE_DIR}/src/sparse_linear_algebra") include_directories("${PROJECT_SOURCE_DIR}/src/tools") include_directories("${PROJECT_SOURCE_DIR}/src") +include_directories("${LIBROM_PATH}/lib") +link_libraries(${LIBROM_LIB}) + # add subdirectories for source files, tests add_subdirectory(src) diff --git a/cmake_modules/FindlibROM.cmake b/cmake_modules/FindlibROM.cmake new file mode 100644 index 00000000..9fed8fb6 --- /dev/null +++ b/cmake_modules/FindlibROM.cmake @@ -0,0 +1,11 @@ +if(NOT LIBROM_PATH) + message(FATAL_ERROR "LIBROM_PATH not specified.") +endif(NOT LIBROM_PATH) + +find_library(LIBROM_LIB libROM.so HINTS "${LIBROM_PATH}/build/lib") +find_path(LIBROM_INCLUDES librom.h HINTS "${LIBROM_PATH}/lib") + +mark_as_advanced(LIBROM_LIB LIBROM_INCLUDES) + +include(FindPackageHandleStandardArgs) +FIND_PACKAGE_HANDLE_STANDARD_ARGS(libROM REQUIRED_VARS LIBROM_LIB LIBROM_INCLUDES) \ No newline at end of file diff --git a/cmake_toolchains/quartz.default.cmake b/cmake_toolchains/quartz.default.cmake new file mode 100644 index 00000000..9901dcd6 --- /dev/null +++ b/cmake_toolchains/quartz.default.cmake @@ -0,0 +1,6 @@ +set(CMAKE_C_COMPILER mpicc) +set(CMAKE_CXX_COMPILER mpicxx) +set(CMAKE_Fortran_COMPILER mpif90) + +set(SCALAPACK_ROOT $ENV{MKLROOT}) +set(SCALAPACK_BLACS_LIBRARY $ENV{MKLROOT}/lib/intel64/libmkl_blacs_intelmpi_lp64.so) diff --git a/docker/Dockerfile b/docker/Dockerfile new file mode 100644 index 00000000..8b1ffa94 --- /dev/null +++ b/docker/Dockerfile @@ -0,0 +1,50 @@ +FROM ubuntu:22.04 + +ENV ENVDIR=env + +# install sudo +RUN apt-get -yq update && apt-get -yq install sudo + +WORKDIR /$ENVDIR + +# install packages +RUN sudo apt-get install -yq git +RUN sudo apt-get install --no-install-recommends -yq make gcc gfortran libssl-dev cmake +RUN sudo apt-get install -yq libopenblas-dev libmpich-dev libblas-dev liblapack-dev libscalapack-mpi-dev libhdf5-mpi-dev +RUN sudo apt-get install -yq libboost-all-dev +RUN sudo apt-get install -yq vim +RUN sudo apt-get install -yq git-lfs +RUN sudo apt-get install -yq valgrind hdf5-tools +RUN sudo apt-get install -yq wget +### clang-format seems to be updated to 14.0. Not using it for now. +# RUN sudo apt-get install -yq clang-format + +# install lldb and gdb for debugging +RUN sudo apt-get install -yq lldb gdb + +RUN sudo apt-get clean -q + +ENV LIB_DIR=/$ENVDIR/dependencies +WORKDIR $LIB_DIR + +# cmake toolchain file for librom +RUN echo 'set(CMAKE_C_COMPILER mpicc)\n\ +set(CMAKE_CXX_COMPILER mpicxx)\n\ +set(CMAKE_Fortran_COMPILER mpif90)' > ./librom_env.cmake +ENV TOOLCHAIN_FILE=$LIB_DIR/librom_env.cmake + +# install libROM for scaleupROM +RUN sudo git clone https://github.com/LLNL/libROM.git +WORKDIR ./libROM/build +# libROM without MFEM. +RUN sudo cmake .. -DCMAKE_TOOLCHAIN_FILE=${TOOLCHAIN_FILE} -DCMAKE_BUILD_TYPE=Optimized -DUSE_MFEM=OFF +RUN sudo make -j 16 + +# create and switch to a user +ENV USERNAME=test +RUN echo "$USERNAME ALL=(ALL) NOPASSWD:ALL" >> /etc/sudoers +RUN useradd --no-log-init -u 1001 --create-home --shell /bin/bash $USERNAME +RUN adduser $USERNAME sudo +USER $USERNAME +WORKDIR /home/$USERNAME + diff --git a/examples/Carbyne/carbyne.rom.cfg b/examples/Carbyne/carbyne.rom.cfg new file mode 100644 index 00000000..cb0cd295 --- /dev/null +++ b/examples/Carbyne/carbyne.rom.cfg @@ -0,0 +1,34 @@ +verbosity=2 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx= 96 +ny= 96 +nz= 192 +[Domain] +ox= -10. +oy= -10. +oz= -20. +lx= 20. +ly= 20. +lz= 40. +[Potentials] +pseudopotential=pseudo.H_ONCV_PBE_SG15 +pseudopotential=pseudo.C_ONCV_PBE_SG15 +[Run] +type=QUENCH +[Quench] +max_steps=5 +atol=1.e-8 +[Orbitals] +initial_type=Fourier +[Restart] +output_level=4 +input_level=4 +input_filename=snapshot0_000 + +[ROM.offline] +restart_filefmt=snapshot0_%03d +restart_min_idx=0 +restart_max_idx=1 +basis_file=carom diff --git a/examples/PinnedH2O/coords.in b/examples/PinnedH2O/coords_test1.in similarity index 100% rename from examples/PinnedH2O/coords.in rename to examples/PinnedH2O/coords_test1.in diff --git a/examples/PinnedH2O/coords_test2.in b/examples/PinnedH2O/coords_test2.in new file mode 100644 index 00000000..e1e31f8b --- /dev/null +++ b/examples/PinnedH2O/coords_test2.in @@ -0,0 +1,3 @@ +O1 1 0.00 0.00 0.00 0 +H1 2 -0.45 1.57 -1.07 1 +H2 2 -0.45 -1.48 -0.97 1 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/get_result.sh b/examples/PinnedH2O/get_result.sh new file mode 100644 index 00000000..b737960d --- /dev/null +++ b/examples/PinnedH2O/get_result.sh @@ -0,0 +1,37 @@ +#filename="offline_PinnedH2O" # FOM +#filename="rom39_PinnedH2O" # ROM compare MD +#filename="39_force_PinnedH2O" # ROM compare force + +#filename="PinnedH2O_test2_ref" # FOM +filename="PinnedH2O_rom_3DOF_test2_2_2_34" # ROM PinnedH2O 3DOF MD + +# Extracting kinetic energy, total energy, temperature from MGmgol output log +awk '/Kinetic/ {print $3}' $filename.out > ke_$filename.txt +awk '/Kinetic/ {print $5}' $filename.out >temp_$filename.txt +awk '/Total/ {print $3}' $filename.out > te_$filename.txt + +# 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 '/O1 / {print $4, $5, $6}' $filename.out > O1_$filename.txt +awk '/H1 / {print $3, $4, $5}' $filename.out > H1_$filename.txt +awk '/H2 / {print $3, $4, $5}' $filename.out > H2_$filename.txt +awk '/O1 / {print $7, $8, $9}' $filename.out > f_O1_$filename.txt +awk '/H1 / {print $6, $7, $8}' $filename.out > f_H1_$filename.txt +awk '/H2 / {print $6, $7, $8}' $filename.out > f_H2_$filename.txt + +# 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.out > H1_rom$filename.txt + sed -n '1~2p' H2_$filename.out > H2_rom$filename.txt + sed -n '1~2p' f_H1_$filename.out > f_H1_rom$filename.txt + sed -n '1~2p' f_H2_$filename.out > f_H2_rom$filename.txt + + sed -n '2~2p' H1_$filename.out > H1_fom$filename.txt + sed -n '2~2p' H2_$filename.out > H2_fom$filename.txt + sed -n '2~2p' f_H1_$filename.out > f_H1_fom$filename.txt + sed -n '2~2p' f_H2_$filename.out > f_H2_fom$filename.txt +fi + +rm -rf snapshot_* diff --git a/examples/PinnedH2O/job.basis_1_50 b/examples/PinnedH2O/job.basis_1_50 new file mode 100644 index 00000000..5291a7b7 --- /dev/null +++ b/examples/PinnedH2O/job.basis_1_50 @@ -0,0 +1,32 @@ +#!/bin/tcsh +#SBATCH -N 1 +#SBATCH -t 0:10:00 +#SBATCH -p pdebug + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 8 + +set maindir = /p/lustre2/cheung26/mgmol + +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 = 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_${increment_md_steps}_${num_md_steps}" + +srun -n $ncpus $exe -f $basis_file $snapshot_files > basis_${increment_md_steps}_${num_md_steps}_PinnedH2O.out + +date diff --git a/examples/PinnedH2O/job.offline b/examples/PinnedH2O/job.offline new file mode 100644 index 00000000..9d7150f8 --- /dev/null +++ b/examples/PinnedH2O/job.offline @@ -0,0 +1,30 @@ +#!/bin/tcsh +#SBATCH -N 1 +#SBATCH -t 1:00:00 +#SBATCH -p pdebug + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 8 + +set maindir = /p/lustre2/cheung26/mgmol + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = mgmol-opt + +cp $maindir/install_quartz/bin/$exe . + +set cfg_offline = mgmol_offline.cfg + +ln -s -f $maindir/potentials/pseudo.O_ONCV_PBE_SG15 . +ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 . + +source $maindir/scripts/modules.quartz + +srun -n $ncpus $exe -c $cfg_offline -i coords.in > offline_PinnedH2O.out + +date diff --git a/examples/PinnedH2O/job.ref b/examples/PinnedH2O/job.ref new file mode 100644 index 00000000..f9826e8e --- /dev/null +++ b/examples/PinnedH2O/job.ref @@ -0,0 +1,36 @@ +#!/bin/tcsh +#SBATCH -N 1 +#SBATCH -t 1:00:00 +#SBATCH -p pdebug + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 1 +set case = 2 + +set maindir = /p/lustre2/cheung26/mgmol + +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 cfg = mgmol_ref_test${case}.cfg +cp $datadir/$cfg . + +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 -i coords_test${case}.in > PinnedH2O_test${case}_ref.out + +date diff --git a/examples/PinnedH2O/job.rom b/examples/PinnedH2O/job.rom new file mode 100644 index 00000000..45285136 --- /dev/null +++ b/examples/PinnedH2O/job.rom @@ -0,0 +1,34 @@ +#!/bin/tcsh +#SBATCH -N 1 +#SBATCH -t 1:00:00 +#SBATCH -p pdebug + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 8 + +set maindir = /p/lustre2/cheung26/mgmol + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = mgmol-opt + +cp $maindir/install_quartz/bin/$exe . + +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 + +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 diff --git a/examples/PinnedH2O/job.rom_3DOF b/examples/PinnedH2O/job.rom_3DOF new file mode 100644 index 00000000..b7dc5d6c --- /dev/null +++ b/examples/PinnedH2O/job.rom_3DOF @@ -0,0 +1,35 @@ +#!/bin/tcsh +#SBATCH -N 1 +#SBATCH -t 1:00:00 +#SBATCH -p pdebug + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 1 +set case = 2 + +set maindir = /p/lustre2/cheung26/mgmol + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = mgmol-opt + +cp $maindir/install_quartz/bin/$exe . + +set bondlength_num_increments = 2 +set bondangle_num_increments = 2 +set num_basis = 34 + +set cfg_rom = mgmol_rom_3DOF_test${case}.cfg + +ln -s -f $maindir/potentials/pseudo.O_ONCV_PBE_SG15 . +ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 . + +source $maindir/scripts/modules.quartz + +srun -n $ncpus $exe -c $cfg_rom -i coords_test${case}.in > PinnedH2O_rom_3DOF_test${case}_${bondlength_num_increments}_${bondangle_num_increments}_${num_basis}.out + +date diff --git a/examples/PinnedH2O/mgmol_offline.cfg b/examples/PinnedH2O/mgmol_offline.cfg new file mode 100644 index 00000000..f535cff2 --- /dev/null +++ b/examples/PinnedH2O/mgmol_offline.cfg @@ -0,0 +1,38 @@ +verbosity=1 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=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] +save_librom_snapshot=true +librom_snapshot_freq=1 diff --git a/examples/PinnedH2O/mgmol.cfg b/examples/PinnedH2O/mgmol_ref_test1.cfg similarity index 85% rename from examples/PinnedH2O/mgmol.cfg rename to examples/PinnedH2O/mgmol_ref_test1.cfg index 32325e47..fd456f92 100644 --- a/examples/PinnedH2O/mgmol.cfg +++ b/examples/PinnedH2O/mgmol_ref_test1.cfg @@ -1,6 +1,6 @@ verbosity=1 xcFunctional=PBE -FDtype=Mehrstellen +FDtype=4th [Mesh] nx=64 ny=64 @@ -18,12 +18,12 @@ pseudopotential=pseudo.H_ONCV_PBE_SG15 [Run] type=MD [MD] -num_steps=50 +num_steps=500 dt=40. thermostat=ON [Thermostat] type=Berendsen -temperature=1000. +temperature=300. relax_time=800. [Quench] max_steps=100 @@ -32,4 +32,4 @@ atol=1.e-8 initial_type=Random initial_width=2. [Restart] -output_level=0 +output_level=4 diff --git a/examples/PinnedH2O/mgmol_ref_test2.cfg b/examples/PinnedH2O/mgmol_ref_test2.cfg new file mode 100644 index 00000000..67aa443c --- /dev/null +++ b/examples/PinnedH2O/mgmol_ref_test2.cfg @@ -0,0 +1,31 @@ +verbosity=1 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=MD +[MD] +num_steps=500 +dt=40. +thermostat=OFF +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +[Restart] +output_level=4 diff --git a/examples/PinnedH2O/mgmol_rom_1_50.cfg b/examples/PinnedH2O/mgmol_rom_1_50.cfg new file mode 100644 index 00000000..13cc3921 --- /dev/null +++ b/examples/PinnedH2O/mgmol_rom_1_50.cfg @@ -0,0 +1,42 @@ +verbosity=1 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=MD +[MD] +num_steps=500 +dt=40. +thermostat=ON +[Thermostat] +type=Berendsen +temperature=300. +relax_time=800. +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +[Restart] +output_level=4 +[ROM] +stage=test_orbital +[ROM.offline] +basis_file=PinnedH2O_orbitals_basis_1_50 +[ROM.basis] +compare_md=false +number_of_orbital_basis=39 diff --git a/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg b/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg new file mode 100644 index 00000000..b7a78a0d --- /dev/null +++ b/examples/PinnedH2O/mgmol_rom_3DOF_test1.cfg @@ -0,0 +1,48 @@ +verbosity=2 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=MD +[MD] +num_steps=500 +dt=40. +thermostat=ON +[Thermostat] +type=Berendsen +temperature=300. +relax_time=800. +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +nempty=30 +[DensityMatrix] +solver=MVP +nb_inner_it=20 +mixing=0.8 +tol=1.e-8 +[Restart] +output_level=4 +[ROM] +stage=online_pinned_H2O_3dof +[ROM.offline] +basis_file=/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_8/PinnedH2O_3DOF_orbitals_basis_2_2 +[ROM.basis] +compare_md=false +number_of_orbital_basis=34 diff --git a/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg b/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg new file mode 100644 index 00000000..bb2912d9 --- /dev/null +++ b/examples/PinnedH2O/mgmol_rom_3DOF_test2.cfg @@ -0,0 +1,44 @@ +verbosity=2 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=MD +[MD] +num_steps=500 +dt=40. +thermostat=OFF +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +nempty=30 +[DensityMatrix] +solver=MVP +nb_inner_it=20 +mixing=0.8 +tol=1.e-8 +[Restart] +output_level=4 +[ROM] +stage=online_pinned_H2O_3dof +[ROM.offline] +basis_file=/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_1/PinnedH2O_3DOF_orbitals_basis_2_2 +[ROM.basis] +compare_md=false +number_of_orbital_basis=34 diff --git a/examples/PinnedH2O/plot_PinnedH2O_force.m b/examples/PinnedH2O/plot_PinnedH2O_force.m new file mode 100644 index 00000000..a41457a5 --- /dev/null +++ b/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 diff --git a/examples/PinnedH2O/plot_PinnedH2O_md.m b/examples/PinnedH2O/plot_PinnedH2O_md.m new file mode 100644 index 00000000..47c309b2 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/examples/PinnedH2O/postprocess_PinnedH2O.py b/examples/PinnedH2O/postprocess_PinnedH2O.py new file mode 100644 index 00000000..142a34ac --- /dev/null +++ b/examples/PinnedH2O/postprocess_PinnedH2O.py @@ -0,0 +1,27 @@ +import subprocess +import re + +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") + +pattern = r"For energy fraction: \d+\.\d+, take first (\d+) of \d+ basis vectors" +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/scripts/build_quartz_libROM.sh b/scripts/build_quartz_libROM.sh new file mode 100644 index 00000000..2b835a73 --- /dev/null +++ b/scripts/build_quartz_libROM.sh @@ -0,0 +1,62 @@ +#!/bin/bash +## An example script to build on LLNL Peloton systems. +## For now, this script assumes intel/ mkl libraries are being used. + +## load some modules +source scripts/modules.quartz + +## set some environment variables. Set them explicitly or use loaded module path (preferred) +## Here we use an explicit path for scalapack to be consistent with the path for the blas libraries and avoid +## benign cmake warnings +##setenv SCALAPACK_ROOT /usr/tce/packages/mkl/mkl-2020.0/lib +#setenv SCALAPACK_ROOT ${MKLROOT} +#setenv HDF5_ROOT /usr/tce/packages/hdf5/hdf5-1.14.0-mvapich2-2.3.6-intel-2022.1.0 +# +## We need to define the cmake blas vendor option here to find the right one. +#set BLAS_VENDOR = Intel10_64lp +# +## manually set the location of BLACS libraries for scalapack +#set BLACS_LIB = ${SCALAPACK_ROOT}/lib/intel64 + +MGMOL_ROOT="$(pwd)" + +INSTALL_DIR=${MGMOL_ROOT}/install_quartz +mkdir -p ${INSTALL_DIR} + +BUILD_DIR=${MGMOL_ROOT}/build_quartz +mkdir -p ${BUILD_DIR} +cd ${BUILD_DIR} + +# clone the libROM GitHub repo in BUILD_DIR +USE_LIBROM="On" +LIBROM_PATH=${BUILD_DIR}/libROM +git clone https://github.com/LLNL/libROM +cd libROM +git checkout 321d18f4d5adfa29f0a3de9be2699fee9732f2bf +#./scripts/compile.sh -t ./cmake/toolchains/default-toss_4_x86_64_ib-librom-dev.cmake +./scripts/compile.sh +cd ${BUILD_DIR} + +# call cmake +cmake -DCMAKE_TOOLCHAIN_FILE=${MGMOL_ROOT}/cmake_toolchains/quartz.default.cmake \ + -DCMAKE_INSTALL_PREFIX=${INSTALL_DIR} \ + -DUSE_LIBROM=${USE_LIBROM} \ + -DLIBROM_PATH=${LIBROM_PATH} \ + .. + +# -DCMAKE_CXX_COMPILER=mpic++ \ +# -DCMAKE_Fortran_COMPILER=mpif77 \ +# -DMPIEXEC_NUMPROC_FLAG="-n" \ +# -DBLA_VENDOR=${BLAS_VENDOR} \ +# -DSCALAPACK_BLACS_LIBRARY=${BLACS_LIB}/libmkl_blacs_intelmpi_lp64.so \ +# -DCMAKE_BUILD_TYPE=DEBUG \ + +# call make install +make -j 16 +### Currently libROM does not have the installation procedure, +### so copying binary file to installation directory will disrupt the relative path to libROM.so, +### causing a run-time error. +make install + +# -DBLAS_LIBRARIES=/usr/tce/packages/mkl/mkl-2022.1.0/mkl/2022.1.0/lib/intel64/lib \ +# -DLAPACK_LIBRARIES=/usr/tce/packages/mkl/mkl-2022.1.0/mkl/2022.1.0/lib/intel64/lib \ diff --git a/scripts/modules.quartz b/scripts/modules.quartz new file mode 100644 index 00000000..82433acb --- /dev/null +++ b/scripts/modules.quartz @@ -0,0 +1,21 @@ +### choose either gcc or intel +#module load intel/2022.1.0 +module load gcc/11.2.1 + +module load cmake +module load hdf5-parallel/1.14.0 +module load boost + +### choose either one +module load mkl-interfaces +#module load mkl + +module load python + +### manually add boost path +#setenv LD_LIBRARY_PATH /usr/tce/packages/boost/boost-1.80.0-mvapich2-2.3.6-gcc-10.3.1/lib:$LD_LIBRARY_PATH + +#setenv MKLROOT $LIBRARY_PATH +#setenv MKLROOT /usr/tce/packages/mkl/mkl-2022.1.0/mkl/2022.1.0 +#setenv HDF5ROOT $LD_LIBRARY_PATH +#setenv HDF5ROOT diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 9a1a02bd..8a2a2021 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,7 @@ +configure_file(${CMAKE_CURRENT_SOURCE_DIR}/mgmol_config.h.in + ${CMAKE_CURRENT_SOURCE_DIR}/mgmol_config.h @ONLY) + add_subdirectory(DistMatrix) add_subdirectory(linear_algebra) add_subdirectory(local_matrices) @@ -8,7 +11,8 @@ add_subdirectory(radial) add_subdirectory(sparse_linear_algebra) add_subdirectory(tools) -set(link_libs mgmol_distmatrix +set(link_libs + mgmol_distmatrix mgmol_linear_algebra mgmol_local_matrices mgmol_numerical_kernels @@ -161,13 +165,23 @@ set(SOURCES magma_singleton.cc ChebyshevApproximation.cc ChebyshevApproximationInterface.cc + rom.cc ) +if(USE_LIBROM) + list(APPEND SOURCES + rom_workflows.cc) +endif(USE_LIBROM) + add_library(mgmol_src ${SOURCES}) target_include_directories(mgmol_src PRIVATE ${HDF5_INCLUDE_DIRS}) target_include_directories(mgmol_src PRIVATE ${Boost_INCLUDE_DIRS}) -target_include_directories (mgmol_src PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) +target_include_directories(mgmol_src PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) +if(USE_LIBROM) + target_include_directories (mgmol_src PRIVATE ${LIBROM_INCLUDES}) + target_link_libraries(mgmol_src ${LIBROM_LIB}) +endif(USE_LIBROM) target_link_libraries(mgmol_src ${link_libs}) if(${MGMOL_WITH_MAGMA}) @@ -200,3 +214,25 @@ if(${MGMOL_WITH_LIBXC}) endif (${MGMOL_WITH_LIBXC}) install(TARGETS mgmol-opt DESTINATION bin) + +# build ROM executable +if(USE_LIBROM) + add_executable(mgmol-rom rom_main.cc) + target_include_directories (mgmol-rom PRIVATE ${Boost_INCLUDE_DIRS}) + + target_link_libraries(mgmol-rom mgmol_src ${link_libs}) + target_link_libraries(mgmol-rom ${SCALAPACK_LIBRARIES}) + target_link_libraries(mgmol-rom ${HDF5_LIBRARIES}) + target_link_libraries(mgmol-rom ${HDF5_HL_LIBRARIES}) + target_link_libraries(mgmol-rom ${BLAS_LIBRARIES}) + target_link_libraries(mgmol-rom ${LAPACK_LIBRARIES}) + target_link_libraries(mgmol-rom ${Boost_LIBRARIES}) + target_link_libraries(mgmol-rom ${LIBROM_LIB}) + if (${OPENMP_CXX_FOUND}) + target_link_libraries(mgmol-rom OpenMP::OpenMP_CXX) + endif() + if(${MGMOL_WITH_LIBXC}) + target_link_libraries(mgmol-rom ${LIBXC_DIR}/lib/libxc.a) + endif (${MGMOL_WITH_LIBXC}) + install(TARGETS mgmol-rom DESTINATION bin) +endif(USE_LIBROM) diff --git a/src/Control.cc b/src/Control.cc index 428e7a4f..94e2ff30 100644 --- a/src/Control.cc +++ b/src/Control.cc @@ -1808,6 +1808,10 @@ void Control::setOptions(const boost::program_options::variables_map& vm) // synchronize all processors sync(); + +#ifdef MGMOL_HAS_LIBROM + setROMOptions(vm); +#endif } int Control::checkOptions() @@ -1950,3 +1954,110 @@ void Control::printPoissonOptions(std::ostream& os) } os << std::endl; } + +void Control::setROMOptions(const boost::program_options::variables_map& vm) +{ + printWithTimeStamp("Control::setROMOptions()...", std::cout); + + if (onpe0) + { + std::string str = vm["ROM.stage"].as(); + if (str.compare("offline") == 0) + rom_pri_option.rom_stage = ROMStage::OFFLINE; + else if (str.compare("online") == 0) + rom_pri_option.rom_stage = ROMStage::ONLINE; + else if (str.compare("build") == 0) + rom_pri_option.rom_stage = ROMStage::BUILD; + else if (str.compare("online_pinned_H2O_3dof") == 0) + rom_pri_option.rom_stage = ROMStage::ONLINE_PINNED_H2O_3DOF; + else if (str.compare("test_orbital") == 0) + rom_pri_option.rom_stage = ROMStage::TEST_ORBITAL; + else if (str.compare("online_poisson") == 0) + rom_pri_option.rom_stage = ROMStage::ONLINE_POISSON; + else if (str.compare("test_poisson") == 0) + rom_pri_option.rom_stage = ROMStage::TEST_POISSON; + else if (str.compare("test_rho") == 0) + rom_pri_option.rom_stage = ROMStage::TEST_RHO; + else if (str.compare("test_ion") == 0) + rom_pri_option.rom_stage = ROMStage::TEST_ION; + else if (str.compare("none") == 0) + rom_pri_option.rom_stage = ROMStage::UNSUPPORTED; + + rom_pri_option.restart_file_fmt = vm["ROM.offline.restart_filefmt"].as(); + rom_pri_option.restart_file_minidx = vm["ROM.offline.restart_min_idx"].as(); + rom_pri_option.restart_file_maxidx = vm["ROM.offline.restart_max_idx"].as(); + rom_pri_option.basis_file = vm["ROM.offline.basis_file"].as(); + + str = vm["ROM.offline.variable"].as(); + if (str.compare("orbitals") == 0) + rom_pri_option.variable = ROMVariable::ORBITALS; + else if (str.compare("potential") == 0) + rom_pri_option.variable = ROMVariable::POTENTIAL; + else + rom_pri_option.variable = ROMVariable::NONE; + + rom_pri_option.save_librom_snapshot = vm["ROM.offline.save_librom_snapshot"].as(); + rom_pri_option.librom_snapshot_freq = vm["ROM.offline.librom_snapshot_freq"].as(); + + rom_pri_option.compare_md = vm["ROM.basis.compare_md"].as(); + rom_pri_option.num_orbbasis = vm["ROM.basis.number_of_orbital_basis"].as(); + rom_pri_option.num_potbasis = vm["ROM.basis.number_of_potential_basis"].as(); + rom_pri_option.pot_rom_file = vm["ROM.potential_rom_file"].as(); + } // onpe0 + + // synchronize all processors + syncROMOptions(); +} + +void Control::syncROMOptions() +{ + if (onpe0 && verbose > 0) + (*MPIdata::sout) << "Control::syncROMOptions()" << std::endl; + + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + + mmpi.bcast(rom_pri_option.restart_file_fmt, comm_global_); + mmpi.bcast(rom_pri_option.basis_file, comm_global_); + mmpi.bcast(rom_pri_option.pot_rom_file, comm_global_); + + auto bcast_check = [](int mpirc) { + if (mpirc != MPI_SUCCESS) + { + (*MPIdata::sout) << "MPI Bcast of Control failed!!!" << std::endl; + MPI_Abort(comm_global_, 2); + } + }; + + short rom_stage = (short)static_cast(rom_pri_option.rom_stage); + int mpirc; + mpirc = MPI_Bcast(&rom_stage, 1, MPI_SHORT, 0, comm_global_); + bcast_check(mpirc); + + mpirc = MPI_Bcast(&rom_pri_option.restart_file_minidx, 1, MPI_INT, 0, comm_global_); + bcast_check(mpirc); + + mpirc = MPI_Bcast(&rom_pri_option.restart_file_maxidx, 1, MPI_INT, 0, comm_global_); + bcast_check(mpirc); + + mpirc = MPI_Bcast(&rom_pri_option.save_librom_snapshot, 1, MPI_C_BOOL, 0, comm_global_); + bcast_check(mpirc); + + mpirc = MPI_Bcast(&rom_pri_option.librom_snapshot_freq, 1, MPI_INT, 0, comm_global_); + bcast_check(mpirc); + + short rom_var = (short)static_cast(rom_pri_option.variable); + mpirc = MPI_Bcast(&rom_var, 1, MPI_SHORT, 0, comm_global_); + bcast_check(mpirc); + + rom_pri_option.rom_stage = static_cast(rom_stage); + rom_pri_option.variable = static_cast(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); +} diff --git a/src/Control.h b/src/Control.h index 115828e9..0956eb07 100644 --- a/src/Control.h +++ b/src/Control.h @@ -13,6 +13,10 @@ #include "Species.h" #include "Timeout.h" +/* enumeration and option variables for libROM */ +#include "mgmol_config.h" +#include "rom_Control.h" + #include #include #include @@ -220,6 +224,9 @@ class Control void printRestartLink(); + /* libROM related options */ + ROMPrivateOptions rom_pri_option; + public: static Control* instance() { @@ -723,6 +730,11 @@ class Control } bool AtomsMove() { return (atoms_dyn_ != 0); } + + /* ROM-related options */ + void setROMOptions(const boost::program_options::variables_map& vm); + void syncROMOptions(); + const ROMPrivateOptions getROMOptions() { return rom_pri_option; } }; #endif diff --git a/src/Electrostatic.h b/src/Electrostatic.h index 06a9e17c..619da18a 100644 --- a/src/Electrostatic.h +++ b/src/Electrostatic.h @@ -47,6 +47,7 @@ class Electrostatic ~Electrostatic(); static Timer solve_tm() { return solve_tm_; } + const bool isDielectric() { return diel_flag_; } pb::GridFunc* getRhoc() { return grhoc_; } Poisson* getPoissonSolver() { return poisson_solver_; } diff --git a/src/ExtendedGridOrbitals.cc b/src/ExtendedGridOrbitals.cc index 589e1f5b..48eafce6 100644 --- a/src/ExtendedGridOrbitals.cc +++ b/src/ExtendedGridOrbitals.cc @@ -30,6 +30,10 @@ #include #include +#ifdef MGMOL_HAS_LIBROM +#include "librom.h" +#endif + #define ORBITAL_OCCUPATION 2. std::string getDatasetName(const std::string& name, const int color); @@ -1810,6 +1814,28 @@ void ExtendedGridOrbitals::initWF( #endif } +#ifdef MGMOL_HAS_LIBROM +template +void ExtendedGridOrbitals::set(std::string file_path, int rdim) +{ + const int dim = getLocNumpt(); + + CAROM::BasisReader reader(file_path); + CAROM::Matrix* orbital_basis = reader.getSpatialBasis(rdim); + + Control& ct = *(Control::instance()); + Mesh* mymesh = Mesh::instance(); + pb::GridFunc gf_psi(mymesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]); + CAROM::Vector psi; + for (int i = 0; i < rdim; ++i) + { + orbital_basis->getColumn(i, psi); + gf_psi.assign(psi.getData()); + setPsi(gf_psi, i); + } +} +#endif + template void ExtendedGridOrbitals::axpy( const ORBDTYPE alpha, const ExtendedGridOrbitals&); diff --git a/src/ExtendedGridOrbitals.h b/src/ExtendedGridOrbitals.h index c641459d..e1fb667c 100644 --- a/src/ExtendedGridOrbitals.h +++ b/src/ExtendedGridOrbitals.h @@ -402,6 +402,11 @@ class ExtendedGridOrbitals : public Orbitals const pb::Grid& mygrid = mymesh->grid(); return mygrid.maxDomainSize(); } + +#ifdef MGMOL_HAS_LIBROM + void set(std::string file_path, int rdim); +#endif + }; #endif diff --git a/src/Ions.cc b/src/Ions.cc index 7dd39229..248323ae 100644 --- a/src/Ions.cc +++ b/src/Ions.cc @@ -200,6 +200,8 @@ void Ions::setup() ions_setup_tm.start(); + Ion::resetIndexCount(); + updateListIons(); //#ifndef NDEBUG diff --git a/src/MGmol.cc b/src/MGmol.cc index cf57a746..0dcd2092 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -108,6 +108,7 @@ extern Timer md_moveVnuc_tm; extern Timer md_updateMasks_tm; extern Timer md_extrapolateOrbitals_tm; extern Timer md_updateRhoAndPot_tm; +extern Timer md_updateDMandEnergy_tm; extern Timer quench_tm; extern Timer ions_setupInteractingIons_tm; extern Timer ions_setup_tm; @@ -934,6 +935,7 @@ void MGmol::printTimers() init_nuc_tm_.print(os_); md_updateMasks_tm.print(os_); md_extrapolateOrbitals_tm.print(os_); + md_updateDMandEnergy_tm.print(os_); quench_tm.print(os_); evnl_tm_.print(os_); ions_setupInteractingIons_tm.print(os_); @@ -1130,6 +1132,18 @@ void MGmol::dumpRestart() if (ierr < 0) os_ << "WARNING: writing restart data failed!!!" << std::endl; + +#ifdef MGMOL_HAS_LIBROM + // Save orbital snapshots + if (ct.getROMOptions().save_librom_snapshot > 0 && ct.AtomsDynamic() == AtomsDynamicType::Quench) + { + ierr = save_orbital_snapshot( + filename, *current_orbitals_); + + if (ierr < 0) + os_ << "WARNING: writing ROM snapshot data failed!!!" << std::endl; + } +#endif } } @@ -1442,6 +1456,33 @@ void MGmol::getAtomicNumbers(std::vector& an) ions_->getAtomicNumbers(an); } +template +void MGmol::updateDMandEnergy(OrbitalsType& orbitals, Ions& ions, double& eks) +{ + // initialize electronic density + rho_->update(orbitals); + + // initialize potential + update_pot(ions); + + // initialize projected matrices + updateHmatrix(orbitals, ions); + proj_matrices_->updateThetaAndHB(); + + // compute DM + std::shared_ptr> dm_strategy( + DMStrategyFactory>::create(comm_, os_, ions, + rho_.get(), energy_.get(), electrostat_.get(), + hamiltonian_.get(), this, proj_matrices_.get(), &orbitals)); + + dm_strategy->update(orbitals); + + // evaluate energy and forces + double ts = 0.; + eks = energy_->evaluateTotal(ts, proj_matrices_.get(), ions, orbitals, 2, os_); +} + template double MGmol::evaluateEnergyAndForces( const std::vector& tau, const std::vector& atnumbers, diff --git a/src/MGmol.h b/src/MGmol.h index 2975fe42..6738c58f 100644 --- a/src/MGmol.h +++ b/src/MGmol.h @@ -10,6 +10,8 @@ #ifndef MGMOL_H #define MGMOL_H +#include "mgmol_config.h" + #include "Energy.h" #include "GridFuncVector.h" #include "Hamiltonian.h" @@ -42,6 +44,7 @@ class IonicAlgorithm; #include "AOMMprojector.h" #include "ClusterOrbitals.h" #include "DMStrategy.h" +#include "Energy.h" #include "ExtendedGridOrbitals.h" #include "Forces.h" #include "Ions.h" @@ -179,6 +182,7 @@ class MGmol : public MGmolInterface return hamiltonian_; } std::shared_ptr> getRho() { return rho_; } + std::shared_ptr getIons() { return ions_; } void run() override; @@ -336,6 +340,11 @@ class MGmol : public MGmolInterface { forces_->force(orbitals, ions); } + void setPositions(const std::vector& positions, + const std::vector& atnumbers) + { + ions_->setPositions(positions, atnumbers); + } /* * simply dump current state @@ -348,6 +357,12 @@ class MGmol : public MGmolInterface { return proj_matrices_; } + +#ifdef MGMOL_HAS_LIBROM + int save_orbital_snapshot(std::string file_path, OrbitalsType& orbitals); + void project_orbital(std::string file_path, int rdim, OrbitalsType& orbitals); +#endif + void updateDMandEnergy(OrbitalsType& orbitals, Ions& ions, double& eks); }; // Instantiate static variables here to avoid clang warnings template diff --git a/src/MGmolInterface.h b/src/MGmolInterface.h index ec913944..3c2384fd 100644 --- a/src/MGmolInterface.h +++ b/src/MGmolInterface.h @@ -41,6 +41,9 @@ class MGmolInterface virtual void getAtomicPositions(std::vector& tau) = 0; virtual void getAtomicNumbers(std::vector& an) = 0; + virtual void setPositions(const std::vector& positions, + const std::vector& atnumbers) + = 0; virtual std::shared_ptr getProjectedMatrices() = 0; virtual void dumpRestart() = 0; diff --git a/src/MGmol_prototypes.h b/src/MGmol_prototypes.h index cedd514a..02e37890 100644 --- a/src/MGmol_prototypes.h +++ b/src/MGmol_prototypes.h @@ -10,9 +10,11 @@ #ifndef MGMOL_PROTOTYPES_H #define MGMOL_PROTOTYPES_H +#include "mgmol_config.h" #include "global.h" #include +namespace po = boost::program_options; class Ions; class KBPsiMatrixSparse; @@ -24,4 +26,8 @@ int read_config(int argc, char** argv, boost::program_options::variables_map& vm, std::string& input_file, std::string& lrs_filename, std::string& constraints_filename, float& total_spin, bool& with_spin); +#ifdef MGMOL_HAS_LIBROM +void setupROMConfigOption(po::options_description &rom_cfg); +#endif + #endif diff --git a/src/Poisson.h b/src/Poisson.h index b64c0301..f7af6468 100644 --- a/src/Poisson.h +++ b/src/Poisson.h @@ -88,6 +88,12 @@ class Poisson : public PoissonInterface Int_vhrho_ = vel * vh_->gdot(rho); Int_vhrhoc_ = vel * vh_->gdot(rhoc); } + + virtual void applyOperator(pb::GridFunc &vh, pb::GridFunc &lhs) + { + std::cerr << "ERROR: Abstract method Poisson::applyOperator()" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + } }; #endif diff --git a/src/Potentials.cc b/src/Potentials.cc index 2fd556f4..960fd74d 100644 --- a/src/Potentials.cc +++ b/src/Potentials.cc @@ -895,6 +895,97 @@ void Potentials::initBackground() if (fabs(background_charge_) < 1.e-10) background_charge_ = 0.; } +void Potentials::evalIonDensityOnSamplePts( + Ions& ions, const std::vector &local_idx, std::vector &sampled_rhoc) +{ + Mesh* mymesh = Mesh::instance(); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const pb::Grid& mygrid = mymesh->grid(); + + const char flag_filter = pot_type(0); + if (flag_filter == 's') + { + if (onpe0) + { + std::cout << "Potentials::evalIonDensityOnSamplePts - flag_filter s is not supported" + << std::endl; + } + mmpi.abort(); + } + + // initialize output vector + sampled_rhoc.resize(local_idx.size()); + for (int d = 0; d < sampled_rhoc.size(); d++) + sampled_rhoc[d] = 0.0; + + // Loop over ions + for (auto& ion : ions.overlappingVL_ions()) + { + const Species& sp(ion->getSpecies()); + + Vector3D position(ion->position(0), ion->position(1), ion->position(2)); + + initializeRadialDataOnSampledPts(position, sp, local_idx, sampled_rhoc); + } + + return; +} + +void Potentials::initializeRadialDataOnSampledPts( + const Vector3D& position, const Species& sp, const std::vector &local_idx, std::vector &sampled_rhoc) +{ + assert(local_idx.size() == sampled_rhoc.size()); + + Control& ct = *(Control::instance()); + + Mesh* mymesh = Mesh::instance(); + const pb::Grid& mygrid = mymesh->grid(); + + const int dim0 = mygrid.dim(0); + const int dim1 = mygrid.dim(1); + const int dim2 = mygrid.dim(2); + + const double start0 = mygrid.start(0); + const double start1 = mygrid.start(1); + const double start2 = mygrid.start(2); + + const double h0 = mygrid.hgrid(0); + const double h1 = mygrid.hgrid(1); + const double h2 = mygrid.hgrid(2); + + Vector3D point(0., 0., 0.); + + const double lrad = sp.lradius(); + + const RadialInter& lpot(sp.local_pot()); + const Vector3D lattice(mygrid.ll(0), mygrid.ll(1), mygrid.ll(2)); + + for(int k = 0; k < local_idx.size(); k++) + { + /* + local_idx provides offset. + offset = iz + iy * dim2 + ix * dim1 * dim2; + evaluate x,y,z indices backward. + */ + int iz = local_idx[k] % dim2; + int ix = local_idx[k] / dim2; + int iy = ix % dim1; + ix /= dim1; + + /* compute grid point position */ + point[0] = start0 + ix * h0; + point[1] = start1 + iy * h1; + point[2] = start2 + iz * h2; + + /* accumulate ion species density */ + const double r = position.minimage(point, lattice, ct.bcPoisson); + if (r < lrad) + sampled_rhoc[k] += sp.getRhoComp(r); + } + + return; +} + int Potentials::read(HDFrestart& h5f_file) { Control& ct = *(Control::instance()); diff --git a/src/Potentials.h b/src/Potentials.h index f5dc4636..3eb41fd3 100644 --- a/src/Potentials.h +++ b/src/Potentials.h @@ -100,6 +100,9 @@ class Potentials void initializeSupersampledRadialDataOnMesh( const Vector3D& position, const Species& sp); + void initializeRadialDataOnSampledPts( + const Vector3D& position, const Species& sp, const std::vector &local_idx, std::vector &sampled_rhoc); + void rescaleRhoComp(); void addBackgroundToRhoComp(); @@ -168,6 +171,8 @@ class Potentials double getChargeInCell() const { return charge_in_cell_; } + const double getBackgroundCharge() const { return background_charge_; } + void getVofRho(std::vector& vrho) const; /*! @@ -203,6 +208,8 @@ class Potentials void resetVhRho2Backup() { vh_rho_ = vh_rho_backup_; } + void evalIonDensityOnSamplePts(Ions& ions, const std::vector &local_idx, std::vector &sampled_rhoc); + #ifdef HAVE_TRICUBIC void readExternalPot(const string filename, const char type); void setupVextTricubic(); diff --git a/src/Rho.h b/src/Rho.h index ea5069da..11063914 100644 --- a/src/Rho.h +++ b/src/Rho.h @@ -71,6 +71,8 @@ class Rho Rho(); ~Rho(){}; + const OrthoType getOrthoType() { return orbitals_type_; } + void setup( const OrthoType orbitals_type, const std::vector>&); void setVerbosityLevel(const int vlevel) { verbosity_level_ = vlevel; } diff --git a/src/md.cc b/src/md.cc index 2e5fe0f9..cc5d1d01 100644 --- a/src/md.cc +++ b/src/md.cc @@ -34,6 +34,12 @@ #include "mgmol_Signal.h" #include "tools.h" +#ifdef MGMOL_HAS_LIBROM +#include "KBPsiMatrixSparse.h" +#include "PinnedH2O.h" +#include "rom_workflows.h" +#endif + #include #include #include @@ -43,6 +49,7 @@ Timer md_tau_tm("md_tau"); Timer md_moveVnuc_tm("md_moveVnuc"); Timer md_updateMasks_tm("md_updateMasks"); Timer md_extrapolateOrbitals_tm("md_extrapolateOrbitals"); +Timer md_updateDMandEnergy_tm("md_updateDMandEnergy"); #define DUMP_MAX_NUM_TRY 5 @@ -398,8 +405,27 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) h5f_file_.reset(); } + bool ROM_MVP = (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF); +#ifdef MGMOL_HAS_LIBROM + // ROM - initialize orbitals and density matrix + if (ROM_MVP) + { + if (onpe0) os_ << "Setup ROM MVP solver..." << std::endl; + ExtendedGridOrbitals** extended_orbitals = reinterpret_cast**>(orbitals); + (*extended_orbitals)->set(ct.getROMOptions().basis_file, ct.numst); + (*extended_orbitals)->orthonormalizeLoewdin(); + (*extended_orbitals)->setDataWithGhosts(true); + (*extended_orbitals)->setIterativeIndex(10); + + std::shared_ptr projmatrices + = getProjectedMatrices(); + projmatrices->setDMuniform(ct.getNelSpin()); + projmatrices->printDM(os_); + } +#endif + // additional SC steps to compensate random start - if (ct.restart_info < 3) + if (ct.restart_info < 3 && !ROM_MVP) { double eks = 0.; quench(**orbitals, ions, ct.max_electronic_steps, 20, eks); @@ -410,6 +436,9 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) bool extrapolated_flag = true; if (ct.dt <= 0.) extrapolated_flag = false; + int librom_snapshot_freq = ct.getROMOptions().librom_snapshot_freq; + if (librom_snapshot_freq == -1) librom_snapshot_freq = ct.md_print_freq; + MDfiles md_files; // main MD iteration loop @@ -428,10 +457,45 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) int retval = 0; bool small_move = true; bool last_move_is_small = true; - do + + bool force_on_ions = true; +#ifdef MGMOL_HAS_LIBROM + // Ions for ROM + Mesh* mymesh = Mesh::instance(); + const pb::Grid& mygrid = mymesh->grid(); + const double lattice[3] = { mygrid.ll(0), mygrid.ll(1), mygrid.ll(2) }; + std::vector positions; + std::vector anumbers; + getAtomicPositions(positions); + getAtomicNumbers(anumbers); + Ions* ROM_ions; + + // Pinned H2O 3 DOF + PinnedH2O H2O_molecule; + if (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF) + { + if (onpe0) os_ << "Rotate Pinned H2O molecule in timestep " << mdstep << std::endl; + H2O_molecule.rotate(positions, anumbers); + if (onpe0) H2O_molecule.print(os_); + ROM_ions = new Ions(positions, anumbers, lattice, ions_->getSpecies()); + setupPotentials(*ROM_ions); + force_on_ions = false; + } +#endif + + if (ROM_MVP) + { + md_updateDMandEnergy_tm.start(); + updateDMandEnergy(**orbitals, *ROM_ions, eks); + md_updateDMandEnergy_tm.stop(); + } + else { retval = quench(**orbitals, ions, ct.max_electronic_steps, 0, eks); + } + do + { // update localization regions if (ct.adaptiveLRs()) { @@ -491,7 +555,58 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) << std::endl; // Compute forces - force(**orbitals, ions); + if (force_on_ions) + force(**orbitals, ions); + +#ifdef MGMOL_HAS_LIBROM + if (ct.getROMOptions().rom_stage == ROMStage::ONLINE_PINNED_H2O_3DOF) + { + force(**orbitals, *ROM_ions); + // Pinned H2O 3 DOF + if (onpe0) os_ << "Transpose rotate the PinnedH2O molecule" << std::endl; + std::vector forces; + ROM_ions->getForces(forces); + H2O_molecule.transpose_rotate(positions, anumbers, forces); + ions.setLocalForces(forces, positions); + delete ROM_ions; + } + + if (ct.getROMOptions().rom_stage == ROMStage::TEST_ORBITAL) + { + if (onpe0) + { + os_ << "Projecting orbitals onto ROM subspaces to compare " + << ((ct.getROMOptions().compare_md) ? "MD dynamics" : "force") << std::endl; + os_ << "Loading ROM basis " << ct.getROMOptions().basis_file << std::endl; + os_ << "ROM basis dimension = " << ct.getROMOptions().num_orbbasis << std::endl; + } + + // Project orbitals to ROM subspace + project_orbital(ct.getROMOptions().basis_file, ct.getROMOptions().num_orbbasis, **orbitals); + if (ct.getROMOptions().compare_md) + { + // overwrite ions force for MD time stepping + force(**orbitals, ions); + } + else + { + // write ROM_ions force for one-step comparison + ROM_ions = new Ions(positions, anumbers, lattice, ions_->getSpecies()); + force(**orbitals, *ROM_ions); + std::string zero = "0"; + if (ions_->getNumIons() < 256 || ct.verbose > 2) + { + if (ct.verbose > 0) ROM_ions->printForcesGlobal(os_); + + } + else if (zero.compare(ct.md_print_filename) == 0) + { + ROM_ions->printForcesLocal(os_); + } + delete ROM_ions; + } + } +#endif // set fion ions.getLocalForces(fion); @@ -639,6 +754,19 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) printWithTimeStamp("dumped restart file...", std::cout); } +#ifdef MGMOL_HAS_LIBROM + // Save orbital snapshots + if (ct.getROMOptions().save_librom_snapshot > 0 && + md_iteration_ % librom_snapshot_freq == 0) + { + int ierr = save_orbital_snapshot( + ct.md_print_filename + "_mdstep" + std::to_string(mdstep), **orbitals); + + if (ierr < 0) + os_ << "WARNING md(): writing ROM snapshot data failed!!!" << std::endl; + } +#endif + md_iterations_tm.stop(); } // md loop diff --git a/src/mgmol_config.h.in b/src/mgmol_config.h.in new file mode 100644 index 00000000..0b099f6c --- /dev/null +++ b/src/mgmol_config.h.in @@ -0,0 +1,19 @@ +// 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 + +// Description: A configuration header for mgmol. +// + +#ifndef INCLUDED_MGMOL_CONFIG_H +#define INCLUDED_MGMOL_CONFIG_H + +/* Have google test library. */ +#cmakedefine MGMOL_HAS_LIBROM + +#endif diff --git a/src/read_config.cc b/src/read_config.cc index d02aecb9..e3272b3b 100644 --- a/src/read_config.cc +++ b/src/read_config.cc @@ -14,6 +14,7 @@ #include #include #include +#include "MGmol_prototypes.h" #include namespace po = boost::program_options; @@ -205,6 +206,10 @@ int read_config(int argc, char** argv, po::variables_map& vm, po::value()->default_value(""), "Output file for dumping cluster information in vtk format"); +#ifdef MGMOL_HAS_LIBROM + setupROMConfigOption(config); +#endif + // Hidden options, will be allowed in config file, but will not be // shown to the user. po::options_description hidden("Hidden options"); @@ -410,3 +415,34 @@ int read_config(int argc, char** argv, po::variables_map& vm, return 0; } + +#ifdef MGMOL_HAS_LIBROM +void setupROMConfigOption(po::options_description &rom_cfg) +{ + rom_cfg.add_options() + ("ROM.stage", po::value()->default_value("none"), + "ROM workflow stage: offline; build; online; none.") + ("ROM.offline.restart_filefmt", po::value()->default_value(""), + "File name format to read for snapshots.") + ("ROM.offline.restart_min_idx", po::value()->default_value(-1), + "Minimum index for snapshot file format.") + ("ROM.offline.restart_max_idx", po::value()->default_value(-1), + "Maximum index for snapshot file format.") + ("ROM.offline.basis_file", po::value()->default_value(""), + "File name for libROM snapshot/POD matrices.") + ("ROM.offline.save_librom_snapshot", po::value()->default_value(false), + "Save libROM snapshot file at FOM simulation.") + ("ROM.offline.librom_snapshot_freq", po::value()->default_value(-1), + "Frequency of saving libROM snapshot file at FOM simulation.") + ("ROM.offline.variable", po::value()->default_value(""), + "FOM variable to perform POD: either orbitals or potential.") + ("ROM.basis.compare_md", po::value()->default_value(false), + "Compare MD or single-step force.") + ("ROM.basis.number_of_orbital_basis", po::value()->default_value(-1), + "Number of orbital POD basis.") + ("ROM.basis.number_of_potential_basis", po::value()->default_value(-1), + "Number of potential POD basis to build Hartree potential ROM operator.") + ("ROM.potential_rom_file", po::value()->default_value(""), + "File name to save/load potential ROM operators."); +} +#endif diff --git a/src/rom.cc b/src/rom.cc new file mode 100644 index 00000000..863d76e9 --- /dev/null +++ b/src/rom.cc @@ -0,0 +1,118 @@ +// 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 "mgmol_config.h" +#ifdef MGMOL_HAS_LIBROM + +#include "LocGridOrbitals.h" +#include "ExtendedGridOrbitals.h" +#include "MGmol.h" + +#include "librom.h" + +#include +#include +#include +#include + +// Save the wavefunction snapshots +template +int MGmol::save_orbital_snapshot(std::string file_path, OrbitalsType& orbitals) +{ + std::string snapshot_filename = file_path; + struct stat s; + if (stat(file_path.c_str(), &s) == 0) + { + if (s.st_mode & S_IFDIR) + { + snapshot_filename = file_path + "/orbital"; + } + else if (s.st_mode & S_IFREG) + { + snapshot_filename = file_path + "_orbital"; + } + else + { + std::cout << file_path << " exists but is not a directory or a file." << std::endl; + return 1; + } + } + + 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, snapshot_filename); + + for (int i = 0; i < totalSamples; ++i) + basis_generator.takeSample(orbitals.getPsi(i)); + + basis_generator.writeSnapshot(); + + return 0; +} + +template +void MGmol::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(); + + 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); + + Control& ct = *(Control::instance()); + Mesh* mesh = Mesh::instance(); + pb::GridFunc 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; + } +} + +//template +//void ExtendedGridOrbitals::set(std::string file_path, int rdim) +//{ +// const int dim = getLocNumpt(); + +// CAROM::BasisReader reader(file_path); +// CAROM::Matrix* orbital_basis = reader.getSpatialBasis(rdim); + +// Control& ct = *(Control::instance()); +// Mesh* mymesh = Mesh::instance(); +// pb::GridFunc gf_psi(mymesh->grid(), ct.bcWF[0], ct.bcWF[1], ct.bcWF[2]); +// CAROM::Vector psi; +// for (int i = 0; i < rdim; ++i) +// { +// orbital_basis->getColumn(i, psi); +// gf_psi.assign(psi.getData()); +// setPsi(gf_psi, i); +// } +//} + +template class MGmol>; +template class MGmol>; + +#endif // MGMOL_HAS_LIBROM diff --git a/src/rom_Control.h b/src/rom_Control.h new file mode 100644 index 00000000..2d143d84 --- /dev/null +++ b/src/rom_Control.h @@ -0,0 +1,63 @@ +// 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 + +#ifndef ROM_CONTROL_H +#define ROM_CONTROL_H + +#include +#include +#include +#include +#include + +enum class ROMStage +{ + OFFLINE, + ONLINE, + RESTORE, // TODO(kevin): what stage is this? + BUILD, + ONLINE_PINNED_H2O_3DOF, + TEST_ORBITAL, + ONLINE_POISSON, + TEST_POISSON, + TEST_RHO, + TEST_ION, + UNSUPPORTED +}; + +enum class ROMVariable +{ + ORBITALS, + POTENTIAL, + NONE +}; + +/* Stored as a private member variable of Control class */ +struct ROMPrivateOptions +{ + ROMStage rom_stage = ROMStage::UNSUPPORTED; + + std::string restart_file_fmt = ""; + int restart_file_minidx = -1; + int restart_file_maxidx = -1; + std::string basis_file = ""; + ROMVariable variable = ROMVariable::NONE; + + /* save librom orbital snapshot matrix at FOM simulation. */ + bool save_librom_snapshot = false; + int librom_snapshot_freq = -1; + + /* options for ROM building */ + bool compare_md = false; + int num_orbbasis = -1; + int num_potbasis = -1; + std::string pot_rom_file = ""; +}; + +#endif // ROM_CONTROL_H diff --git a/src/rom_main.cc b/src/rom_main.cc new file mode 100644 index 00000000..d29419b9 --- /dev/null +++ b/src/rom_main.cc @@ -0,0 +1,178 @@ +// 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 + +// +// main.cc +// +// Description: +// Real grid, finite difference, molecular dynamics program +// for with nonorthogonal localized orbitals. +// +// Uses Mehrstellen operators, multigrid accelerations, and +// non-local pseudopotentials. +// +// Includes LDA and PBE exchange and correlation functionals. +// +// Units: +// Potentials, eigenvalues and operators in Rydberg +// Energies in Hartree +// +#include "rom_workflows.h" + +// A helper function +template +std::ostream& operator<<(std::ostream& os, const std::vector& v) +{ + copy(v.begin(), v.end(), std::ostream_iterator(std::cout, " ")); + return os; +} + +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; + + 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; + + // use configure file if it can be found + // std::string config_filename("mgmol.cfg"); + + // read options 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(); + + Control::setup(global_comm, with_spin, total_spin); + Control& ct = *(Control::instance()); + + ct.setOptions(vm); + + int ret = ct.checkOptions(); + if (ret < 0) return ret; + + unsigned ngpts[3] = { ct.ngpts_[0], ct.ngpts_[1], ct.ngpts_[2] }; + double origin[3] = { ct.ox_, ct.oy_, ct.oz_ }; + const double cell[3] = { ct.lx_, ct.ly_, ct.lz_ }; + Mesh::setup(mmpi.commSpin(), ngpts, origin, cell, ct.lap_type); + + mmpi.bcastGlobal(input_filename); + mmpi.bcastGlobal(lrs_filename); + + /* Get ROM driver mode */ + ROMPrivateOptions rom_options = ct.getROMOptions(); + ROMStage rom_stage = rom_options.rom_stage; + + // Enter main scope + { + MGmolInterface* mgmol; + if (ct.isLocMode()) + mgmol = new MGmol>(global_comm, *MPIdata::sout, + input_filename, lrs_filename, constraints_filename); + else + mgmol = new MGmol>(global_comm, *MPIdata::sout, + input_filename, lrs_filename, constraints_filename); + + mgmol->setup(); + + switch (rom_stage) + { + case (ROMStage::OFFLINE): + if (ct.isLocMode()) + readRestartFiles>(mgmol); + else + readRestartFiles>(mgmol); + break; + + case (ROMStage::BUILD): + if (ct.isLocMode()) + buildROMPoissonOperator>(mgmol); + else + buildROMPoissonOperator>(mgmol); + break; + + case (ROMStage::ONLINE_POISSON): + if (ct.isLocMode()) + runPoissonROM>(mgmol); + else + runPoissonROM>(mgmol); + break; + + case (ROMStage::TEST_POISSON): + if (ct.isLocMode()) + testROMPoissonOperator>(mgmol); + else + testROMPoissonOperator>(mgmol); + break; + + case (ROMStage::TEST_RHO): + if (ct.isLocMode()) + testROMRhoOperator>(mgmol); + else + testROMRhoOperator>(mgmol); + + case (ROMStage::TEST_ION): + if (ct.isLocMode()) + testROMIonDensity>(mgmol); + else + testROMIonDensity>(mgmol); + + break; + + default: + std::cerr << "rom_main error: Unknown ROM stage" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + break; + } + + + + 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; + + // MemTrack::TrackDumpBlocks(); + + // MemTrack::TrackListMemoryUsage(); + + return 0; +} diff --git a/src/rom_workflows.cc b/src/rom_workflows.cc new file mode 100644 index 00000000..00ffefa5 --- /dev/null +++ b/src/rom_workflows.cc @@ -0,0 +1,961 @@ +// 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 "rom_workflows.h" +#include "Electrostatic.h" +#include +#include +#include + +namespace CAROM { + +/* + This is not implemented in libROM yet. + A temporary placeholder until it is merged into libROM. +*/ +void mult(const Matrix &A, const std::vector &sample_row_A, const Matrix &B, Matrix &AB) +{ + CAROM_VERIFY(!B.distributed()); + CAROM_VERIFY(A.distributed() == AB.distributed()); + CAROM_VERIFY(A.numColumns() == B.numRows()); + CAROM_VERIFY(sample_row_A.size() <= A.numRows()); + + const int num_rows = sample_row_A.size(); + const int num_cols = B.numColumns(); + AB.setSize(num_rows, num_cols); + + // Do the multiplication. + const int Acol = A.numColumns(); + for (int r = 0; r < num_rows; ++r) { + double *d_Arow = A.getData() + sample_row_A[r] * Acol; + for (int c = 0; c < num_cols; ++c) { + double *d_Aptr = d_Arow; + double *d_Bptr = B.getData() + c; + double result_val = 0.0; + for (int ac = 0; ac < Acol; ++ac, d_Aptr++) { + result_val += (*d_Aptr) * (*d_Bptr); + d_Bptr += num_cols; + } + AB.item(r, c) = result_val; + } + } +} + +} + +template +std::string string_format( const std::string& format, Args ... args ) +{ + int size_s = std::snprintf( nullptr, 0, format.c_str(), args ... ) + 1; // Extra space for '\0' + if( size_s <= 0 ){ throw std::runtime_error( "Error during formatting." ); } + auto size = static_cast( size_s ); + std::unique_ptr buf( new char[ size ] ); + std::snprintf( buf.get(), size, format.c_str(), args ... ); + return std::string( buf.get(), buf.get() + size - 1 ); // We don't want the '\0' inside +} + +template +void readRestartFiles(MGmolInterface *mgmol_) +{ + Control& ct = *(Control::instance()); + Mesh* mymesh = Mesh::instance(); + const pb::PEenv& myPEenv = mymesh->peenv(); + + ROMPrivateOptions rom_options = ct.getROMOptions(); + /* type of variable we intend to run POD */ + ROMVariable rom_var = rom_options.variable; + + /* number of restart files, start/end indices */ + assert(rom_options.restart_file_minidx >= 0); + assert(rom_options.restart_file_maxidx >= 0); + const int minidx = rom_options.restart_file_minidx; + const int maxidx = rom_options.restart_file_maxidx; + const int num_restart = maxidx - minidx + 1; + + MGmol *mgmol = static_cast *>(mgmol_); + OrbitalsType *orbitals = mgmol->getOrbitals(); + Potentials& pot = mgmol->getHamiltonian()->potential(); + std::string filename; + + /* Determine basis prefix, dimension, and sample size */ + std::string basis_prefix = rom_options.basis_file; + int dim; + int totalSamples = num_restart; + const int chrom_num = orbitals->chromatic_number(); + switch (rom_var) + { + case ROMVariable::ORBITALS: + basis_prefix += "_orbitals"; + dim = orbitals->getLocNumpt(); + /* if orbitals, each sample have chromatic number of wave functions */ + totalSamples *= orbitals->chromatic_number(); + break; + + case ROMVariable::POTENTIAL: + basis_prefix += "_potential"; + dim = pot.size(); + break; + + default: + ct.global_exit(); + break; + } + + /* Initialize libROM classes */ + CAROM::Options svd_options(dim, totalSamples, 1); + CAROM::BasisGenerator basis_generator(svd_options, false, basis_prefix); + + /* Collect the restart files */ + for (int k = minidx; k <= maxidx; k++) + { + filename = string_format(rom_options.restart_file_fmt, k); + mgmol->loadRestartFile(filename); + assert(dim == orbitals->getLocNumpt()); + assert(chrom_num == orbitals->chromatic_number()); + + switch (rom_var) + { + case ROMVariable::ORBITALS: + for (int i = 0; i < chrom_num; ++i) + basis_generator.takeSample(orbitals->getPsi(i)); + break; + + case ROMVariable::POTENTIAL: + basis_prefix += "_potential"; + /* we save hartree potential */ + // TODO: consider create overloaded takeSample with const + basis_generator.takeSample(const_cast(pot.vh_rho().data())); + break; + } + } + basis_generator.writeSnapshot(); + basis_generator.endSamples(); +} + +template +void buildROMPoissonOperator(MGmolInterface *mgmol_) +{ + Control& ct = *(Control::instance()); + Mesh* mymesh = Mesh::instance(); + const pb::PEenv& myPEenv = mymesh->peenv(); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const int rank = mmpi.mypeGlobal(); + const int nprocs = mmpi.size(); + + ROMPrivateOptions rom_options = ct.getROMOptions(); + /* type of variable we intend to run POD */ + ROMVariable rom_var = rom_options.variable; + if (rom_var != ROMVariable::POTENTIAL) + { + std::cerr << "buildROMPoissonOperator error: ROM variable must be POTENTIAL to run this stage!\n" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + } + + /* Load Hartree potential basis matrix */ + std::string basis_file = rom_options.basis_file; + const int num_pot_basis = rom_options.num_potbasis; + CAROM::BasisReader basis_reader(basis_file); + CAROM::Matrix *pot_basis = basis_reader.getSpatialBasis(num_pot_basis); + + /* Load PoissonSolver pointer */ + MGmol *mgmol = static_cast *>(mgmol_); + Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); + + /* GridFunc initialization inputs */ + const pb::Grid &grid(poisson->vh().grid()); + short bc[3]; + for (int d = 0; d < 3; d++) + bc[d] = poisson->vh().bc(d); + + /* Initialize ROM matrix (undistributed) */ + CAROM::Matrix pot_rom(num_pot_basis, num_pot_basis, false); + + pb::GridFunc col_gf(grid, bc[0], bc[1], bc[2]); + pb::GridFunc opcol_gf(grid, bc[0], bc[1], bc[2]); + CAROM::Vector op_col(pot_basis->numRows(), true); + CAROM::Vector rom_col(num_pot_basis, false); + for (int c = 0; c < num_pot_basis; c++) + { + /* copy c-th column librom vector to GridFunc gf_col */ + CAROM::Vector *col = pot_basis->getColumn(c); + col_gf.assign(col->getData(), 'd'); + + /* apply Laplace operator */ + poisson->applyOperator(col_gf, opcol_gf); + + /* get librom view-vector of gf_opcol */ + opcol_gf.init_vect(op_col.getData(), 'd'); + + /* Compute basis projection of the column */ + /* Resulting vector is undistributed */ + pot_basis->transposeMult(op_col, rom_col); + + /* libROM matrix is row-major, so data copy is necessary */ + for (int r = 0; r < num_pot_basis; r++) + pot_rom(r, c) = rom_col(r); + + delete col; + } // for (int c = 0; c < num_pot_basis; c++) + + /* DEIM hyperreduction */ + CAROM::Matrix pot_rhs_rom(num_pot_basis, num_pot_basis, false); + std::vector global_sampled_row(num_pot_basis), sampled_rows_per_proc(nprocs); + DEIM(pot_basis, num_pot_basis, global_sampled_row, sampled_rows_per_proc, + pot_rhs_rom, rank, nprocs); + + /* ROM rescaleTotalCharge operator */ + CAROM::Vector fom_ones(pot_basis->numRows(), true); + CAROM::Vector rom_ones(num_pot_basis, false); + fom_ones = mymesh->grid().vel(); // volume element + pot_basis->transposeMult(fom_ones, rom_ones); + + /* Save ROM operator */ + // write the file from PE0 only + if (MPIdata::onpe0) + { + std::string rom_oper = rom_options.pot_rom_file; + CAROM::HDFDatabase h5_helper; + h5_helper.create(rom_oper); + h5_helper.putInteger("number_of_potential_basis", num_pot_basis); + h5_helper.putDoubleArray("potential_rom_operator", pot_rom.getData(), + num_pot_basis * num_pot_basis, false); + + /* save the inverse as well */ + pot_rom.inverse(); + h5_helper.putDoubleArray("potential_rom_inverse", pot_rom.getData(), + num_pot_basis * num_pot_basis, false); + + /* save right-hand side hyper-reduction operator */ + h5_helper.putDoubleArray("potential_rhs_rom_inverse", pot_rhs_rom.getData(), + num_pot_basis * num_pot_basis, false); + + /* save right-hand side rescaling operator */ + h5_helper.putDoubleArray("potential_rhs_rescaler", rom_ones.getData(), + num_pot_basis, false); + + h5_helper.close(); + } +} + +template +void runPoissonROM(MGmolInterface *mgmol_) +{ + Control& ct = *(Control::instance()); + Mesh* mymesh = Mesh::instance(); + const pb::PEenv& myPEenv = mymesh->peenv(); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const int rank = mmpi.mypeGlobal(); + const int nprocs = mmpi.size(); + + ROMPrivateOptions rom_options = ct.getROMOptions(); + /* type of variable we intend to run POD */ + ROMVariable rom_var = rom_options.variable; + if (rom_var != ROMVariable::POTENTIAL) + { + std::cerr << "runPoissonROM error: ROM variable must be POTENTIAL to run this stage!\n" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + } + + /* Load Hartree potential basis matrix */ + std::string basis_file = rom_options.basis_file; + const int num_pot_basis = rom_options.num_potbasis; + CAROM::BasisReader basis_reader(basis_file); + CAROM::Matrix *pot_basis = basis_reader.getSpatialBasis(num_pot_basis); + + /* initialize rom operator variables */ + CAROM::Matrix pot_rom(num_pot_basis, num_pot_basis, false); + CAROM::Matrix pot_rom_inv(num_pot_basis, num_pot_basis, false); + CAROM::Matrix pot_rhs_rom(num_pot_basis, num_pot_basis, false); + CAROM::Vector pot_rhs_rescaler(num_pot_basis, false); + + /* Load ROM operator */ + // read the file from PE0 only + if (MPIdata::onpe0) + { + std::string rom_oper = rom_options.pot_rom_file; + CAROM::HDFDatabase h5_helper; + h5_helper.open(rom_oper, "r"); + int num_pot_basis_file = -1; + h5_helper.getInteger("number_of_potential_basis", num_pot_basis_file); + CAROM_VERIFY(num_pot_basis_file == num_pot_basis); + + h5_helper.getDoubleArray("potential_rom_operator", pot_rom.getData(), + num_pot_basis * num_pot_basis, false); + + /* load the inverse as well */ + h5_helper.getDoubleArray("potential_rom_inverse", pot_rom_inv.getData(), + num_pot_basis * num_pot_basis, false); + + /* load right-hand side hyper-reduction operator */ + h5_helper.getDoubleArray("potential_rhs_rom_inverse", pot_rhs_rom.getData(), + num_pot_basis * num_pot_basis, false); + + /* load right-hand side rescaling operator */ + h5_helper.getDoubleArray("potential_rhs_rescaler", pot_rhs_rescaler.getData(), + num_pot_basis, false); + + h5_helper.close(); + } +} + +/* test routines */ + +template +void testROMPoissonOperator(MGmolInterface *mgmol_) +{ + Control& ct = *(Control::instance()); + Mesh* mymesh = Mesh::instance(); + const pb::PEenv& myPEenv = mymesh->peenv(); + + ROMPrivateOptions rom_options = ct.getROMOptions(); + + /* Load MGmol pointer and Potentials */ + MGmol *mgmol = static_cast *>(mgmol_); + Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); + Potentials& pot = mgmol->getHamiltonian()->potential(); + const int dim = pot.size(); + printf("pot size: %d\n", dim); + + /* GridFunc initialization inputs */ + const pb::Grid &grid(poisson->vh().grid()); + short bc[3]; + for (int d = 0; d < 3; d++) + bc[d] = poisson->vh().bc(d); + + /* fictitious snapshot numbers */ + const int nsnapshot = 3; + + /* Set compensating charges to zero now */ + pb::GridFunc rhoc(grid, bc[0], bc[1], bc[2]); + rhoc = 0.0; + + /* Generate fictitious right-hand sides and snapshots */ + std::vector> rhs(nsnapshot), fom_sol(nsnapshot); + for (int s = 0; s < nsnapshot; s++) + { + rhs[s].resize(dim); + for (int d = 0; d < dim; d++) + rhs[s][d] = ran0(); + + /* average out for periodic bc */ + pb::GridFunc rhs_gf(grid, bc[0], bc[1], bc[2]); + rhs_gf.assign(rhs[s].data(), 'd'); + double avg = rhs_gf.get_average(); + rhs_gf -= avg; + + /* copy back to rhs */ + rhs_gf.init_vect(rhs[s].data(), 'd'); + + poisson->solve(rhs_gf, rhoc); + + fom_sol[s].resize(dim); + poisson->vh().init_vect(fom_sol[s].data(), 'd'); + + /* check if the solution is correct */ + pb::GridFunc res(grid, bc[0], bc[1], bc[2]); + pb::GridFunc sol_gf(grid, bc[0], bc[1], bc[2]); + sol_gf.assign(fom_sol[s].data()); + /* apply Laplace operator */ + poisson->applyOperator(sol_gf, res); + /* FD operator scales rhs by 4pi */ + res.axpy(- 4. * M_PI, rhs_gf); + printf("FOM res norm: %.3e\n", res.norm2()); + } + + /* Initialize libROM classes */ + std::string basis_prefix = "test_poisson"; + CAROM::Options svd_options(dim, nsnapshot, 1); + CAROM::BasisGenerator basis_generator(svd_options, false, basis_prefix); + + /* Collect snapshots and train POD basis */ + for (int s = 0; s < nsnapshot; s++) + basis_generator.takeSample(fom_sol[s].data()); + basis_generator.endSamples(); + + /* Load POD basis. We use the maximum number of basis vectors. */ + const CAROM::Matrix *pot_basis = basis_generator.getSpatialBasis(); + + /* Check if full projection preserves FOM solution */ + for (int c = 0; c < nsnapshot; c++) + { + CAROM::Vector *fom_sol_vec = nullptr; + /* get librom view-vector of fom_sol */ + fom_sol_vec = new CAROM::Vector(fom_sol[c].data(), pot_basis->numRows(), true, false); + + CAROM::Vector *rom_proj = pot_basis->transposeMult(*fom_sol_vec); + CAROM::Vector *reconstruct = pot_basis->mult(*rom_proj); + + /* error on libROM side */ + CAROM::Vector *librom_error = reconstruct->minus(fom_sol_vec); + printf("librom reconstruction error: %.3e\n", librom_error->norm()); + + /* error on mgmol side */ + pb::GridFunc recon_gf(grid, bc[0], bc[1], bc[2]); + pb::GridFunc fom_gf(grid, bc[0], bc[1], bc[2]); + recon_gf.assign(reconstruct->getData(), 'd'); + fom_gf.assign(fom_sol[c].data(), 'd'); + recon_gf -= fom_gf; + printf("mgmol reconstruction error: %.3e\n", recon_gf.norm2()); + + delete fom_sol_vec; + delete rom_proj; + delete reconstruct; + delete librom_error; + } + + /* Check FOM axpy is equivalent to ROM axpy */ + for (int s = 0; s < nsnapshot; s++) + { + CAROM::Vector fom_res(pot_basis->numRows(), true); + CAROM::Vector rom_res(nsnapshot, false); + CAROM::Vector fom_rhs(pot_basis->numRows(), true); + CAROM::Vector rom_rhs(nsnapshot, false); + + pb::GridFunc res(grid, bc[0], bc[1], bc[2]); + pb::GridFunc fomsol_gf(grid, bc[0], bc[1], bc[2]); + pb::GridFunc mgmol_rhs(grid, bc[0], bc[1], bc[2]); + fomsol_gf.assign(fom_sol[s].data(), 'd'); + mgmol_rhs.assign(rhs[s].data(), 'd'); + + /* apply Laplace operator */ + poisson->applyOperator(fomsol_gf, res); + + /* get librom view-vector of fom_res */ + res.init_vect(fom_res.getData(), 'd'); + pot_basis->transposeMult(fom_res, rom_res); + + /* get librom view-vector of fom_rhs */ + mgmol_rhs.init_vect(fom_rhs.getData(), 'd'); + pot_basis->transposeMult(fom_rhs, rom_rhs); + + /* ROM residual: FD operator scales rhs by 4pi */ + rom_rhs *= 4. * M_PI; + rom_res -= rom_rhs; + printf("ROM res norm: %.3e\n", rom_res.norm()); + + /* FOM residual: FD operator scales rhs by 4pi */ + res.axpy(- 4. * M_PI, mgmol_rhs); + printf("FOM res norm: %.3e\n", res.norm2()); + + /* projection of the residual */ + res.init_vect(fom_res.getData(), 'd'); + CAROM::Vector *res_proj = pot_basis->transposeMult(fom_res); + printf("FOM res projection norm: %.3e\n", res_proj->norm()); + + delete res_proj; + } + + /* Initialize Projection ROM matrix (undistributed) */ + CAROM::Matrix pot_rom(nsnapshot, nsnapshot, false); + + /* Build Projection of Poisson operator */ + for (int c = 0; c < nsnapshot; c++) + { + pb::GridFunc col_gf(grid, bc[0], bc[1], bc[2]); + pb::GridFunc opcol_gf(grid, bc[0], bc[1], bc[2]); + CAROM::Vector op_col(pot_basis->numRows(), true); + + /* copy c-th column librom vector to GridFunc gf_col */ + CAROM::Vector *col = pot_basis->getColumn(c); + col_gf.assign(col->getData(), 'd'); + + /* apply Laplace operator */ + poisson->applyOperator(col_gf, opcol_gf); + + /* get librom view-vector of gf_opcol */ + opcol_gf.init_vect(op_col.getData(), 'd'); + + /* Compute basis projection of the column */ + /* Resulting vector is undistributed */ + CAROM::Vector *rom_col = pot_basis->transposeMult(op_col); + + /* libROM matrix is row-major, so data copy is necessary */ + for (int r = 0; r < nsnapshot; r++) + pot_rom(r, c) = (*rom_col)(r); + + delete col; + delete rom_col; + } // for (int c = 0; c < num_pot_basis; c++) + + /* Inverse of the projection ROM matrix */ + CAROM::Matrix pot_rom_inv(pot_rom); + pot_rom_inv.inverse(); + + /* Check the inverse */ + CAROM::Matrix *identity = pot_rom_inv.mult(pot_rom); + printf("pot_rom_inv * pot_rom = identity\n"); + for (int i = 0; i < nsnapshot; i++) + { + for (int j = 0; j < nsnapshot; j++) + printf("%.3e\t", identity->item(i, j)); + printf("\n"); + } + delete identity; + + /* Test with sample RHS. ROM must be able to 100% reproduce the FOM solution. */ + std::vector rom_sol(0), rom_rhs(0); + std::vector> test_sol(nsnapshot); + for (int s = 0; s < nsnapshot; s++) + { + /* get librom view-vector of rhs[s] */ + CAROM::Vector fom_rhs(rhs[s].data(), dim, true, false); + + /* project onto POD basis */ + rom_rhs.push_back(pot_basis->transposeMult(fom_rhs)); + + /* FOM FD operator scales rhs by 4pi */ + *rom_rhs.back() *= 4. * M_PI; + + /* solve ROM */ + rom_sol.push_back(pot_rom_inv.mult(*rom_rhs.back())); + + /* check ROM solution */ + CAROM::Vector &res(*pot_rom.mult(*rom_sol.back())); + res -= *rom_rhs.back(); + printf("rom res norm: %.3e\n", res.norm()); + + /* initialize lift-up FOM solution */ + test_sol[s].resize(dim); + /* get librom view-vector of test_sol[s] */ + CAROM::Vector test_sol_vec(test_sol[s].data(), dim, true, false); + pot_basis->mult(*rom_sol.back(), test_sol_vec); + } + + /* Compute relative errors */ + for (int s = 0; s < nsnapshot; s++) + { + pb::GridFunc testsol_gf(grid, bc[0], bc[1], bc[2]); + pb::GridFunc fomsol_gf(grid, bc[0], bc[1], bc[2]); + + testsol_gf.assign(test_sol[s].data(), 'd'); + fomsol_gf.assign(fom_sol[s].data(), 'd'); + + testsol_gf -= fomsol_gf; + double rel_error = testsol_gf.norm2() / fomsol_gf.norm2(); + printf("%d-th sample relative error: %.3e\n", s, rel_error); + + if (rel_error > 1.0e-9) + abort(); + } + + /* clean up pointers */ + for (int s = 0; s < nsnapshot; s++) + { + delete rom_sol[s]; + delete rom_rhs[s]; + } +} + +template +void testROMRhoOperator(MGmolInterface *mgmol_) +{ + Control& ct = *(Control::instance()); + Mesh* mymesh = Mesh::instance(); + const int subdivx = mymesh->subdivx(); + const pb::PEenv& myPEenv = mymesh->peenv(); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const int rank = mmpi.mypeGlobal(); + const int nprocs = mmpi.size(); + + // if (ct.isLocMode()) + // printf("LocMode is On!\n"); + // else + // printf("LocMode is Off!\n"); + + ROMPrivateOptions rom_options = ct.getROMOptions(); + + /* Load MGmol pointer and Potentials */ + MGmol *mgmol = static_cast *>(mgmol_); + Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); + Potentials& pot = mgmol->getHamiltonian()->potential(); + std::shared_ptr> rho = NULL; // mgmol->getRho(); + const OrthoType ortho_type = rho->getOrthoType(); + assert(ortho_type == OrthoType::Nonorthogonal); + + /* potential should have the same size as rho */ + const int dim = pot.size(); + + /* number of restart files, start/end indices */ + assert(rom_options.restart_file_minidx >= 0); + assert(rom_options.restart_file_maxidx >= 0); + const int minidx = rom_options.restart_file_minidx; + const int maxidx = rom_options.restart_file_maxidx; + const int num_restart = maxidx - minidx + 1; + + /* Initialize libROM classes */ + /* + In practice, we do not produce rho POD basis. + This rho POD basis is for the sake of verification. + */ + CAROM::Options svd_options(dim, num_restart, 1); + svd_options.static_svd_preserve_snapshot = true; + CAROM::BasisGenerator basis_generator(svd_options, false, "test_rho"); + + /* Collect the restart files */ + std::string filename; + for (int k = minidx; k <= maxidx; k++) + { + filename = string_format(rom_options.restart_file_fmt, k); + mgmol->loadRestartFile(filename); + basis_generator.takeSample(&rho->rho_[0][0]); + } + // basis_generator.writeSnapshot(); + const CAROM::Matrix rho_snapshots(*basis_generator.getSnapshotMatrix()); + basis_generator.endSamples(); + + const CAROM::Matrix *rho_basis = basis_generator.getSpatialBasis(); + CAROM::Matrix *proj_rho = rho_basis->transposeMult(rho_snapshots); + + /* DEIM hyperreduction */ + CAROM::Matrix rho_basis_inv(num_restart, num_restart, false); + std::vector global_sampled_row(num_restart), sampled_rows_per_proc(nprocs); + DEIM(rho_basis, num_restart, global_sampled_row, sampled_rows_per_proc, + rho_basis_inv, rank, nprocs); + if (rank == 0) + { + int num_sample_rows = 0; + for (int k = 0; k < sampled_rows_per_proc.size(); k++) + num_sample_rows += sampled_rows_per_proc[k]; + printf("number of sampled row: %d\n", num_sample_rows); + } + + /* get local sampled row */ + std::vector offsets, sampled_row(sampled_rows_per_proc[rank]); + int num_global_sample = CAROM::get_global_offsets(sampled_rows_per_proc[rank], offsets); + for (int s = 0, gs = offsets[rank]; gs < offsets[rank+1]; gs++, s++) + sampled_row[s] = global_sampled_row[gs]; + + /* load only the first restart file for now */ + const int test_idx = 2; + + filename = string_format(rom_options.restart_file_fmt, test_idx + minidx); + /* + currently, this does not update rho. + computeRhoOnSamplePts computes with the new density matrix, + while mgmol rho remains the same as the initial condition. + Commenting line 516 gives a consistent result, both for the initial condition. + */ + mgmol->loadRestartFile(filename); + + const int nrows = mymesh->locNumpt(); + // printf("mesh::locNumpt: %d\n", nrows); + + OrbitalsType *orbitals = mgmol->getOrbitals(); + // printf("orbitals::locNumpt: %d\n", orbitals->getLocNumpt()); + + /* NOTE(kevin): we assume we only use ProjectedMatrices class */ + // ProjectedMatrices> *proj_matrices = + // static_cast> *>(orbitals->getProjMatrices()); + ProjectedMatricesInterface *proj_matrices = orbitals->getProjMatrices(); + proj_matrices->updateSubMatX(); + SquareLocalMatrices& localX(proj_matrices->getLocalX()); + + // printf("localX nmat: %d\n", localX.nmat()); + // printf("localX n: %d\n", localX.n()); + // printf("localX m: %d\n", localX.m()); + + bool dm_distributed = (localX.nmat() > 1); + assert(!dm_distributed); + + /* copy density matrix */ + CAROM::Matrix dm(localX.getRawPtr(), localX.m(), localX.n(), dm_distributed, true); + + // /* random density matrix */ + // /* NOTE(kevin): Due to rescaleTotalCharge, the result is slightly inconsistent. */ + // CAROM::Matrix dm(localX.m(), localX.n(), dm_distributed, true); + // dm = 0.0; + // dist_matrix::DistMatrix dm_mgmol(proj_matrices->dm()); + // for (int i = 0; i < dm_mgmol.m(); i++) + // { + // for (int j = 0; j < dm_mgmol.m(); j++) + // { + // if (dm_mgmol.getVal(i, j) == 0.0) + // continue; + + // dm(i, j) = dm_mgmol.getVal(i, j) * (0.975 + 0.05 * ran0()); + // dm_mgmol.setVal(i, j, dm(i, j)); + // } + // } + + /* update rho first */ + // rho->computeRho(*orbitals, dm_mgmol); + // rho->update(*orbitals); + + const int chrom_num = orbitals->chromatic_number(); + CAROM::Matrix psi(dim, chrom_num, true); + for (int c = 0; c < chrom_num; c++) + { + ORBDTYPE *d_psi = orbitals->getPsi(c); + for (int d = 0; d < dim ; d++) + psi.item(d, c) = *(d_psi + d); + } + + CAROM::Matrix rom_psi(chrom_num, chrom_num, false); + for (int i = 0; i < chrom_num; i++) + for (int j = 0; j < chrom_num; j++) + rom_psi(i, j) = (i == j) ? 1 : 0; + + /* this will be resized in computeRhoOnSamplePts */ + CAROM::Vector sample_rho(1, true); + + computeRhoOnSamplePts(dm, psi, rom_psi, sampled_row, sample_rho); + + for (int s = 0; s < sampled_row.size(); s++) + { + const double error = abs(rho->rho_[0][sampled_row[s]] - sample_rho(s)); + if (error > 1.0e-4) + printf("rank %d, rho[%d]: %.5e, sample_rho: %.5e, librom_snapshot: %.5e\n", + rank, sampled_row[s], rho->rho_[0][sampled_row[s]], sample_rho(s), rho_snapshots(sampled_row[s], test_idx)); + CAROM_VERIFY(error < 1.0e-4); + } + + sample_rho.gather(); + + CAROM::Vector *rom_rho = rho_basis_inv.mult(sample_rho); + for (int d = 0; d < rom_rho->dim(); d++) + { + if ((rank == 0) && (abs(proj_rho->item(d, test_idx) - rom_rho->item(d)) > 1.0e-3)) + printf("rom_rho error: %.3e\n", abs(proj_rho->item(d, test_idx) - rom_rho->item(d))); + CAROM_VERIFY(abs(proj_rho->item(d, test_idx) - rom_rho->item(d)) < 1.0e-3); + } + + CAROM::Vector *fom_rho = rho_basis->mult(*rom_rho); + + CAROM_VERIFY(fom_rho->dim() == rho->rho_[0].size()); + for (int d = 0; d < fom_rho->dim(); d++) + CAROM_VERIFY(abs(fom_rho->item(d) - rho->rho_[0][d]) < 1.0e-4); + + delete rom_rho; + delete fom_rho; +} + +/* + dm: density matrix converted to CAROM::Matrix + phi_basis: POD basis matrix of orbitals, or orbitals themselves + rom_phi: ROM coefficients of POD Basis. If phi_basis is orbitals themselves, then simply an identity. + local_idx: sampled local grid indices on this processor. + sampled_rho: FOM density on sampled grid points. +*/ +void computeRhoOnSamplePts(const CAROM::Matrix &dm, + const CAROM::Matrix &phi_basis, const CAROM::Matrix &rom_phi, + const std::vector &local_idx, CAROM::Vector &sampled_rho) +{ + assert(!dm.distributed()); + assert(sampled_rho.distributed() == phi_basis.distributed()); + + // this will be resized. + CAROM::Matrix sampled_phi(1, 1, phi_basis.distributed()); + + CAROM::mult(phi_basis, local_idx, rom_phi, sampled_phi); + + /* same product as in computeRhoSubdomainUsingBlas3 */ + CAROM::Matrix product(1, 1, sampled_phi.distributed()); + sampled_phi.mult(dm, product); + + sampled_rho.setSize(sampled_phi.numRows()); + double *d_product = product.getData(); + double *d_phi = sampled_phi.getData(); + for (int d = 0; d < sampled_rho.dim(); d++) + { + double val = 0.0; + /* CAROM Matrices are row-major */ + for (int c = 0; c < sampled_phi.numColumns(); c++, d_product++, d_phi++) + val += (*d_product) * (*d_phi); + + sampled_rho(d) = val; + } + + /* + TODO(kevin): need to figure out what these functions do. + and probably need to make ROM-equivalent functions with another hyper-reduction? + */ + // gatherSpin(); + + /* + rescaleTotalCharge is handled after hyperreduction. + */ + // rescaleTotalCharge(); +} + +template +void testROMIonDensity(MGmolInterface *mgmol_) +{ + /* random number generator */ + static std::random_device rd; // Will be used to obtain a seed for the random number engine + static std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd(){} + static std::uniform_real_distribution<> dis(0.0, 1.0); + + Control& ct = *(Control::instance()); + Mesh* mymesh = Mesh::instance(); + const pb::Grid& mygrid = mymesh->grid(); + const int subdivx = mymesh->subdivx(); + const pb::PEenv& myPEenv = mymesh->peenv(); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + const int rank = mmpi.mypeGlobal(); + const int nprocs = mmpi.size(); + + ROMPrivateOptions rom_options = ct.getROMOptions(); + + /* Load MGmol pointer and Potentials */ + MGmol *mgmol = static_cast *>(mgmol_); + Poisson *poisson = mgmol->electrostat_->getPoissonSolver(); + Potentials& pot = mgmol->getHamiltonian()->potential(); + const int dim = pot.size(); + std::shared_ptr ions = mgmol->getIons(); + + /* get the extent of global domain */ + const double origin[3] = { mygrid.origin(0), mygrid.origin(1), mygrid.origin(2) }; + const double lattice[3] = { mygrid.ll(0), mygrid.ll(1), mygrid.ll(2) }; + if (rank == 0) + { + printf("origin: (%.3e, %.3e, %.3e)\n", origin[0], origin[1], origin[2]); + printf("lattice: (%.3e, %.3e, %.3e)\n", lattice[0], lattice[1], lattice[2]); + } + + /* get global atomic numbers */ + const int num_ions = ions->getNumIons(); + std::vector atnumbers(num_ions); + ions->getAtomicNumbers(atnumbers); + + assert(!(mgmol->electrostat_->isDielectric())); + assert(pot.getBackgroundCharge() <= 0.0); + + const int num_snap = 3; + + /* 3 fictitious ion configurations */ + std::vector> cfgs(num_snap); + for (int idx = 0; idx < num_snap; idx++) + { + cfgs[idx].resize(3 * num_ions); + if (rank == 0) + for (int k = 0; k < num_ions; k++) + for (int d = 0; d < 3; d++) + cfgs[idx][3 * k + d] = origin[d] + lattice[d] * dis(gen); + + mmpi.bcastGlobal(cfgs[idx].data(), 3 * num_ions, 0); + } + + /* Artificial ions object to avoid repeated setPositions */ + Ions* new_ions; + /* Collect fictitious ion density based on each configuration */ + std::vector> fom_rhoc(num_snap); + /* Sanity check for overlappingVL_ions */ + std::vector>> fom_overlap_ions(num_snap); + for (int idx = 0; idx < num_snap; idx++) + { + /* set ion positions */ + new_ions = new Ions(cfgs[idx], atnumbers, lattice, ions->getSpecies()); + + /* save overlapping ions for sanity check */ + fom_overlap_ions[idx].resize(new_ions->overlappingVL_ions().size()); + for (int k = 0; k < new_ions->overlappingVL_ions().size(); k++) + { + fom_overlap_ions[idx][k].resize(3); + for (int d = 0; d < 3; d++) + fom_overlap_ions[idx][k][d] = new_ions->overlappingVL_ions()[k]->position(d); + } + + /* compute resulting ion density */ + /* NOTE: we exclude rescaling for the sake of verification */ + mgmol->setupPotentials(*new_ions); + + fom_rhoc[idx].resize(dim); + mgmol->electrostat_->getRhoc()->init_vect(fom_rhoc[idx].data(), 'd'); + + delete new_ions; + } + + /* Initialize libROM classes */ + /* + In practice, we do not produce rhoc POD basis. + This rhoc POD basis is for the sake of verification. + */ + CAROM::Options svd_options(dim, num_snap, 1); + svd_options.static_svd_preserve_snapshot = true; + CAROM::BasisGenerator basis_generator(svd_options, false, "test_rhoc"); + for (int idx = 0; idx < num_snap; idx++) + basis_generator.takeSample(&fom_rhoc[idx][0]); + + const CAROM::Matrix rhoc_snapshots(*basis_generator.getSnapshotMatrix()); + basis_generator.endSamples(); + + const CAROM::Matrix *rhoc_basis = basis_generator.getSpatialBasis(); + CAROM::Matrix *proj_rhoc = rhoc_basis->transposeMult(rhoc_snapshots); + + /* DEIM hyperreduction */ + CAROM::Matrix rhoc_basis_inv(num_snap, num_snap, false); + std::vector global_sampled_row(num_snap), sampled_rows_per_proc(nprocs); + DEIM(rhoc_basis, num_snap, global_sampled_row, sampled_rows_per_proc, + rhoc_basis_inv, rank, nprocs); + if (rank == 0) + { + int num_sample_rows = 0; + for (int k = 0; k < sampled_rows_per_proc.size(); k++) + num_sample_rows += sampled_rows_per_proc[k]; + printf("number of sampled row: %d\n", num_sample_rows); + } + + /* get local sampled row */ + std::vector offsets, sampled_row(sampled_rows_per_proc[rank]); + int num_global_sample = CAROM::get_global_offsets(sampled_rows_per_proc[rank], offsets); + for (int s = 0, gs = offsets[rank]; gs < offsets[rank+1]; gs++, s++) + sampled_row[s] = global_sampled_row[gs]; + + /* test one solution */ + std::uniform_int_distribution<> distrib(0, num_snap-1); + int test_idx = distrib(gen); + mmpi.bcastGlobal(&test_idx); + if (rank == 0) printf("test index: %d\n", test_idx); + + /* set ion positions */ + new_ions = new Ions(cfgs[test_idx], atnumbers, lattice, ions->getSpecies()); + + /* Sanity check for overlapping ions */ + CAROM_VERIFY(fom_overlap_ions[test_idx].size() == new_ions->overlappingVL_ions().size()); + for (int k = 0; k < new_ions->overlappingVL_ions().size(); k++) + for (int d = 0; d < 3; d++) + CAROM_VERIFY(abs(fom_overlap_ions[test_idx][k][d] - new_ions->overlappingVL_ions()[k]->position(d)) < 1.0e-12); + + /* set up potentials */ + mgmol->setupPotentials(*new_ions); + + /* eval ion density on sample grid points */ + std::vector sampled_rhoc(sampled_row.size()); + pot.evalIonDensityOnSamplePts(*new_ions, sampled_row, sampled_rhoc); + + // For now, we relax the threshold to allow the slight difference of the values at the sampled indices. + // This is because the rescaling procedure after getting the radial data on mesh requires global information, + // thus cannot be done in the ROM level only with sampled values. + // After we have the DEIM reconstruction, we will do rescaling and revisit this. + for (int d = 0; d < sampled_row.size(); d++) + { + printf("rank %d, fom rhoc[%d]: %.3e, rom rhoc: %.3e\n", rank, sampled_row[d], fom_rhoc[test_idx][sampled_row[d]], sampled_rhoc[d]); + CAROM_VERIFY(abs(fom_rhoc[test_idx][sampled_row[d]] - sampled_rhoc[d]) < 1.0e-2); + } + + delete new_ions; +} + +template void readRestartFiles>(MGmolInterface *mgmol_); +template void readRestartFiles>(MGmolInterface *mgmol_); + +template void buildROMPoissonOperator>(MGmolInterface *mgmol_); +template void buildROMPoissonOperator>(MGmolInterface *mgmol_); + +template void runPoissonROM>(MGmolInterface *mgmol_); +template void runPoissonROM>(MGmolInterface *mgmol_); + +template void testROMPoissonOperator>(MGmolInterface *mgmol_); +template void testROMPoissonOperator>(MGmolInterface *mgmol_); + +template void testROMRhoOperator>(MGmolInterface *mgmol_); +template void testROMRhoOperator>(MGmolInterface *mgmol_); + +template void testROMIonDensity>(MGmolInterface *mgmol_); +template void testROMIonDensity>(MGmolInterface *mgmol_); diff --git a/src/rom_workflows.h b/src/rom_workflows.h new file mode 100644 index 00000000..9222249f --- /dev/null +++ b/src/rom_workflows.h @@ -0,0 +1,64 @@ +// 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 + +#ifndef ROM_WORKFLOWS_H +#define ROM_WORKFLOWS_H + +#include "Control.h" +#include "ExtendedGridOrbitals.h" +#include "ProjectedMatrices.h" +#include "LocGridOrbitals.h" +#include "Potentials.h" +#include "MGmol.h" +#include "MGmol_MPI.h" +#include "MPIdata.h" +#include "Mesh.h" +#include "mgmol_run.h" +#include "tools.h" + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +namespace po = boost::program_options; + +#include "librom.h" +#include "utils/HDFDatabase.h" +#include "utils/mpi_utils.h" + +template +void readRestartFiles(MGmolInterface *mgmol_); + +template +void buildROMPoissonOperator(MGmolInterface *mgmol_); + +template +void runPoissonROM(MGmolInterface *mgmol_); + +template +void testROMPoissonOperator(MGmolInterface *mgmol_); + +template +void testROMRhoOperator(MGmolInterface *mgmol_); + +template +void testROMIonDensity(MGmolInterface *mgmol_); + +void computeRhoOnSamplePts(const CAROM::Matrix &dm, + const CAROM::Matrix &phi_basis, const CAROM::Matrix &rom_phi, + const std::vector &local_idx, CAROM::Vector &sampled_rho); + +#endif // ROM_WORKFLOWS_H diff --git a/src/tools/CMakeLists.txt b/src/tools/CMakeLists.txt index baae6cda..ca187a0e 100644 --- a/src/tools/CMakeLists.txt +++ b/src/tools/CMakeLists.txt @@ -8,6 +8,7 @@ set(SOURCES mgmol_mpi_tools.cc random.cc coloring.cc SymmetricMatrix.cc + PinnedH2O.cc ) add_library(mgmol_tools ${SOURCES}) diff --git a/src/tools/PinnedH2O.cc b/src/tools/PinnedH2O.cc new file mode 100644 index 00000000..22f4677a --- /dev/null +++ b/src/tools/PinnedH2O.cc @@ -0,0 +1,237 @@ +// 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 "PinnedH2O.h" +#include + +PinnedH2O::PinnedH2O() + : flipped_bond(false), + O1_idx(-1), + H1_idx(-1), + H2_idx(-1), + planar_rotation_angle(0.0) +{ + for (int i = 0; i < 3; ++i) { + for (int j = 0; j < 3; ++j) { + out_of_plane_rotation_matrix[i][j] = 0.0; + } + } +} + +double PinnedH2O::calculate_bondlength(const double atom1[3], const double atom2[3]) const +{ + return sqrt(pow(atom1[0] - atom2[0], 2) + pow(atom1[1] - atom2[1], 2) + pow(atom1[2] - atom2[2], 2)); +} + +double PinnedH2O::calculate_bondangle(const double atom1[3], const double atom2[3], const double atom3[3]) const +{ + 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); + + return angle; +} + +void PinnedH2O::rotation_matrix(const double axis[3], double angle, double matrix[3][3]) const +{ + 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 PinnedH2O::normalize(double vec[3]) const +{ + double norm = sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]); + vec[0] /= norm; + vec[1] /= norm; + vec[2] /= norm; +} + +void PinnedH2O::cross(const double a[3], const double b[3], double result[3]) const +{ + 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 PinnedH2O::apply_rotation(const double matrix[3][3], const double vec[3], double result[3]) const +{ + result[0] = matrix[0][0] * vec[0] + matrix[0][1] * vec[1] + matrix[0][2] * vec[2]; + result[1] = matrix[1][0] * vec[0] + matrix[1][1] * vec[1] + matrix[1][2] * vec[2]; + result[2] = matrix[2][0] * vec[0] + matrix[2][1] * vec[1] + matrix[2][2] * vec[2]; +} + +void PinnedH2O::apply_transpose_rotation(const double matrix[3][3], const double vec[3], double result[3]) const +{ + 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]; +} + +void PinnedH2O::rotate(std::vector& positions, std::vector& anumbers) +{ + for (int i = 0; i < 3; i++) { + if (positions[3 * i] == 0.0 && positions[3 * i + 1] == 0.0 && positions[3 * i + 2] == 0.0) { + O1_idx = i; + break; + } + } + if (O1_idx == -1) return; + + H1_idx = (O1_idx + 1) % 3; + H2_idx = (O1_idx + 2) % 3; + + double O1[3] = {positions[3 * O1_idx], positions[3 * O1_idx + 1], positions[3 * O1_idx + 2]}; + double H1[3] = {positions[3 * H1_idx], positions[3 * H1_idx + 1], positions[3 * H1_idx + 2]}; + double H2[3] = {positions[3 * H2_idx], positions[3 * H2_idx + 1], positions[3 * H2_idx + 2]}; + + bondlength1 = calculate_bondlength(H1, O1); + bondlength2 = calculate_bondlength(H2, O1); + bondangle = calculate_bondangle(H1, O1, H2); + + double H1_temp[3], H2_temp[3]; + double H1_rotated[3], H2_rotated[3]; + + double plane_normal[3]; + cross(H2, H1, plane_normal); + normalize(plane_normal); + double target_plane_normal[3] = {0, 0, 1}; + 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 = std::acos(std::min(std::max(dot_product, -1.0), 1.0)); + double axis_to_align[3]; + if (abs(dot_product) > 1.0 - 1e-8) + { + axis_to_align[0] = 1.0; + axis_to_align[1] = 0.0; + axis_to_align[2] = 0.0; + } + else + { + cross(plane_normal, target_plane_normal, axis_to_align); + normalize(axis_to_align); + } + + rotation_matrix(axis_to_align, angle_to_align, out_of_plane_rotation_matrix); + apply_rotation(out_of_plane_rotation_matrix, H1, H1_temp); + apply_rotation(out_of_plane_rotation_matrix, H2, H2_temp); + + double theta1 = std::atan2(H1_temp[1], H1_temp[0]); + planar_rotation_angle = -theta1 + bondangle / 2.0; + double planar_rotation_matrix[3][3]; + rotation_matrix(target_plane_normal, planar_rotation_angle, planar_rotation_matrix); + apply_rotation(planar_rotation_matrix, H1_temp, H1_rotated); + apply_rotation(planar_rotation_matrix, H2_temp, H2_rotated); + + flipped_bond = (bondlength1 < bondlength2); + if (flipped_bond) + { + H1_rotated[1] *= -1.0; + H2_rotated[1] *= -1.0; + std::swap(H1_rotated, H2_rotated); + std::swap(bondlength1, bondlength2); + } + + positions[0] = H2_rotated[0]; + positions[1] = H2_rotated[1]; + positions[2] = H2_rotated[2]; + positions[3] = 0.0; + positions[4] = 0.0; + positions[5] = 0.0; + positions[6] = H1_rotated[0]; + positions[7] = H1_rotated[1]; + positions[8] = H1_rotated[2]; + + anumbers[0] = 1; + anumbers[1] = 8; + anumbers[2] = 1; +} + +void PinnedH2O::transpose_rotate(std::vector& positions, std::vector& anumbers, std::vector& forces) +{ + double H2_rotated[3] = {positions[0], positions[1], positions[2]}; + double O1_rotated[3] = {positions[3], positions[4], positions[5]}; + double H1_rotated[3] = {positions[6], positions[7], positions[8]}; + + double f_H2_rotated[3] = {forces[0], forces[1], forces[2]}; + double f_O1_rotated[3] = {forces[3], forces[4], forces[5]}; + double f_H1_rotated[3] = {forces[6], forces[7], forces[8]}; + + if (flipped_bond) + { + H1_rotated[1] *= -1.0; + H2_rotated[1] *= -1.0; + f_O1_rotated[1] *= -1.0; + f_H1_rotated[1] *= -1.0; + f_H2_rotated[1] *= -1.0; + } + + double planar_rotation_matrix[3][3]; + double target_plane_normal[3] = {0, 0, 1}; + rotation_matrix(target_plane_normal, planar_rotation_angle, planar_rotation_matrix); + + double H1_temp[3], H2_temp[3]; + apply_transpose_rotation(planar_rotation_matrix, H1_rotated, H1_temp); + apply_transpose_rotation(planar_rotation_matrix, H2_rotated, H2_temp); + + double f_O1_temp[3], f_H1_temp[3], f_H2_temp[3]; + apply_transpose_rotation(planar_rotation_matrix, f_O1_rotated, f_O1_temp); + apply_transpose_rotation(planar_rotation_matrix, f_H1_rotated, f_H1_temp); + apply_transpose_rotation(planar_rotation_matrix, f_H2_rotated, f_H2_temp); + + double H1_restored[3], H2_restored[3]; + apply_transpose_rotation(out_of_plane_rotation_matrix, H1_temp, H1_restored); + apply_transpose_rotation(out_of_plane_rotation_matrix, H2_temp, H2_restored); + + double f_O1_restored[3], f_H1_restored[3], f_H2_restored[3]; + apply_transpose_rotation(out_of_plane_rotation_matrix, f_O1_temp, f_O1_restored); + apply_transpose_rotation(out_of_plane_rotation_matrix, f_H1_temp, f_H1_restored); + apply_transpose_rotation(out_of_plane_rotation_matrix, f_H2_temp, f_H2_restored); + + positions[3*H2_idx] = H2_restored[0]; + positions[3*H2_idx+1] = H2_restored[1]; + positions[3*H2_idx+2] = H2_restored[2]; + positions[3*O1_idx] = 0.0; + positions[3*O1_idx+1] = 0.0; + positions[3*O1_idx+2] = 0.0; + positions[3*H1_idx] = H1_restored[0]; + positions[3*H1_idx+1] = H1_restored[1]; + positions[3*H1_idx+2] = H1_restored[2]; + + anumbers[H2_idx] = 1; + anumbers[O1_idx] = 8; + anumbers[H1_idx] = 1; + + forces[3*H2_idx] = f_H2_restored[0]; + forces[3*H2_idx+1] = f_H2_restored[1]; + forces[3*H2_idx+2] = f_H2_restored[2]; + forces[3*O1_idx] = f_O1_restored[0]; + forces[3*O1_idx+1] = f_O1_restored[1]; + forces[3*O1_idx+2] = f_O1_restored[2]; + forces[3*H1_idx] = f_H1_restored[0]; + forces[3*H1_idx+1] = f_H1_restored[1]; + forces[3*H1_idx+2] = f_H1_restored[2]; +} diff --git a/src/tools/PinnedH2O.h b/src/tools/PinnedH2O.h new file mode 100644 index 00000000..bee880d8 --- /dev/null +++ b/src/tools/PinnedH2O.h @@ -0,0 +1,54 @@ +// 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 + +#ifndef PINNED_H2O_H +#define PINNED_H2O_H + +#include +#include +#include +#include + +class PinnedH2O +{ + +public: + PinnedH2O(); + ~PinnedH2O() = default; + + void rotate(std::vector& positions, std::vector& anumbers); + void transpose_rotate(std::vector& positions, std::vector& anumbers, std::vector& forces); + void print(std::ostream& os) + { + os << "Bondlengths = " << bondlength1 << ", " << bondlength2 << " Bohrs; " + << "Bondangle = " << bondangle * 180.0 / M_PI << " degrees." << std::endl; + } + +private: + double calculate_bondlength(const double atom1[3], const double atom2[3]) const; + double calculate_bondangle(const double atom1[3], const double atom2[3], const double atom3[3]) const; + void rotation_matrix(const double axis[3], double angle, double matrix[3][3]) const; + void normalize(double vec[3]) const; + void cross(const double a[3], const double b[3], double result[3]) const; + void apply_rotation(const double matrix[3][3], const double vec[3], double result[3]) const; + void apply_transpose_rotation(const double matrix[3][3], const double vec[3], double result[3]) const; + + double bondlength1; + double bondlength2; + double bondangle; + + bool flipped_bond; + int O1_idx; + int H1_idx; + int H2_idx; + double planar_rotation_angle; + double out_of_plane_rotation_matrix[3][3]; +}; + +#endif // PINNED_H2O_H diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 03dffc89..5779198d 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -292,6 +292,9 @@ add_executable(testDMandEnergyAndForces ${CMAKE_SOURCE_DIR}/tests/DMandEnergyAndForces/testDMandEnergyAndForces.cc) add_executable(testRestartEnergyAndForces ${CMAKE_SOURCE_DIR}/tests/RestartEnergyAndForces/testRestartEnergyAndForces.cc) +add_executable(testPinnedH2O_3DOF + ${CMAKE_SOURCE_DIR}/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc + ${CMAKE_SOURCE_DIR}/src/tools/PinnedH2O.cc) if(${MAGMA_FOUND}) add_executable(testOpenmpOffload @@ -416,6 +419,14 @@ add_test(NAME testRestartEnergyAndForces ${CMAKE_CURRENT_SOURCE_DIR}/RestartEnergyAndForces/restart.cfg ${CMAKE_CURRENT_SOURCE_DIR}/RestartEnergyAndForces/h2o.xyz ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) +add_test(NAME testPinnedH2O_3DOF + COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/PinnedH2O_3DOF/test.py + ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} + ${CMAKE_CURRENT_BINARY_DIR}/testPinnedH2O_3DOF + ${CMAKE_CURRENT_SOURCE_DIR}/PinnedH2O_3DOF/mgmol.cfg + ${CMAKE_CURRENT_SOURCE_DIR}/PinnedH2O_3DOF/coords.in + ${CMAKE_CURRENT_SOURCE_DIR}/PinnedH2O_3DOF/lrs.in + ${CMAKE_CURRENT_SOURCE_DIR}/../potentials) add_test(NAME testRhoVhRestart COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/RhoVhRestart/test.py ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS} @@ -664,6 +675,7 @@ target_link_libraries(testEnergyAndForces PRIVATE mgmol_src) target_link_libraries(testWFEnergyAndForces PRIVATE mgmol_src) target_link_libraries(testDMandEnergyAndForces PRIVATE mgmol_src) target_link_libraries(testRestartEnergyAndForces PRIVATE mgmol_src) +target_link_libraries(testPinnedH2O_3DOF PRIVATE mgmol_src) target_link_libraries(testIons PRIVATE mgmol_src) target_link_libraries(testDensityMatrix PRIVATE ${HDF5_LIBRARIES}) target_link_libraries(testRhoVhRestart mgmol_src) diff --git a/tests/PinnedH2O_3DOF/README.md b/tests/PinnedH2O_3DOF/README.md new file mode 100644 index 00000000..90551b47 --- /dev/null +++ b/tests/PinnedH2O_3DOF/README.md @@ -0,0 +1 @@ +python3 test.py srun -p debug -n 64 ../../build_quartz/tests/testPinnedH2O_3DOF mgmol_online.cfg coords_1.00_1.00_0.0.in ../../potentials diff --git a/tests/PinnedH2O_3DOF/assemble_results.sh b/tests/PinnedH2O_3DOF/assemble_results.sh new file mode 100644 index 00000000..6a9a4d0d --- /dev/null +++ b/tests/PinnedH2O_3DOF/assemble_results.sh @@ -0,0 +1,31 @@ +#!/bin/bash + +rm -rf PinnedH2O_3DOF_assembled_results +mkdir PinnedH2O_3DOF_assembled_results + +for d in results_*; do + echo "Processing $d" + rm $d/energy_and_forces_*.out + rm $d/*.png + + suffix=$(echo "$d" | cut -d'_' -f2,3,4) + N_l=$(echo "$suffix" | cut -d'_' -f1) + N_theta=$(echo "$suffix" | cut -d'_' -f2) + rdim=$(echo "$suffix" | cut -d'_' -f3) + python3 compare_energy_and_forces.py --N_l "$N_l" --N_theta "$N_theta" --rdim "$rdim" + + cp "$d/Eks_difference_histogram_reproductive.png" "PinnedH2O_3DOF_assembled_results/Eks_difference_histogram_${suffix}_reproductive.png" + cp "$d/f_H1_difference_histogram_reproductive.png" "PinnedH2O_3DOF_assembled_results/f_H1_difference_histogram_${suffix}_reproductive.png" + cp "$d/f_H2_difference_histogram_reproductive.png" "PinnedH2O_3DOF_assembled_results/f_H2_difference_histogram_${suffix}_reproductive.png" + cp "$d/f_O1_difference_histogram_reproductive.png" "PinnedH2O_3DOF_assembled_results/f_O1_difference_histogram_${suffix}_reproductive.png" + + cp "$d/Eks_difference_histogram_predictive.png" "PinnedH2O_3DOF_assembled_results/Eks_difference_histogram_${suffix}_predictive.png" + cp "$d/f_H1_difference_histogram_predictive.png" "PinnedH2O_3DOF_assembled_results/f_H1_difference_histogram_${suffix}_predictive.png" + cp "$d/f_H2_difference_histogram_predictive.png" "PinnedH2O_3DOF_assembled_results/f_H2_difference_histogram_${suffix}_predictive.png" + cp "$d/f_O1_difference_histogram_predictive.png" "PinnedH2O_3DOF_assembled_results/f_O1_difference_histogram_${suffix}_predictive.png" +done + +tar -cvf PinnedH2O_3DOF_assembled_results.tar PinnedH2O_3DOF_assembled_results +mv PinnedH2O_3DOF_assembled_results.tar /p/lustre2/cheung26/scp_local + +echo "Finished assembling results." diff --git a/tests/PinnedH2O_3DOF/aux/PinnedH2O_3DOF.png b/tests/PinnedH2O_3DOF/aux/PinnedH2O_3DOF.png new file mode 100644 index 00000000..b484f4c9 Binary files /dev/null and b/tests/PinnedH2O_3DOF/aux/PinnedH2O_3DOF.png differ diff --git a/tests/PinnedH2O_3DOF/aux/coords_1.00_1.00_0.0.in b/tests/PinnedH2O_3DOF/aux/coords_1.00_1.00_0.0.in new file mode 100644 index 00000000..46f5681d --- /dev/null +++ b/tests/PinnedH2O_3DOF/aux/coords_1.00_1.00_0.0.in @@ -0,0 +1,3 @@ +O1 1 0.00 0.00 0.00 0 +H1 2 1.12 1.45 0.00 1 +H2 2 1.12 -1.45 0.00 1 diff --git a/tests/PinnedH2O_3DOF/aux/coords_1.02_0.98_2.0.in b/tests/PinnedH2O_3DOF/aux/coords_1.02_0.98_2.0.in new file mode 100644 index 00000000..4c440443 --- /dev/null +++ b/tests/PinnedH2O_3DOF/aux/coords_1.02_0.98_2.0.in @@ -0,0 +1,3 @@ +O1 1 0.00 0.00 0.00 0 +H1 2 1.12 1.50 0.00 1 +H2 2 1.07 -1.44 0.00 1 diff --git a/tests/PinnedH2O_3DOF/aux/get_ROM_table.py b/tests/PinnedH2O_3DOF/aux/get_ROM_table.py new file mode 100644 index 00000000..5fea1e38 --- /dev/null +++ b/tests/PinnedH2O_3DOF/aux/get_ROM_table.py @@ -0,0 +1,31 @@ +import subprocess +import re + +bondlength_num_increments = (2, 5, 10) +bondangle_num_increments = (2, 5, 10) + +pattern = r"For energy fraction: \d+\.\d+, take first (\d+) of \d+ basis vectors" + +print("\\begin{tabular}{|c|c||c|c|c|c|c|c|c|}") +print("\\hline") +print("$N_L$ & $N_\\theta$ & $\\varepsilon = 10^{-1}$ & $\\varepsilon = 10^{-2}$ & $\\varepsilon = 10^{-3}$ & $\\varepsilon = 10^{-4}$ & $\\varepsilon = 10^{-5}$ & Snapshots \\\\") +print("\\hline") + +for _, N_L in enumerate(bondlength_num_increments): + for _, N_theta in enumerate(bondangle_num_increments): + snapshots = 2*(N_L+1)*(N_L+2)*(N_theta+1) + grep_command = f"grep 'take first' basis_{N_L}_{N_theta}_PinnedH2O_3DOF.out" + result = subprocess.run(grep_command, shell=True, capture_output=True, text=True) + matches = re.findall(pattern, result.stdout) + energy_fractions = { + "0.9": matches[0], + "0.99": matches[1], + "0.999": matches[2], + "0.9999": matches[3], + "0.99999": matches[4], + } + line = f"{N_L} & {N_theta} & {energy_fractions['0.9']} & {energy_fractions['0.99']} & {energy_fractions['0.999']} & {energy_fractions['0.9999']} & {energy_fractions['0.99999']} & {snapshots} \\\\" + print(line) + +print("\\hline") +print("\\end{tabular}") diff --git a/tests/PinnedH2O_3DOF/aux/plot_PinnedH2O_3DOF.py b/tests/PinnedH2O_3DOF/aux/plot_PinnedH2O_3DOF.py new file mode 100644 index 00000000..e679c691 --- /dev/null +++ b/tests/PinnedH2O_3DOF/aux/plot_PinnedH2O_3DOF.py @@ -0,0 +1,69 @@ +import matplotlib +import matplotlib.pyplot as plt +import numpy as np + +ref_bondlength = 1.83 +ref_bondangle = 104.5 + +# factors and increments for bond lengths and bond angle +bondlength_factor = np.linspace(0.95, 1.05, 3) +bondangle_increment = np.linspace(-5, 5, 3) + +fig, ax = plt.subplots() + +for d_bondangle in bondangle_increment: + bondangle = ref_bondangle + d_bondangle + x = ref_bondlength * np.cos(np.radians(bondangle / 2)) + y = ref_bondlength * np.sin(np.radians(bondangle / 2)) + for i, f_bondlength1 in enumerate(bondlength_factor): + H1_x = f_bondlength1*x + H1_y = f_bondlength1*y + ax.plot(H1_x, H1_y, 'ro', markersize=8, alpha=0.1) + for f_bondlength2 in bondlength_factor[:(i+1)]: + H2_x = f_bondlength2*x + H2_y = -f_bondlength2*y + ax.plot(H2_x, H2_y, 'bo', markersize=8, alpha=0.1) + +ref_bondangle = np.radians(ref_bondangle) +x = ref_bondlength * np.cos(ref_bondangle / 2) +y = ref_bondlength * np.sin(ref_bondangle / 2) + +ax.plot([0, x], [0, y], 'r-', lw=2) +ax.plot([0, x], [0, -y], 'b-', lw=2) +ax.plot(0, 0, 'ko', markersize=8) +ax.plot(x, y, 'ro', markersize=8) +ax.plot(x, -y, 'bo', markersize=8) + +arc1 = matplotlib.patches.Arc(xy=(0, 0), + width=0.5, + height=0.5, + angle=0, + theta1=0, + theta2=ref_bondangle/2, + color='red', + linewidth=2) +ax.add_patch(arc1) +ax.text(0.4, 0.25, r"$\frac{\theta}{2}$", ha='center', va='center', fontsize=12, color='red') + +arc2 = matplotlib.patches.Arc(xy=(0, 0), + width=0.5, + height=0.5, + angle=0, + theta1=-ref_bondangle/2, + theta2=0, + color='blue', + linewidth=2) +ax.add_patch(arc2) +ax.text(0.4, -0.25, r"$\frac{\theta}{2}$", ha='center', va='center', fontsize=12, color='blue') + +ax.text(x / 2 - 0.1, y / 2, r"$L_1$", ha='right', color='red') +ax.text(x / 2 - 0.1, -y / 2, r"$L_2$", ha='right', color='blue') + +ax.plot([-2, 2], [0, 0], 'k--', lw=1) +ax.set_xlim(-2.0, 2.0) +ax.set_ylim(-2.0, 2.0) +ax.set_aspect('equal') +ax.set_xlabel("x (Bohr)") +ax.set_ylabel("y (Bohr)") +ax.grid(True) +plt.savefig("PinnedH2O_3DOF.png") diff --git a/tests/PinnedH2O_3DOF/aux/rotation_test.cc b/tests/PinnedH2O_3DOF/aux/rotation_test.cc new file mode 100644 index 00000000..21cfd93a --- /dev/null +++ b/tests/PinnedH2O_3DOF/aux/rotation_test.cc @@ -0,0 +1,165 @@ +#include +#include +#include +#include + +using namespace std; + +double calculate_bondlength(const double atom1[3], const double atom2[3]) +{ + return sqrt(pow(atom1[0] - atom2[0], 2) + pow(atom1[1] - atom2[1], 2) + pow(atom1[2] - atom2[2], 2)); +} + +double calculate_bondangle(const double atom1[3], const double atom2[3], const double atom3[3], bool radian) +{ + double vector1[3] = {atom1[0] - atom2[0], atom1[1] - atom2[1], atom1[2] - atom2[2]}; + double vector2[3] = {atom3[0] - atom2[0], atom3[1] - atom2[1], atom3[2] - atom2[2]}; + + double dot_product = vector1[0] * vector2[0] + vector1[1] * vector2[1] + vector1[2] * vector2[2]; + double magnitude_product = sqrt(pow(vector1[0], 2) + pow(vector1[1], 2) + pow(vector1[2], 2)) * + sqrt(pow(vector2[0], 2) + pow(vector2[1], 2) + pow(vector2[2], 2)); + double angle = acos(dot_product / magnitude_product); + + if (!radian) { + angle = angle * 180.0 / M_PI; + } + return angle; +} + +void rotation_matrix(const double axis[3], double angle, double matrix[3][3]) +{ + double cos_theta = cos(angle); + double sin_theta = sin(angle); + double ux = axis[0], uy = axis[1], uz = axis[2]; + + matrix[0][0] = cos_theta + ux * ux * (1 - cos_theta); + matrix[0][1] = ux * uy * (1 - cos_theta) - uz * sin_theta; + matrix[0][2] = ux * uz * (1 - cos_theta) + uy * sin_theta; + + matrix[1][0] = uy * ux * (1 - cos_theta) + uz * sin_theta; + matrix[1][1] = cos_theta + uy * uy * (1 - cos_theta); + matrix[1][2] = uy * uz * (1 - cos_theta) - ux * sin_theta; + + matrix[2][0] = uz * ux * (1 - cos_theta) - uy * sin_theta; + matrix[2][1] = uz * uy * (1 - cos_theta) + ux * sin_theta; + matrix[2][2] = cos_theta + uz * uz * (1 - cos_theta); +} + +void normalize(double vec[3]) +{ + double norm = sqrt(vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2]); + vec[0] /= norm; + vec[1] /= norm; + vec[2] /= norm; +} + +void cross(const double a[3], const double b[3], double result[3]) +{ + result[0] = a[1] * b[2] - a[2] * b[1]; + result[1] = a[2] * b[0] - a[0] * b[2]; + result[2] = a[0] * b[1] - a[1] * b[0]; +} + +void apply_rotation(const double matrix[3][3], const double vec[3], double result[3]) +{ + result[0] = matrix[0][0] * vec[0] + matrix[0][1] * vec[1] + matrix[0][2] * vec[2]; + result[1] = matrix[1][0] * vec[0] + matrix[1][1] * vec[1] + matrix[1][2] * vec[2]; + result[2] = matrix[2][0] * vec[0] + matrix[2][1] * vec[1] + matrix[2][2] * vec[2]; +} + +void apply_transpose_rotation(const double matrix[3][3], const double vec[3], double result[3]) +{ + result[0] = matrix[0][0] * vec[0] + matrix[1][0] * vec[1] + matrix[2][0] * vec[2]; + result[1] = matrix[0][1] * vec[0] + matrix[1][1] * vec[1] + matrix[2][1] * vec[2]; + result[2] = matrix[0][2] * vec[0] + matrix[1][2] * vec[1] + matrix[2][2] * vec[2]; +} + +int main() +{ + double O1[3] = {0.00, 0.00, 0.00}; + double H1[3] = {-0.45, -1.48, -0.97}; + double H2[3] = {-0.45, 1.42, -1.07}; + + double 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_temp[3], H2_temp[3]; + double H1_rotated[3], H2_rotated[3]; + + 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 out_of_plane_rotation_matrix[3][3]; + rotation_matrix(axis_to_align, angle_to_align, out_of_plane_rotation_matrix); + apply_rotation(out_of_plane_rotation_matrix, H1, H1_temp); + apply_rotation(out_of_plane_rotation_matrix, H2, H2_temp); + + double theta1 = atan2(H1_temp[1], H1_temp[0]); + double planar_rotation_angle = -theta1 + bondangle / 360.0 * M_PI; + double planar_rotation_matrix[3][3]; + rotation_matrix(target_plane_normal, planar_rotation_angle, planar_rotation_matrix); + apply_rotation(planar_rotation_matrix, H1_temp, H1_rotated); + apply_rotation(planar_rotation_matrix, H2_temp, H2_rotated); + + if (bondlength1 < bondlength2) + { + swap(H1_rotated, H2_rotated); + H1_rotated[1] *= -1.0; + H2_rotated[1] *= -1.0; + } + + double bondlength1_rotated = calculate_bondlength(H1_rotated, O1); + double bondlength2_rotated = calculate_bondlength(H2_rotated, O1); + double bondangle_rotated = calculate_bondangle(H1_rotated, O1, H2_rotated, false); + + cout << "Reference system (in z=0 plane, symmetric about x=0 axis, with longer bondlength in Q1)" << 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_rotated << endl; + cout << "Bondlength of O1-H2 = " << bondlength2_rotated << endl; + cout << "Angle between O1-H1 and O1-H2 = " << bondangle_rotated << endl; + + double H1_restored[3], H2_restored[3]; + + if (bondlength1 < bondlength2) + { + swap(H1_rotated, H2_rotated); + H1_rotated[1] *= -1.0; + H2_rotated[1] *= -1.0; + } + + apply_transpose_rotation(planar_rotation_matrix, H1_rotated, H1_temp); + apply_transpose_rotation(planar_rotation_matrix, H2_rotated, H2_temp); + + apply_transpose_rotation(out_of_plane_rotation_matrix, H1_temp, H1_restored); + apply_transpose_rotation(out_of_plane_rotation_matrix, H2_temp, H2_restored); + + double bondlength1_restored = calculate_bondlength(H1_restored, O1); + double bondlength2_restored = calculate_bondlength(H2_restored, O1); + double bondangle_restored = 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_restored << endl; + cout << "Bondlength of O1-H2 = " << bondlength2_restored << endl; + cout << "Angle between O1-H1 and O1-H2 = " << bondangle_restored << endl; + + return 0; +} diff --git a/tests/PinnedH2O_3DOF/aux/rotation_test.py b/tests/PinnedH2O_3DOF/aux/rotation_test.py new file mode 100644 index 00000000..41e20296 --- /dev/null +++ b/tests/PinnedH2O_3DOF/aux/rotation_test.py @@ -0,0 +1,99 @@ +import numpy as np + +O1 = np.array([0.00, 0.00, 0.00]) +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) + +def calculate_bondangle(atom1, atom2, atom3, radian): + vector1 = atom1 - atom2 + vector2 = atom3 - atom2 + dot_product = np.dot(vector1, vector2) + magnitude_product = np.linalg.norm(vector1) * np.linalg.norm(vector2) + angle = np.arccos(dot_product / magnitude_product) + if not radian: + angle = np.degrees(angle) + return angle + +bondlength1 = calculate_bondlength(H1, O1) +bondlength2 = calculate_bondlength(H2, O1) +bondangle = calculate_bondangle(H1, O1, H2, False) + +print('Original system') +print(f'H1 = {H1}') +print(f'H2 = {H2}') +print(f'Bondlength of O1-H1 = {bondlength1}') +print(f'Bondlength of O1-H2 = {bondlength2}') +print(f'Angle between O1-H1 and O1-H2 = {bondangle}') + +def rotation_matrix(axis, angle): + cos_theta = np.cos(angle) + sin_theta = np.sin(angle) + ux, uy, uz = axis + return np.array([ + [cos_theta + ux**2 * (1 - cos_theta), ux * uy * (1 - cos_theta) - uz * sin_theta, ux * uz * (1 - cos_theta) + uy * sin_theta], + [uy * ux * (1 - cos_theta) + uz * sin_theta, cos_theta + uy**2 * (1 - cos_theta), uy * uz * (1 - cos_theta) - ux * sin_theta], + [uz * ux * (1 - cos_theta) - uy * sin_theta, uz * uy * (1 - cos_theta) + ux * sin_theta, cos_theta + uz**2 * (1 - cos_theta)] + ]) + +H1_rotated, H2_rotated = H1, H2 + +plane_normal = np.cross(H2, H1) +plane_normal = plane_normal / np.linalg.norm(plane_normal) + +target_plane_normal = np.array([0, 0, 1]) +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)) +out_of_plane_rotation_matrix = rotation_matrix(axis_to_align, angle_to_align) +H1_rotated = np.dot(out_of_plane_rotation_matrix, H1) +H2_rotated = np.dot(out_of_plane_rotation_matrix, H2) + +theta1 = np.arctan2(H1_rotated[1], H1_rotated[0]) +theta2 = np.arctan2(H2_rotated[1], H2_rotated[0]) +planar_rotation_angle = -theta1 + np.radians(bondangle) / 2 +planar_rotation_matrix = rotation_matrix(target_plane_normal, planar_rotation_angle) +H1_rotated = np.dot(planar_rotation_matrix, H1_rotated) +H2_rotated = np.dot(planar_rotation_matrix, H2_rotated) + +if bondlength1 < bondlength2: + H1_rotated, H2_rotated = H2_rotated, H1_rotated + H1_rotated[1] *= -1 + H2_rotated[1] *= -1 + +bondlength1_rotated = calculate_bondlength(H1_rotated, O1) +bondlength2_rotated = calculate_bondlength(H2_rotated, O1) +bondangle_rotated = calculate_bondangle(H1_rotated, O1, H2_rotated, False) + +print('Reference system (z=0 plane, symmetric about x=0 axis, with longer bondlength in Q1)') +print(f'H1 = {H1_rotated}') +print(f'H2 = {H2_rotated}') +print(f'Bondlength of O1-H1 = {bondlength1_rotated}') +print(f'Bondlength of O1-H2 = {bondlength2_rotated}') +print(f'Angle between O1-H1 and O1-H2 = {bondangle_rotated}') + +H1_restored, H2_restored = H1_rotated, H2_rotated + +if bondlength1 < bondlength2: + H1_rotated, H2_rotated = H2_rotated, H1_rotated + H1_rotated[1] *= -1 + H2_rotated[1] *= -1 + +H1_restored = np.dot(planar_rotation_matrix.T, H1_restored) +H2_restored = np.dot(planar_rotation_matrix.T, H2_restored) + +H1_restored = np.dot(out_of_plane_rotation_matrix.T, H1_restored) +H2_restored = np.dot(out_of_plane_rotation_matrix.T, H2_restored) + +bondlength1_restored = calculate_bondlength(H1_restored, O1) +bondlength2_restored = calculate_bondlength(H2_restored, O1) +bondangle_restored = calculate_bondangle(H1_restored, O1, H2_restored, False) + +print('Restored system') +print(f'H1 = {H1_restored}') +print(f'H2 = {H2_restored}') +print(f'Bondlength of O1-H1 = {bondlength1_restored}') +print(f'Bondlength of O1-H2 = {bondlength2_restored}') +print(f'Angle between O1-H1 and O1-H2 = {bondangle_restored}') diff --git a/tests/PinnedH2O_3DOF/compare_energy_and_forces.py b/tests/PinnedH2O_3DOF/compare_energy_and_forces.py new file mode 100644 index 00000000..b11dbe91 --- /dev/null +++ b/tests/PinnedH2O_3DOF/compare_energy_and_forces.py @@ -0,0 +1,233 @@ +import numpy as np +import os +import matplotlib.pyplot as plt +import matplotlib.ticker as ticker +import argparse + +parser = argparse.ArgumentParser(description="Calculate and plot differences in energies and forces.") + +parser.add_argument("--N_l", type=int, default=2, help="Sampling frequency of bond length points") +parser.add_argument("--N_theta", type=int, default=2, help="Sampling frequency of bond angle points") +parser.add_argument("--rdim", type=int, default=18, help="Dimension of the ROM basis") + +args = parser.parse_args() + +N_l = args.N_l +N_theta = args.N_theta +rdim = args.rdim + +bondlength_min = 0.95 +bondlength_max = 1.05 +bondlength_num_increments = 10 + +bondangle_min = -5.0 +bondangle_max = 5.0 +bondangle_num_increments = 10 + +output_dir = f'results_{N_l}_{N_theta}_{rdim}' + +def process_data(s_l1, s_l2, s_theta, N_l, N_theta, rdim): + try: + data_dir = f'data/{s_l1}_{s_l2}_{s_theta}' + offline_file = os.path.join(data_dir, 'offline_PinnedH2O.out') + log_file = f'{output_dir}/{s_l1}_{s_l2}_{s_theta}.log' + + f_O1_fom = None + f_H1_fom = None + f_H2_fom = None + + if not os.path.exists(offline_file): + print(f"Error: File not found: {offline_file}") + return None + + with open(offline_file, 'r') as f: + eks_lines = [line for line in f if " SC ENERGY [Ha]" in line] + if eks_lines: + Eks_fom = float(eks_lines[-1].split()[-1]) + else: + print(f"Error: Could not find 'SC ENERGY [Ha]' line in {offline_file}") + return None + + with open(offline_file, 'r') as f: + for line in f: + if "O1" in line: + f_O1_fom = np.array(line.split()[6:9], dtype=float) + elif "H1" in line: + f_H1_fom = np.array(line.split()[5:8], dtype=float) + elif "H2" in line: + f_H2_fom = np.array(line.split()[5:8], dtype=float) + + if f_O1_fom is not None and f_H1_fom is not None and f_H2_fom is not None: + break # Exit the loop once we found all lines + + if f_O1_fom is None or f_H1_fom is None or f_H2_fom is None: + print(f"Error: Could not find O1, H1, or H2 lines in {offline_file}") + return None + + f_H2_rom = None + f_O1_rom = None + f_H1_rom = None + + if not os.path.exists(log_file): + print(f"Error: File not found: {log_file}") + return None + + with open(log_file, 'r') as f: + lines = f.readlines() + + eks_index = -1 + for i in range(len(lines) - 1, -1, -1): # Search from the end + if "Eks:" in lines[i]: + eks_index = i + break + + if eks_index == -1: + print(f"Error: Could not find 'Eks:' line in {log_file}") + return None + + Eks_rom = float(lines[eks_index].split()[-1]) # Get the last element + + force_index = -1 + for i in range(len(lines) - 1, -1, -1): # Search from the end + if "Forces:" in lines[i]: + force_index = i + break + + if force_index == -1: + print(f"Error: Could not find 'Forces:' line in {log_file}") + return None + + if force_index + 3 >= len(lines): + print(f"Error: Not enough lines after 'Forces:' in {log_file}") + return None + + f_H2_rom = np.array(lines[force_index + 1].split(), dtype=float) + f_O1_rom = np.array(lines[force_index + 2].split(), dtype=float) + f_H1_rom = np.array(lines[force_index + 3].split(), dtype=float) + + return Eks_fom, f_O1_fom, f_H1_fom, f_H2_fom, Eks_rom, f_O1_rom, f_H1_rom, f_H2_rom + + except Exception as e: + print(f"An error occurred: {e}") + return None + +Eks_diff_reproductive = [] +f_O1_diff_reproductive = [] +f_H1_diff_reproductive = [] +f_H2_diff_reproductive = [] + +Eks_diff_predictive = [] +f_O1_diff_predictive = [] +f_H1_diff_predictive = [] +f_H2_diff_predictive = [] + +for i in range(bondlength_num_increments + 1): + bondlength_one = round(bondlength_min + i * (bondlength_max - bondlength_min) / bondlength_num_increments, 2) + for j in range(i + 1): + bondlength_two = round(bondlength_min + j * (bondlength_max - bondlength_min) / bondlength_num_increments, 2) + for k in range(bondangle_num_increments + 1): + bondangle = round(bondangle_min + k * (bondangle_max - bondangle_min) / bondangle_num_increments, 1) + + s_l1 = f"{bondlength_one:.2f}" + s_l2 = f"{bondlength_two:.2f}" + s_theta = f"{bondangle:.1f}" + tag = f'{s_l1}_{s_l2}_{s_theta}' + + output_filename = f'{output_dir}/energy_and_forces_{tag}.out' + + with open(output_filename, 'w') as outfile: + print("s_l1:", s_l1, file=outfile) + print("s_l2:", s_l2, file=outfile) + print("s_theta:", s_theta, file=outfile) + print("N_l:", N_l, file=outfile) + print("N_theta:", N_theta, file=outfile) + print("rdim:", rdim, file=outfile) + + results = process_data(s_l1, s_l2, s_theta, N_l, N_theta, rdim) + + if results: + Eks_fom, f_O1_fom, f_H1_fom, f_H2_fom, Eks_rom, f_O1_rom, f_H1_rom, f_H2_rom = results + + print("Eks_fom:", Eks_fom, file=outfile) + print("Eks_rom:", Eks_rom, file=outfile) + + print("f_O1_fom:", f_O1_fom, file=outfile) + print("f_O1_rom:", f_O1_rom, file=outfile) + + print("f_H1_fom:", f_H1_fom, file=outfile) + print("f_H1_rom:", f_H1_rom, file=outfile) + + print("f_H2_fom:", f_H2_fom, file=outfile) + print("f_H2_rom:", f_H2_rom, file=outfile) + + def calculate_differences(fom, rom, name): + abs_diff = np.linalg.norm(fom - rom) + rel_diff = abs_diff / np.linalg.norm(fom) if np.linalg.norm(fom)!= 0 else float('inf') + print(f"Absolute difference in {name}:", abs_diff, file=outfile) + print(f"Relative difference in {name}:", rel_diff, file=outfile) + return abs_diff + + Eks_diff = calculate_differences(Eks_fom, Eks_rom, "Eks") + f_O1_diff = calculate_differences(f_O1_fom, f_O1_rom, "f_O1") + f_H1_diff = calculate_differences(f_H1_fom, f_H1_rom, "f_H1") + f_H2_diff = calculate_differences(f_H2_fom, f_H2_rom, "f_H2") + + if i * N_l % bondlength_num_increments == 0 and j * N_l % bondlength_num_increments == 0 and k * N_theta % bondangle_num_increments == 0: + Eks_diff_reproductive.append(Eks_diff) + f_O1_diff_reproductive.append(f_O1_diff) + f_H1_diff_reproductive.append(f_H1_diff) + f_H2_diff_reproductive.append(f_H2_diff) + else: + Eks_diff_predictive.append(Eks_diff) + f_O1_diff_predictive.append(f_O1_diff) + f_H1_diff_predictive.append(f_H1_diff) + f_H2_diff_predictive.append(f_H2_diff) + else: + print("Error occurred during data processing. Differences cannot be calculated.", file=outfile) + +def plot_histogram(data, quantity, test_case): + plt.figure(figsize=(8, 6)) + plt.hist(data, bins=20, color='skyblue', edgecolor='black') + + if quantity == "Eks": + quantity_name = "absolute difference in total energy" + elif quantity.startswith("f_"): + quantity_name = f"magnitude of difference in force on {quantity[2:]}" + else: + raise ValueError("Invalid input quantity") + + plt.title(f'Histogram of {quantity_name}') + plt.xlabel('Difference') + plt.ylabel('Frequency') + + min_val, max_val = np.min(data), np.max(data) + plt.xlim(min_val, max_val) + num_ticks = 8 + xticks = np.linspace(min_val, max_val, num_ticks) + plt.xticks(xticks) + formatter = ticker.ScalarFormatter(useMathText=True) + formatter.set_scientific(True) + formatter.set_powerlimits((0, 0)) + plt.gca().xaxis.set_major_formatter(formatter) + + total_count = len(data) + mean_val = np.mean(data) + max_val = np.max(data) + stats_text = (f"Total {test_case} cases: {total_count}\n" + f"Mean: {mean_val:.3e}") + plt.text(0.95, 0.95, stats_text, transform=plt.gca().transAxes, + fontsize=12, verticalalignment='top', horizontalalignment='right', + bbox=dict(facecolor='white', alpha=0.7, edgecolor='black')) + + plt.tight_layout() + plt.savefig(f"{output_dir}/{quantity}_difference_histogram_{test_case}.png") + +plot_histogram(Eks_diff_reproductive, "Eks", "reproductive") +plot_histogram(f_O1_diff_reproductive, "f_O1", "reproductive") +plot_histogram(f_H1_diff_reproductive, "f_H1", "reproductive") +plot_histogram(f_H2_diff_reproductive, "f_H2", "reproductive") + +plot_histogram(Eks_diff_predictive, "Eks", "predictive") +plot_histogram(f_O1_diff_predictive, "f_O1", "predictive") +plot_histogram(f_H1_diff_predictive, "f_H1", "predictive") +plot_histogram(f_H2_diff_predictive, "f_H2", "predictive") diff --git a/tests/PinnedH2O_3DOF/coords_rotate1.in b/tests/PinnedH2O_3DOF/coords_rotate1.in new file mode 100644 index 00000000..7498ff55 --- /dev/null +++ b/tests/PinnedH2O_3DOF/coords_rotate1.in @@ -0,0 +1,3 @@ +O1 1 0.00000 0.00000 0.00 0 +H1 2 -1.08248 -1.47039 0.00 1 +H2 2 -1.14768 1.43060 0.00 1 diff --git a/tests/PinnedH2O_3DOF/coords_rotate2.in b/tests/PinnedH2O_3DOF/coords_rotate2.in new file mode 100644 index 00000000..307e170c --- /dev/null +++ b/tests/PinnedH2O_3DOF/coords_rotate2.in @@ -0,0 +1,3 @@ +O1 1 0.00 0.00 0.00 0 +H1 2 -0.45 1.42 -1.07 1 +H2 2 -0.45 -1.48 -0.97 1 diff --git a/tests/PinnedH2O_3DOF/generate_coord.py b/tests/PinnedH2O_3DOF/generate_coord.py new file mode 100644 index 00000000..194c2df7 --- /dev/null +++ b/tests/PinnedH2O_3DOF/generate_coord.py @@ -0,0 +1,31 @@ +import numpy as np +import os + +O1 = np.array([0.00, 0.00, 0.00]) + +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, 11) +bondangle_increment = np.linspace(-5, 5, 11) + +# output directory +output_dir = "data" + +# generation +os.makedirs(output_dir, exist_ok=True) + +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): + for f_bondlength2 in bondlength_factor[:(i+1)]: + H1 = np.array([f_bondlength1*x, f_bondlength1*y, 0.0]) + H2 = np.array([f_bondlength2*x, -f_bondlength2*y, 0.0]) + filename = f"{output_dir}/coords_{f_bondlength1:.2f}_{f_bondlength2:.2f}_{d_bondangle}.in" + with open(filename, "w") as file: + file.write(f"O1 1 {O1[0]:.2f} {O1[1]:.2f} {O1[2]:.2f} 0\n") + file.write(f"H1 2 {H1[0]:.2f} {H1[1]:.2f} {H1[2]:.2f} 1\n") + file.write(f"H2 2 {H2[0]:.2f} {H2[1]:.2f} {H2[2]:.2f} 1\n") diff --git a/tests/PinnedH2O_3DOF/job.basis b/tests/PinnedH2O_3DOF/job.basis new file mode 100644 index 00000000..35f564a3 --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.basis @@ -0,0 +1,47 @@ +#!/bin/tcsh +#SBATCH -N 1 +#SBATCH -t 0:10:00 +#SBATCH -p pdebug + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 8 + +set maindir = /p/lustre2/cheung26/mgmol + +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 = 2 +set bondangle_min = -5.0 +set bondangle_max = 5.0 +set bondangle_num_increments = 2 + +foreach i (`seq 0 $bondlength_num_increments`) + set bondlength_one = `echo "scale=2; $bondlength_min + $i * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_one = `printf "%.2f" $bondlength_one` + foreach j (`seq 0 $i`) + set bondlength_two = `echo "scale=2; $bondlength_min + $j * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_two = `printf "%.2f" $bondlength_two` + foreach k (`seq 0 $bondangle_num_increments`) + set bondangle = `echo "scale=1; $bondangle_min + $k * ($bondangle_max - $bondangle_min) / ($bondangle_num_increments)" | bc` + set bondangle = `printf "%.1f" $bondangle` + echo "bondlength1: $bondlength_one, bondlength2: $bondlength_two, bondangle: $bondangle" + set tag = ${bondlength_one}_${bondlength_two}_${bondangle} + set snapshot_files = "$snapshot_files data/${tag}/orbital_snapshot" + end + end +end +echo "Snapshot files: $snapshot_files" + +set basis_file = "data/PinnedH2O_3DOF_orbitals_basis_${bondlength_num_increments}_${bondangle_num_increments}" +srun -n $ncpus $exe -f $basis_file $snapshot_files > data/basis_${bondlength_num_increments}_${bondangle_num_increments}_PinnedH2O_3DOF.out + +date diff --git a/tests/PinnedH2O_3DOF/job.offline b/tests/PinnedH2O_3DOF/job.offline new file mode 100644 index 00000000..5b4224f0 --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.offline @@ -0,0 +1,59 @@ +#!/bin/tcsh +#SBATCH -N 1 +#SBATCH -t 4:00:00 +#SBATCH -p pbatch + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 8 + +set maindir = /p/lustre2/cheung26/mgmol + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = mgmol-opt + +cp $maindir/install_quartz/bin/$exe . + +set cfg_offline = mgmol_offline.cfg + +ln -s -f $maindir/potentials/pseudo.O_ONCV_PBE_SG15 . +ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 . + +source $maindir/scripts/modules.quartz + +set bondlength_min = 0.95 +set bondlength_max = 1.05 +set bondlength_num_increments = 10 +set bondangle_min = -5.0 +set bondangle_max = 5.0 +set bondangle_num_increments = 10 + +foreach i (`seq 0 $bondlength_num_increments`) + set bondlength_one = `echo "scale=2; $bondlength_min + $i * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_one = `printf "%.2f" $bondlength_one` + foreach j (`seq 0 $i`) + set bondlength_two = `echo "scale=2; $bondlength_min + $j * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_two = `printf "%.2f" $bondlength_two` + foreach k (`seq 0 $bondangle_num_increments`) + set bondangle = `echo "scale=1; $bondangle_min + $k * ($bondangle_max - $bondangle_min) / ($bondangle_num_increments)" | bc` + set bondangle = `printf "%.1f" $bondangle` + echo "bondlength1: $bondlength_one, bondlength2: $bondlength_two, bondangle: $bondangle" + set tag = ${bondlength_one}_${bondlength_two}_${bondangle} + if (-e data/${tag}) then + continue + else + mkdir data/${tag} + endif + cp $cfg_offline data/${tag}/${cfg_offline} + mv data/coords_${tag}.in data/${tag}/coords.in + sed -i "s/CUSTOM_TAG/${tag}/g" data/${tag}/${cfg_offline} + srun -n $ncpus $exe -c data/${tag}/${cfg_offline} -i data/${tag}/coords.in > data/${tag}/offline_PinnedH2O.out + end + end +end + +date diff --git a/tests/PinnedH2O_3DOF/job.online b/tests/PinnedH2O_3DOF/job.online new file mode 100644 index 00000000..ea75e59f --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.online @@ -0,0 +1,67 @@ +#!/bin/tcsh +#SBATCH -N 1 +#SBATCH -t 1:00:00 +#SBATCH -p pdebug + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 8 + +set maindir = /p/lustre2/cheung26/mgmol + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = testPinnedH2O_3DOF + +cp $maindir/build_quartz/tests/$exe . + +set cfg_online = mgmol_online.cfg + +ln -s -f $maindir/potentials/pseudo.O_ONCV_PBE_SG15 . +ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 . + +source $maindir/scripts/modules.quartz + +set sample_bondlength_num_increments = 2 +set sample_bondangle_num_increments = 2 +set num_orbital_basis = 34 + +set bondlength_min = 0.95 +set bondlength_max = 1.05 +set bondlength_num_increments = 10 +set bondangle_min = -5.0 +set bondangle_max = 5.0 +set bondangle_num_increments = 10 + +set case = ${sample_bondlength_num_increments}_${sample_bondangle_num_increments}_${num_orbital_basis} +mkdir results_${case} +cp ${cfg_online} ${case}_${cfg_online} + +@ num_empty_orbital = $num_orbital_basis - 4 +sed -i "s/CUSTOM_NE/${num_empty_orbital}/g" ${case}_${cfg_online} +sed -i "s/CUSTOM_NL/${sample_bondlength_num_increments}/g" ${case}_${cfg_online} +sed -i "s/CUSTOM_NT/${sample_bondangle_num_increments}/g" ${case}_${cfg_online} +sed -i "s/CUSTOM_NB/${num_orbital_basis}/g" ${case}_${cfg_online} + +foreach i (`seq 0 $bondlength_num_increments`) + set bondlength_one = `echo "scale=2; $bondlength_min + $i * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_one = `printf "%.2f" $bondlength_one` + foreach j (`seq 0 $i`) + set bondlength_two = `echo "scale=2; $bondlength_min + $j * ($bondlength_max - $bondlength_min) / ($bondlength_num_increments)" | bc` + set bondlength_two = `printf "%.2f" $bondlength_two` + foreach k (`seq 0 $bondangle_num_increments`) + set bondangle = `echo "scale=1; $bondangle_min + $k * ($bondangle_max - $bondangle_min) / ($bondangle_num_increments)" | bc` + set bondangle = `printf "%.1f" $bondangle` + echo "bondlength1: $bondlength_one, bondlength2: $bondlength_two, bondangle: $bondangle" + set tag = ${bondlength_one}_${bondlength_two}_${bondangle} + srun -p pdebug -n $ncpus $exe -c ${case}_${cfg_online} -i /usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_${ncpus}/${tag}/coords.in >& results_${case}/${tag}.log + end + end +end + +mv ${case}_${cfg_online} results_${case} + +date diff --git a/tests/PinnedH2O_3DOF/job.online_rotate b/tests/PinnedH2O_3DOF/job.online_rotate new file mode 100644 index 00000000..a2065b72 --- /dev/null +++ b/tests/PinnedH2O_3DOF/job.online_rotate @@ -0,0 +1,36 @@ +#!/bin/tcsh +#SBATCH -N 1 +#SBATCH -t 0:10:00 +#SBATCH -p pdebug + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 8 + +set maindir = /p/lustre2/cheung26/mgmol + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +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 + +rm -rf snapshot0_* + +set coord_file = coords_rotate2.in + +set exe = mgmol-opt +set cfg_online = mgmol_ref_rotate.cfg +cp $maindir/install_quartz/bin/$exe . +srun -p pdebug -n $ncpus $exe -c $cfg_online -i $coord_file + +set exe = testPinnedH2O_3DOF +set cfg_online = mgmol_online_rotate.cfg +cp $maindir/build_quartz/tests/$exe . +srun -p pdebug -n $ncpus $exe -c $cfg_online -i $coord_file + +date diff --git a/tests/PinnedH2O_3DOF/mgmol_offline.cfg b/tests/PinnedH2O_3DOF/mgmol_offline.cfg new file mode 100644 index 00000000..7b2dbcfa --- /dev/null +++ b/tests/PinnedH2O_3DOF/mgmol_offline.cfg @@ -0,0 +1,30 @@ +verbosity=1 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=quench +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +[Restart] +output_level=4 +output_filename=CUSTOM_TAG +[ROM.offline] +save_librom_snapshot=true diff --git a/tests/PinnedH2O_3DOF/mgmol_online.cfg b/tests/PinnedH2O_3DOF/mgmol_online.cfg new file mode 100644 index 00000000..6f66067d --- /dev/null +++ b/tests/PinnedH2O_3DOF/mgmol_online.cfg @@ -0,0 +1,36 @@ +verbosity=2 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=quench +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +nempty=CUSTOM_NE +[Restart] +output_level=4 +[DensityMatrix] +solver=MVP +nb_inner_it=100 +[ROM.offline] +basis_file=/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_8/PinnedH2O_3DOF_orbitals_basis_CUSTOM_NL_CUSTOM_NT +[ROM.basis] +compare_md=false +number_of_orbital_basis=CUSTOM_NB diff --git a/tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg b/tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg new file mode 100644 index 00000000..8c598248 --- /dev/null +++ b/tests/PinnedH2O_3DOF/mgmol_online_rotate.cfg @@ -0,0 +1,36 @@ +verbosity=2 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=quench +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +nempty=30 +[Restart] +output_level=4 +[DensityMatrix] +solver=MVP +nb_inner_it=100 +[ROM.offline] +basis_file=/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_8/PinnedH2O_3DOF_orbitals_basis_2_2 +[ROM.basis] +compare_md=false +number_of_orbital_basis=34 diff --git a/tests/PinnedH2O_3DOF/mgmol_ref_rotate.cfg b/tests/PinnedH2O_3DOF/mgmol_ref_rotate.cfg new file mode 100644 index 00000000..96be6b8d --- /dev/null +++ b/tests/PinnedH2O_3DOF/mgmol_ref_rotate.cfg @@ -0,0 +1,27 @@ +verbosity=1 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=quench +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +[Restart] +output_level=4 diff --git a/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc new file mode 100644 index 00000000..ad9ef743 --- /dev/null +++ b/tests/PinnedH2O_3DOF/testPinnedH2O_3DOF.cc @@ -0,0 +1,256 @@ +// 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" +#include "PinnedH2O.h" + +#ifdef MGMOL_HAS_LIBROM +#include "librom.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +namespace po = boost::program_options; + +int main(int argc, char** argv) +{ + int mpirc = MPI_Init(&argc, &argv); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Initialization failed!!!" << std::endl; + MPI_Abort(MPI_COMM_WORLD, 0); + } + + MPI_Comm comm = MPI_COMM_WORLD; + + /* + * Initialize general things, like magma, openmp, IO, ... + */ + mgmol_init(comm); + + /* + * read runtime parameters + */ + std::string input_filename(""); + std::string lrs_filename; + std::string constraints_filename(""); + + float total_spin = 0.; + bool with_spin = false; + + po::variables_map vm; + + // read from PE0 only + if (MPIdata::onpe0) + { + read_config(argc, argv, vm, input_filename, lrs_filename, + constraints_filename, total_spin, with_spin); + } + + MGmol_MPI::setup(comm, std::cout, with_spin); + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); + MPI_Comm global_comm = mmpi.commGlobal(); + + /* + * Setup control struct with run time parameters + */ + Control::setup(global_comm, with_spin, total_spin); + Control& ct = *(Control::instance()); + + ct.setOptions(vm); + + int ret = ct.checkOptions(); + if (ret < 0) return ret; + + mmpi.bcastGlobal(input_filename); + mmpi.bcastGlobal(lrs_filename); + + // Enter main scope + { + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "Construct MGmol object..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + + MGmolInterface* mgmol = new MGmol>(global_comm, + *MPIdata::sout, input_filename, lrs_filename, constraints_filename); + + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "MGmol setup..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + mgmol->setup(); + + if (MPIdata::onpe0) + { + std::cout << "-------------------------" << std::endl; + std::cout << "Setup done..." << std::endl; + std::cout << "-------------------------" << std::endl; + } + + // here we just use the atomic positions read in and used + // to initialize MGmol + std::vector positions; + mgmol->getAtomicPositions(positions); + std::vector anumbers; + mgmol->getAtomicNumbers(anumbers); + if (MPIdata::onpe0) + { + std::cout << "Positions:" << std::endl; + std::vector::iterator ita = anumbers.begin(); + for (std::vector::iterator it = positions.begin(); + it != positions.end(); it += 3) + { + std::cout << *ita; + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + ita++; + } + } + + // rotate the molecule to the reference coordinate system + PinnedH2O H2O_molecule; + H2O_molecule.rotate(positions, anumbers); + + Mesh* mymesh = Mesh::instance(); + const pb::Grid& mygrid = mymesh->grid(); + const pb::PEenv& myPEenv = mymesh->peenv(); + + // compute energy and forces again with projected problem onto ROM subspace + const int rdim = ct.getROMOptions().num_orbbasis; + if (rdim != ct.numst) + { + std::cerr << "The number of functions in the ROM basis file, " + << rdim << " is not equal to ct.numst, " << ct.numst + << std::endl; + MPI_Abort(mmpi.commSameSpin(), 0); + } + + std::shared_ptr projmatrices + = mgmol->getProjectedMatrices(); + + ExtendedGridOrbitals orbitals("new_orbitals", mygrid, mymesh->subdivx(), + ct.numst, ct.bcWF, projmatrices.get(), nullptr, nullptr, nullptr, + nullptr); + + orbitals.set(ct.getROMOptions().basis_file, ct.numst); + orbitals.orthonormalizeLoewdin(); + orbitals.setDataWithGhosts(true); + + // set the iterative index to 1 to differentiate it from first instance + // in MGmol initial() function. This is not very clean and could be + // better designed, but works for now + orbitals.setIterativeIndex(10); + + // set initial DM with uniform occupations + projmatrices->setDMuniform(ct.getNelSpin()); + projmatrices->printDM(std::cout); + + // + // evaluate energy and forces with ROM bases just read + // + std::vector forces; + double eks = mgmol->evaluateDMandEnergyAndForces( + &orbitals, positions, anumbers, forces); + + // print out results + if (MPIdata::onpe0) + { + std::cout << "Eks: " << eks << std::endl; + std::cout << "Positions in the reference domain:" << std::endl; + std::vector::iterator ita = anumbers.begin(); + for (std::vector::iterator it = positions.begin(); + it != positions.end(); it += 3) + { + std::cout << *ita; + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + ita++; + } + + std::cout << "Forces in the reference domain:" << std::endl; + ita = anumbers.begin(); + for (std::vector::iterator it = forces.begin(); + it != forces.end(); it += 3) + { + std::cout << *ita; + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + ita++; + } + } + + // rotate the forces to the original coordinate system + H2O_molecule.transpose_rotate(positions, anumbers, forces); + + // print out results + if (MPIdata::onpe0) + { + std::cout << "Positions:" << std::endl; + std::vector::iterator ita = anumbers.begin(); + for (std::vector::iterator it = positions.begin(); + it != positions.end(); it += 3) + { + std::cout << *ita; + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + ita++; + } + std::cout << "Forces:" << std::endl; + ita = anumbers.begin(); + for (std::vector::iterator it = forces.begin(); + it != forces.end(); it += 3) + { + std::cout << *ita; + for (int i = 0; i < 3; i++) + std::cout << " " << *(it + i); + std::cout << std::endl; + ita++; + } + } + + delete mgmol; + + } // close main scope + + mgmol_finalize(); + + mpirc = MPI_Finalize(); + if (mpirc != MPI_SUCCESS) + { + std::cerr << "MPI Finalize failed!!!" << std::endl; + } + + time_t tt; + time(&tt); + if (onpe0) std::cout << " Run ended at " << ctime(&tt) << std::endl; + + return 0; +} +#endif // MGMOL_HAS_LIBROM diff --git a/tests/PinnedH2O_3DOF/transfer_nlrom.sh b/tests/PinnedH2O_3DOF/transfer_nlrom.sh new file mode 100644 index 00000000..62824de2 --- /dev/null +++ b/tests/PinnedH2O_3DOF/transfer_nlrom.sh @@ -0,0 +1,20 @@ +#!/bin/bash + +ncpus=8 +target_dir="/usr/workspace/nlrom/MGmol/PinnedH2O_3DOF/data_${ncpus}" + +if [ ! -d $target_dir ]; then + mkdir -p $target_dir +fi + +for dir in data/*/; do + dirname=$(basename "$dir") + new_dir="$target_dir/${dirname}" + if [ ! -d $new_dir ]; then + mkdir -p $new_dir + fi + mv -f $dir/* $new_dir/ + rm -rf $dir +done + +rm -rf data diff --git a/tests/ROM/test_rom_poisson/carbyne.in b/tests/ROM/test_rom_poisson/carbyne.in new file mode 100644 index 00000000..9f4e975a --- /dev/null +++ b/tests/ROM/test_rom_poisson/carbyne.in @@ -0,0 +1,14 @@ +H00 1 -0.0000 -0.0000 15.2674 +C01 2 -0.0000 0.0000 13.2519 +C02 2 -0.0000 0.0000 10.9495 +C03 2 -0.0000 -0.0000 8.4221 +C04 2 0.0000 0.0000 6.0897 +C05 2 -0.0000 0.0000 3.5892 +C06 2 -0.0000 -0.0000 1.2470 +C07 2 0.0000 -0.0000 -1.2469 +C08 2 0.0000 -0.0000 -3.5891 +C09 2 -0.0000 -0.0000 -6.0897 +C10 2 -0.0000 0.0000 -8.4221 +C11 2 0.0000 -0.0000 -10.9495 +C12 2 0.0000 0.0000 -13.2520 +H13 1 0.0000 0.0000 -15.2675 diff --git a/tests/ROM/test_rom_poisson/carbyne.ion.cfg b/tests/ROM/test_rom_poisson/carbyne.ion.cfg new file mode 100644 index 00000000..d98232ad --- /dev/null +++ b/tests/ROM/test_rom_poisson/carbyne.ion.cfg @@ -0,0 +1,63 @@ +verbosity=2 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx= 96 +ny= 96 +nz= 192 +[Domain] +ox= -10. +oy= -10. +oz= -20. +lx= 20. +ly= 20. +lz= 40. +[Poisson] +FDtype=4th +#max_steps_initial=99 +#max_steps=99 +[Potentials] +pseudopotential=pseudo.H_ONCV_PBE_SG15 +pseudopotential=pseudo.C_ONCV_PBE_SG15 +[Run] +#type=QUENCH +type=MD +[Quench] +#solver=PR +max_steps=300 +atol=1.e-8 +[MD] +num_steps=3000 +dt=40. +print_interval=1 +[Orbitals] +initial_type=Fourier +nempty=10 +temperature=300. +[ProjectedMatrices] +solver=exact +[DensityMatrix] +solver=MVP +nb_inner_it=1 + +[Thermostat] +type=Langevin +temperature=300. +relax_time=1000. + +[Restart] +#input_filename=init_cond_144proc +#input_level=4 +output_level=4 +interval=1 + +[ROM] +stage=test_ion + +[ROM.offline] +restart_filefmt=case-300K/snapshot%05d +restart_min_idx=800 +restart_max_idx=1999 +basis_file=basis_300K_2/test_300K +variable=potential + diff --git a/tests/ROM/test_rom_poisson/carbyne.poisson.cfg b/tests/ROM/test_rom_poisson/carbyne.poisson.cfg new file mode 100644 index 00000000..6aac926c --- /dev/null +++ b/tests/ROM/test_rom_poisson/carbyne.poisson.cfg @@ -0,0 +1,63 @@ +verbosity=2 +xcFunctional=PBE +FDtype=4th +[Mesh] +nx= 96 +ny= 96 +nz= 192 +[Domain] +ox= -10. +oy= -10. +oz= -20. +lx= 20. +ly= 20. +lz= 40. +[Poisson] +FDtype=4th +#max_steps_initial=99 +#max_steps=99 +[Potentials] +pseudopotential=pseudo.H_ONCV_PBE_SG15 +pseudopotential=pseudo.C_ONCV_PBE_SG15 +[Run] +#type=QUENCH +type=MD +[Quench] +#solver=PR +max_steps=300 +atol=1.e-8 +[MD] +num_steps=3000 +dt=40. +print_interval=1 +[Orbitals] +initial_type=Fourier +nempty=10 +temperature=300. +[ProjectedMatrices] +solver=exact +[DensityMatrix] +solver=MVP +nb_inner_it=1 + +[Thermostat] +type=Langevin +temperature=300. +relax_time=1000. + +[Restart] +#input_filename=init_cond_144proc +#input_level=4 +output_level=4 +interval=1 + +[ROM] +stage=test_poisson + +[ROM.offline] +restart_filefmt=case-300K/snapshot%05d +restart_min_idx=800 +restart_max_idx=1999 +basis_file=basis_300K_2/test_300K +variable=potential + diff --git a/tests/RestartEnergyAndForces/job.test b/tests/RestartEnergyAndForces/job.test new file mode 100644 index 00000000..89b4865f --- /dev/null +++ b/tests/RestartEnergyAndForces/job.test @@ -0,0 +1,28 @@ +#!/bin/tcsh +#SBATCH -N 2 +#SBATCH -t 1:00:00 +#SBATCH -p pbatch + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 2 + +set maindir = /p/lustre2/cheung26/mgmol-20241220 + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = testRestartEnergyAndForces + +cp $maindir/build_quartz/tests/$exe . + +ln -s -f $maindir/potentials/pseudo.O_ONCV_PBE_SG15 . +ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 . + +source $maindir/scripts/modules.quartz + +python3 test.py srun -n $ncpus ${maindir}/install_quartz/bin/mgmol-opt ${exe} mgmol.cfg restart.cfg h2o.xyz . + +date diff --git a/tests/WFEnergyAndForces/job.test b/tests/WFEnergyAndForces/job.test new file mode 100644 index 00000000..7a9e49ee --- /dev/null +++ b/tests/WFEnergyAndForces/job.test @@ -0,0 +1,27 @@ +#!/bin/tcsh +#SBATCH -N 2 +#SBATCH -t 1:00:00 +#SBATCH -p pbatch + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 2 + +set maindir = /p/lustre2/cheung26/mgmol-20241220 + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = testWFEnergyAndForces + +cp $maindir/build_quartz/tests/$exe . + +set cfg_online = mgmol.cfg + +source $maindir/scripts/modules.quartz + +python3 test.py srun -p pdebug -n $ncpus $exe $cfg_online sih4.xyz $maindir/potentials + +date