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 753a455d..6c48a2e9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -126,6 +126,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) + set(MGMOL_HAS_LIBROM 1) + 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) @@ -249,6 +265,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/generate_coord.py b/examples/PinnedH2O/generate_coord.py new file mode 100644 index 00000000..c239eb87 --- /dev/null +++ b/examples/PinnedH2O/generate_coord.py @@ -0,0 +1,63 @@ +import numpy as np +import os + +# coords.in +O1 = np.array([0.00, 0.00, 0.00]) +ref_H1 = np.array([-0.45, 1.42, -1.07]) +ref_H2 = np.array([-0.45, -1.48, -0.97]) + +# factors and increments for bond lengths and bond angle +bondlength1_factor = np.linspace(0.95, 1.05, 11) +bondlength2_factor = np.linspace(0.95, 1.05, 11) +bondangle_increment = np.linspace(-5, 5, 11) + +# output directory +output_dir = "PinnedH2O_3dof_coords" + +# utilities +def calculate_bondlength(atom1, atom2): + return np.linalg.norm(atom1 - atom2) + +def calculate_bondangle(atom1, atom2, atom3): + vector1 = atom1 - atom2 + vector2 = atom3 - atom2 + dot_product = np.dot(vector1, vector2) + magnitude_product = np.linalg.norm(vector1) * np.linalg.norm(vector2) + angle_rad = np.arccos(dot_product / magnitude_product) + angle_deg = np.degrees(angle_rad) + return angle_deg + +# Rodrigues' rotation formula +def rotation_matrix(axis, angle_degrees): + angle = np.radians(angle_degrees) + 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)] + ]) + +# generation +os.makedirs(output_dir, exist_ok=True) + +ref_bondlength1 = calculate_bondlength(ref_H1, O1) +ref_bondlength2 = calculate_bondlength(ref_H2, O1) +ref_bondangle = calculate_bondangle(ref_H1, O1, ref_H2) + +normal_vector = np.cross(ref_H1, ref_H2) +normal_unit_vector = normal_vector / np.linalg.norm(normal_vector) + +for d_bondangle in bondangle_increment: + Q = rotation_matrix(normal_unit_vector, d_bondangle) + Q_ref_H2 = np.dot(Q, ref_H2) + for f_bondlength1 in bondlength1_factor: + for f_bondlength2 in bondlength2_factor: + H1 = f_bondlength1 * ref_H1 + H2 = f_bondlength2 * Q_ref_H2 + 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/examples/PinnedH2O/generate_coord_simple.py b/examples/PinnedH2O/generate_coord_simple.py new file mode 100644 index 00000000..7e799bdd --- /dev/null +++ b/examples/PinnedH2O/generate_coord_simple.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 = "PinnedH2O_3dof_coords" + +# 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/examples/PinnedH2O/get_result.sh b/examples/PinnedH2O/get_result.sh new file mode 100644 index 00000000..443e4cdd --- /dev/null +++ b/examples/PinnedH2O/get_result.sh @@ -0,0 +1,27 @@ +filename="offline_PinnedH2O.out" # FOM +#filename="rom39_PinnedH2O.out" # compare MD +#filename="39_force_PinnedH2O.out" # compare force + +# Extracting H1, H2, F1, F2 from MGmgol output log +# if FOM, these files contain the FOM results +# if compare MD, these files contain the results with projected orbitals +awk '/H1 / {print $3, $4, $5}' $filename > H1_$filename +awk '/H2 / {print $3, $4, $5}' $filename > H2_$filename +awk '/F1 / {print $6, $7, $8}' $filename > F1_$filename +awk '/F2 / {print $6, $7, $8}' $filename > F2_$filename + +# if compare force, files with "_fom" contain the FOM results +# files with "_rom" contain the results with projected orbitals +if [[ "$filename" == *"force_"* ]]; then + sed -n '1~2p' H1_$filename > H1_rom$filename + sed -n '1~2p' H2_$filename > H2_rom$filename + sed -n '1~2p' F1_$filename > F1_rom$filename + sed -n '1~2p' F2_$filename > F2_rom$filename + + sed -n '2~2p' H1_$filename > H1_fom$filename + sed -n '2~2p' H2_$filename > H2_fom$filename + sed -n '2~2p' F1_$filename > F1_fom$filename + sed -n '2~2p' F2_$filename > F2_fom$filename +fi + +rm -rf snapshot0_* diff --git a/examples/PinnedH2O/job.basis_1_50 b/examples/PinnedH2O/job.basis_1_50 new file mode 100644 index 00000000..a5069570 --- /dev/null +++ b/examples/PinnedH2O/job.basis_1_50 @@ -0,0 +1,32 @@ +#!/bin/tcsh +#SBATCH -N 2 +#SBATCH -t 0:10:00 +#SBATCH -p pbatch + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 64 + +set maindir = /p/lustre2/cheung26/mgmol-20241012 + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = ${maindir}/build_quartz/libROM/build/examples/misc/combine_samples + +set snapshot_files = "" +set 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..0f8ef90e --- /dev/null +++ b/examples/PinnedH2O/job.offline @@ -0,0 +1,35 @@ +#!/bin/tcsh +#SBATCH -N 2 +#SBATCH -t 1:00:00 +#SBATCH -p pbatch + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 64 + +set maindir = /p/lustre2/cheung26/mgmol-20241012 + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = mgmol-opt + +cp $maindir/install_quartz/bin/$exe . + +set datadir = $maindir/examples/PinnedH2O + +set cfg_offline = mgmol_offline.cfg +cp $datadir/$cfg_offline . + +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_offline -i coords.in > offline_PinnedH2O.out + +date diff --git a/examples/PinnedH2O/job.rom_1_50 b/examples/PinnedH2O/job.rom_1_50 new file mode 100644 index 00000000..aa6df335 --- /dev/null +++ b/examples/PinnedH2O/job.rom_1_50 @@ -0,0 +1,39 @@ +#!/bin/tcsh +#SBATCH -N 2 +#SBATCH -t 1:00:00 +#SBATCH -p pbatch + +date + +setenv OMP_NUM_THREADS 1 +#setenv KMP_DETERMINISTIC_REDUCTION 1 + +set ncpus = 64 + +set maindir = /p/lustre2/cheung26/mgmol-20241012 + +setenv LD_LIBRARY_PATH ${maindir}/build_quartz/libROM/build/lib:$LD_LIBRARY_PATH + +set exe = mgmol-opt + +cp $maindir/install_quartz/bin/$exe . + +set datadir = $maindir/examples/PinnedH2O + +set increment_md_steps = 1 +set num_md_steps = 50 +set basis_file = PinnedH2O_orbitals_basis_${increment_md_steps}_${num_md_steps} + +set cfg_rom = mgmol_rom_${increment_md_steps}_${num_md_steps}.cfg +cp $datadir/$cfg_rom . + +cp $datadir/coords.in . + +ln -s -f $maindir/potentials/pseudo.O_ONCV_PBE_SG15 . +ln -s -f $maindir/potentials/pseudo.H_ONCV_PBE_SG15 . + +source $maindir/scripts/modules.quartz + +srun -n $ncpus $exe -c $cfg_rom -i coords.in > rom_${increment_md_steps}_${num_md_steps}_PinnedH2O.out + +date diff --git a/examples/PinnedH2O/mgmol_offline.cfg b/examples/PinnedH2O/mgmol_offline.cfg new file mode 100644 index 00000000..332dd77d --- /dev/null +++ b/examples/PinnedH2O/mgmol_offline.cfg @@ -0,0 +1,38 @@ +verbosity=1 +xcFunctional=PBE +FDtype=Mehrstellen +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=MD +[MD] +num_steps=500 +dt=40. +thermostat=ON +[Thermostat] +type=Berendsen +temperature=1000. +relax_time=800. +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +[Restart] +output_level=4 +[ROM.offline] +save_librom_snapshot=true +librom_snapshot_freq=1 diff --git a/examples/PinnedH2O/mgmol_rom_1_50.cfg b/examples/PinnedH2O/mgmol_rom_1_50.cfg new file mode 100644 index 00000000..337a86eb --- /dev/null +++ b/examples/PinnedH2O/mgmol_rom_1_50.cfg @@ -0,0 +1,40 @@ +verbosity=1 +xcFunctional=PBE +FDtype=Mehrstellen +[Mesh] +nx=64 +ny=64 +nz=64 +[Domain] +ox=-6. +oy=-6. +oz=-6. +lx=12. +ly=12. +lz=12. +[Potentials] +pseudopotential=pseudo.O_ONCV_PBE_SG15 +pseudopotential=pseudo.H_ONCV_PBE_SG15 +[Run] +type=MD +[MD] +num_steps=500 +dt=40. +thermostat=ON +[Thermostat] +type=Berendsen +temperature=1000. +relax_time=800. +[Quench] +max_steps=100 +atol=1.e-8 +[Orbitals] +initial_type=Random +initial_width=2. +[Restart] +output_level=4 +[ROM.offline] +basis_file=PinnedH2O_orbitals_basis_1_50 +[ROM.basis] +compare_md=false +number_of_orbital_basis=39 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/examples/PinnedH2O/rotation_test.py b/examples/PinnedH2O/rotation_test.py new file mode 100644 index 00000000..79f1aaef --- /dev/null +++ b/examples/PinnedH2O/rotation_test.py @@ -0,0 +1,64 @@ +import numpy as np + +O1 = np.array([0.00, 0.00, 0.00]) +H1 = np.array([-0.45, 1.42, -1.07]) +H2 = np.array([-0.45, -1.48, -0.97]) + +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 + +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)] + ]) + +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)) + +rot_matrix_align_plane = rotation_matrix(axis_to_align, angle_to_align) +H1_rotated = np.dot(rot_matrix_align_plane, H1) +H2_rotated = np.dot(rot_matrix_align_plane, H2) + +bondlength1 = calculate_bondlength(H1, O1) +bondlength2 = calculate_bondlength(H2, O1) +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}') + +bondlength1 = calculate_bondlength(H1_rotated, O1) +bondlength2 = calculate_bondlength(H2_rotated, O1) +bondangle = calculate_bondangle(H1_rotated, O1, H2_rotated, False) +if bondlength1 < bondlength2: + H1_rotated, H2_rotated, bondlength1, bondlength2 = H2_rotated, H1_rotated, bondlength2, bondlength1 + +print('Rotated system in z=0 plane about x=0 axis, with longer bondlength in H1') +print(f'H1 = {H1_rotated}') +print(f'H2 = {H2_rotated}') +print(f'Bondlength of O1-H1 = {bondlength1}') +print(f'Bondlength of O1-H2 = {bondlength2}') +print(f'Angle between O1-H1 and O1-H2 = {bondangle}') diff --git a/scripts/build_quartz_libROM.sh b/scripts/build_quartz_libROM.sh new file mode 100644 index 00000000..7092f3d5 --- /dev/null +++ b/scripts/build_quartz_libROM.sh @@ -0,0 +1,61 @@ +#!/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 +#./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/build_ubuntu22_openmpi.sh b/scripts/build_ubuntu22_openmpi.sh index 23f510b6..a0ca470a 100755 --- a/scripts/build_ubuntu22_openmpi.sh +++ b/scripts/build_ubuntu22_openmpi.sh @@ -25,7 +25,7 @@ cmake -DCMAKE_INSTALL_PREFIX=${INSTALL_DIR} \ -DMPIEXEC_PREFLAGS="--oversubscribe" \ -DMGMOL_WITH_CLANG_FORMAT=ON \ -DCMAKE_PREFIX_PATH=${HOME}/bin \ - -D CMAKE_CXX_FLAGS="-Wall -pedantic -Wextra" \ + -DCMAKE_CXX_FLAGS="-Wall -pedantic -Wextra" \ .. # call make install 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 3e1188fa..e602101f 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 @@ -155,13 +159,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}) @@ -194,3 +208,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 dd786678..3f8a11cc 100644 --- a/src/Control.cc +++ b/src/Control.cc @@ -1882,6 +1882,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() @@ -2016,3 +2020,106 @@ 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_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 5134c718..43da0be7 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() { @@ -706,6 +713,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/DensityMatrix.cc b/src/DensityMatrix.cc index 78d7008e..139aa02d 100644 --- a/src/DensityMatrix.cc +++ b/src/DensityMatrix.cc @@ -43,6 +43,12 @@ DensityMatrix::DensityMatrix(const int ndim) { assert(ndim > 0); + dim_ = ndim; + + occ_uptodate_ = false; + stripped_ = false; + uniform_occ_ = false; + MGmol_MPI& mmpi = *(MGmol_MPI::instance()); orbital_occupation_ = mmpi.nspin() > 1 ? 1. : 2.; diff --git a/src/Electrostatic.h b/src/Electrostatic.h index 072ce0cc..a8eb4709 100644 --- a/src/Electrostatic.h +++ b/src/Electrostatic.h @@ -47,6 +47,11 @@ 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_; } + void setup(const short max_sweeps); void setupPB(const double e0, const double rho0, const double drho0, Potentials& pot); diff --git a/src/Hartree_CG.h b/src/Hartree_CG.h index 78e9da65..efa50919 100644 --- a/src/Hartree_CG.h +++ b/src/Hartree_CG.h @@ -40,6 +40,13 @@ class Hartree_CG : public Poisson void solve(const pb::GridFunc& rho, const pb::GridFunc& rhoc) override; + + void applyOperator(pb::GridFunc &vh, + pb::GridFunc &lhs) override + { + T *oper = poisson_solver_->getOperator(); + oper->apply(vh, lhs); + } }; #endif diff --git a/src/MGmol.cc b/src/MGmol.cc index 4aa1fa26..acc2b944 100644 --- a/src/MGmol.cc +++ b/src/MGmol.cc @@ -1121,6 +1121,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 } } diff --git a/src/MGmol.h b/src/MGmol.h index 4d82dd64..1eb8b2d3 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" @@ -183,6 +186,8 @@ class MGmol : public MGmolInterface { return hamiltonian_; } + std::shared_ptr> getRho() { return rho_; } + std::shared_ptr getIons() { return ions_; } void run() override; @@ -352,6 +357,11 @@ class MGmol : public MGmolInterface { return proj_matrices_; } + +#ifdef MGMOL_HAS_LIBROM + int save_orbital_snapshot(std::string snapshot_dir, OrbitalsType& orbitals); + void project_orbital(std::string snapshot_dir, int rdim, OrbitalsType& orbitals); +#endif }; // Instantiate static variables here to avoid clang warnings template 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/PCGSolver.h b/src/PCGSolver.h index 04e1028f..f6db7cc6 100644 --- a/src/PCGSolver.h +++ b/src/PCGSolver.h @@ -95,6 +95,8 @@ class PCGSolver double getFinalResidual() const { return final_residual_; } double getResidualReduction() const { return residual_reduction_; } + T* getOperator() { return &oper_; } + // Destructor ~PCGSolver() { clear(); } }; diff --git a/src/Poisson.h b/src/Poisson.h index 87107542..14716cda 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 380f05a0..a54f0c47 100644 --- a/src/Potentials.cc +++ b/src/Potentials.cc @@ -899,6 +899,97 @@ void Potentials::initBackground(Ions& ions) 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) + { + cout << "Potentials::evalIonDensityOnSamplePts - flag_filter s is not supported" + << 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; +} + template void Potentials::setVxc( const double* const vxc, const int iterativeIndex); template void Potentials::setVxc( diff --git a/src/Potentials.h b/src/Potentials.h index 205fb538..0ddd10fa 100644 --- a/src/Potentials.h +++ b/src/Potentials.h @@ -95,6 +95,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); + public: Potentials(const bool vh_frozen = false); @@ -159,6 +162,8 @@ class Potentials double getChargeInCell() const { return charge_in_cell_; } + const double getBackgroundCharge() const { return background_charge_; } + /*! * initialize total potential as local pseudopotential */ @@ -196,6 +201,8 @@ class Potentials void initBackground(Ions& ions); void addBackgroundToRhoComp(); + 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 b50618c8..36236f85 100644 --- a/src/md.cc +++ b/src/md.cc @@ -405,6 +405,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 @@ -488,6 +491,41 @@ void MGmol::md(OrbitalsType** orbitals, Ions& ions) // Compute forces force(**orbitals, ions); +#ifdef MGMOL_HAS_LIBROM + // TODO: cleanup + if (ct.getROMOptions().num_orbbasis > 0) + { + 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_orbital(ct.getROMOptions().basis_file, ct.getROMOptions().num_orbbasis, **orbitals); + if (ct.getROMOptions().compare_md) + { + force(**orbitals, ions); + } + else + { + double shift[3]; + for (short i = 0; i < 3; i++) shift[i] = 0.; + Ions ROM_ions(ions, shift); + 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_); + } + } + } +#endif + // set fion ions.getLocalForces(fion); @@ -634,6 +672,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 642b1b99..884aa8d3 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; @@ -200,6 +201,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"); @@ -404,3 +409,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..19ff66a3 --- /dev/null +++ b/src/rom.cc @@ -0,0 +1,97 @@ +// 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 "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 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..200061cb --- /dev/null +++ b/src/rom_Control.h @@ -0,0 +1,61 @@ +// 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_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..7eb323a7 --- /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..52ad6e47 --- /dev/null +++ b/src/rom_workflows.cc @@ -0,0 +1,948 @@ +// 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(-1); + 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 */ + basis_generator.takeSample(pot.vh_rho()); + 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); + } + + /* 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 */ + ions->setPositions(cfgs[idx], atnumbers); + + /* save overlapping ions for sanity check */ + fom_overlap_ions[idx].resize(ions->overlappingVL_ions().size()); + for (int k = 0; k < 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] = ions->overlappingVL_ions()[k]->position(d); + } + + /* compute resulting ion density */ + /* NOTE: we exclude rescaling for the sake of verification */ + pot.initialize(*ions); + + mgmol->electrostat_->setupRhoc(pot.rho_comp()); + fom_rhoc[idx].resize(dim); + mgmol->electrostat_->getRhoc()->init_vect(fom_rhoc[idx].data(), 'd'); + } + + /* 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 */ + ions->setPositions(cfgs[test_idx], atnumbers); + + /* Sanity check for overlapping ions */ + CAROM_VERIFY(fom_overlap_ions[test_idx].size() == ions->overlappingVL_ions().size()); + for (int k = 0; k < ions->overlappingVL_ions().size(); k++) + for (int d = 0; d < 3; d++) + CAROM_VERIFY(abs(fom_overlap_ions[test_idx][k][d] - ions->overlappingVL_ions()[k]->position(d)) < 1.0e-12); + + /* eval ion density on sample grid points */ + std::vector sampled_rhoc(sampled_row.size()); + pot.evalIonDensityOnSamplePts(*ions, sampled_row, sampled_rhoc); + + 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-12); + } +} + +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/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 +