diff --git a/INSTALL.rst b/INSTALL.rst index 229bf40a..e37d31eb 100644 --- a/INSTALL.rst +++ b/INSTALL.rst @@ -2,9 +2,9 @@ Installation ============ -The |qp| can be downloaded on GitHub as an `archive -`_ or as a `git -repository `_. +|qp| can be downloaded on GitHub as an `archive +`_ or as a `git +repository `_. .. code:: bash @@ -19,16 +19,16 @@ Before anything, go into your :file:`quantum_package` directory and run This script will create the :file:`quantum_package.rc` bash script, which -sets all the environment variables required for the normal operation of the -*Quantum Package*. It will also initialize the git submodules that are +sets all the environment variables required for the normal operation of +|qp|. It will also initialize the git submodules that are required, and tell you which external dependencies are missing and need to be installed. The required dependencies are located in the -`external/qp2-dependencies` directory, such that once QP is configured the +`external/qp2-dependencies` directory, such that once |qp| is configured the internet connection is not needed any more. When all dependencies have been installed, (the :command:`configure` will -inform you) source the :file:`quantum_package.rc` in order to load all -environment variables and compile the |QP|. +inform you what is missing) source the :file:`quantum_package.rc` in order to +load all environment variables and compile |QP|. Now all the requirements are met, you can compile the programs using @@ -37,6 +37,15 @@ Now all the requirements are met, you can compile the programs using make +Installation of dependencies via a Conda environment +==================================================== + +.. code:: bash + + conda env create -f qp2.yml + + + Requirements ============ @@ -64,8 +73,8 @@ architecture. Modify it if needed, and run :command:`configure` with .. code:: bash - cp ./config/gfortran.example config/gfortran.cfg - ./configure -c config/gfortran.cfg + cp ./config/gfortran.example config/gfortran_avx.cfg + ./configure -c config/gfortran_avx.cfg .. note:: @@ -86,45 +95,33 @@ The command is to be used as follows: .. code:: bash - ./configure --install= + ./configure -i The following packages are supported by the :command:`configure` installer: * ninja -* irpf90 * zeromq * f77zmq * gmp * ocaml (:math:`\approx` 5 minutes) -* ezfio * docopt * resultsFile * bats +* zlib Example: .. code:: bash - ./configure -i ezfio + ./configure -i ninja -.. note:: - - When installing the ocaml package, you will be asked the location of where - it should be installed. A safe option is to enter the path proposed by the - |QP|: - - QP>> Please install it here: /your_quantum_package_directory/bin - - So just enter the proposition of the |QP| and press enter. If the :command:`configure` executable fails to install a specific dependency ----------------------------------------------------------------------------- -If the :command:`configure` executable does not succeed to install a specific -dependency, there are some proposition of how to download and install the -minimal dependencies to compile and use the |QP|. - +If the :command:`configure` executable does not succeed in installing a specific +dependency, you should try to install the dependency on your system by yourself. Before doing anything below, try to install the packages with your package manager (:command:`apt`, :command:`yum`, etc). @@ -149,11 +146,11 @@ IRPF90 *IRPF90* is a Fortran code generator for programming using the Implicit Reference to Parameters (IRP) method. -If you have *pip* for Python2, you can do +If you have *pip* for Python2, you can do .. code:: bash - python2 -m pip install --user irpf90 + python3 -m pip install --user irpf90 Otherwise, @@ -209,7 +206,7 @@ ZeroMQ and its Fortran binding .. code:: bash - cp f77_zmq_free.h ${QP_ROOT}/src/ZMQ/f77_zmq.h + cp f77_zmq_free.h ${QP_ROOT}/src/zmq/f77_zmq.h Zlib @@ -262,53 +259,6 @@ With Debian or Ubuntu, you can use sudo apt install libgmp-dev -libcap ------- - -Libcap is a library for getting and setting POSIX.1e draft 15 capabilities. - -* Download the latest version of libcap here: - ``_ - and move it in the :file:`${QP_ROOT}/external` directory - -* Extract the archive, go into the :file:`libcap-*/libcap` directory and run - the following command - -.. code:: bash - - prefix=$QP_ROOT make install - -With Debian or Ubuntu, you can use - -.. code:: bash - - sudo apt install libcap-dev - - -Bubblewrap ----------- - -Bubblewrap is an unprivileged sandboxing tool. - -* Download Bubblewrap here: - ``_ - and move it in the :file:`${QP_ROOT}/external` directory - -* Extract the archive, go into the :file:`bubblewrap-*` directory and run - the following commands - -.. code:: bash - - ./configure --prefix=$QP_ROOT && make -j 8 - make install-exec-am - - -With Debian or Ubuntu, you can use - -.. code:: bash - - sudo apt install bubblewrap - OCaml @@ -327,7 +277,7 @@ OCaml ``_ and move it in the :file:`${QP_ROOT}/external` directory -* If you use OCaml only with the |qp|, you can install the OPAM directory +* If you use OCaml only with |qp|, you can install the OPAM directory containing the compiler and all the installed libraries in the :file:`${QP_ROOT}/external` directory as @@ -352,14 +302,14 @@ OCaml .. code:: bash - opam init --comp=4.07.1 + opam init --comp=4.11.1 eval `${QP_ROOT}/bin/opam env` If the installation fails because of bwrap, you can initialize opam using: .. code:: bash - opam init --disable-sandboxing --comp=4.07.1 + opam init --disable-sandboxing --comp=4.11.1 eval `${QP_ROOT}/bin/opam env` * Install the required external OCaml libraries @@ -369,17 +319,6 @@ OCaml opam install ocamlbuild cryptokit zmq sexplib ppx_sexp_conv ppx_deriving getopt -EZFIO ------ - -*EZFIO* is the Easy Fortran Input/Output library generator. - -* Download EZFIO here : ``_ and move - the downloaded archive in the :file:`${QP_ROOT}/external` directory - -* Extract the archive, and rename it as :file:`${QP_ROOT}/external/ezfio` - - Docopt ------ @@ -406,7 +345,7 @@ resultsFile *resultsFile* is a Python package to extract data from output files of quantum chemistry codes. -If you have *pip* for Python3, you can do +If you have *pip* for Python3, you can do .. code:: bash @@ -414,3 +353,4 @@ If you have *pip* for Python3, you can do + diff --git a/RELEASE_NOTES.org b/RELEASE_NOTES.org index 7178f0a0..7724d1d1 100644 --- a/RELEASE_NOTES.org +++ b/RELEASE_NOTES.org @@ -30,7 +30,8 @@ - Fixed bug in DIIS - Fixed bug in molden (Au -> Angs) - Fixed bug with non-contiguous MOs in active space and deleter MOs - - Network-free installation + - Complete network-free installation + - Fixed bug in selection when computing full PT2 - Updated version of f77-zmq *** User interface @@ -59,6 +60,7 @@ symmetry in matrices - qp_export_as_tgz exports also plugin codes - Added a basis module containing basis set information + - Added qp_run truncate_wf *** Code diff --git a/bin/qp_convert_output_to_ezfio b/bin/qp_convert_output_to_ezfio index 1ce78c73..091423e4 100755 --- a/bin/qp_convert_output_to_ezfio +++ b/bin/qp_convert_output_to_ezfio @@ -195,48 +195,52 @@ def write_ezfio(res, filename): # P a r s i n g # # ~#~#~#~#~#~#~ # + inucl = {} + for i, a in enumerate(res.geometry): + inucl[a.coord] = i + nbasis = 0 - nucl_center = [] + nucl_index = [] curr_center = -1 nucl_shell_num = [] ang_mom = [] nshell = 0 - shell_prim_index = [1] + nshell_tot = 0 + shell_index = [] shell_prim_num = [] for b in res.basis: s = b.sym if str.count(s, "y") + str.count(s, "x") == 0: - c = b.center + c = inucl[b.center] nshell += 1 + nshell_tot += 1 if c != curr_center: curr_center = c - nucl_center.append(nbasis+1) nucl_shell_num.append(nshell) nshell = 0 nbasis += 1 + nucl_index.append(c+1) coefficient += b.coef[:len(b.prim)] exponent += [p.expo for p in b.prim] ang_mom.append(str.count(s, "z")) - shell_prim_index.append(len(exponent)+1) shell_prim_num.append(len(b.prim)) - - nucl_shell_num.append(nshell+1) - nucl_shell_num = nucl_shell_num[1:] + shell_index += [nshell_tot+1] * len(b.prim) # ~#~#~#~#~ # # W r i t e # # ~#~#~#~#~ # ezfio.set_basis_basis("Read from ResultsFile") - ezfio.set_basis_basis_nucleus_index(nucl_center) - ezfio.set_basis_prim_num(len(coefficient)) ezfio.set_basis_shell_num(len(ang_mom)) + ezfio.set_basis_basis_nucleus_index(nucl_index) + ezfio.set_basis_prim_num(len(coefficient)) + ezfio.set_basis_nucleus_shell_num(nucl_shell_num) ezfio.set_basis_prim_coef(coefficient) ezfio.set_basis_prim_expo(exponent) ezfio.set_basis_shell_ang_mom(ang_mom) ezfio.set_basis_shell_prim_num(shell_prim_num) - ezfio.set_basis_shell_prim_index(shell_prim_index) + ezfio.set_basis_shell_index(shell_index) print("OK") diff --git a/bin/qp_plugins b/bin/qp_plugins index ef0f5a45..c9158422 100755 --- a/bin/qp_plugins +++ b/bin/qp_plugins @@ -6,6 +6,7 @@ Usage: qp_plugins download [-n ] qp_plugins install ... qp_plugins uninstall + qp_plugins remove qp_plugins update [-r ] qp_plugins create -n [-r ] [...] @@ -24,6 +25,8 @@ Options: uninstall Uninstall a plugin + remove Uninstall a plugin + update Update the repository create @@ -274,7 +277,7 @@ def main(arguments): subprocess.check_call(["qp_create_ninja", "update"]) print("[ OK ]") - elif arguments["uninstall"]: + elif arguments["uninstall"] or arguments["remove"]: m_instance = ModuleHandler([QP_SRC]) d_descendant = m_instance.dict_descendant diff --git a/bin/qp_set_frozen_core b/bin/qp_set_frozen_core index 879c71de..bc6f6834 100755 --- a/bin/qp_set_frozen_core +++ b/bin/qp_set_frozen_core @@ -7,12 +7,13 @@ setting all MOs as Active, except the n/2 first ones which are set as Core. If pseudo-potentials are used, all the MOs are set as Active. Usage: - qp_set_frozen_core [-q|--query] [(-l|-s|--large|--small)] EZFIO_DIR + qp_set_frozen_core [-q|--query] [(-l|-s|-u|--large|--small|--unset)] EZFIO_DIR Options: -q --query Prints in the standard output the number of frozen MOs -l --large Use a small core -s --small Use a large core + -u --unset Unset frozen core Default numbers of frozen electrons: @@ -88,7 +89,9 @@ def main(arguments): elif charge <= 54: n_frozen += 9 elif charge <= 86: n_frozen += 18 elif charge <= 118: n_frozen += 27 + elif arguments["--unset"]: + n_frozen = 0 else: # default for charge in ezfio.nuclei_nucl_charge: if charge <= 4: pass diff --git a/config/gfortran.cfg b/config/gfortran.cfg index 342acae9..56bb6ba4 100644 --- a/config/gfortran.cfg +++ b/config/gfortran.cfg @@ -13,7 +13,7 @@ FC : gfortran -g -ffree-line-length-none -I . -fPIC LAPACK_LIB : -lblas -llapack IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 --assert +IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED # Global options ################ @@ -35,14 +35,14 @@ OPENMP : 1 ; Append OpenMP flags # -ffast-math and the Fortran-specific # -fno-protect-parens and -fstack-arrays. [OPT] -FCFLAGS : -Ofast +FCFLAGS : -Ofast # Profiling flags ################# # [PROFILE] FC : -p -g -FCFLAGS : -Ofast +FCFLAGS : -Ofast # Debugging flags ################# @@ -58,5 +58,5 @@ FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurpris # [OPENMP] FC : -fopenmp -IRPF90_FLAGS : --openmp +IRPF90_FLAGS : --openmp diff --git a/config/gfortran_avx.cfg b/config/gfortran_avx.cfg index 4f45e3a1..747dff67 100644 --- a/config/gfortran_avx.cfg +++ b/config/gfortran_avx.cfg @@ -13,7 +13,7 @@ FC : gfortran -ffree-line-length-none -I . -mavx -g -fPIC LAPACK_LIB : -llapack -lblas IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 +IRPF90_FLAGS : --ninja --align=32 -DSET_NESTED # Global options ################ @@ -42,7 +42,7 @@ FCFLAGS : -Ofast -mavx # [PROFILE] FC : -p -g -FCFLAGS : -Ofast +FCFLAGS : -Ofast # Debugging flags ################# @@ -51,12 +51,12 @@ FCFLAGS : -Ofast # -g : Extra debugging information # [DEBUG] -FCFLAGS : -fcheck=all -g +FCFLAGS : -fcheck=all -g # OpenMP flags ################# # [OPENMP] FC : -fopenmp -IRPF90_FLAGS : --openmp +IRPF90_FLAGS : --openmp diff --git a/config/gfortran_debug.cfg b/config/gfortran_debug.cfg index 926255e0..51e5a500 100644 --- a/config/gfortran_debug.cfg +++ b/config/gfortran_debug.cfg @@ -13,7 +13,7 @@ FC : gfortran -g -ffree-line-length-none -I . -fPIC LAPACK_LIB : -lblas -llapack IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 --assert +IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED # Global options ################ @@ -35,14 +35,14 @@ OPENMP : 1 ; Append OpenMP flags # -ffast-math and the Fortran-specific # -fno-protect-parens and -fstack-arrays. [OPT] -FCFLAGS : -Ofast +FCFLAGS : -Ofast # Profiling flags ################# # [PROFILE] FC : -p -g -FCFLAGS : -Ofast +FCFLAGS : -Ofast # Debugging flags ################# @@ -59,5 +59,5 @@ FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurpris # [OPENMP] FC : -fopenmp -IRPF90_FLAGS : --openmp +IRPF90_FLAGS : --openmp diff --git a/config/gfortran_mpi.cfg b/config/gfortran_mpi.cfg index d72160c1..1af3ca45 100644 --- a/config/gfortran_mpi.cfg +++ b/config/gfortran_mpi.cfg @@ -13,7 +13,7 @@ FC : mpif90 -ffree-line-length-none -I . -g -fPIC LAPACK_LIB : -lblas -llapack IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 -DMPI +IRPF90_FLAGS : --ninja --align=32 -DMPI -DSET_NESTED # Global options ################ @@ -35,14 +35,14 @@ OPENMP : 1 ; Append OpenMP flags # -ffast-math and the Fortran-specific # -fno-protect-parens and -fstack-arrays. [OPT] -FCFLAGS : -Ofast -msse4.2 +FCFLAGS : -Ofast -msse4.2 # Profiling flags ################# # [PROFILE] FC : -p -g -FCFLAGS : -Ofast -msse4.2 +FCFLAGS : -Ofast -msse4.2 # Debugging flags ################# @@ -51,7 +51,7 @@ FCFLAGS : -Ofast -msse4.2 # -g : Extra debugging information # [DEBUG] -FCFLAGS : -fcheck=all -g +FCFLAGS : -fcheck=all -g # OpenMP flags ################# diff --git a/config/ifort_2019_avx.cfg b/config/ifort_2019_avx.cfg new file mode 100644 index 00000000..661a0e8f --- /dev/null +++ b/config/ifort_2019_avx.cfg @@ -0,0 +1,63 @@ +# Common flags +############## +# +# -mkl=[parallel|sequential] : Use the MKL library +# --ninja : Allow the utilisation of ninja. It is mandatory ! +# --align=32 : Align all provided arrays on a 32-byte boundary +# +[COMMON] +FC : ifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED + +# Global options +################ +# +# 1 : Activate +# 0 : Deactivate +# +[OPTION] +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +CACHE : 0 ; Enable cache_compile.py +OPENMP : 1 ; Append OpenMP flags + +# Optimization flags +#################### +# +# -xHost : Compile a binary optimized for the current architecture +# -O2 : O3 not better than O2. +# -ip : Inter-procedural optimizations +# -ftz : Flushes denormal results to zero +# +[OPT] +FC : -traceback +FCFLAGS : -xAVX -O2 -ip -ftz -g + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -xSSE4.2 -O2 -ip -ftz + +# Debugging flags +################# +# +# -traceback : Activate backtrace on runtime +# -fpe0 : All floating point exaceptions +# -C : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# -xSSE2 : Valgrind needs a very simple x86 executable +# +[DEBUG] +FC : -g -traceback +FCFLAGS : -xSSE2 -C -fpe0 -implicitnone + +# OpenMP flags +################# +# +[OPENMP] +FC : -qopenmp +IRPF90_FLAGS : --openmp + diff --git a/config/ifort_2019_avx_mpi.cfg b/config/ifort_2019_avx_mpi.cfg new file mode 100644 index 00000000..2d212db5 --- /dev/null +++ b/config/ifort_2019_avx_mpi.cfg @@ -0,0 +1,64 @@ +# Common flags +############## +# +# -mkl=[parallel|sequential] : Use the MKL library +# --ninja : Allow the utilisation of ninja. It is mandatory ! +# --align=32 : Align all provided arrays on a 32-byte boundary +# +[COMMON] +FC : mpiifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL -DSET_NESTED + +# Global options +################ +# +# 1 : Activate +# 0 : Deactivate +# +[OPTION] +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +CACHE : 0 ; Enable cache_compile.py +OPENMP : 1 ; Append OpenMP flags + +# Optimization flags +#################### +# +# -xHost : Compile a binary optimized for the current architecture +# -O2 : O3 not better than O2. +# -ip : Inter-procedural optimizations +# -ftz : Flushes denormal results to zero +# +[OPT] +FCFLAGS : -mavx -axAVX -O2 -ip -ftz -g -traceback + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -march=corei7 -O2 -ip -ftz + + +# Debugging flags +################# +# +# -traceback : Activate backtrace on runtime +# -fpe0 : All floating point exaceptions +# -C : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# -xSSE2 : Valgrind needs a very simple x86 executable +# +[DEBUG] +FC : -g -traceback +FCFLAGS : -xSSE2 -C -fpe0 -implicitnone + + +# OpenMP flags +################# +# +[OPENMP] +FC : -qopenmp +IRPF90_FLAGS : --openmp + diff --git a/config/ifort_2019_mpi_rome.cfg b/config/ifort_2019_mpi_rome.cfg new file mode 100644 index 00000000..171219e6 --- /dev/null +++ b/config/ifort_2019_mpi_rome.cfg @@ -0,0 +1,63 @@ +# Common flags +############## +# +# -mkl=[parallel|sequential] : Use the MKL library +# --ninja : Allow the utilisation of ninja. It is mandatory ! +# --align=32 : Align all provided arrays on a 32-byte boundary +# +[COMMON] +FC : mpiifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED + +# Global options +################ +# +# 1 : Activate +# 0 : Deactivate +# +[OPTION] +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +CACHE : 0 ; Enable cache_compile.py +OPENMP : 1 ; Append OpenMP flags + +# Optimization flags +#################### +# +# -xHost : Compile a binary optimized for the current architecture +# -O2 : O3 not better than O2. +# -ip : Inter-procedural optimizations +# -ftz : Flushes denormal results to zero +# +[OPT] +FC : -traceback -shared-intel +FCFLAGS : -O2 -ip -g -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -xSSE4.2 -O2 -ip -ftz + +# Debugging flags +################# +# +# -traceback : Activate backtrace on runtime +# -fpe0 : All floating point exaceptions +# -C : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# -xSSE2 : Valgrind needs a very simple x86 executable +# +[DEBUG] +FC : -g -traceback +FCFLAGS : -xSSE2 -C -fpe0 -implicitnone + +# OpenMP flags +################# +# +[OPENMP] +FC : -qopenmp +IRPF90_FLAGS : --openmp + diff --git a/config/ifort_2019_rome.cfg b/config/ifort_2019_rome.cfg new file mode 100644 index 00000000..e923a1dd --- /dev/null +++ b/config/ifort_2019_rome.cfg @@ -0,0 +1,63 @@ +# Common flags +############## +# +# -mkl=[parallel|sequential] : Use the MKL library +# --ninja : Allow the utilisation of ninja. It is mandatory ! +# --align=32 : Align all provided arrays on a 32-byte boundary +# +[COMMON] +FC : ifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED + +# Global options +################ +# +# 1 : Activate +# 0 : Deactivate +# +[OPTION] +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +CACHE : 0 ; Enable cache_compile.py +OPENMP : 1 ; Append OpenMP flags + +# Optimization flags +#################### +# +# -xHost : Compile a binary optimized for the current architecture +# -O2 : O3 not better than O2. +# -ip : Inter-procedural optimizations +# -ftz : Flushes denormal results to zero +# +[OPT] +FC : -traceback -shared-intel +FCFLAGS : -O2 -ip -g -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -xSSE4.2 -O2 -ip -ftz + +# Debugging flags +################# +# +# -traceback : Activate backtrace on runtime +# -fpe0 : All floating point exaceptions +# -C : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# -xSSE2 : Valgrind needs a very simple x86 executable +# +[DEBUG] +FC : -g -traceback +FCFLAGS : -xSSE2 -C -fpe0 -implicitnone + +# OpenMP flags +################# +# +[OPENMP] +FC : -qopenmp +IRPF90_FLAGS : --openmp + diff --git a/config/ifort_2019_sse4.cfg b/config/ifort_2019_sse4.cfg new file mode 100644 index 00000000..a3aa7cbd --- /dev/null +++ b/config/ifort_2019_sse4.cfg @@ -0,0 +1,63 @@ +# Common flags +############## +# +# -mkl=[parallel|sequential] : Use the MKL library +# --ninja : Allow the utilisation of ninja. It is mandatory ! +# --align=32 : Align all provided arrays on a 32-byte boundary +# +[COMMON] +FC : ifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=32 -DINTEL -DSET_NESTED + +# Global options +################ +# +# 1 : Activate +# 0 : Deactivate +# +[OPTION] +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +CACHE : 0 ; Enable cache_compile.py +OPENMP : 1 ; Append OpenMP flags + +# Optimization flags +#################### +# +# -xHost : Compile a binary optimized for the current architecture +# -O2 : O3 not better than O2. +# -ip : Inter-procedural optimizations +# -ftz : Flushes denormal results to zero +# +[OPT] +FC : -traceback +FCFLAGS : -xSSE4.2 -O2 -ip -ftz -g + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -xSSE4.2 -O2 -ip -ftz + +# Debugging flags +################# +# +# -traceback : Activate backtrace on runtime +# -fpe0 : All floating point exaceptions +# -C : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# -xSSE2 : Valgrind needs a very simple x86 executable +# +[DEBUG] +FC : -g -traceback +FCFLAGS : -xSSE2 -C -fpe0 -implicitnone + +# OpenMP flags +################# +# +[OPENMP] +FC : -qopenmp +IRPF90_FLAGS : --openmp + diff --git a/config/ifort_2019_sse4_mpi.cfg b/config/ifort_2019_sse4_mpi.cfg new file mode 100644 index 00000000..6959d176 --- /dev/null +++ b/config/ifort_2019_sse4_mpi.cfg @@ -0,0 +1,64 @@ +# Common flags +############## +# +# -mkl=[parallel|sequential] : Use the MKL library +# --ninja : Allow the utilisation of ninja. It is mandatory ! +# --align=32 : Align all provided arrays on a 32-byte boundary +# +[COMMON] +FC : mpiifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL -DSET_NESTED + +# Global options +################ +# +# 1 : Activate +# 0 : Deactivate +# +[OPTION] +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +CACHE : 0 ; Enable cache_compile.py +OPENMP : 1 ; Append OpenMP flags + +# Optimization flags +#################### +# +# -xHost : Compile a binary optimized for the current architecture +# -O2 : O3 not better than O2. +# -ip : Inter-procedural optimizations +# -ftz : Flushes denormal results to zero +# +[OPT] +FCFLAGS : -msse4.2 -O2 -ip -ftz -g -traceback + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -msse4.2 -O2 -ip -ftz + + +# Debugging flags +################# +# +# -traceback : Activate backtrace on runtime +# -fpe0 : All floating point exaceptions +# -C : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# -xSSE2 : Valgrind needs a very simple x86 executable +# +[DEBUG] +FC : -g -traceback +FCFLAGS : -xSSE2 -C -fpe0 -implicitnone + + +# OpenMP flags +################# +# +[OPENMP] +FC : -qopenmp +IRPF90_FLAGS : --openmp + diff --git a/config/ifort_2019_xHost.cfg b/config/ifort_2019_xHost.cfg new file mode 100644 index 00000000..22d28803 --- /dev/null +++ b/config/ifort_2019_xHost.cfg @@ -0,0 +1,63 @@ +# Common flags +############## +# +# -mkl=[parallel|sequential] : Use the MKL library +# --ninja : Allow the utilisation of ninja. It is mandatory ! +# --align=32 : Align all provided arrays on a 32-byte boundary +# +[COMMON] +FC : ifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=64 -DINTEL -DSET_NESTED + +# Global options +################ +# +# 1 : Activate +# 0 : Deactivate +# +[OPTION] +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +CACHE : 0 ; Enable cache_compile.py +OPENMP : 1 ; Append OpenMP flags + +# Optimization flags +#################### +# +# -xHost : Compile a binary optimized for the current architecture +# -O2 : O3 not better than O2. +# -ip : Inter-procedural optimizations +# -ftz : Flushes denormal results to zero +# +[OPT] +FC : -traceback +FCFLAGS : -xHost -O2 -ip -ftz -g + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -xSSE4.2 -O2 -ip -ftz + +# Debugging flags +################# +# +# -traceback : Activate backtrace on runtime +# -fpe0 : All floating point exaceptions +# -C : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# -xSSE2 : Valgrind needs a very simple x86 executable +# +[DEBUG] +FC : -g -traceback +FCFLAGS : -xSSE2 -C -fpe0 -implicitnone + +# OpenMP flags +################# +# +[OPENMP] +FC : -qopenmp +IRPF90_FLAGS : --openmp + diff --git a/config/ifort_avx.cfg b/config/ifort_2021_avx.cfg similarity index 97% rename from config/ifort_avx.cfg rename to config/ifort_2021_avx.cfg index a2cb4c8a..6f657052 100644 --- a/config/ifort_avx.cfg +++ b/config/ifort_2021_avx.cfg @@ -9,7 +9,7 @@ FC : ifort -fpic LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 -DINTEL +IRPF90_FLAGS : --ninja --align=32 -DINTEL # Global options ################ diff --git a/config/ifort_avx_mpi.cfg b/config/ifort_2021_avx_mpi.cfg similarity index 96% rename from config/ifort_avx_mpi.cfg rename to config/ifort_2021_avx_mpi.cfg index f2bb8889..c991a4a9 100644 --- a/config/ifort_avx_mpi.cfg +++ b/config/ifort_2021_avx_mpi.cfg @@ -9,7 +9,7 @@ FC : mpiifort -fpic LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL +IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL # Global options ################ diff --git a/config/ifort_2021_mpi_rome.cfg b/config/ifort_2021_mpi_rome.cfg new file mode 100644 index 00000000..8413d23d --- /dev/null +++ b/config/ifort_2021_mpi_rome.cfg @@ -0,0 +1,63 @@ +# Common flags +############## +# +# -mkl=[parallel|sequential] : Use the MKL library +# --ninja : Allow the utilisation of ninja. It is mandatory ! +# --align=32 : Align all provided arrays on a 32-byte boundary +# +[COMMON] +FC : mpiifort -fpic +LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps +IRPF90 : irpf90 +IRPF90_FLAGS : --ninja --align=32 -DINTEL + +# Global options +################ +# +# 1 : Activate +# 0 : Deactivate +# +[OPTION] +MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +CACHE : 0 ; Enable cache_compile.py +OPENMP : 1 ; Append OpenMP flags + +# Optimization flags +#################### +# +# -xHost : Compile a binary optimized for the current architecture +# -O2 : O3 not better than O2. +# -ip : Inter-procedural optimizations +# -ftz : Flushes denormal results to zero +# +[OPT] +FC : -traceback -shared-intel +FCFLAGS : -O2 -ip -g -march=core-avx2 -align array64byte -fma -ftz -fomit-frame-pointer + +# Profiling flags +################# +# +[PROFILE] +FC : -p -g +FCFLAGS : -xSSE4.2 -O2 -ip -ftz + +# Debugging flags +################# +# +# -traceback : Activate backtrace on runtime +# -fpe0 : All floating point exaceptions +# -C : Checks uninitialized variables, array subscripts, etc... +# -g : Extra debugging information +# -xSSE2 : Valgrind needs a very simple x86 executable +# +[DEBUG] +FC : -g -traceback +FCFLAGS : -xSSE2 -C -fpe0 -implicitnone + +# OpenMP flags +################# +# +[OPENMP] +FC : -qopenmp +IRPF90_FLAGS : --openmp + diff --git a/config/ifort_rome.cfg b/config/ifort_2021_rome.cfg similarity index 97% rename from config/ifort_rome.cfg rename to config/ifort_2021_rome.cfg index 5ed01227..b3023186 100644 --- a/config/ifort_rome.cfg +++ b/config/ifort_2021_rome.cfg @@ -9,7 +9,7 @@ FC : ifort -fpic LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 -DINTEL +IRPF90_FLAGS : --ninja --align=32 -DINTEL # Global options ################ diff --git a/config/ifort.cfg b/config/ifort_2021_sse4.cfg similarity index 97% rename from config/ifort.cfg rename to config/ifort_2021_sse4.cfg index 714c4b10..a6299665 100644 --- a/config/ifort.cfg +++ b/config/ifort_2021_sse4.cfg @@ -9,7 +9,7 @@ FC : ifort -fpic LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 -DINTEL +IRPF90_FLAGS : --ninja --align=32 -DINTEL # Global options ################ diff --git a/config/ifort_mpi.cfg b/config/ifort_2021_sse4_mpi.cfg similarity index 96% rename from config/ifort_mpi.cfg rename to config/ifort_2021_sse4_mpi.cfg index e0d489a0..6ae56d2a 100644 --- a/config/ifort_mpi.cfg +++ b/config/ifort_2021_sse4_mpi.cfg @@ -9,7 +9,7 @@ FC : mpiifort -fpic LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL +IRPF90_FLAGS : --ninja --align=32 -DMPI -DINTEL # Global options ################ diff --git a/config/ifort_xHost.cfg b/config/ifort_2021_xHost.cfg similarity index 97% rename from config/ifort_xHost.cfg rename to config/ifort_2021_xHost.cfg index ddb4aa2d..1e76a69d 100644 --- a/config/ifort_xHost.cfg +++ b/config/ifort_2021_xHost.cfg @@ -9,7 +9,7 @@ FC : ifort -fpic LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=64 -DINTEL +IRPF90_FLAGS : --ninja --align=64 -DINTEL # Global options ################ diff --git a/config/ifort_debug.cfg b/config/ifort_debug.cfg deleted file mode 100644 index 9b718380..00000000 --- a/config/ifort_debug.cfg +++ /dev/null @@ -1,66 +0,0 @@ -# Common flags -############## -# -# -mkl=[parallel|sequential] : Use the MKL library -# --ninja : Allow the utilisation of ninja. It is mandatory ! -# --align=32 : Align all provided arrays on a 32-byte boundary -# -[COMMON] -FC : ifort -fpic -LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps -IRPF90 : irpf90 -IRPF90_FLAGS : --ninja --align=32 --assert -DINTEL - -# Global options -################ -# -# 1 : Activate -# 0 : Deactivate -# -[OPTION] -MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below -CACHE : 0 ; Enable cache_compile.py -OPENMP : 1 ; Append OpenMP flags - -# Optimization flags -#################### -# -# -xHost : Compile a binary optimized for the current architecture -# -O2 : O3 not better than O2. -# -ip : Inter-procedural optimizations -# -ftz : Flushes denormal results to zero -# -[OPT] -FC : -traceback -FCFLAGS : -msse4.2 -O2 -ip -ftz -g - - -# Profiling flags -################# -# -[PROFILE] -FC : -p -g -FCFLAGS : -msse4.2 -O2 -ip -ftz - - -# Debugging flags -################# -# -# -traceback : Activate backtrace on runtime -# -fpe0 : All floating point exaceptions -# -C : Checks uninitialized variables, array subscripts, etc... -# -g : Extra debugging information -# -msse4.2 : Valgrind needs a very simple x86 executable -# -[DEBUG] -FC : -g -traceback -FCFLAGS : -msse4.2 -check all -debug all -fpe-all=0 -implicitnone - - -# OpenMP flags -################# -# -[OPENMP] -FC : -qopenmp -IRPF90_FLAGS : --openmp - diff --git a/etc/qp.rc b/etc/qp.rc index 7661aadf..c56661c7 100644 --- a/etc/qp.rc +++ b/etc/qp.rc @@ -204,6 +204,9 @@ _qp_Complete() uninstall) COMPREPLY=( $(compgen -W "$(qp_plugins list -i)" -- $cur ) ) return 0;; + remove) + COMPREPLY=( $(compgen -W "$(qp_plugins list -i)" -- $cur ) ) + return 0;; create) COMPREPLY=( $(compgen -W "-n " -- $cur ) ) return 0;; diff --git a/ocaml/qp_create_ezfio.ml b/ocaml/qp_create_ezfio.ml index a4865e2b..be6c305b 100644 --- a/ocaml/qp_create_ezfio.ml +++ b/ocaml/qp_create_ezfio.ml @@ -91,7 +91,7 @@ let run ?o b au c d m p cart xyz_file = | Element e -> Element.to_string e | Int_elem (i,e) -> Printf.sprintf "%d,%s" (Nucl_number.to_int i) (Element.to_string e) in - Hashtbl.find basis_table key + Hashtbl.find basis_table key in let temp_filename = @@ -132,7 +132,7 @@ let run ?o b au c d m p cart xyz_file = Element.to_string elem.Atom.element in Hashtbl.add basis_table key new_channel - ) nuclei + ) nuclei end | Some (key, basis) -> (*Aux basis *) begin @@ -277,6 +277,16 @@ let run ?o b au c d m p cart xyz_file = ) nuclei in + let z_core = + List.map (fun x -> + Positive_int.to_int x.Pseudo.n_elec + |> float_of_int + ) pseudo + in + let nucl_num = (List.length z_core) in + Ezfio.set_pseudo_nucl_charge_remove (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| nucl_num |] ~data:z_core); + let molecule = let n_elec_to_remove = List.fold_left (fun accu x -> @@ -293,13 +303,13 @@ let run ?o b au c d m p cart xyz_file = Molecule.nuclei = let charges = list_map (fun x -> Positive_int.to_int x.Pseudo.n_elec - |> Float.of_int) pseudo + |> Float.of_int) pseudo |> Array.of_list in List.mapi (fun i x -> { x with Atom.charge = (Charge.to_float x.Atom.charge) -. charges.(i) |> Charge.of_float } - ) molecule.Molecule.nuclei + ) molecule.Molecule.nuclei } in let nuclei = @@ -356,11 +366,11 @@ let run ?o b au c d m p cart xyz_file = in if (x > accu) then x else accu - ) 0 x.Pseudo.non_local + ) 0 x.Pseudo.non_local in if (x > accu) then x else accu - ) 0 pseudo + ) 0 pseudo in let kmax = @@ -368,10 +378,10 @@ let run ?o b au c d m p cart xyz_file = list_map (fun x -> List.filter (fun (y,_) -> (Positive_int.to_int y.Pseudo.GaussianPrimitive_non_local.proj) = i) - x.Pseudo.non_local - |> List.length ) pseudo + x.Pseudo.non_local + |> List.length ) pseudo |> List.fold_left (fun accu x -> - if accu > x then accu else x) 0 + if accu > x then accu else x) 0 ) |> Array.fold_left (fun accu i -> if i > accu then i else accu) 0 @@ -396,11 +406,11 @@ let run ?o b au c d m p cart xyz_file = in tmp_array_dz_k.(i).(j) <- y; tmp_array_n_k.(i).(j) <- z; - ) x.Pseudo.local + ) x.Pseudo.local ) pseudo ; let concat_2d tmp_array = let data = - Array.map Array.to_list tmp_array + Array.map Array.to_list tmp_array |> Array.to_list |> List.concat in @@ -438,14 +448,14 @@ let run ?o b au c d m p cart xyz_file = tmp_array_dz_kl.(k).(i).(j) <- y; tmp_array_n_kl.(k).(i).(j) <- z; last_idx.(k) <- i+1; - ) x.Pseudo.non_local + ) x.Pseudo.non_local ) pseudo ; let concat_3d tmp_array = let data = Array.map (fun x -> Array.map Array.to_list x |> Array.to_list - |> List.concat) tmp_array + |> List.concat) tmp_array |> Array.to_list |> List.concat in @@ -513,8 +523,8 @@ let run ?o b au c d m p cart xyz_file = Ezfio.set_ao_basis_ao_num ao_num; Ezfio.set_ao_basis_ao_basis b; Ezfio.set_basis_basis b; - let ao_prim_num = list_map (fun (_,g,_) -> List.length g.Gto.lc) long_basis - and ao_nucl = list_map (fun (_,_,n) -> Nucl_number.to_int n) long_basis + let ao_prim_num = list_map (fun (_,g,_) -> List.length g.Gto.lc) long_basis + and ao_nucl = list_map (fun (_,_,n) -> Nucl_number.to_int n) long_basis and ao_power= let l = list_map (fun (x,_,_) -> x) long_basis in (list_map (fun t -> Positive_int.to_int Angmom.Xyz.(t.x)) l)@ @@ -526,7 +536,7 @@ let run ?o b au c d m p cart xyz_file = else s) 0 ao_prim_num in let gtos = - list_map (fun (_,x,_) -> x) long_basis + list_map (fun (_,x,_) -> x) long_basis in let create_expo_coef ec = @@ -534,10 +544,10 @@ let run ?o b au c d m p cart xyz_file = begin match ec with | `Coefs -> list_map (fun x-> list_map (fun (_,coef) -> - AO_coef.to_float coef) x.Gto.lc) gtos + AO_coef.to_float coef) x.Gto.lc) gtos | `Expos -> list_map (fun x-> list_map (fun (prim,_) -> AO_expo.to_float - prim.GaussianPrimitive.expo) x.Gto.lc) gtos + prim.GaussianPrimitive.expo) x.Gto.lc) gtos end in let rec get_n n accu = function @@ -567,7 +577,7 @@ let run ?o b au c d m p cart xyz_file = list_map ( fun (g,_) -> g.Gto.lc ) basis in let ang_mom = - list_map (fun (l : (GaussianPrimitive.t * Qptypes.AO_coef.t) list) -> + list_map (fun (l : (GaussianPrimitive.t * Qptypes.AO_coef.t) list) -> let x, _ = List.hd l in Angmom.to_l x.GaussianPrimitive.sym |> Qptypes.Positive_int.to_int ) lc @@ -577,7 +587,7 @@ let run ?o b au c d m p cart xyz_file = |> List.concat in let coef = - list_map (fun l -> + list_map (fun l -> list_map (fun (_,x) -> Qptypes.AO_coef.to_float x) l ) lc |> List.concat @@ -585,12 +595,16 @@ let run ?o b au c d m p cart xyz_file = let shell_prim_num = list_map List.length lc in - let shell_prim_idx = + let shell_idx = + let rec make_list n accu = function + | 0 -> accu + | i -> make_list n (n :: accu) (i-1) + in let rec aux count accu = function | [] -> List.rev accu | l::rest -> - let newcount = count+(List.length l) in - aux newcount (count::accu) rest + let new_l = make_list count accu (List.length l) in + aux (count+1) new_l rest in aux 1 [] lc in @@ -602,26 +616,18 @@ let run ?o b au c d m p cart xyz_file = ~rank:1 ~dim:[| shell_num |] ~data:shell_prim_num); Ezfio.set_basis_shell_ang_mom (Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| shell_num |] ~data:ang_mom ) ; - Ezfio.set_basis_shell_prim_index (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| shell_num |] ~data:shell_prim_idx) ; + Ezfio.set_basis_shell_index (Ezfio.ezfio_array_of_list + ~rank:1 ~dim:[| prim_num |] ~data:shell_idx) ; Ezfio.set_basis_basis_nucleus_index (Ezfio.ezfio_array_of_list - ~rank:1 ~dim:[| nucl_num |] - ~data:( - list_map (fun (_,n) -> Nucl_number.to_int n) basis - |> List.fold_left (fun accu i -> - match accu with - | [] -> [] - | (h,j) :: rest -> if j == i then ((h+1,j)::rest) else ((h+1,i)::(h+1,j)::rest) - ) [(0,0)] - |> List.rev - |> List.map fst - )) ; + ~rank:1 ~dim:[| shell_num |] + ~data:( list_map (fun (_,n) -> Nucl_number.to_int n) basis) + ) ; Ezfio.set_basis_nucleus_shell_num(Ezfio.ezfio_array_of_list ~rank:1 ~dim:[| nucl_num |] ~data:( list_map (fun (_,n) -> Nucl_number.to_int n) basis - |> List.fold_left (fun accu i -> - match accu with + |> List.fold_left (fun accu i -> + match accu with | [] -> [(1,i)] | (h,j) :: rest -> if j == i then ((h+1,j)::rest) else ((1,i)::(h,j)::rest) ) [] @@ -717,7 +723,7 @@ If a file with the same name as the basis set exists, this file will be read. O anonymous "FILE" Mandatory "Input file in xyz format or z-matrix."; ] - |> set_specs + |> set_specs end; @@ -741,7 +747,7 @@ If a file with the same name as the basis set exists, this file will be read. O | None -> 0 | Some x -> ( if x.[0] = 'm' then ~- (int_of_string (String.sub x 1 (String.length x - 1))) - else + else int_of_string x ) in diff --git a/scripts/cipsi_save.sh b/scripts/cipsi_save.sh new file mode 100644 index 00000000..a4d9b65e --- /dev/null +++ b/scripts/cipsi_save.sh @@ -0,0 +1,27 @@ +#!/bin/bash +# +# This script runs a CIPSI calculation as a sequence of single CIPSI iterations. +# After each iteration, the EZFIO directory is saved. +# +# Usage: cipsi_save [EZFIO_FILE] [NDET] +# +# Example: cipsi_save file.ezfio 10000 + +EZ=$1 +NDETMAX=$2 + +qp set_file ${EZ} +qp reset -d +qp set determinants read_wf true +declare -i NDET +NDET=1 +while [[ ${NDET} -lt ${NDETMAX} ]] +do + NDET=$(($NDET + $NDET)) + qp set determinants n_det_max $NDET + qp run fci > ${EZ}.out + NDET=$(qp get determinants n_det) + mv ${EZ}.out ${EZ}.${NDET}.out + cp -r ${EZ} ${EZ}.${NDET} +done + diff --git a/scripts/module/module_handler.py b/scripts/module/module_handler.py index d66918e2..6dd74f34 100755 --- a/scripts/module/module_handler.py +++ b/scripts/module/module_handler.py @@ -116,6 +116,7 @@ def get_l_module_descendant(d_child, l_module): print("Error: ", file=sys.stderr) print("`{0}` is not a submodule".format(module), file=sys.stderr) print("Check the typo (spelling, case, '/', etc.) ", file=sys.stderr) +# pass sys.exit(1) return list(set(l)) diff --git a/scripts/verif_omp/check_omp.f90 b/scripts/verif_omp/check_omp.f90 new file mode 100644 index 00000000..ca6af8bd --- /dev/null +++ b/scripts/verif_omp/check_omp.f90 @@ -0,0 +1,175 @@ +program check_omp_v2 + + use omp_lib + + implicit none + + integer :: accu, accu2 + integer :: s, n_setting + logical :: verbose, test_versions + logical, allocatable :: is_working(:) + + verbose = .False. + test_versions = .True. + n_setting = 4 + + allocate(is_working(n_setting)) + + is_working = .False. + + ! set the number of threads + call omp_set_num_threads(2) + + do s = 1, n_setting + + accu = 0 + accu2 = 0 + + call omp_set_max_active_levels(1) + call omp_set_nested(.False.) + + if (s==1) then + !call set_multiple_levels_omp() + cycle + elseif (s==2) then + call omp_set_max_active_levels(5) + elseif (s==3) then + call omp_set_nested(.True.) + else + call omp_set_nested(.True.) + call omp_set_max_active_levels(5) + endif + + ! Level 1 + !$OMP PARALLEL + if (verbose) then + print*,'Num threads level 1:',omp_get_num_threads() + endif + + ! Level 2 + !$OMP PARALLEL + if (verbose) then + print*,'Num threads level 2:',omp_get_num_threads() + endif + + ! Level 3 + !$OMP PARALLEL + if (verbose) then + print*,'Num threads level 3:',omp_get_num_threads() + endif + + call check_omp_in_subroutine(accu2) + + ! Level 4 + !$OMP PARALLEL + + if (verbose) then + print*,'Num threads level 4:',omp_get_num_threads() + endif + + !$OMP ATOMIC + accu = accu + 1 + !$OMP END ATOMIC + + !$OMP END PARALLEL + + + !$OMP END PARALLEL + + + !$OMP END PARALLEL + + + !$OMP END PARALLEL + + if (verbose) then + print*,'Setting:',s,'accu=',accu + print*,'Setting:',s,'accu2=',accu2 + endif + + if (accu == 16 .and. accu2 == 16) then + is_working(s) = .True. + endif + + enddo + + if (verbose) then + if (is_working(2)) then + print*,'The parallelization works on 4 levels with:' + print*,'call omp_set_max_active_levels(5)' + print*,'' + print*,'Please use the irpf90 flags -DSET_MAX_ACT in qp2/config/${compiler_name}.cfg' + elseif (is_working(3)) then + print*,'The parallelization works on 4 levels with:' + print*,'call omp_set_nested(.True.)' + print*,'' + print*,'Please use the irpf90 flag -DSET_NESTED in qp2/config/${compiler_name}.cfg' + elseif (is_working(4)) then + print*,'The parallelization works on 4 levels with:' + print*,'call omp_set_nested(.True.)' + print*,'+' + print*,'call omp_set_max_active_levels(5)' + print*,'' + print*,'Please use the irpf90 flags -DSET_NESTED -DSET_MAX_ACT in qp2/config/${compiler_name}.cfg' + else + print*,'The parallelization on multiple levels does not work with:' + print*,'call omp_set_max_active_levels(5)' + print*,'or' + print*,'call omp_set_nested(.True.)' + print*,'or' + print*,'call omp_set_nested(.True.)' + print*,'+' + print*,'call omp_set_max_active_levels(5)' + print*,'' + print*,'Try an other compiler and good luck...' + endif + + ! if (is_working(1)) then + ! print*,'' + ! print*,'==========================================================' + ! print*,'Your actual set up works for parallelization with 4 levels' + ! print*,'==========================================================' + ! print*,'' + ! else + ! print*,'' + ! print*,'===================================================================' + ! print*,'Your actual set up does not work for parallelization with 4 levels' + ! print*,'Please look at the previous messages to understand the requirements' + ! print*,'===================================================================' + ! print*,'' + ! endif + endif + + ! List of working flags + if (test_versions) then + print*,'Tests:',is_working(2:4) + endif + + ! IRPF90_FLAGS + if (is_working(2)) then + print*,'-DSET_MAX_ACT' + elseif (is_working(3)) then + print*,'-DSET_NESTED' + elseif (is_working(4)) then + print*,'-DSET_MAX_ACT -DSET_NESTED' + else + print*,'ERROR' + endif + +end + +subroutine check_omp_in_subroutine(accu2) + + implicit none + + integer, intent(inout) :: accu2 + + !$OMP PARALLEL + + !$OMP ATOMIC + accu2 = accu2 + 1 + !$OMP END ATOMIC + + !$OMP END PARALLEL + +end diff --git a/scripts/verif_omp/check_required_setup.sh b/scripts/verif_omp/check_required_setup.sh new file mode 100755 index 00000000..367530b6 --- /dev/null +++ b/scripts/verif_omp/check_required_setup.sh @@ -0,0 +1,19 @@ +#!/bin/sh + +# take one argument which is the compiler used +# return the required IRPF90_FLAGS for the $1 compiler + +if [ -z "$1" ] +then + echo "Give the compiler in argument" +else + +$1 --version > /dev/null \ +&& $1 -O0 -fopenmp check_omp.f90 \ +&& ./a.out | tail -n 1 + + +# if there is an error or if the compiler is not found +$1 --version > /dev/null || echo 'compiler not found' + +fi diff --git a/scripts/verif_omp/study_omp.sh b/scripts/verif_omp/study_omp.sh new file mode 100755 index 00000000..900d04e1 --- /dev/null +++ b/scripts/verif_omp/study_omp.sh @@ -0,0 +1,30 @@ +#!/bin/sh + +# list of compilers +list_comp="ifort gfortran-7 gfortran-8 gfortran-9" + +# file to store the results +FILE=results.dat + +touch $FILE +rm $FILE + +# Comments +echo "1: omp_set_max_active_levels(5)" >> $FILE +echo "2: omp_set_nested(.True.)" >> $FILE +echo "3: 1 + 2" >> $FILE +echo "" >> $FILE +echo "1 2 3" >> $FILE + +# loop on the comp +for comp in $list_comp +do + $comp --version > /dev/null \ + && $comp -O0 -fopenmp check_omp.f90 \ + && echo $(./a.out | grep "Tests:" | cut -d ":" -f2- ) $(echo " : ") $($comp --version | head -n 1) >> $FILE + +done + +# Display +cat $FILE + diff --git a/scripts/verif_omp/update_comp.sh b/scripts/verif_omp/update_comp.sh new file mode 100755 index 00000000..14b644de --- /dev/null +++ b/scripts/verif_omp/update_comp.sh @@ -0,0 +1,49 @@ +#!/bin/bash + +# Compiler +COMP=$1 + +# Path to file.cfg +config_PATH="../../config/" +END="*.cfg" +CONFIG="/config/" + +#LIST=${config_PATH}${COMP}${END} # without ${QP_ROOT} +LIST=${QP_ROOT}${CONFIG}${COMP}${END} + +if [ -z "$1" ] +then + echo "Give the compiler in argument" +else + + # List of the config files for the compiler + #list_files=$(ls ../../config/$comp*.cfg) #does not give the right list + list_files=${LIST} + echo "Files that will be modified:" + echo $list_files + + # Flags that must be added + FLAGS=$(./check_required_setup.sh $COMP) + + # Add the flags + for file in $list_files + do + echo $file + BASE="IRPF90_FLAGS : --ninja" + ACTUAL=$(grep "$BASE" $file) + + # To have only one time each flag + grep " -DSET_MAX_ACT" $file && ${ACTUAL/" -DSET_MAX"/""} + grep " -DSET_NESTED" $file && ${ACTUAL/" -DSET_NESTED"/""} + SPACE=" " + + NEW=${ACTUAL}${SPACE}${FLAGS} + + # Debug + #echo ${NEW} + + sed "s/${ACTUAL}/${NEW}/" $file + # -i # to change the files + done + +fi diff --git a/src/ao_basis/aos_in_r.irp.f b/src/ao_basis/aos_in_r.irp.f index 7fcb980a..902827eb 100644 --- a/src/ao_basis/aos_in_r.irp.f +++ b/src/ao_basis/aos_in_r.irp.f @@ -12,21 +12,21 @@ double precision function ao_value(i,r) integer :: power_ao(3) double precision :: accu,dx,dy,dz,r2 num_ao = ao_nucl(i) - power_ao(1:3)= ao_power(i,1:3) - center_ao(1:3) = nucl_coord(num_ao,1:3) - dx = (r(1) - center_ao(1)) - dy = (r(2) - center_ao(2)) - dz = (r(3) - center_ao(3)) - r2 = dx*dx + dy*dy + dz*dz - dx = dx**power_ao(1) - dy = dy**power_ao(2) - dz = dz**power_ao(3) +! power_ao(1:3)= ao_power(i,1:3) +! center_ao(1:3) = nucl_coord(num_ao,1:3) +! dx = (r(1) - center_ao(1)) +! dy = (r(2) - center_ao(2)) +! dz = (r(3) - center_ao(3)) +! r2 = dx*dx + dy*dy + dz*dz +! dx = dx**power_ao(1) +! dy = dy**power_ao(2) +! dz = dz**power_ao(3) accu = 0.d0 - do m=1,ao_prim_num(i) - beta = ao_expo_ordered_transp(m,i) - accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2) - enddo +! do m=1,ao_prim_num(i) +! beta = ao_expo_ordered_transp(m,i) +! accu += ao_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2) +! enddo ao_value = accu * dx * dy * dz end diff --git a/src/ao_basis/spherical_to_cartesian.irp.f b/src/ao_basis/spherical_to_cartesian.irp.f index 33a3bc89..336161f8 100644 --- a/src/ao_basis/spherical_to_cartesian.irp.f +++ b/src/ao_basis/spherical_to_cartesian.irp.f @@ -1,7 +1,7 @@ ! Spherical to cartesian transformation matrix obtained with ! Horton (http://theochem.github.com/horton/, 2015) -! First index is the index of the carteisan AO, obtained by ao_power_index +! First index is the index of the cartesian AO, obtained by ao_power_index ! Second index is the index of the spherical AO BEGIN_PROVIDER [ double precision, cart_to_sphe_0, (1,1) ] diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg index a6df8e7a..a6864418 100644 --- a/src/basis/EZFIO.cfg +++ b/src/basis/EZFIO.cfg @@ -37,16 +37,16 @@ doc: Number of primitives in a shell size: (basis.shell_num) interface: ezfio, provider -[shell_prim_index] +[shell_index] type: integer -doc: Index of the first primitive of the shell -size: (basis.shell_num) +doc: Index of the shell for each primitive +size: (basis.prim_num) interface: ezfio, provider [basis_nucleus_index] type: integer -doc: Index of the nucleus on which the shell is centered -size: (nuclei.nucl_num) +doc: Nucleus on which the shell is centered +size: (basis.shell_num) interface: ezfio, provider [prim_normalization_factor] diff --git a/src/basis/basis.irp.f b/src/basis/basis.irp.f index 6a406e28..b750d75a 100644 --- a/src/basis/basis.irp.f +++ b/src/basis/basis.irp.f @@ -30,8 +30,10 @@ BEGIN_PROVIDER [ double precision, shell_normalization_factor , (shell_num) ] powA(3) = 0 norm = 0.d0 - do k=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1 - do j=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1 + do k=1, prim_num + if (shell_index(k) /= i) cycle + do j=1, prim_num + if (shell_index(j) /= i) cycle call overlap_gaussian_xyz(C_A,C_A,prim_expo(j),prim_expo(k), & powA,powA,overlap_x,overlap_y,overlap_z,c,nz) norm = norm+c*prim_coef(j)*prim_coef(k) * prim_normalization_factor(j) * prim_normalization_factor(k) @@ -91,7 +93,8 @@ BEGIN_PROVIDER [ double precision, prim_normalization_factor , (prim_num) ] powA(2) = 0 powA(3) = 0 - do k=shell_prim_index(i),shell_prim_index(i)+shell_prim_num(i)-1 + do k=1, prim_num + if (shell_index(k) /= i) cycle call overlap_gaussian_xyz(C_A,C_A,prim_expo(k),prim_expo(k), & powA,powA,overlap_x,overlap_y,overlap_z,norm,nz) prim_normalization_factor(k) = 1.d0/dsqrt(norm) diff --git a/src/becke_numerical_grid/grid_becke_vector.irp.f b/src/becke_numerical_grid/grid_becke_vector.irp.f index a72200f7..343bd054 100644 --- a/src/becke_numerical_grid/grid_becke_vector.irp.f +++ b/src/becke_numerical_grid/grid_becke_vector.irp.f @@ -58,3 +58,17 @@ END_PROVIDER enddo END_PROVIDER + +BEGIN_PROVIDER [double precision, final_grid_points_transp, (n_points_final_grid,3)] + implicit none + BEGIN_DOC +! Transposed final_grid_points + END_DOC + + integer :: i,j + do j=1,3 + do i=1,n_points_final_grid + final_grid_points_transp(i,j) = final_grid_points(j,i) + enddo + enddo +END_PROVIDER diff --git a/src/bitmask/bitmasks_routines.irp.f b/src/bitmask/bitmasks_routines.irp.f index c34d54dc..9c6f4f0c 100644 --- a/src/bitmask/bitmasks_routines.irp.f +++ b/src/bitmask/bitmasks_routines.irp.f @@ -268,6 +268,21 @@ subroutine print_spindet(string,Nint) end +subroutine print_det_one_dimension(string,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Subroutine to print the content of a determinant using the '+-' notation + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: string(Nint) + character*(2048) :: output(1) + + call bitstring_to_str( output(1), string, Nint ) + print *, trim(output(1)) + +end + logical function is_integer_in_string(bite,string,Nint) use bitmasks implicit none diff --git a/src/cipsi/EZFIO.cfg b/src/cipsi/EZFIO.cfg index 19b45ac1..e01359c5 100644 --- a/src/cipsi/EZFIO.cfg +++ b/src/cipsi/EZFIO.cfg @@ -1,9 +1,3 @@ -[pert_2rdm] -type: logical -doc: If true, computes the one- and two-body rdms with perturbation theory -interface: ezfio,provider,ocaml -default: False - [save_wf_after_selection] type: logical doc: If true, saves the wave function after the selection, before the diagonalization @@ -40,3 +34,9 @@ doc: Maximum number of excitation for beta determinants with respect to the Hart interface: ezfio,ocaml,provider default: -1 +[twice_hierarchy_max] +type: integer +doc: Twice the maximum hierarchy parameter (excitation degree plus half the seniority number). Using -1 selects all determinants +interface: ezfio,ocaml,provider +default: -1 + diff --git a/src/cipsi/NEED b/src/cipsi/NEED index bfbc559a..85d01f79 100644 --- a/src/cipsi/NEED +++ b/src/cipsi/NEED @@ -2,5 +2,4 @@ perturbation zmq mpi iterations -two_body_rdm csf diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index 6e715531..da77a527 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -70,8 +70,8 @@ subroutine run_cipsi do while ( & (N_det < N_det_max) .and. & - (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. & - (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and. & + (sum(abs(pt2_data % pt2(1:N_states)) * state_average_weight(1:N_states)) > pt2_max) .and. & + (sum(abs(pt2_data % variance(1:N_states)) * state_average_weight(1:N_states)) > variance_max) .and. & (correlation_energy_ratio <= correlation_energy_ratio_max) & ) write(*,'(A)') '--------------------------------------------------------------------------------' diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f deleted file mode 100644 index eca8decc..00000000 --- a/src/cipsi/pert_rdm_providers.irp.f +++ /dev/null @@ -1,183 +0,0 @@ - -use bitmasks -use omp_lib - -BEGIN_PROVIDER [ integer(omp_lock_kind), pert_2rdm_lock] - use f77_zmq - implicit none - call omp_init_lock(pert_2rdm_lock) -END_PROVIDER - -BEGIN_PROVIDER [integer, n_orb_pert_rdm] - implicit none - n_orb_pert_rdm = n_act_orb -END_PROVIDER - -BEGIN_PROVIDER [integer, list_orb_reverse_pert_rdm, (mo_num)] - implicit none - list_orb_reverse_pert_rdm = list_act_reverse - -END_PROVIDER - -BEGIN_PROVIDER [integer, list_orb_pert_rdm, (n_orb_pert_rdm)] - implicit none - list_orb_pert_rdm = list_act - -END_PROVIDER - -BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm)] - implicit none - pert_2rdm_provider = 0.d0 - -END_PROVIDER - -subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, psi_det_connection, psi_coef_connection_reverse, n_det_connection) - use bitmasks - use selection_types - implicit none - - integer, intent(in) :: n_det_connection - double precision, intent(in) :: psi_coef_connection_reverse(N_states,n_det_connection) - integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) - integer, intent(in) :: i_generator, sp, h1, h2 - double precision, intent(in) :: mat(N_states, mo_num, mo_num) - logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num) - double precision, intent(in) :: fock_diag_tmp(mo_num) - double precision, intent(in) :: E0(N_states) - type(pt2_type), intent(inout) :: pt2_data - type(selection_buffer), intent(inout) :: buf - logical :: ok - integer :: s1, s2, p1, p2, ib, j, istate, jstate - integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) - double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp, alpha_h_psi, coef(N_states) - double precision, external :: diag_H_mat_elem_fock - double precision :: E_shift - - logical, external :: detEq - double precision, allocatable :: values(:) - integer, allocatable :: keys(:,:) - integer :: nkeys - integer :: sze_buff - sze_buff = 5 * mo_num ** 2 - allocate(keys(4,sze_buff),values(sze_buff)) - nkeys = 0 - if(sp == 3) then - s1 = 1 - s2 = 2 - else - s1 = sp - s2 = sp - end if - call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int) - E_shift = 0.d0 - - if (h0_type == 'CFG') then - j = det_to_configuration(i_generator) - E_shift = psi_det_Hii(i_generator) - psi_configuration_Hii(j) - endif - - do p1=1,mo_num - if(bannedOrb(p1, s1)) cycle - ib = 1 - if(sp /= 3) ib = p1+1 - - do p2=ib,mo_num - -! ----- -! /!\ Generating only single excited determinants doesn't work because a -! determinant generated by a single excitation may be doubly excited wrt -! to a determinant of the future. In that case, the determinant will be -! detected as already generated when generating in the future with a -! double excitation. -! -! if (.not.do_singles) then -! if ((h1 == p1) .or. (h2 == p2)) then -! cycle -! endif -! endif -! -! if (.not.do_doubles) then -! if ((h1 /= p1).and.(h2 /= p2)) then -! cycle -! endif -! endif -! ----- - - if(bannedOrb(p2, s2)) cycle - if(banned(p1,p2)) cycle - - - if( sum(abs(mat(1:N_states, p1, p2))) == 0d0) cycle - call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) - - if (do_only_cas) then - integer, external :: number_of_holes, number_of_particles - if (number_of_particles(det)>0) then - cycle - endif - if (number_of_holes(det)>0) then - cycle - endif - endif - - if (do_ddci) then - logical, external :: is_a_two_holes_two_particles - if (is_a_two_holes_two_particles(det)) then - cycle - endif - endif - - if (do_only_1h1p) then - logical, external :: is_a_1h1p - if (.not.is_a_1h1p(det)) cycle - endif - - - Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) - - sum_e_pert = 0d0 - integer :: degree - call get_excitation_degree(det,HF_bitmask,degree,N_int) - if(degree == 2)cycle - do istate=1,N_states - delta_E = E0(istate) - Hii + E_shift - alpha_h_psi = mat(istate, p1, p2) - val = alpha_h_psi + alpha_h_psi - tmp = dsqrt(delta_E * delta_E + val * val) - if (delta_E < 0.d0) then - tmp = -tmp - endif - e_pert = 0.5d0 * (tmp - delta_E) - coef(istate) = e_pert / alpha_h_psi - print*,e_pert,coef,alpha_h_psi - pt2_data % pt2(istate) += e_pert - pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi - enddo - - do istate=1,N_states - alpha_h_psi = mat(istate, p1, p2) - e_pert = coef(istate) * alpha_h_psi - do jstate=1,N_states - pt2_data % overlap(jstate,jstate) = coef(istate) * coef(jstate) - enddo - - if (weight_selection /= 5) then - ! Energy selection - sum_e_pert = sum_e_pert + e_pert * selection_weight(istate) - - else - ! Variance selection - sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate) - endif - end do - call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff) - - if(sum_e_pert <= buf%mini) then - call add_to_selection_buffer(buf, det, sum_e_pert) - end if - end do - end do - call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) -end - - diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index b366a268..c7cee1ac 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -117,7 +117,6 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) integer(ZMQ_PTR) :: zmq_to_qp_run_socket, zmq_socket_pull integer, intent(in) :: N_in -! integer, intent(inout) :: N_in double precision, intent(in) :: relative_error, E(N_states) type(pt2_type), intent(inout) :: pt2_data, pt2_data_err ! @@ -132,8 +131,8 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted PROVIDE psi_det_hii selection_weight pseudo_sym - PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max + PROVIDE list_act list_inact list_core list_virt list_del seniority_max + PROVIDE excitation_beta_max excitation_alpha_max excitation_max if (h0_type == 'CFG') then PROVIDE psi_configuration_hii det_to_configuration @@ -288,11 +287,11 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) call write_int(6,nproc_target,'Number of threads for PT2') call write_double(6,mem,'Memory (Gb)') - call omp_set_max_active_levels(1) + call set_multiple_levels_omp(.False.) print '(A)', '========== ======================= ===================== ===================== ===========' - print '(A)', ' Samples Energy Variance Norm^2 Seconds' + print '(A)', ' Samples Energy Variance Norm^2 Seconds' print '(A)', '========== ======================= ===================== ===================== ===========' PROVIDE global_selection_buffer @@ -315,14 +314,14 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in) endif !$OMP END PARALLEL call end_parallel_job(zmq_to_qp_run_socket, zmq_socket_pull, 'pt2') - call omp_set_max_active_levels(8) + call set_multiple_levels_omp(.True.) print '(A)', '========== ======================= ===================== ===================== ===========' - do k=1,N_states - pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) - enddo - SOFT_TOUCH pt2_overlap + do k=1,N_states + pt2_overlap(pt2_stoch_istate,k) = pt2_data % overlap(k,pt2_stoch_istate) + enddo + SOFT_TOUCH pt2_overlap enddo FREE pt2_stoch_istate @@ -524,21 +523,21 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ ! 1/(N-1.5) : see Brugger, The American Statistician (23) 4 p. 32 (1969) if(c > 2) then eqt = dabs((pt2_data_S2(t) % pt2(pt2_stoch_istate) / c) - (pt2_data_S(t) % pt2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqt = sqrt(eqt / (dble(c) - 1.5d0)) + eqt = dsqrt(eqt / (dble(c) - 1.5d0)) pt2_data_err % pt2(pt2_stoch_istate) = eqt eqt = dabs((pt2_data_S2(t) % variance(pt2_stoch_istate) / c) - (pt2_data_S(t) % variance(pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqt = sqrt(eqt / (dble(c) - 1.5d0)) + eqt = dsqrt(eqt / (dble(c) - 1.5d0)) pt2_data_err % variance(pt2_stoch_istate) = eqt eqta(:) = dabs((pt2_data_S2(t) % overlap(:,pt2_stoch_istate) / c) - (pt2_data_S(t) % overlap(:,pt2_stoch_istate)/c)**2) ! dabs for numerical stability - eqta(:) = sqrt(eqta(:) / (dble(c) - 1.5d0)) + eqta(:) = dsqrt(eqta(:) / (dble(c) - 1.5d0)) pt2_data_err % overlap(:,pt2_stoch_istate) = eqta(:) if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then time1 = time - print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.4)', c, & + print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1)', c, & pt2_data % pt2(pt2_stoch_istate) +E, & pt2_data_err % pt2(pt2_stoch_istate), & pt2_data % variance(pt2_stoch_istate), & @@ -576,11 +575,11 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_ endif do i=1,n_tasks if(index(i).gt.size(pt2_data_I,1).or.index(i).lt.1)then - print*,'PB !!!' - print*,'If you see this, send a bug report with the following content' - print*,irp_here - print*,'i,index(i),size(pt2_data_I,1) = ',i,index(i),size(pt2_data_I,1) - stop -1 + print*,'PB !!!' + print*,'If you see this, send a bug report with the following content' + print*,irp_here + print*,'i,index(i),size(pt2_data_I,1) = ',i,index(i),size(pt2_data_I,1) + stop -1 endif call pt2_add(pt2_data_I(index(i)),1.d0,pt2_data_task(i)) f(index(i)) -= 1 @@ -843,9 +842,8 @@ END_PROVIDER do t=1, pt2_N_teeth tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t)) if (tooth_width == 0.d0) then - tooth_width = sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1))) + tooth_width = max(1.d-15,sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1)))) endif - ASSERT(tooth_width > 0.d0) do i=pt2_n_0(t)+1, pt2_n_0(t+1) pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width end do diff --git a/src/cipsi/run_pt2_slave.irp.f b/src/cipsi/run_pt2_slave.irp.f index a72d3dbb..30fc7ce0 100644 --- a/src/cipsi/run_pt2_slave.irp.f +++ b/src/cipsi/run_pt2_slave.irp.f @@ -31,12 +31,11 @@ subroutine run_pt2_slave(thread,iproc,energy) double precision, intent(in) :: energy(N_states_diag) integer, intent(in) :: thread, iproc - call run_pt2_slave_large(thread,iproc,energy) -! if (N_det > nproc*(elec_alpha_num * (mo_num-elec_alpha_num))**2) then -! call run_pt2_slave_large(thread,iproc,energy) -! else -! call run_pt2_slave_small(thread,iproc,energy) -! endif + if (N_det > 100000 ) then + call run_pt2_slave_large(thread,iproc,energy) + else + call run_pt2_slave_small(thread,iproc,energy) + endif end subroutine run_pt2_slave_small(thread,iproc,energy) @@ -67,7 +66,6 @@ subroutine run_pt2_slave_small(thread,iproc,energy) double precision, external :: memory_of_double, memory_of_int integer :: bsize ! Size of selection buffers -! logical :: sending allocate(task_id(pt2_n_tasks_max), task(pt2_n_tasks_max)) allocate(pt2_data(pt2_n_tasks_max), i_generator(pt2_n_tasks_max), subset(pt2_n_tasks_max)) @@ -85,7 +83,6 @@ subroutine run_pt2_slave_small(thread,iproc,energy) buffer_ready = .False. n_tasks = 1 -! sending = .False. done = .False. do while (.not.done) @@ -119,14 +116,13 @@ subroutine run_pt2_slave_small(thread,iproc,energy) do k=1,n_tasks call pt2_alloc(pt2_data(k),N_states) b%cur = 0 -!double precision :: time2 -!call wall_time(time2) +! double precision :: time2 +! call wall_time(time2) call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k))) -!call wall_time(time1) -!print *, i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1)) +! call wall_time(time1) +! print *, i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1)) enddo call wall_time(time1) -!print *, '-->', i_generator(1), time1-time0, n_tasks integer, external :: tasks_done_to_taskserver if (tasks_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_tasks) == -1) then @@ -164,6 +160,11 @@ end subroutine subroutine run_pt2_slave_large(thread,iproc,energy) use selection_types use f77_zmq + BEGIN_DOC +! This subroutine can miss important determinants when the PT2 is completely +! computed. It should be called only for large workloads where the PT2 is +! interrupted before the end + END_DOC implicit none double precision, intent(in) :: energy(N_states_diag) @@ -189,8 +190,12 @@ subroutine run_pt2_slave_large(thread,iproc,energy) integer :: bsize ! Size of selection buffers logical :: sending + double precision :: time_shift + PROVIDE global_selection_buffer global_selection_buffer_lock + call random_number(time_shift) + time_shift = time_shift*15.d0 zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() @@ -208,6 +213,9 @@ subroutine run_pt2_slave_large(thread,iproc,energy) sending = .False. done = .False. + double precision :: time0, time1 + call wall_time(time0) + time0 = time0+time_shift do while (.not.done) integer, external :: get_tasks_from_taskserver @@ -234,25 +242,28 @@ subroutine run_pt2_slave_large(thread,iproc,energy) ASSERT (b%N == bsize) endif - double precision :: time0, time1 - call wall_time(time0) call pt2_alloc(pt2_data,N_states) b%cur = 0 call select_connected(i_generator,energy,pt2_data,b,subset,pt2_F(i_generator)) - call wall_time(time1) integer, external :: tasks_done_to_taskserver if (tasks_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_tasks) == -1) then done = .true. endif call sort_selection_buffer(b) + + call wall_time(time1) +! if (time1-time0 > 15.d0) then + call omp_set_lock(global_selection_buffer_lock) + global_selection_buffer%mini = b%mini + call merge_selection_buffers(b,global_selection_buffer) + b%cur=0 + call omp_unset_lock(global_selection_buffer_lock) + call wall_time(time0) +! endif + call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending) - call omp_set_lock(global_selection_buffer_lock) - global_selection_buffer%mini = b%mini - call merge_selection_buffers(b,global_selection_buffer) - b%cur=0 - call omp_unset_lock(global_selection_buffer_lock) - if ( iproc == 1 ) then + if ( iproc == 1 .or. i_generator < 100 .or. done) then call omp_set_lock(global_selection_buffer_lock) call push_pt2_results_async_send(zmq_socket_push, (/i_generator/), (/pt2_data/), global_selection_buffer, (/task_id/), 1,sending) global_selection_buffer%cur = 0 diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index eda9642c..d4f184f3 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -195,7 +195,10 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d integer :: l_a, nmax, idx integer, allocatable :: indices(:), exc_degree(:), iorder(:) - double precision, parameter :: norm_thr = 1.d-16 + + ! Removed to avoid introducing determinants already presents in the wf + !double precision, parameter :: norm_thr = 1.d-16 + allocate (indices(N_det), & exc_degree(max(N_det_alpha_unique,N_det_beta_unique))) @@ -215,10 +218,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d i = psi_bilinear_matrix_rows(l_a) if (nt + exc_degree(i) <= 4) then idx = psi_det_sorted_order(psi_bilinear_matrix_order(l_a)) - if (psi_average_norm_contrib_sorted(idx) > norm_thr) then + ! Removed to avoid introducing determinants already presents in the wf + !if (psi_average_norm_contrib_sorted(idx) > norm_thr) then indices(k) = idx k=k+1 - endif + !endif endif enddo enddo @@ -242,10 +246,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d idx = psi_det_sorted_order( & psi_bilinear_matrix_order( & psi_bilinear_matrix_transp_order(l_a))) - if (psi_average_norm_contrib_sorted(idx) > norm_thr) then + ! Removed to avoid introducing determinants already presents in the wf + !if (psi_average_norm_contrib_sorted(idx) > norm_thr) then indices(k) = idx k=k+1 - endif + !endif endif enddo enddo @@ -464,27 +469,21 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d allocate (fullminilist (N_int, 2, fullinteresting(0)), & minilist (N_int, 2, interesting(0)) ) - if(pert_2rdm)then - allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) - do i=1,fullinteresting(0) - do j = 1, N_states - coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) - enddo - enddo - endif +! if(pert_2rdm)then +! allocate(coef_fullminilist_rev(N_states,fullinteresting(0))) +! do i=1,fullinteresting(0) +! do j = 1, N_states +! coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j) +! enddo +! enddo +! endif do i=1,fullinteresting(0) - do k=1,N_int - fullminilist(k,1,i) = psi_det_sorted(k,1,fullinteresting(i)) - fullminilist(k,2,i) = psi_det_sorted(k,2,fullinteresting(i)) - enddo + fullminilist(:,:,i) = psi_det_sorted(:,:,fullinteresting(i)) enddo do i=1,interesting(0) - do k=1,N_int - minilist(k,1,i) = psi_det_sorted(k,1,interesting(i)) - minilist(k,2,i) = psi_det_sorted(k,2,interesting(i)) - enddo + minilist(:,:,i) = psi_det_sorted(:,:,interesting(i)) enddo do s2=s1,2 @@ -531,19 +530,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting) - if(.not.pert_2rdm)then +! if(.not.pert_2rdm)then call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf) - else - call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) - endif +! else +! call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0)) +! endif end if enddo if(s1 /= s2) monoBdo = .false. enddo deallocate(fullminilist,minilist) - if(pert_2rdm)then - deallocate(coef_fullminilist_rev) - endif +! if(pert_2rdm)then +! deallocate(coef_fullminilist_rev) +! endif enddo enddo deallocate(preinteresting, prefullinteresting, interesting, fullinteresting) @@ -572,6 +571,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d double precision, external :: diag_H_mat_elem_fock double precision :: E_shift double precision :: s_weight(N_states,N_states) + logical, external :: is_in_wavefunction PROVIDE dominant_dets_of_cfgs N_dominant_dets_of_cfgs do jstate=1,N_states do istate=1,N_states @@ -713,6 +713,25 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if (do_cycle) cycle endif + if (twice_hierarchy_max >= 0) then + s = 0 + do k=1,N_int + s = s + popcnt(ieor(det(k,1),det(k,2))) + enddo + if ( mod(s,2)>0 ) stop 'For now, hierarchy CI is defined only for an even number of electrons' + if (excitation_ref == 1) then + call get_excitation_degree(HF_bitmask,det(1,1),degree,N_int) + else if (excitation_ref == 2) then + stop 'For now, hierarchy CI is defined only for a single reference determinant' +! do k=1,N_dominant_dets_of_cfgs +! call get_excitation_degree(dominant_dets_of_cfgs(1,1,k),det(1,1),degree,N_int) +! enddo + endif + integer :: twice_hierarchy + twice_hierarchy = degree + s/2 + if (twice_hierarchy > twice_hierarchy_max) cycle + endif + Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int) w = 0d0 @@ -834,8 +853,27 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d endif end select + + ! To force the inclusion of determinants with a positive pt2 contribution + if (e_pert(istate) > 1d-8) then + w = -huge(1.0) + endif + end do +!!!BEGIN_DEBUG +! ! To check if the pt2 is taking determinants already in the wf +! if (is_in_wavefunction(det(N_int,1),N_int)) then +! print*, 'A determinant contributing to the pt2 is already in' +! print*, 'the wave function:' +! call print_det(det(N_int,1),N_int) +! print*,'contribution to the pt2 for the states:', e_pert(:) +! print*,'error in the filtering in' +! print*, 'cipsi/selection.irp.f sub: selecte_singles_and_doubles' +! print*, 'abort' +! call abort +! endif +!!!END_DEBUG integer(bit_kind) :: occ(N_int,2), n if (h0_type == 'CFG') then @@ -1556,7 +1594,7 @@ subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint) use bitmasks implicit none BEGIN_DOC - ! Gives the inidices(+1) of the bits set to 1 in the bit string + ! Gives the indices(+1) of the bits set to 1 in the bit string END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: string(Nint) diff --git a/src/cipsi/selection_buffer.irp.f b/src/cipsi/selection_buffer.irp.f index 10132086..1f743e0e 100644 --- a/src/cipsi/selection_buffer.irp.f +++ b/src/cipsi/selection_buffer.irp.f @@ -60,6 +60,7 @@ subroutine add_to_selection_buffer(b, det, val) b%val(b%cur) = val if(b%cur == size(b%val)) then call sort_selection_buffer(b) + b%cur = b%cur-1 end if end if end subroutine @@ -86,43 +87,56 @@ subroutine merge_selection_buffers(b1, b2) double precision :: rss double precision, external :: memory_of_double sze = max(size(b1%val), size(b2%val)) - rss = memory_of_double(sze) + 2*N_int*memory_of_double(sze) - call check_mem(rss,irp_here) +! rss = memory_of_double(sze) + 2*N_int*memory_of_double(sze) +! call check_mem(rss,irp_here) allocate(val(sze), detmp(N_int, 2, sze)) i1=1 i2=1 - do i=1,nmwen - if ( (i1 > b1%cur).and.(i2 > b2%cur) ) then - exit - else if (i1 > b1%cur) then - val(i) = b2%val(i2) - detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2) - detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2) - i2=i2+1 - else if (i2 > b2%cur) then - val(i) = b1%val(i1) - detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1) - detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1) - i1=i1+1 - else - if (b1%val(i1) <= b2%val(i2)) then - val(i) = b1%val(i1) - detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1) - detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1) - i1=i1+1 + + select case (N_int) +BEGIN_TEMPLATE + case $case + do i=1,nmwen + if ( (i1 > b1%cur).and.(i2 > b2%cur) ) then + exit + else if (i1 > b1%cur) then + val(i) = b2%val(i2) + detmp(1:$N_int,1,i) = b2%det(1:$N_int,1,i2) + detmp(1:$N_int,2,i) = b2%det(1:$N_int,2,i2) + i2=i2+1 + else if (i2 > b2%cur) then + val(i) = b1%val(i1) + detmp(1:$N_int,1,i) = b1%det(1:$N_int,1,i1) + detmp(1:$N_int,2,i) = b1%det(1:$N_int,2,i1) + i1=i1+1 else - val(i) = b2%val(i2) - detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2) - detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2) - i2=i2+1 + if (b1%val(i1) <= b2%val(i2)) then + val(i) = b1%val(i1) + detmp(1:$N_int,1,i) = b1%det(1:$N_int,1,i1) + detmp(1:$N_int,2,i) = b1%det(1:$N_int,2,i1) + i1=i1+1 + else + val(i) = b2%val(i2) + detmp(1:$N_int,1,i) = b2%det(1:$N_int,1,i2) + detmp(1:$N_int,2,i) = b2%det(1:$N_int,2,i2) + i2=i2+1 + endif endif - endif - enddo + enddo + do i=nmwen+1,b2%N + val(i) = 0.d0 +! detmp(1:$N_int,1,i) = 0_bit_kind +! detmp(1:$N_int,2,i) = 0_bit_kind + enddo +SUBST [ case, N_int ] +(1); 1;; +(2); 2;; +(3); 3;; +(4); 4;; +default; N_int;; +END_TEMPLATE + end select deallocate(b2%det, b2%val) - do i=nmwen+1,b2%N - val(i) = 0.d0 - detmp(1:N_int,1:2,i) = 0_bit_kind - enddo b2%det => detmp b2%val => val b2%mini = min(b2%mini,b2%val(b2%N)) @@ -144,8 +158,8 @@ subroutine sort_selection_buffer(b) double precision :: rss double precision, external :: memory_of_double, memory_of_int - rss = memory_of_int(b%cur) + 2*N_int*memory_of_double(size(b%det,3)) - call check_mem(rss,irp_here) +! rss = memory_of_int(b%cur) + 2*N_int*memory_of_double(size(b%det,3)) +! call check_mem(rss,irp_here) allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3))) do i=1,b%cur iorder(i) = i @@ -225,14 +239,14 @@ subroutine make_selection_buffer_s2(b) endif dup = .True. do k=1,N_int - if ( (tmp_array(k,1,i) /= tmp_array(k,1,j)) & - .or. (tmp_array(k,2,i) /= tmp_array(k,2,j)) ) then + if ( (tmp_array(k,1,i) /= tmp_array(k,1,j)) .or. & + (tmp_array(k,2,i) /= tmp_array(k,2,j)) ) then dup = .False. exit endif enddo if (dup) then - val(i) = max(val(i), val(j)) + val(i) = min(val(i), val(j)) duplicate(j) = .True. endif j+=1 @@ -282,9 +296,6 @@ subroutine make_selection_buffer_s2(b) call configuration_to_dets_size(o(1,1,i),sze,elec_alpha_num,N_int) n_d = n_d + sze if (n_d > b%cur) then -! if (n_d - b%cur > b%cur - n_d + sze) then -! n_d = n_d - sze -! endif exit endif enddo @@ -329,10 +340,11 @@ subroutine remove_duplicates_in_selection_buffer(b) integer(bit_kind), allocatable :: tmp_array(:,:,:) logical, allocatable :: duplicate(:) - n_d = b%cur logical :: found_duplicates double precision :: rss double precision, external :: memory_of_double + + n_d = b%cur rss = (4*N_int+4)*memory_of_double(n_d) call check_mem(rss,irp_here) diff --git a/src/cipsi/selection_weight.irp.f b/src/cipsi/selection_weight.irp.f index 3c09e59a..756c65a1 100644 --- a/src/cipsi/selection_weight.irp.f +++ b/src/cipsi/selection_weight.irp.f @@ -38,11 +38,11 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st) avg = sum(pt2(1:N_st)) / dble(N_st) + 1.d-32 ! Avoid future division by zero - dt = 8.d0 !* selection_factor + dt = 4.d0 !* selection_factor do k=1,N_st - element = exp(dt*(pt2(k)/avg - 1.d0)) - element = min(2.0d0 , element) - element = max(0.5d0 , element) + element = pt2(k) !exp(dt*(pt2(k)/avg - 1.d0)) +! element = min(2.0d0 , element) +! element = max(0.5d0 , element) pt2_match_weight(k) *= element enddo @@ -50,9 +50,9 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st) avg = sum(variance(1:N_st)) / dble(N_st) + 1.d-32 ! Avoid future division by zero do k=1,N_st - element = exp(dt*(variance(k)/avg -1.d0)) - element = min(2.0d0 , element) - element = max(0.5d0 , element) + element = variance(k) ! exp(dt*(variance(k)/avg -1.d0)) +! element = min(2.0d0 , element) +! element = max(0.5d0 , element) variance_match_weight(k) *= element enddo @@ -62,6 +62,9 @@ subroutine update_pt2_and_variance_weights(pt2_data, N_st) variance_match_weight(:) = 1.d0 endif + pt2_match_weight(:) = pt2_match_weight(:)/sum(pt2_match_weight(:)) + variance_match_weight(:) = variance_match_weight(:)/sum(variance_match_weight(:)) + threshold_davidson_pt2 = min(1.d-6, & max(threshold_davidson, 1.e-1 * PT2_relative_error * minval(abs(pt2(1:N_states)))) ) @@ -87,7 +90,7 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] selection_weight(1:N_states) = c0_weight(1:N_states) case (2) - print *, 'Using pt2-matching weight in selection' + print *, 'Using PT2-matching weight in selection' selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) print *, '# PT2 weight ', real(pt2_match_weight(:),4) @@ -97,7 +100,7 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] print *, '# var weight ', real(variance_match_weight(:),4) case (4) - print *, 'Using variance- and pt2-matching weights in selection' + print *, 'Using variance- and PT2-matching weights in selection' selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) print *, '# PT2 weight ', real(pt2_match_weight(:),4) print *, '# var weight ', real(variance_match_weight(:),4) @@ -112,7 +115,7 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] selection_weight(1:N_states) = c0_weight(1:N_states) case (7) - print *, 'Input weights multiplied by variance- and pt2-matching' + print *, 'Input weights multiplied by variance- and PT2-matching' selection_weight(1:N_states) = c0_weight(1:N_states) * sqrt(variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)) * state_average_weight(1:N_states) print *, '# PT2 weight ', real(pt2_match_weight(:),4) print *, '# var weight ', real(variance_match_weight(:),4) @@ -128,6 +131,7 @@ BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] print *, '# var weight ', real(variance_match_weight(:),4) end select + selection_weight(:) = selection_weight(:)/sum(selection_weight(:)) print *, '# Total weight ', real(selection_weight(:),4) END_PROVIDER diff --git a/src/cipsi/slave_cipsi.irp.f b/src/cipsi/slave_cipsi.irp.f index 510c667b..f96aaa6a 100644 --- a/src/cipsi/slave_cipsi.irp.f +++ b/src/cipsi/slave_cipsi.irp.f @@ -4,7 +4,7 @@ subroutine run_slave_cipsi ! Helper program for distributed parallelism END_DOC - call omp_set_max_active_levels(1) + call set_multiple_levels_omp(.False.) distributed_davidson = .False. read_wf = .False. SOFT_TOUCH read_wf distributed_davidson @@ -171,9 +171,9 @@ subroutine run_slave_main call write_double(6,(t1-t0),'Broadcast time') !--- - call omp_set_max_active_levels(8) + call set_multiple_levels_omp(.True.) call davidson_slave_tcp(0) - call omp_set_max_active_levels(1) + call set_multiple_levels_omp(.False.) print *, mpi_rank, ': Davidson done' !--- @@ -311,7 +311,7 @@ subroutine run_slave_main if (mpi_master) then print *, 'Running PT2' endif - !$OMP PARALLEL PRIVATE(i) NUM_THREADS(nproc_target+1) + !$OMP PARALLEL PRIVATE(i) NUM_THREADS(nproc_target) i = omp_get_thread_num() call run_pt2_slave(0,i,pt2_e0_denominator) !$OMP END PARALLEL diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 781fcda6..5fc9db0f 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -69,8 +69,8 @@ subroutine run_stochastic_cipsi do while ( & (N_det < N_det_max) .and. & - (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and. & - (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and. & + (sum(abs(pt2_data % pt2(1:N_states)) * state_average_weight(1:N_states)) > pt2_max) .and. & + (sum(abs(pt2_data % variance(1:N_states)) * state_average_weight(1:N_states)) > variance_max) .and. & (correlation_energy_ratio <= correlation_energy_ratio_max) & ) write(*,'(A)') '--------------------------------------------------------------------------------' diff --git a/src/cipsi/update_2rdm.irp.f b/src/cipsi/update_2rdm.irp.f deleted file mode 100644 index 260c48fd..00000000 --- a/src/cipsi/update_2rdm.irp.f +++ /dev/null @@ -1,223 +0,0 @@ -use bitmasks - -subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: n_det_connection,sze_buff - double precision, intent(in) :: coef(N_states) - integer(bit_kind), intent(in) :: det(N_int,2) - integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) - double precision, intent(in) :: psi_coef_connection_reverse(N_states,n_det_connection) - integer, intent(inout) :: keys(4,sze_buff),nkeys - double precision, intent(inout) :: values(sze_buff) - integer :: i,j - integer :: exc(0:2,2,2) - integer :: degree - double precision :: phase, contrib - do i = 1, n_det_connection - call get_excitation(det,psi_det_connection(1,1,i),exc,degree,phase,N_int) - if(degree.gt.2)cycle - contrib = 0.d0 - do j = 1, N_states - contrib += state_average_weight(j) * psi_coef_connection_reverse(j,i) * phase * coef(j) - enddo - ! case of single excitations - if(degree == 1)then - if (nkeys + 6 * elec_alpha_num .ge. sze_buff)then - call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) - nkeys = 0 - endif - call update_buffer_single_exc_rdm(det,psi_det_connection(1,1,i),exc,phase,contrib,nkeys,keys,values,sze_buff) - else - !! case of double excitations - ! if (nkeys + 4 .ge. sze_buff)then - ! call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) - ! nkeys = 0 - ! endif - ! call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) - endif - enddo -!call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock) -!nkeys = 0 - -end - -subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: sze_buff - integer(bit_kind), intent(in) :: det1(N_int,2) - integer(bit_kind), intent(in) :: det2(N_int,2) - integer,intent(in) :: exc(0:2,2,2) - double precision,intent(in) :: phase, contrib - integer, intent(inout) :: nkeys, keys(4,sze_buff) - double precision, intent(inout):: values(sze_buff) - - integer :: occ(N_int*bit_kind_size,2) - integer :: n_occ_ab(2),ispin,other_spin - integer :: h1,h2,p1,p2,i - call bitstring_to_list_ab(det1, occ, n_occ_ab, N_int) - - if (exc(0,1,1) == 1) then - ! Mono alpha - h1 = exc(1,1,1) - p1 = exc(1,2,1) - ispin = 1 - other_spin = 2 - else - ! Mono beta - h1 = exc(1,1,2) - p1 = exc(1,2,2) - ispin = 2 - other_spin = 1 - endif - if(list_orb_reverse_pert_rdm(h1).lt.0)return - h1 = list_orb_reverse_pert_rdm(h1) - if(list_orb_reverse_pert_rdm(p1).lt.0)return - p1 = list_orb_reverse_pert_rdm(p1) - !update the alpha/beta part - do i = 1, n_occ_ab(other_spin) - h2 = occ(i,other_spin) - if(list_orb_reverse_pert_rdm(h2).lt.0)return - h2 = list_orb_reverse_pert_rdm(h2) - - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = h2 - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = h2 - keys(4,nkeys) = p1 - enddo - !update the same spin part -!do i = 1, n_occ_ab(ispin) -! h2 = occ(i,ispin) -! if(list_orb_reverse_pert_rdm(h2).lt.0)return -! h2 = list_orb_reverse_pert_rdm(h2) - -! nkeys += 1 -! values(nkeys) = 0.5d0 * contrib * phase -! keys(1,nkeys) = h1 -! keys(2,nkeys) = h2 -! keys(3,nkeys) = p1 -! keys(4,nkeys) = h2 - -! nkeys += 1 -! values(nkeys) = - 0.5d0 * contrib * phase -! keys(1,nkeys) = h1 -! keys(2,nkeys) = h2 -! keys(3,nkeys) = h2 -! keys(4,nkeys) = p1 -! -! nkeys += 1 -! values(nkeys) = 0.5d0 * contrib * phase -! keys(1,nkeys) = h2 -! keys(2,nkeys) = h1 -! keys(3,nkeys) = h2 -! keys(4,nkeys) = p1 - -! nkeys += 1 -! values(nkeys) = - 0.5d0 * contrib * phase -! keys(1,nkeys) = h2 -! keys(2,nkeys) = h1 -! keys(3,nkeys) = p1 -! keys(4,nkeys) = h2 -!enddo - -end - -subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff) - implicit none - integer, intent(in) :: sze_buff - integer,intent(in) :: exc(0:2,2,2) - double precision,intent(in) :: phase, contrib - integer, intent(inout) :: nkeys, keys(4,sze_buff) - double precision, intent(inout):: values(sze_buff) - integer :: h1,h2,p1,p2 - - if (exc(0,1,1) == 1) then - ! Double alpha/beta - h1 = exc(1,1,1) - h2 = exc(1,1,2) - p1 = exc(1,2,1) - p2 = exc(1,2,2) - ! check if the orbitals involved are within the orbital range - if(list_orb_reverse_pert_rdm(h1).lt.0)return - h1 = list_orb_reverse_pert_rdm(h1) - if(list_orb_reverse_pert_rdm(h2).lt.0)return - h2 = list_orb_reverse_pert_rdm(h2) - if(list_orb_reverse_pert_rdm(p1).lt.0)return - p1 = list_orb_reverse_pert_rdm(p1) - if(list_orb_reverse_pert_rdm(p2).lt.0)return - p2 = list_orb_reverse_pert_rdm(p2) - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = p1 - keys(2,nkeys) = p2 - keys(3,nkeys) = h1 - keys(4,nkeys) = h2 - - else - if (exc(0,1,1) == 2) then - ! Double alpha/alpha - h1 = exc(1,1,1) - h2 = exc(2,1,1) - p1 = exc(1,2,1) - p2 = exc(2,2,1) - else if (exc(0,1,2) == 2) then - ! Double beta - h1 = exc(1,1,2) - h2 = exc(2,1,2) - p1 = exc(1,2,2) - p2 = exc(2,2,2) - endif - ! check if the orbitals involved are within the orbital range - if(list_orb_reverse_pert_rdm(h1).lt.0)return - h1 = list_orb_reverse_pert_rdm(h1) - if(list_orb_reverse_pert_rdm(h2).lt.0)return - h2 = list_orb_reverse_pert_rdm(h2) - if(list_orb_reverse_pert_rdm(p1).lt.0)return - p1 = list_orb_reverse_pert_rdm(p1) - if(list_orb_reverse_pert_rdm(p2).lt.0)return - p2 = list_orb_reverse_pert_rdm(p2) - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - - nkeys += 1 - values(nkeys) = - 0.5d0 * contrib * phase - keys(1,nkeys) = h1 - keys(2,nkeys) = h2 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = 0.5d0 * contrib * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p2 - keys(4,nkeys) = p1 - - nkeys += 1 - values(nkeys) = - 0.5d0 * contrib * phase - keys(1,nkeys) = h2 - keys(2,nkeys) = h1 - keys(3,nkeys) = p1 - keys(4,nkeys) = p2 - endif - -end - - diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index 58630709..1bfe87c0 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -22,7 +22,7 @@ subroutine ZMQ_selection(N_in, pt2_data) PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max - PROVIDE pert_2rdm excitation_beta_max excitation_alpha_max excitation_max + PROVIDE excitation_beta_max excitation_alpha_max excitation_max call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection') diff --git a/src/cis/EZFIO.cfg b/src/cis/EZFIO.cfg index 7e0eeb03..955d1bef 100644 --- a/src/cis/EZFIO.cfg +++ b/src/cis/EZFIO.cfg @@ -5,4 +5,3 @@ interface: ezfio size: (determinants.n_states) - diff --git a/src/cis/cis.irp.f b/src/cis/cis.irp.f index acec29c2..f72197c2 100644 --- a/src/cis/cis.irp.f +++ b/src/cis/cis.irp.f @@ -62,6 +62,7 @@ subroutine run else call H_apply_cis endif + print*,'' print *, 'N_det = ', N_det print*,'******************************' print *, 'Energies of the states:' @@ -69,16 +70,18 @@ subroutine run print *, i, CI_energy(i) enddo if (N_states > 1) then - print*,'******************************' - print*,'Excitation energies ' + print*,'' + print*,'******************************************************' + print*,'Excitation energies (au) (eV)' do i = 2, N_states - print*, i ,CI_energy(i) - CI_energy(1) + print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1))/0.0367502d0 enddo + print*,'' endif call ezfio_set_cis_energy(CI_energy) psi_coef = ci_eigenvectors SOFT_TOUCH psi_coef - call save_wavefunction_truncated(1.d-12) + call save_wavefunction_truncated(save_threshold) end diff --git a/src/cisd/cisd.irp.f b/src/cisd/cisd.irp.f index 6c55e2ff..fca3b10e 100644 --- a/src/cisd/cisd.irp.f +++ b/src/cisd/cisd.irp.f @@ -63,7 +63,7 @@ subroutine run endif psi_coef = ci_eigenvectors SOFT_TOUCH psi_coef - call save_wavefunction + call save_wavefunction_truncated(save_threshold) call ezfio_set_cisd_energy(CI_energy) do i = 1,N_states diff --git a/src/csf/configurations.irp.f b/src/csf/configurations.irp.f index 8e2a513c..ce5d48ab 100644 --- a/src/csf/configurations.irp.f +++ b/src/csf/configurations.irp.f @@ -779,6 +779,7 @@ subroutine binary_search_cfg(cfgInp,addcfg) end subroutine BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ] +&BEGIN_PROVIDER [ integer, psi_configuration_n_det, (N_configuration) ] &BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ] implicit none @@ -867,6 +868,29 @@ end subroutine enddo deallocate(dets, old_order) + integer :: ndet_conf + do i = 1, N_configuration + ndet_conf = psi_configuration_to_psi_det(2,i) - psi_configuration_to_psi_det(1,i) + 1 + psi_configuration_n_det(i) = ndet_conf + enddo END_PROVIDER + +BEGIN_PROVIDER [ integer, n_elec_alpha_for_psi_configuration, (N_configuration)] + implicit none + integer :: i,j,k,l + integer(bit_kind) :: det_tmp(N_int,2),det_alpha(N_int) + n_elec_alpha_for_psi_configuration = 0 + do i = 1, N_configuration + j = psi_configuration_to_psi_det(2,i) + det_tmp(:,:) = psi_det(:,:,j) + k = 0 + do l = 1, N_int + det_alpha(N_int) = iand(det_tmp(l,1),psi_configuration(l,1,i)) + k += popcnt(det_alpha(l)) + enddo + n_elec_alpha_for_psi_configuration(i) = k + enddo + +END_PROVIDER diff --git a/src/csf/conversion.irp.f b/src/csf/conversion.irp.f index fecc6123..75f6e539 100644 --- a/src/csf/conversion.irp.f +++ b/src/csf/conversion.irp.f @@ -1,3 +1,15 @@ +BEGIN_PROVIDER [ double precision, psi_csf_coef, (N_csf, N_states) ] + implicit none + BEGIN_DOC + ! Wafe function in CSF basis + END_DOC + + double precision, allocatable :: buffer(:,:) + allocate ( buffer(N_det, N_states) ) + buffer(1:N_det, 1:N_states) = psi_coef(1:N_det, 1:N_states) + call convertWFfromDETtoCSF(N_states, buffer, psi_csf_coef) +END_PROVIDER + subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) use cfunctions use bitmasks @@ -26,6 +38,8 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) integer s, bfIcfg integer countcsf + integer MS + MS = elec_alpha_num-elec_beta_num countcsf = 0 phasedet = 1.0d0 do i = 1,N_configuration @@ -44,12 +58,17 @@ subroutine convertWFfromDETtoCSF(N_st,psi_coef_det_in, psi_coef_cfg_out) enddo enddo - s = 0 + s = 0 ! s == total number of SOMOs do k=1,N_int if (psi_configuration(k,1,i) == 0_bit_kind) cycle s = s + popcnt(psi_configuration(k,1,i)) enddo - bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1)))) + + if(iand(s,1) .EQ. 0) then + bfIcfg = max(1,nint((binom(s,s/2)-binom(s,(s/2)+1)))) + else + bfIcfg = max(1,nint((binom(s,(s+1)/2)-binom(s,((s+1)/2)+1)))) + endif ! perhaps blocking with CFGs of same seniority ! can be more efficient diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index edf6c400..5aaba9a3 100644 --- a/src/csf/sigma_vector.irp.f +++ b/src/csf/sigma_vector.irp.f @@ -1,9 +1,12 @@ real*8 function logabsgamma(x) implicit none real*8, intent(in) :: x - logabsgamma = log(abs(gamma(x))) + logabsgamma = 1.d32 ! Avoid floating point exception + if (x>0.d0) then + logabsgamma = log(abs(gamma(x))) + endif end function logabsgamma - + BEGIN_PROVIDER [ integer, NSOMOMax] &BEGIN_PROVIDER [ integer, NCSFMax] &BEGIN_PROVIDER [ integer*8, NMO] @@ -48,11 +51,24 @@ if(cfg_seniority_index(i+2) > ncfgpersomo) then ncfgpersomo = cfg_seniority_index(i+2) else - k = 0 - do while(cfg_seniority_index(i+2+k) < ncfgpersomo) - k = k + 2 - ncfgpersomo = cfg_seniority_index(i+2+k) + ! l = i+k+2 + ! Loop over l with a constraint to ensure that l <= size(cfg_seniority_index,1)-1 + ! Old version commented just below + do l = min(size(cfg_seniority_index,1)-1, i+2), size(cfg_seniority_index,1)-1, 2 + if (cfg_seniority_index(l) >= ncfgpersomo) then + ncfgpersomo = cfg_seniority_index(l) + endif enddo + !k = 0 + !if ((i+2+k) < size(cfg_seniority_index,1)) then + ! do while(cfg_seniority_index(i+2+k) < ncfgpersomo) + ! k = k + 2 + ! if ((i+2+k) >= size(cfg_seniority_index,1)) then + ! exit + ! endif + ! ncfgpersomo = cfg_seniority_index(i+2+k) + ! enddo + !endif endif endif ncfg = ncfgpersomo - ncfgprev @@ -62,34 +78,33 @@ dimcsfpercfg = 2 else if(iand(MS,1) .EQ. 0) then - !dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1)))) - binom1 = dexp(logabsgamma(1.0d0*(i+1)) & - - logabsgamma(1.0d0*((i/2)+1)) & - - logabsgamma(1.0d0*(i-((i/2))+1))); - binom2 = dexp(logabsgamma(1.0d0*(i+1)) & - - logabsgamma(1.0d0*(((i/2)+1)+1)) & - - logabsgamma(1.0d0*(i-((i/2)+1)+1))); - dimcsfpercfg = max(1,nint(binom1 - binom2)) + dimcsfpercfg = max(1,nint((binom(i,i/2)-binom(i,i/2+1)))) else - !dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2)))) - binom1 = dexp(logabsgamma(1.0d0*(i+1)) & - - logabsgamma(1.0d0*(((i+1)/2)+1)) & - - logabsgamma(1.0d0*(i-(((i+1)/2))+1))); - binom2 = dexp(logabsgamma(1.0d0*(i+1)) & - - logabsgamma(1.0d0*((((i+3)/2)+1)+1)) & - - logabsgamma(1.0d0*(i-(((i+3)/2)+1)+1))); - dimcsfpercfg = max(1,nint(binom1 - binom2)) + dimcsfpercfg = max(1,nint((binom(i,(i+1)/2)-binom(i,(i+3)/2)))) endif endif n_CSF += ncfg * dimcsfpercfg if(cfg_seniority_index(i+2) > ncfgprev) then ncfgprev = cfg_seniority_index(i+2) else - k = 0 - do while(cfg_seniority_index(i+2+k) < ncfgprev) - k = k + 2 - ncfgprev = cfg_seniority_index(i+2+k) + ! l = i+k+2 + ! Loop over l with a constraint to ensure that l <= size(cfg_seniority_index,1)-1 + ! Old version commented just below + do l = min(size(cfg_seniority_index,1)-1, i+2), size(cfg_seniority_index,1)-1, 2 + if (cfg_seniority_index(l) >= ncfgprev) then + ncfgprev = cfg_seniority_index(l) + endif enddo + !k = 0 + !if ((i+2+k) < size(cfg_seniority_index,1)) then + ! do while(cfg_seniority_index(i+2+k) < ncfgprev) + ! k = k + 2 + ! if ((i+2+k) >= size(cfg_seniority_index,1)) then + ! exit + ! endif + ! ncfgprev = cfg_seniority_index(i+2+k) + ! enddo + !endif endif enddo END_PROVIDER diff --git a/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f b/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f new file mode 100644 index 00000000..243e9995 --- /dev/null +++ b/src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f @@ -0,0 +1,451 @@ + +subroutine davidson_general_ext_rout_diag_dressed(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc) + use mmap_module + implicit none + BEGIN_DOC + ! Generic Davidson diagonalization with ONE DIAGONAL DRESSING OPERATOR + ! + ! Dress_jj : DIAGONAL DRESSING of the Hamiltonian + ! + ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson + ! + ! u_in : guess coefficients on the various states. Overwritten on exit + ! + ! sze : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! N_st_diag_in : Number of states in which H is diagonalized. Assumed > sze + ! + ! Initial guess vectors are not necessarily orthonormal + ! + ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output) + END_DOC + integer, intent(in) :: sze, N_st, N_st_diag_in + double precision, intent(in) :: H_jj(sze),Dress_jj(sze) + double precision, intent(inout) :: u_in(sze,N_st_diag_in) + double precision, intent(out) :: energies(N_st) + external hcalc + + integer :: iter, N_st_diag + integer :: i,j,k,l,m + logical, intent(inout) :: converged + + double precision, external :: u_dot_v, u_dot_u + + integer :: k_pairs, kl + + integer :: iter2, itertot + double precision, allocatable :: y(:,:), h(:,:), lambda(:) + double precision, allocatable :: residual_norm(:) + character*(16384) :: write_buffer + double precision :: to_print(2,N_st) + double precision :: cpu, wall + integer :: shift, shift2, itermax, istate + double precision :: r1, r2, alpha + integer :: nproc_target + integer :: order(N_st_diag_in) + double precision :: cmax + double precision, allocatable :: U(:,:), overlap(:,:)!, S_d(:,:) + double precision, pointer :: W(:,:) + logical :: disk_based + double precision :: energy_shift(N_st_diag_in*davidson_sze_max) + + include 'constants.include.F' + + N_st_diag = N_st_diag_in + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda + if (N_st_diag*3 > sze) then + print *, 'error in Davidson :' + print *, 'Increase n_det_max_full to ', N_st_diag*3 + stop -1 + endif + + itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1 + itertot = 0 + + if (state_following) then + allocate(overlap(N_st_diag*itermax, N_st_diag*itermax)) + else + allocate(overlap(1,1)) ! avoid 'if' for deallocate + endif + overlap = 0.d0 + + provide threshold_davidson !nthreads_davidson + call write_time(6) + write(6,'(A)') '' + write(6,'(A)') 'Davidson Diagonalization' + write(6,'(A)') '------------------------' + write(6,'(A)') '' + + ! Find max number of cores to fit in memory + ! ----------------------------------------- + + nproc_target = nproc + double precision :: rss + integer :: maxab + maxab = sze + + m=1 + disk_based = .False. + call resident_memory(rss) + do + r1 = 8.d0 * &! bytes + ( dble(sze)*(N_st_diag*itermax) &! U + + 1.d0*dble(sze*m)*(N_st_diag*itermax) &! W + + 2.0d0*(N_st_diag*itermax)**2 &! h,y + + 2.d0*(N_st_diag*itermax) &! s2,lambda + + 1.d0*(N_st_diag) &! residual_norm + ! In H_S2_u_0_nstates_zmq + + 3.d0*(N_st_diag*N_det) &! u_t, v_t, s_t on collector + + 3.d0*(N_st_diag*N_det) &! u_t, v_t, s_t on slave + + 0.5d0*maxab &! idx0 in H_S2_u_0_nstates_openmp_work_* + + nproc_target * &! In OMP section + ( 1.d0*(N_int*maxab) &! buffer + + 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx + ) / 1024.d0**3 + + if (nproc_target == 0) then + call check_mem(r1,irp_here) + nproc_target = 1 + exit + endif + + if (r1+rss < qp_max_mem) then + exit + endif + + if (itermax > 4) then + itermax = itermax - 1 + else if (m==1.and.disk_based_davidson) then + m=0 + disk_based = .True. + itermax = 6 + else + nproc_target = nproc_target - 1 + endif + + enddo + nthreads_davidson = nproc_target + TOUCH nthreads_davidson + call write_int(6,N_st,'Number of states') + call write_int(6,N_st_diag,'Number of states in diagonalization') + call write_int(6,sze,'Number of basis functions') + call write_int(6,nproc_target,'Number of threads for diagonalization') + call write_double(6, r1, 'Memory(Gb)') + if (disk_based) then + print *, 'Using swap space to reduce RAM' + endif + + double precision, allocatable :: H_jj_tmp(:) + ASSERT (N_st > 0) + ASSERT (sze > 0) + allocate(H_jj_tmp(sze)) + + do i=1,sze + H_jj_tmp(i) = H_jj(i) + Dress_jj(i) + enddo + + !--------------- + + write(6,'(A)') '' + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = 'Iter' + do i=1,N_st + write_buffer = trim(write_buffer)//' Energy Residual ' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + + + allocate(W(sze,N_st_diag*itermax)) + + allocate( & + ! Large + U(sze,N_st_diag*itermax), & + ! Small + h(N_st_diag*itermax,N_st_diag*itermax), & + y(N_st_diag*itermax,N_st_diag*itermax), & + residual_norm(N_st_diag), & + lambda(N_st_diag*itermax)) + + h = 0.d0 + U = 0.d0 + y = 0.d0 + + + ASSERT (N_st > 0) + ASSERT (N_st_diag >= N_st) + ASSERT (sze > 0) + + ! Davidson iterations + ! =================== + + converged = .False. + + ! Initialize from N_st to N_st_diat with gaussian random numbers + ! to be sure to have overlap with any eigenvectors + do k=N_st+1,N_st_diag + u_in(k,k) = 10.d0 + do i=1,sze + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + u_in(i,k) = r1*dcos(r2) + enddo + enddo + ! Normalize all states + do k=1,N_st_diag + call normalize(u_in(1,k),sze) + enddo + + ! Copy from the guess input "u_in" to the working vectors "U" + do k=1,N_st_diag + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + + do while (.not.converged) + itertot = itertot+1 + if (itertot == 8) then + exit + endif + + do iter=1,itermax-1 + + shift = N_st_diag*(iter-1) + shift2 = N_st_diag*iter + + if ((iter > 1).or.(itertot == 1)) then + ! Compute |W_k> = \sum_i |i> + ! ----------------------------------- + + ! Gram-Schmidt to orthogonalize all new guess with the previous vectors + call ortho_qr(U,size(U,1),sze,shift2) + call ortho_qr(U,size(U,1),sze,shift2) + ! it does W = H U with W(sze,N_st_diag),U(sze,N_st_diag) + ! where sze is the size of the vector, N_st_diag is the number of states + call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag,sze) + ! Compute then the DIAGONAL PART OF THE DRESSING + ! += Dress_jj(i) * + call dressing_diag_uv(W(1,shift+1),U(1,shift+1),Dress_jj,N_st_diag_in,sze) + else + ! Already computed in update below + continue + endif + + ! Compute h_kl = = + ! ------------------------------------------- + + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), W, size(W,1), & + 0.d0, h, size(h,1)) + + ! Diagonalize h y = lambda y + ! --------------- + + call lapack_diag(lambda,y,h,size(h,1),shift2) + + if (state_following) then + + overlap = -1.d0 + do k=1,shift2 + do i=1,shift2 + overlap(k,i) = dabs(y(k,i)) + enddo + enddo + do k=1,N_st + cmax = -1.d0 + do i=1,N_st + if (overlap(i,k) > cmax) then + cmax = overlap(i,k) + order(k) = i + endif + enddo + do i=1,N_st_diag + overlap(order(k),i) = -1.d0 + enddo + enddo + overlap = y + do k=1,N_st + l = order(k) + if (k /= l) then + y(1:shift2,k) = overlap(1:shift2,l) + endif + enddo + do k=1,N_st + overlap(k,1) = lambda(k) + enddo + do k=1,N_st + l = order(k) + if (k /= l) then + lambda(k) = overlap(l,1) + endif + enddo + + endif + + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1)) + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1)) + + ! Compute residual vector and davidson step + ! ----------------------------------------- + + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) + do k=1,N_st_diag + do i=1,sze + U(i,shift2+k) = & + (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & + /max(H_jj_tmp(i) - lambda (k),1.d-2) + enddo + + if (k <= N_st) then + residual_norm(k) = u_dot_u(U(1,shift2+k),sze) + to_print(1,k) = lambda(k) + to_print(2,k) = residual_norm(k) + endif + enddo + !$OMP END PARALLEL DO + + + if ((itertot>1).and.(iter == 1)) then + !don't print + continue + else + write(*,'(1X,I3,1X,100(1X,F16.10,1X,F11.6,1X,E11.3))') iter-1, to_print(1:2,1:N_st) + endif + + ! Check convergence + if (iter > 1) then + converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson + endif + + + do k=1,N_st + if (residual_norm(k) > 1.e8) then + print *, 'Davidson failed' + stop -1 + endif + enddo + if (converged) then + exit + endif + + logical, external :: qp_stop + if (qp_stop()) then + converged = .True. + exit + endif + + + enddo + + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & + W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag + do i=1,sze + W(i,k) = u_in(i,k) + enddo + enddo + + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & + U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + call ortho_qr(U,size(U,1),sze,N_st_diag) + call ortho_qr(U,size(U,1),sze,N_st_diag) + do j=1,N_st_diag + k=1 + do while ((k sze + ! + ! Initial guess vectors are not necessarily orthonormal + ! + ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output) + END_DOC + implicit none + integer, intent(in) :: sze, N_st, N_st_diag, idx_dress + double precision, intent(in) :: H_jj(sze),Dress_jj(sze),Dressing_vec(sze,N_st) + double precision, intent(inout) :: u_in(sze,N_st_diag) + double precision, intent(out) :: energies(N_st_diag) + logical, intent(out) :: converged + external hcalc + + double precision, allocatable :: H_jj_tmp(:) + ASSERT (N_st > 0) + ASSERT (sze > 0) + allocate(H_jj_tmp(sze)) + + do i=1,sze + H_jj_tmp(i) = H_jj(i) + Dress_jj(i) + enddo + do k=1,N_st + do i=1,sze + H_jj_tmp(i) += u_in(i,k) * Dressing_vec(i,k) + enddo + enddo + + integer :: iter + integer :: i,j,k,l,m + + double precision, external :: u_dot_v, u_dot_u + + integer :: k_pairs, kl + + integer :: iter2, itertot + double precision, allocatable :: y(:,:), h(:,:), lambda(:) + double precision, allocatable :: s_tmp(:,:) + double precision, allocatable :: residual_norm(:),inv_c_idx_dress_vec(:) + character*(16384) :: write_buffer + double precision :: to_print(2,N_st),inv_c_idx_dress + double precision :: cpu, wall + integer :: shift, shift2, itermax, istate + double precision :: r1, r2, alpha + logical :: state_ok(N_st_diag*davidson_sze_max) + integer :: nproc_target + integer :: order(N_st_diag) + double precision :: cmax + double precision, allocatable :: U(:,:), overlap(:,:) + double precision, pointer :: W(:,:) + logical :: disk_based + double precision :: energy_shift(N_st_diag*davidson_sze_max) + + + allocate(inv_c_idx_dress_vec(N_st)) + inv_c_idx_dress = 1.d0/u_in(idx_dress,1) + do i = 1, N_st + inv_c_idx_dress_vec(i) = 1.d0/u_in(idx_dress,i) + enddo + include 'constants.include.F' + + integer :: N_st_diag_in + N_st_diag_in = N_st_diag + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda + if (N_st_diag_in*3 > sze) then + print *, 'error in Davidson :' + print *, 'Increase n_det_max_full to ', N_st_diag_in*3 + stop -1 + endif + + itermax = max(2,min(davidson_sze_max, sze/N_st_diag_in))+1 + itertot = 0 + + if (state_following) then + allocate(overlap(N_st_diag_in*itermax, N_st_diag_in*itermax)) + else + allocate(overlap(1,1)) ! avoid 'if' for deallocate + endif + overlap = 0.d0 + + call write_time(6) + write(6,'(A)') '' + write(6,'(A)') 'Davidson Diagonalization' + write(6,'(A)') '------------------------' + write(6,'(A)') '' + + ! Find max number of cores to fit in memory + ! ----------------------------------------- + + nproc_target = nproc + double precision :: rss + integer :: maxab + maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 + + m=1 + disk_based = .False. + call resident_memory(rss) + do + r1 = 8.d0 * &! bytes + ( dble(sze)*(N_st_diag_in*itermax) &! U + + 1.0d0*dble(sze*m)*(N_st_diag_in*itermax) &! W + + 3.0d0*(N_st_diag_in*itermax)**2 &! h,y,s_tmp + + 1.d0*(N_st_diag_in*itermax) &! lambda + + 1.d0*(N_st_diag_in) &! residual_norm + ! In H_u_0_nstates_zmq + + 2.d0*(N_st_diag_in*N_det) &! u_t, v_t, on collector + + 2.d0*(N_st_diag_in*N_det) &! u_t, v_t, on slave + + 0.5d0*maxab &! idx0 in H_u_0_nstates_openmp_work_* + + nproc_target * &! In OMP section + ( 1.d0*(N_int*maxab) &! buffer + + 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx + ) / 1024.d0**3 + + if (nproc_target == 0) then + call check_mem(r1,irp_here) + nproc_target = 1 + exit + endif + + if (r1+rss < qp_max_mem) then + exit + endif + + if (itermax > 4) then + itermax = itermax - 1 + else if (m==1.and.disk_based_davidson) then + m=0 + disk_based = .True. + itermax = 6 + else + nproc_target = nproc_target - 1 + endif + + enddo + nthreads_davidson = nproc_target + TOUCH nthreads_davidson + call write_int(6,N_st,'Number of states') + call write_int(6,N_st_diag_in,'Number of states in diagonalization') + call write_int(6,sze,'Number of basis functions ') + call write_int(6,nproc_target,'Number of threads for diagonalization') + call write_double(6, r1, 'Memory(Gb)') + if (disk_based) then + print *, 'Using swap space to reduce RAM' + endif + + !--------------- + + write(6,'(A)') '' + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = 'Iter' + do i=1,N_st + write_buffer = trim(write_buffer)//' Energy Residual ' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + + + allocate(W(sze,N_st_diag_in*itermax)) + + allocate( & + ! Large + U(sze,N_st_diag_in*itermax), & + + ! Small + h(N_st_diag_in*itermax,N_st_diag_in*itermax), & + y(N_st_diag_in*itermax,N_st_diag_in*itermax), & + s_tmp(N_st_diag_in*itermax,N_st_diag_in*itermax), & + residual_norm(N_st_diag_in), & + lambda(N_st_diag_in*itermax)) + + h = 0.d0 + U = 0.d0 + y = 0.d0 + s_tmp = 0.d0 + + + ASSERT (N_st > 0) + ASSERT (N_st_diag_in >= N_st) + ASSERT (sze > 0) + + ! Davidson iterations + ! =================== + + converged = .False. + + do k=N_st+1,N_st_diag_in + do i=1,sze + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st) + enddo + u_in(k,k) = u_in(k,k) + 10.d0 + enddo + do k=1,N_st_diag_in + call normalize(u_in(1,k),sze) + enddo + + do k=1,N_st_diag_in + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + + do while (.not.converged) + itertot = itertot+1 + if (itertot == 2) then + exit + endif + + do iter=1,itermax-1 + + shift = N_st_diag_in*(iter-1) + shift2 = N_st_diag_in*iter + + if ((iter > 1).or.(itertot == 1)) then + ! Compute |W_k> = \sum_i |i> + ! ----------------------------------- + call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag_in,sze) + ! Compute then the DIAGONAL PART OF THE DRESSING + ! += Dress_jj(i) * + call dressing_diag_uv(W(1,shift+1),U(1,shift+1),Dress_jj,N_st_diag_in,sze) + else + ! Already computed in update below + continue + endif + + + if (N_st == 1) then + + l = idx_dress + double precision :: f + f = inv_c_idx_dress + do istate=1,N_st_diag_in + do i=1,sze + W(i,shift+istate) += Dressing_vec(i,1) *f * U(l,shift+istate) + W(l,shift+istate) += Dressing_vec(i,1) *f * U(i,shift+istate) + enddo + enddo + + else + print*,'dav_double_dressed routine not yet implemented for N_st > 1' +! +! call dgemm('T','N', N_st, N_st_diag_in, sze, 1.d0, & +! psi_coef, size(psi_coef,1), & +! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) +! +! call dgemm('N','N', sze, N_st_diag_in, N_st, 1.0d0, & +! Dressing_vec, size(Dressing_vec,1), s_tmp, size(s_tmp,1), & +! 1.d0, W(1,shift+1), size(W,1)) +! +! +! call dgemm('T','N', N_st, N_st_diag_in, sze, 1.d0, & +! Dressing_vec, size(Dressing_vec,1), & +! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) +! +! call dgemm('N','N', sze, N_st_diag_in, N_st, 1.0d0, & +! psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), & +! 1.d0, W(1,shift+1), size(W,1)) +! + endif + + ! Compute h_kl = = + ! ------------------------------------------- + + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), W, size(W,1), & + 0.d0, h, size(h,1)) + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), U, size(U,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + ! Diagonalize h + ! --------------- + + integer :: lwork, info + double precision, allocatable :: work(:) + + y = h + lwork = -1 + allocate(work(1)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + lwork = int(work(1)) + deallocate(work) + allocate(work(lwork)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + deallocate(work) + if (info /= 0) then + stop 'DSYGV Diagonalization failed' + endif + + ! Compute Energy for each eigenvector + ! ----------------------------------- + + call dgemm('N','N',shift2,shift2,shift2, & + 1.d0, h, size(h,1), y, size(y,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + call dgemm('T','N',shift2,shift2,shift2, & + 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & + 0.d0, h, size(h,1)) + + do k=1,shift2 + lambda(k) = h(k,k) + enddo + + if (state_following) then + + overlap = -1.d0 + do k=1,shift2 + do i=1,shift2 + overlap(k,i) = dabs(y(k,i)) + enddo + enddo + do k=1,N_st + cmax = -1.d0 + do i=1,N_st + if (overlap(i,k) > cmax) then + cmax = overlap(i,k) + order(k) = i + endif + enddo + do i=1,N_st_diag_in + overlap(order(k),i) = -1.d0 + enddo + enddo + overlap = y + do k=1,N_st + l = order(k) + if (k /= l) then + y(1:shift2,k) = overlap(1:shift2,l) + endif + enddo + do k=1,N_st + overlap(k,1) = lambda(k) + enddo + + endif + + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + + call dgemm('N','N', sze, N_st_diag_in, shift2, & + 1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1)) + call dgemm('N','N', sze, N_st_diag_in, shift2, & + 1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1)) + + ! Compute residual vector and davidson step + ! ----------------------------------------- + + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) + do k=1,N_st_diag_in + do i=1,sze + U(i,shift2+k) = & + (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & + /max(H_jj_tmp(i) - lambda (k),1.d-2) + enddo + + if (k <= N_st) then + residual_norm(k) = u_dot_u(U(1,shift2+k),sze) + to_print(1,k) = lambda(k) + to_print(2,k) = residual_norm(k) + endif + enddo + !$OMP END PARALLEL DO + + + if ((itertot>1).and.(iter == 1)) then + !don't print + continue + else + write(*,'(1X,I3,1X,100(1X,F16.10,1X,E11.3))') iter-1, to_print(1:2,1:N_st) + endif + + ! Check convergence + if (iter > 1) then + converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson + endif + + do k=1,N_st + if (residual_norm(k) > 1.d8) then + print *, 'Davidson failed' + stop -1 + endif + enddo + if (converged) then + exit + endif + + logical, external :: qp_stop + if (qp_stop()) then + converged = .True. + exit + endif + + + enddo + + ! Re-contract U and update W + ! -------------------------------- + + call dgemm('N','N', sze, N_st_diag_in, shift2, 1.d0, & + W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag_in + do i=1,sze + W(i,k) = u_in(i,k) + enddo + enddo + + call dgemm('N','N', sze, N_st_diag_in, shift2, 1.d0, & + U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + + do k=1,N_st_diag_in + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + enddo + + + call nullify_small_elements(sze,N_st_diag_in,U,size(U,1),threshold_davidson_pt2) + do k=1,N_st_diag_in + do i=1,sze + u_in(i,k) = U(i,k) + enddo + enddo + + do k=1,N_st_diag_in + energies(k) = lambda(k) + enddo + write_buffer = '======' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') trim(write_buffer) + write(6,'(A)') '' + call write_time(6) + + deallocate(W) + + deallocate ( & + residual_norm, & + U, overlap, & + h, y, s_tmp, & + lambda & + ) + FREE nthreads_davidson +end + + +subroutine dressing_diag_uv(v,u,dress_diag,N_st,sze) + implicit none + BEGIN_DOC + ! Routine that computes the diagonal part of the dressing + ! + ! v(i) += u(i) * dress_diag(i) + ! + ! !!!!!!!! WARNING !!!!!!!! the vector v is not initialized + ! + ! !!!!!!!! SO MAKE SURE THERE ARE SOME MEANINGFUL VALUES IN THERE + END_DOC + integer, intent(in) :: N_st,sze + double precision, intent(in) :: u(sze,N_st),dress_diag(sze) + double precision, intent(inout) :: v(sze,N_st) + integer :: i,istate + do istate = 1, N_st + do i = 1, sze + v(i,istate) += dress_diag(i) * u(i,istate) + enddo + enddo +end + + + + + + + + + + + diff --git a/src/dav_general_mat/dav_dressed_ext_rout.irp.f b/src/dav_general_mat/dav_dressed_ext_rout.irp.f new file mode 100644 index 00000000..c3bfe91a --- /dev/null +++ b/src/dav_general_mat/dav_dressed_ext_rout.irp.f @@ -0,0 +1,485 @@ +subroutine davidson_general_ext_rout_dressed(u_in,H_jj,energies,sze,N_st,N_st_diag,dressing_state,dressing_vec,idress,converged,hcalc) + use mmap_module + implicit none + BEGIN_DOC + ! Davidson diagonalization. + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! sze : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: sze, N_st, N_st_diag,idress + double precision, intent(inout) :: u_in(sze,N_st_diag) + double precision, intent(inout) :: H_jj(sze) + double precision, intent(out) :: energies(N_st_diag) + double precision, intent(in) :: dressing_vec(sze,N_st) + integer, intent(in) :: dressing_state + logical, intent(out) :: converged + external hcalc + + double precision :: f + + integer :: iter + integer :: i,j,k,l,m + + double precision, external :: u_dot_v, u_dot_u + + integer :: k_pairs, kl + + integer :: iter2, itertot + double precision, allocatable :: y(:,:), h(:,:), lambda(:) + double precision, allocatable :: s_tmp(:,:) + double precision :: diag_h_mat_elem + double precision, allocatable :: residual_norm(:) + character*(16384) :: write_buffer + double precision :: to_print(2,N_st) + double precision :: cpu, wall + integer :: shift, shift2, itermax, istate + double precision :: r1, r2, alpha + logical :: state_ok(N_st_diag*davidson_sze_max) + integer :: nproc_target + integer :: order(N_st_diag) + double precision :: cmax + double precision, allocatable :: U(:,:), overlap(:,:) + double precision, pointer :: W(:,:) + logical :: disk_based + double precision :: energy_shift(N_st_diag*davidson_sze_max) + !!!! TO CHANGE !!!! + integer :: idx_dress(1) + idx_dress = idress + + + if (dressing_state > 0) then + do k=1,N_st + do i=1,sze + H_jj(i) += u_in(i,k) * dressing_vec(i,k) + enddo + enddo + endif + + l = idx_dress(1) + f = 1.0d0/u_in(l,1) + + include 'constants.include.F' + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, y, h, lambda + if (N_st_diag*3 > sze) then + print *, 'error in Davidson :' + print *, 'Increase n_det_max_full to ', N_st_diag*3 + stop -1 + endif + + itermax = max(2,min(davidson_sze_max, sze/N_st_diag))+1 + itertot = 0 + + if (state_following) then + allocate(overlap(N_st_diag*itermax, N_st_diag*itermax)) + else + allocate(overlap(1,1)) ! avoid 'if' for deallocate + endif + overlap = 0.d0 + + + provide threshold_davidson !nthreads_davidson + call write_time(6) + write(6,'(A)') '' + write(6,'(A)') 'Davidson Diagonalization' + write(6,'(A)') '------------------------' + write(6,'(A)') '' + + ! Find max number of cores to fit in memory + ! ----------------------------------------- + + nproc_target = nproc + double precision :: rss + integer :: maxab +! maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 + maxab = sze + + m=1 + disk_based = .False. + call resident_memory(rss) + do + r1 = 8.d0 * &! bytes + ( dble(sze)*(N_st_diag*itermax) &! U + + 1.0d0*dble(sze*m)*(N_st_diag*itermax) &! W + + 3.0d0*(N_st_diag*itermax)**2 &! h,y,s_tmp + + 1.d0*(N_st_diag*itermax) &! lambda + + 1.d0*(N_st_diag) &! residual_norm + ! In H_u_0_nstates_zmq + + 2.d0*(N_st_diag*N_det) &! u_t, v_t, on collector + + 2.d0*(N_st_diag*N_det) &! u_t, v_t, on slave + + 0.5d0*maxab &! idx0 in H_u_0_nstates_openmp_work_* + + nproc_target * &! In OMP section + ( 1.d0*(N_int*maxab) &! buffer + + 3.5d0*(maxab) ) &! singles_a, singles_b, doubles, idx + ) / 1024.d0**3 + + if (nproc_target == 0) then + call check_mem(r1,irp_here) + nproc_target = 1 + exit + endif + + if (r1+rss < qp_max_mem) then + exit + endif + + if (itermax > 4) then + itermax = itermax - 1 + else if (m==1.and.disk_based_davidson) then + m=0 + disk_based = .True. + itermax = 6 + else + nproc_target = nproc_target - 1 + endif + + enddo + nthreads_davidson = nproc_target + TOUCH nthreads_davidson + call write_int(6,N_st,'Number of states') + call write_int(6,N_st_diag,'Number of states in diagonalization') + call write_int(6,sze,'Number of basis function') + call write_int(6,nproc_target,'Number of threads for diagonalization') + call write_double(6, r1, 'Memory(Gb)') + if (disk_based) then + print *, 'Using swap space to reduce RAM' + endif + + !--------------- + + write(6,'(A)') '' + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = 'Iter' + do i=1,N_st + write_buffer = trim(write_buffer)//' Energy Residual ' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + write_buffer = '=====' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*N_st) + + + allocate(W(sze,N_st_diag*itermax)) + + allocate( & + ! Large + U(sze,N_st_diag*itermax), & + + ! Small + h(N_st_diag*itermax,N_st_diag*itermax), & + y(N_st_diag*itermax,N_st_diag*itermax), & + s_tmp(N_st_diag*itermax,N_st_diag*itermax), & + residual_norm(N_st_diag), & + lambda(N_st_diag*itermax)) + + h = 0.d0 + U = 0.d0 + y = 0.d0 + s_tmp = 0.d0 + + + ASSERT (N_st > 0) + ASSERT (N_st_diag >= N_st) + ASSERT (sze > 0) + + ! Davidson iterations + ! =================== + + converged = .False. + + do k=N_st+1,N_st_diag + do i=1,sze + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + u_in(i,k) = r1*dcos(r2) * u_in(i,k-N_st) + enddo + u_in(k,k) = u_in(k,k) + 10.d0 + enddo + ! Normalize all states + do k=1,N_st_diag + call normalize(u_in(1,k),sze) + enddo + ! Copy from the guess input "u_in" to the working vectors "U" + + do k=1,N_st_diag + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + + do while (.not.converged) + itertot = itertot+1 + if (itertot == 8) then + exit + endif + + do iter=1,itermax-1 + + shift = N_st_diag*(iter-1) + shift2 = N_st_diag*iter + + if ((iter > 1).or.(itertot == 1)) then + ! Compute |W_k> = \sum_i |i> + ! ----------------------------------- + ! Gram-Schmidt to orthogonalize all new guess with the previous vectors + call ortho_qr(U,size(U,1),sze,shift2) + call ortho_qr(U,size(U,1),sze,shift2) + ! it does W = H U with W(sze,N_st_diag),U(sze,N_st_diag) + ! where sze is the size of the vector, N_st_diag is the number of states + call hcalc(W(1,shift+1),U(1,shift+1),N_st_diag,sze) + else + ! Already computed in update below + continue + endif + + if (dressing_state > 0) then + + if (N_st == 1) then + + + do istate=1,N_st_diag + do i=1,sze + W(i,shift+istate) += dressing_vec(i,1) *f * U(l,shift+istate) + W(l,shift+istate) += dressing_vec(i,1) *f * U(i,shift+istate) + enddo + + enddo + + else + print*,'Not implemented yet for multi state ...' + stop +! call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, & +! psi_coef, size(psi_coef,1), & +! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) +! +! call dgemm('N','N', sze, N_st_diag, N_st, 1.0d0, & +! dressing_vec, size(dressing_vec,1), s_tmp, size(s_tmp,1), & +! 1.d0, W(1,shift+1), size(W,1)) +! +! +! call dgemm('T','N', N_st, N_st_diag, sze, 1.d0, & +! dressing_vec, size(dressing_vec,1), & +! U(1,shift+1), size(U,1), 0.d0, s_tmp, size(s_tmp,1)) +! +! call dgemm('N','N', sze, N_st_diag, N_st, 1.0d0, & +! psi_coef, size(psi_coef,1), s_tmp, size(s_tmp,1), & +! 1.d0, W(1,shift+1), size(W,1)) + + endif + endif + + ! Compute h_kl = = + ! ------------------------------------------- + + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), W, size(W,1), & + 0.d0, h, size(h,1)) + call dgemm('T','N', shift2, shift2, sze, & + 1.d0, U, size(U,1), U, size(U,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + ! Diagonalize h + ! --------------- + + integer :: lwork, info + double precision, allocatable :: work(:) + + y = h + lwork = -1 + allocate(work(1)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + lwork = int(work(1)) + deallocate(work) + allocate(work(lwork)) + call dsygv(1,'V','U',shift2,y,size(y,1), & + s_tmp,size(s_tmp,1), lambda, work,lwork,info) + deallocate(work) + if (info /= 0) then + stop 'DSYGV Diagonalization failed' + endif + + ! Compute Energy for each eigenvector + ! ----------------------------------- + + call dgemm('N','N',shift2,shift2,shift2, & + 1.d0, h, size(h,1), y, size(y,1), & + 0.d0, s_tmp, size(s_tmp,1)) + + call dgemm('T','N',shift2,shift2,shift2, & + 1.d0, y, size(y,1), s_tmp, size(s_tmp,1), & + 0.d0, h, size(h,1)) + + do k=1,shift2 + lambda(k) = h(k,k) + enddo + + if (state_following) then + + overlap = -1.d0 + do k=1,shift2 + do i=1,shift2 + overlap(k,i) = dabs(y(k,i)) + enddo + enddo + do k=1,N_st + cmax = -1.d0 + do i=1,N_st + if (overlap(i,k) > cmax) then + cmax = overlap(i,k) + order(k) = i + endif + enddo + do i=1,N_st_diag + overlap(order(k),i) = -1.d0 + enddo + enddo + overlap = y + do k=1,N_st + l = order(k) + if (k /= l) then + y(1:shift2,k) = overlap(1:shift2,l) + endif + enddo + do k=1,N_st + overlap(k,1) = lambda(k) + enddo + + endif + + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, U, size(U,1), y, size(y,1), 0.d0, U(1,shift2+1), size(U,1)) + call dgemm('N','N', sze, N_st_diag, shift2, & + 1.d0, W, size(W,1), y, size(y,1), 0.d0, W(1,shift2+1), size(W,1)) + + ! Compute residual vector and davidson step + ! ----------------------------------------- + + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) + do k=1,N_st_diag + do i=1,sze + U(i,shift2+k) = & + (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) & + /max(H_jj(i) - lambda (k),1.d-2) + enddo + + if (k <= N_st) then + residual_norm(k) = u_dot_u(U(1,shift2+k),sze) + to_print(1,k) = lambda(k) + to_print(2,k) = residual_norm(k) + endif + enddo + !$OMP END PARALLEL DO + + + if ((itertot>1).and.(iter == 1)) then + !don't print + continue + else + write(*,'(1X,I3,1X,100(1X,F16.10,1X,E11.3))') iter-1, to_print(1:2,1:N_st) + endif + + ! Check convergence + if (iter > 1) then + converged = dabs(maxval(residual_norm(1:N_st))) < threshold_davidson + endif + + do k=1,N_st + if (residual_norm(k) > 1.d8) then + print *, 'Davidson failed' + stop -1 + endif + enddo + if (converged) then + exit + endif + + logical, external :: qp_stop + if (qp_stop()) then + converged = .True. + exit + endif + + + enddo + + ! Re-contract U and update W + ! -------------------------------- + + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & + W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + do k=1,N_st_diag + do i=1,sze + W(i,k) = u_in(i,k) + enddo + enddo + + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & + U, size(U,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) + + do k=1,N_st_diag + do i=1,sze + U(i,k) = u_in(i,k) + enddo + enddo + + enddo + + + call nullify_small_elements(sze,N_st_diag,U,size(U,1),threshold_davidson_pt2) + do k=1,N_st_diag + do i=1,sze + u_in(i,k) = U(i,k) + enddo + enddo + + do k=1,N_st_diag + energies(k) = lambda(k) + enddo + write_buffer = '======' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') trim(write_buffer) + write(6,'(A)') '' + call write_time(6) + + deallocate(W) + + deallocate ( & + residual_norm, & + U, overlap, & + h, y, s_tmp, & + lambda & + ) + FREE nthreads_davidson +end + + + + + + + diff --git a/src/dav_general_mat/dav_ext_rout.irp.f b/src/dav_general_mat/dav_ext_rout.irp.f index a2aad413..868d928b 100644 --- a/src/dav_general_mat/dav_ext_rout.irp.f +++ b/src/dav_general_mat/dav_ext_rout.irp.f @@ -1,15 +1,15 @@ -subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,hcalc) +subroutine davidson_general_ext_rout(u_in,H_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc) use mmap_module implicit none BEGIN_DOC - ! Davidson diagonalization with specific diagonal elements of the H matrix + ! Generic Davidson diagonalization ! ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson ! ! u_in : guess coefficients on the various states. Overwritten on exit ! - ! dim_in : leftmost dimension of u_in + ! sze : leftmost dimension of u_in ! ! sze : Number of determinants ! @@ -21,9 +21,9 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia ! ! hcalc subroutine to compute W = H U (see routine hcalc_template for template of input/output) END_DOC - integer, intent(in) :: dim_in, sze, N_st, N_st_diag_in + integer, intent(in) :: sze, N_st, N_st_diag_in double precision, intent(in) :: H_jj(sze) - double precision, intent(inout) :: u_in(dim_in,N_st_diag_in) + double precision, intent(inout) :: u_in(sze,N_st_diag_in) double precision, intent(out) :: energies(N_st) external hcalc @@ -157,19 +157,7 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia write(6,'(A)') write_buffer(1:6+41*N_st) -! if (disk_based) then -! ! Create memory-mapped files for W and S -! type(c_ptr) :: ptr_w, ptr_s -! integer :: fd_s, fd_w -! call mmap(trim(ezfio_work_dir)//'davidson_w', (/int(sze,8),int(N_st_diag*itermax,8)/),& -! 8, fd_w, .False., ptr_w) -! call mmap(trim(ezfio_work_dir)//'davidson_s', (/int(sze,8),int(N_st_diag*itermax,8)/),& -! 4, fd_s, .False., ptr_s) -! call c_f_pointer(ptr_w, w, (/sze,N_st_diag*itermax/)) -! call c_f_pointer(ptr_s, s, (/sze,N_st_diag*itermax/)) -! else - allocate(W(sze,N_st_diag*itermax)) -! endif + allocate(W(sze,N_st_diag*itermax)) allocate( & ! Large @@ -233,7 +221,6 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia if ((iter > 1).or.(itertot == 1)) then ! Compute |W_k> = \sum_i |i> ! ----------------------------------- - ! Gram-Schmidt to orthogonalize all new guess with the previous vectors call ortho_qr(U,size(U,1),sze,shift2) call ortho_qr(U,size(U,1),sze,shift2) @@ -357,6 +344,9 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia enddo + ! Re-contract U and update W + ! -------------------------------- + call dgemm('N','N', sze, N_st_diag, shift2, 1.d0, & W, size(W,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) do k=1,N_st_diag @@ -372,8 +362,8 @@ subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_dia U(i,k) = u_in(i,k) enddo enddo - call ortho_qr(U,size(U,1),sze,N_st_diag) - call ortho_qr(U,size(U,1),sze,N_st_diag) + call ortho_qr(U,size(U,1),sze,N_st_diag) + call ortho_qr(U,size(U,1),sze,N_st_diag) do j=1,N_st_diag k=1 do while ((k 1).or.(itertot == 1)) then +! if ((iter > 1).or.(itertot == 1)) then ! Compute |W_k> = \sum_i |i> ! ----------------------------------- @@ -309,10 +309,10 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N else call H_u_0_nstates_openmp(W,U,N_st_diag,sze) endif - else - ! Already computed in update below - continue - endif +! else +! ! Already computed in update below +! continue +! endif if (dressing_state > 0) then @@ -508,17 +508,8 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N enddo - ! Re-contract U and update W - ! -------------------------------- - - call dgemm('N','N', sze_csf, N_st_diag, shift2, 1.d0, & - W_csf, size(W_csf,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) - do k=1,N_st_diag - do i=1,sze_csf - W_csf(i,k) = u_in(i,k) - enddo - enddo - call convertWFfromCSFtoDET(N_st_diag,W_csf,W) + ! Re-contract U + ! ------------- call dgemm('N','N', sze_csf, N_st_diag, shift2, 1.d0, & U_csf, size(U_csf,1), y, size(y,1), 0.d0, u_in, size(u_in,1)) diff --git a/src/davidson/diagonalization_hs2_dressed.irp.f b/src/davidson/diagonalization_hs2_dressed.irp.f index 1a27a75e..d37b7386 100644 --- a/src/davidson/diagonalization_hs2_dressed.irp.f +++ b/src/davidson/diagonalization_hs2_dressed.irp.f @@ -349,7 +349,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ shift = N_st_diag*(iter-1) shift2 = N_st_diag*iter - if ((iter > 1).or.(itertot == 1)) then +! if ((iter > 1).or.(itertot == 1)) then ! Compute |W_k> = \sum_i |i> ! ----------------------------------- @@ -359,10 +359,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_ call H_S2_u_0_nstates_openmp(W(1,shift+1),S_d,U(1,shift+1),N_st_diag,sze) endif S(1:sze,shift+1:shift+N_st_diag) = real(S_d(1:sze,1:N_st_diag)) - else - ! Already computed in update below - continue - endif +! else +! ! Already computed in update below +! continue +! endif if (dressing_state > 0) then diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f index fb991b65..46ad8f78 100644 --- a/src/davidson/diagonalize_ci.irp.f +++ b/src/davidson/diagonalize_ci.irp.f @@ -1,9 +1,9 @@ - BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ] implicit none BEGIN_DOC ! :c:data:`n_states` lowest eigenvalues of the |CI| matrix END_DOC + PROVIDE distributed_davidson integer :: j character*(8) :: st @@ -247,6 +247,7 @@ subroutine diagonalize_CI ! eigenstates of the |CI| matrix. END_DOC integer :: i,j + PROVIDE distributed_davidson do j=1,N_states do i=1,N_det psi_coef(i,j) = CI_eigenvectors(i,j) diff --git a/src/davidson/u0_hs2_u0.irp.f b/src/davidson/u0_hs2_u0.irp.f index 8f7bf06b..38fb56bd 100644 --- a/src/davidson/u0_hs2_u0.irp.f +++ b/src/davidson/u0_hs2_u0.irp.f @@ -203,7 +203,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend, integer, allocatable :: doubles(:) integer, allocatable :: singles_a(:) integer, allocatable :: singles_b(:) - integer, allocatable :: idx(:), idx0(:) + integer, allocatable :: idx(:), buffer_lrow(:), idx0(:) integer :: maxab, n_singles_a, n_singles_b, kcol_prev integer*8 :: k8 logical :: compute_singles @@ -253,7 +253,7 @@ compute_singles=.True. !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, & !$OMP lcol, lrow, l_a, l_b, utl, kk, u_is_sparse, & !$OMP buffer, doubles, n_doubles, umax, & - !$OMP tmp_det2, hij, sij, idx, l, kcol_prev, & + !$OMP tmp_det2, hij, sij, idx, buffer_lrow, l, kcol_prev, & !$OMP singles_a, n_singles_a, singles_b, ratio, & !$OMP n_singles_b, k8, last_found,left,right,right_max) @@ -264,7 +264,7 @@ compute_singles=.True. singles_a(maxab), & singles_b(maxab), & doubles(maxab), & - idx(maxab), utl(N_st,block_size)) + idx(maxab), buffer_lrow(maxab), utl(N_st,block_size)) kcol_prev=-1 @@ -332,18 +332,20 @@ compute_singles=.True. l_a = psi_bilinear_matrix_columns_loc(lcol) ASSERT (l_a <= N_det) - !DIR$ UNROLL(8) - !DIR$ LOOP COUNT avg(50000) do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) lrow = psi_bilinear_matrix_rows(l_a) ASSERT (lrow <= N_det_alpha_unique) - buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow) ! hot spot + buffer_lrow(j) = lrow ASSERT (l_a <= N_det) idx(j) = l_a l_a = l_a+1 enddo + + do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol) + buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, buffer_lrow(j)) ! hot spot + enddo j = j-1 call get_all_spin_singles_$N_int( & @@ -789,7 +791,7 @@ compute_singles=.True. end do !$OMP END DO - deallocate(buffer, singles_a, singles_b, doubles, idx, utl) + deallocate(buffer, singles_a, singles_b, doubles, idx, buffer_lrow, utl) !$OMP END PARALLEL end diff --git a/src/davidson_dressed/diagonalize_ci.irp.f b/src/davidson_dressed/diagonalize_ci.irp.f index 7619532a..b58ce9c0 100644 --- a/src/davidson_dressed/diagonalize_ci.irp.f +++ b/src/davidson_dressed/diagonalize_ci.irp.f @@ -12,7 +12,7 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ] enddo do j=1,min(N_det,N_states) write(st,'(I4)') j - call write_double(6,CI_energy_dressed(j),'Energy of state '//trim(st)) + call write_double(6,CI_energy_dressed(j),'Energy dressed of state '//trim(st)) call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st)) enddo @@ -21,133 +21,201 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ] &BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ] - BEGIN_DOC - ! Eigenvectors/values of the CI matrix - END_DOC - implicit none - double precision :: ovrlp,u_dot_v - integer :: i_good_state - integer, allocatable :: index_good_state_array(:) - logical, allocatable :: good_state_array(:) - double precision, allocatable :: s2_values_tmp(:) - integer :: i_other_state - double precision, allocatable :: eigenvectors(:,:), eigenvectors_s2(:,:), eigenvalues(:) - integer :: i_state - double precision :: e_0 - integer :: i,j,k,mrcc_state - double precision, allocatable :: s2_eigvalues(:) - double precision, allocatable :: e_array(:) - integer, allocatable :: iorder(:) - - PROVIDE threshold_davidson nthreads_davidson - ! Guess values for the "N_states" states of the CI_eigenvectors_dressed - do j=1,min(N_states,N_det) - do i=1,N_det - CI_eigenvectors_dressed(i,j) = psi_coef(i,j) - enddo - enddo - - do j=min(N_states,N_det)+1,N_states_diag - do i=1,N_det - CI_eigenvectors_dressed(i,j) = 0.d0 - enddo - enddo - - if (diag_algorithm == "Davidson") then - - do j=1,min(N_states,N_det) - do i=1,N_det - CI_eigenvectors_dressed(i,j) = psi_coef(i,j) - enddo - enddo - logical :: converged - converged = .False. - call davidson_diag_HS2(psi_det,CI_eigenvectors_dressed, CI_eigenvectors_s2_dressed,& - size(CI_eigenvectors_dressed,1), CI_electronic_energy_dressed,& - N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1,converged) - - else if (diag_algorithm == "Lapack") then - - allocate (eigenvectors(size(H_matrix_dressed,1),N_det)) - allocate (eigenvalues(N_det)) - - call lapack_diag(eigenvalues,eigenvectors, & - H_matrix_dressed,size(H_matrix_dressed,1),N_det) - CI_electronic_energy_dressed(:) = 0.d0 - if (s2_eig) then - i_state = 0 - allocate (s2_eigvalues(N_det)) - allocate(index_good_state_array(N_det),good_state_array(N_det)) - good_state_array = .False. - - call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,& - N_det,size(eigenvectors,1)) - do j=1,N_det - ! Select at least n_states states with S^2 values closed to "expected_s2" - if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then - i_state +=1 - index_good_state_array(i_state) = j - good_state_array(j) = .True. - endif - if(i_state.eq.N_states) then - exit - endif + BEGIN_DOC + ! Eigenvectors/values of the CI matrix + END_DOC + implicit none + double precision :: ovrlp,u_dot_v + integer :: i_good_state + integer, allocatable :: index_good_state_array(:) + logical, allocatable :: good_state_array(:) + double precision, allocatable :: s2_values_tmp(:) + integer :: i_other_state + double precision, allocatable :: eigenvectors(:,:), eigenvectors_s2(:,:), eigenvalues(:) + integer :: i_state + double precision :: e_0 + integer :: i,j,k,mrcc_state + double precision, allocatable :: s2_eigvalues(:) + double precision, allocatable :: e_array(:) + integer, allocatable :: iorder(:) + logical :: converged + logical :: do_csf + + PROVIDE threshold_davidson nthreads_davidson + ! Guess values for the "N_states" states of the CI_eigenvectors_dressed + do j=1,min(N_states,N_det) + do i=1,N_det + CI_eigenvectors_dressed(i,j) = psi_coef(i,j) + enddo + enddo + + do j=min(N_states,N_det)+1,N_states_diag + do i=1,N_det + CI_eigenvectors_dressed(i,j) = 0.d0 + enddo + enddo + + do_csf = s2_eig .and. only_expected_s2 .and. csf_based + + if (diag_algorithm == "Davidson") then + + do j=1,min(N_states,N_det) + do i=1,N_det + CI_eigenvectors_dressed(i,j) = psi_coef(i,j) enddo - if(i_state .ne.0)then - ! Fill the first "i_state" states that have a correct S^2 value - do j = 1, i_state - do i=1,N_det - CI_eigenvectors_dressed(i,j) = eigenvectors(i,index_good_state_array(j)) - enddo - CI_electronic_energy_dressed(j) = eigenvalues(index_good_state_array(j)) - CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(j)) - enddo - i_other_state = 0 - do j = 1, N_det - if(good_state_array(j))cycle - i_other_state +=1 - if(i_state+i_other_state.gt.n_states_diag)then - exit - endif - do i=1,N_det - CI_eigenvectors_dressed(i,i_state+i_other_state) = eigenvectors(i,j) - enddo - CI_electronic_energy_dressed(i_state+i_other_state) = eigenvalues(j) - CI_eigenvectors_s2_dressed(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state) - enddo - else - print*,'' - print*,'!!!!!!!! WARNING !!!!!!!!!' - print*,' Within the ',N_det,'determinants selected' - print*,' and the ',N_states_diag,'states requested' - print*,' We did not find any state with S^2 values close to ',expected_s2 - print*,' We will then set the first N_states eigenvectors of the H matrix' - print*,' as the CI_eigenvectors_dressed' - print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space' - print*,'' - do j=1,min(N_states_diag,N_det) - do i=1,N_det - CI_eigenvectors_dressed(i,j) = eigenvectors(i,j) - enddo - CI_electronic_energy_dressed(j) = eigenvalues(j) - CI_eigenvectors_s2_dressed(j) = s2_eigvalues(j) - enddo - endif - deallocate(index_good_state_array,good_state_array) - deallocate(s2_eigvalues) + enddo + converged = .False. + if (do_csf) then + call davidson_diag_H_csf(psi_det,CI_eigenvectors_dressed, & + size(CI_eigenvectors_dressed,1),CI_electronic_energy_dressed, & + N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,1,converged) else - call u_0_S2_u_0(CI_eigenvectors_s2_dressed,eigenvectors,N_det,psi_det,N_int,& - min(N_det,N_states_diag),size(eigenvectors,1)) - ! Select the "N_states_diag" states of lowest energy - do j=1,min(N_det,N_states_diag) - do i=1,N_det - CI_eigenvectors_dressed(i,j) = eigenvectors(i,j) - enddo - CI_electronic_energy_dressed(j) = eigenvalues(j) - enddo + call davidson_diag_HS2(psi_det,CI_eigenvectors_dressed, CI_eigenvectors_s2_dressed,& + size(CI_eigenvectors_dressed,1), CI_electronic_energy_dressed,& + N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1,converged) endif - deallocate(eigenvectors,eigenvalues) - endif + + integer :: N_states_diag_save + N_states_diag_save = N_states_diag + do while (.not.converged) + double precision, allocatable :: CI_electronic_energy_tmp (:) + double precision, allocatable :: CI_eigenvectors_tmp (:,:) + double precision, allocatable :: CI_s2_tmp (:) + + N_states_diag *= 2 + TOUCH N_states_diag + + if (do_csf) then + + allocate (CI_electronic_energy_tmp (N_states_diag) ) + allocate (CI_eigenvectors_tmp (N_det,N_states_diag) ) + + CI_electronic_energy_tmp(1:N_states_diag_save) = CI_electronic_energy_dressed(1:N_states_diag_save) + CI_eigenvectors_tmp(1:N_det,1:N_states_diag_save) = CI_eigenvectors_dressed(1:N_det,1:N_states_diag_save) + + call davidson_diag_H_csf(psi_det,CI_eigenvectors_tmp, & + size(CI_eigenvectors_tmp,1),CI_electronic_energy_tmp, & + N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,1,converged) + + CI_electronic_energy_dressed(1:N_states_diag_save) = CI_electronic_energy_tmp(1:N_states_diag_save) + CI_eigenvectors_dressed(1:N_det,1:N_states_diag_save) = CI_eigenvectors_tmp(1:N_det,1:N_states_diag_save) + + deallocate (CI_electronic_energy_tmp) + deallocate (CI_eigenvectors_tmp) + + else + + allocate (CI_electronic_energy_tmp (N_states_diag) ) + allocate (CI_eigenvectors_tmp (N_det,N_states_diag) ) + allocate (CI_s2_tmp (N_states_diag) ) + + CI_electronic_energy_tmp(1:N_states_diag_save) = CI_electronic_energy_dressed(1:N_states_diag_save) + CI_eigenvectors_tmp(1:N_det,1:N_states_diag_save) = CI_eigenvectors_dressed(1:N_det,1:N_states_diag_save) + CI_s2_tmp(1:N_states_diag_save) = CI_eigenvectors_s2_dressed(1:N_states_diag_save) + + call davidson_diag_HS2(psi_det,CI_eigenvectors_tmp, CI_s2_tmp, & + size(CI_eigenvectors_tmp,1),CI_electronic_energy_tmp, & + N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1,converged) + + CI_electronic_energy_dressed(1:N_states_diag_save) = CI_electronic_energy_tmp(1:N_states_diag_save) + CI_eigenvectors_dressed(1:N_det,1:N_states_diag_save) = CI_eigenvectors_tmp(1:N_det,1:N_states_diag_save) + CI_eigenvectors_s2_dressed(1:N_states_diag_save) = CI_s2_tmp(1:N_states_diag_save) + + deallocate (CI_electronic_energy_tmp) + deallocate (CI_eigenvectors_tmp) + deallocate (CI_s2_tmp) + + endif + + enddo + if (N_states_diag > N_states_diag_save) then + N_states_diag = N_states_diag_save + TOUCH N_states_diag + endif + + else if (diag_algorithm == "Lapack") then + + print *, 'Diagonalization of H using Lapack' + allocate (eigenvectors(size(H_matrix_dressed,1),N_det)) + allocate (eigenvalues(N_det)) + + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_dressed,size(H_matrix_dressed,1),N_det) + CI_electronic_energy_dressed(:) = 0.d0 + if (s2_eig) then + i_state = 0 + allocate (s2_eigvalues(N_det)) + allocate(index_good_state_array(N_det),good_state_array(N_det)) + good_state_array = .False. + + call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,& + N_det,size(eigenvectors,1)) + do j=1,N_det + ! Select at least n_states states with S^2 values closed to "expected_s2" + if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then + i_state +=1 + index_good_state_array(i_state) = j + good_state_array(j) = .True. + endif + if(i_state.eq.N_states) then + exit + endif + enddo + if(i_state .ne.0)then + ! Fill the first "i_state" states that have a correct S^2 value + do j = 1, i_state + do i=1,N_det + CI_eigenvectors_dressed(i,j) = eigenvectors(i,index_good_state_array(j)) + enddo + CI_electronic_energy_dressed(j) = eigenvalues(index_good_state_array(j)) + CI_eigenvectors_s2_dressed(j) = s2_eigvalues(index_good_state_array(j)) + enddo + i_other_state = 0 + do j = 1, N_det + if(good_state_array(j))cycle + i_other_state +=1 + if(i_state+i_other_state.gt.n_states_diag)then + exit + endif + do i=1,N_det + CI_eigenvectors_dressed(i,i_state+i_other_state) = eigenvectors(i,j) + enddo + CI_electronic_energy_dressed(i_state+i_other_state) = eigenvalues(j) + CI_eigenvectors_s2_dressed(i_state+i_other_state) = s2_eigvalues(i_state+i_other_state) + enddo + else + print*,'' + print*,'!!!!!!!! WARNING !!!!!!!!!' + print*,' Within the ',N_det,'determinants selected' + print*,' and the ',N_states_diag,'states requested' + print*,' We did not find any state with S^2 values close to ',expected_s2 + print*,' We will then set the first N_states eigenvectors of the H matrix' + print*,' as the CI_eigenvectors_dressed' + print*,' You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space' + print*,'' + do j=1,min(N_states_diag,N_det) + do i=1,N_det + CI_eigenvectors_dressed(i,j) = eigenvectors(i,j) + enddo + CI_electronic_energy_dressed(j) = eigenvalues(j) + CI_eigenvectors_s2_dressed(j) = s2_eigvalues(j) + enddo + endif + deallocate(index_good_state_array,good_state_array) + deallocate(s2_eigvalues) + else + call u_0_S2_u_0(CI_eigenvectors_s2_dressed,eigenvectors,N_det,psi_det,N_int,& + min(N_det,N_states_diag),size(eigenvectors,1)) + ! Select the "N_states_diag" states of lowest energy + do j=1,min(N_det,N_states_diag) + do i=1,N_det + CI_eigenvectors_dressed(i,j) = eigenvectors(i,j) + enddo + CI_electronic_energy_dressed(j) = eigenvalues(j) + enddo + endif + deallocate(eigenvectors,eigenvalues) + endif END_PROVIDER diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg index 662c6fbb..9eefa66c 100644 --- a/src/determinants/EZFIO.cfg +++ b/src/determinants/EZFIO.cfg @@ -42,13 +42,13 @@ default: 2 [weight_selection] type: integer -doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: rPT2 matching, 3: variance matching, 4: variance and rPT2 matching, 5: variance minimization and matching, 6: CI coefficients 7: input state-average multiplied by variance and rPT2 matching 8: input state-average multiplied by rPT2 matching 9: input state-average multiplied by variance matching +doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: PT2 matching, 3: variance matching, 4: variance and PT2 matching, 5: variance minimization and matching, 6: CI coefficients 7: input state-average multiplied by variance and PT2 matching 8: input state-average multiplied by PT2 matching 9: input state-average multiplied by variance matching interface: ezfio,provider,ocaml default: 1 [threshold_generators] type: Threshold -doc: Thresholds on generators (fraction of the square of the norm) +doc: Thresholds on generators (fraction of the square of the norm) interface: ezfio,provider,ocaml default: 0.999 @@ -80,7 +80,7 @@ type: integer [psi_coef] interface: ezfio doc: Coefficients of the wave function -type: double precision +type: double precision size: (determinants.n_det,determinants.n_states) [psi_det] @@ -92,7 +92,7 @@ size: (determinants.n_int*determinants.bit_kind/8,2,determinants.n_det) [psi_coef_qp_edit] interface: ezfio doc: Coefficients of the wave function -type: double precision +type: double precision size: (determinants.n_det_qp_edit,determinants.n_states) [psi_det_qp_edit] @@ -126,13 +126,18 @@ default: 1. [thresh_sym] type: Threshold -doc: Thresholds to check if a determinant is connected with HF +doc: Thresholds to check if a determinant is connected with HF interface: ezfio,provider,ocaml default: 1.e-15 [pseudo_sym] type: logical -doc: If |true|, discard any Slater determinants with an interaction smaller than thresh_sym with HF. +doc: If |true|, discard any Slater determinants with an interaction smaller than thresh_sym with HF. interface: ezfio,provider,ocaml default: False +[save_threshold] +type: Threshold +doc: Cut-off to apply to the CI coefficients when the wave function is stored +interface: ezfio,provider,ocaml +default: 1.e-14 diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f index 6a474662..1a1d92b5 100644 --- a/src/determinants/density_matrix.irp.f +++ b/src/determinants/density_matrix.irp.f @@ -268,6 +268,44 @@ subroutine set_natural_mos soft_touch mo_occ end + +subroutine save_natural_mos_canon_label + implicit none + BEGIN_DOC + ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in + ! the |MO| basis + END_DOC + call set_natural_mos_canon_label + call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10) + call orthonormalize_mos + call save_mos +end + +subroutine set_natural_mos_canon_label + implicit none + BEGIN_DOC + ! Set natural orbitals, obtained by diagonalization of the one-body density matrix + ! in the |MO| basis + END_DOC + character*(64) :: label + double precision, allocatable :: tmp(:,:) + + label = "Canonical" + integer :: i,j,iorb,jorb + do i = 1, n_virt_orb + iorb = list_virt(i) + do j = 1, n_core_inact_act_orb + jorb = list_core_inact_act(j) + enddo + enddo + call mo_as_svd_vectors_of_mo_matrix_eig(one_e_dm_mo,size(one_e_dm_mo,1),mo_num,mo_num,mo_occ,label) + soft_touch mo_occ +end + + + + + subroutine set_natorb_no_ov_rot implicit none BEGIN_DOC @@ -330,12 +368,12 @@ BEGIN_PROVIDER [ double precision, c0_weight, (N_states) ] c = maxval(psi_coef(:,i) * psi_coef(:,i)) c0_weight(i) = 1.d0/(c+1.d-20) enddo - c = 1.d0/minval(c0_weight(:)) + c = 1.d0/sum(c0_weight(:)) do i=1,N_states c0_weight(i) = c0_weight(i) * c enddo else - c0_weight = 1.d0 + c0_weight(:) = 1.d0 endif END_PROVIDER @@ -352,7 +390,7 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ] if (weight_one_e_dm == 0) then state_average_weight(:) = c0_weight(:) else if (weight_one_e_dm == 1) then - state_average_weight(:) = 1./N_states + state_average_weight(:) = 1.d0/N_states else call ezfio_has_determinants_state_average_weight(exists) if (exists) then diff --git a/src/determinants/determinants.irp.f b/src/determinants/determinants.irp.f index 2a6057de..eeadf779 100644 --- a/src/determinants/determinants.irp.f +++ b/src/determinants/determinants.irp.f @@ -77,29 +77,31 @@ BEGIN_PROVIDER [ integer, psi_det_size ] END_DOC PROVIDE ezfio_filename logical :: exists - if (mpi_master) then - call ezfio_has_determinants_n_det(exists) - if (exists) then - call ezfio_get_determinants_n_det(psi_det_size) - else - psi_det_size = 1 + psi_det_size = 1 + PROVIDE mpi_master + if (read_wf) then + if (mpi_master) then + call ezfio_has_determinants_n_det(exists) + if (exists) then + call ezfio_get_determinants_n_det(psi_det_size) + else + psi_det_size = 1 + endif + call write_int(6,psi_det_size,'Dimension of the psi arrays') endif - psi_det_size = max(psi_det_size,100000) - call write_int(6,psi_det_size,'Dimension of the psi arrays') + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST( psi_det_size, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read psi_det_size with MPI' + endif + IRP_ENDIF endif - IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - IRP_ENDIF - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST( psi_det_size, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read psi_det_size with MPI' - endif - IRP_ENDIF - END_PROVIDER diff --git a/src/determinants/dipole_moments.irp.f b/src/determinants/dipole_moments.irp.f index 8a5f1a2d..06fca0cd 100644 --- a/src/determinants/dipole_moments.irp.f +++ b/src/determinants/dipole_moments.irp.f @@ -9,7 +9,7 @@ double precision :: weight, r(3) double precision :: cpu0,cpu1,nuclei_part_z,nuclei_part_y,nuclei_part_x - call cpu_time(cpu0) +! call cpu_time(cpu0) z_dipole_moment = 0.d0 y_dipole_moment = 0.d0 x_dipole_moment = 0.d0 @@ -26,10 +26,10 @@ enddo enddo - print*,'electron part for z_dipole = ',z_dipole_moment - print*,'electron part for y_dipole = ',y_dipole_moment - print*,'electron part for x_dipole = ',x_dipole_moment - +! print*,'electron part for z_dipole = ',z_dipole_moment +! print*,'electron part for y_dipole = ',y_dipole_moment +! print*,'electron part for x_dipole = ',x_dipole_moment +! nuclei_part_z = 0.d0 nuclei_part_y = 0.d0 nuclei_part_x = 0.d0 @@ -38,28 +38,39 @@ nuclei_part_y += nucl_charge(i) * nucl_coord(i,2) nuclei_part_x += nucl_charge(i) * nucl_coord(i,1) enddo - print*,'nuclei part for z_dipole = ',nuclei_part_z - print*,'nuclei part for y_dipole = ',nuclei_part_y - print*,'nuclei part for x_dipole = ',nuclei_part_x - +! print*,'nuclei part for z_dipole = ',nuclei_part_z +! print*,'nuclei part for y_dipole = ',nuclei_part_y +! print*,'nuclei part for x_dipole = ',nuclei_part_x +! do istate = 1, N_states z_dipole_moment(istate) += nuclei_part_z y_dipole_moment(istate) += nuclei_part_y x_dipole_moment(istate) += nuclei_part_x enddo - call cpu_time(cpu1) - print*,'Time to provide the dipole moment :',cpu1-cpu0 +! call cpu_time(cpu1) +! print*,'Time to provide the dipole moment :',cpu1-cpu0 END_PROVIDER - subroutine print_z_dipole_moment_only + subroutine print_dipole_moments implicit none + integer :: i print*, '' print*, '' print*, '****************************************' - print*, 'z_dipole_moment = ',z_dipole_moment + write(*,'(A10)',advance='no') ' State : ' + do i = 1,N_states + write(*,'(i16)',advance='no') i + end do + write(*,*) '' + write(*,'(A17,100(1pE16.8))') 'x_dipole_moment = ',x_dipole_moment + write(*,'(A17,100(1pE16.8))') 'y_dipole_moment = ',y_dipole_moment + write(*,'(A17,100(1pE16.8))') 'z_dipole_moment = ',z_dipole_moment + !print*, 'x_dipole_moment = ',x_dipole_moment + !print*, 'y_dipole_moment = ',y_dipole_moment + !print*, 'z_dipole_moment = ',z_dipole_moment print*, '****************************************' end diff --git a/src/determinants/h_apply.irp.f b/src/determinants/h_apply.irp.f index 98fafb4a..d01ad1c7 100644 --- a/src/determinants/h_apply.irp.f +++ b/src/determinants/h_apply.irp.f @@ -322,10 +322,7 @@ subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) enddo do i=1,n_selected - do j=1,N_int - H_apply_buffer(iproc)%det(j,1,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,1,i) - H_apply_buffer(iproc)%det(j,2,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,2,i) - enddo + H_apply_buffer(iproc)%det(:,:,i+H_apply_buffer(iproc)%N_det) = det_buffer(:,:,i) ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i+H_apply_buffer(iproc)%N_det)) )== elec_alpha_num) ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i+H_apply_buffer(iproc)%N_det))) == elec_beta_num) enddo diff --git a/src/determinants/s2.irp.f b/src/determinants/s2.irp.f index d73b2dbf..2c1a8757 100644 --- a/src/determinants/s2.irp.f +++ b/src/determinants/s2.irp.f @@ -103,13 +103,17 @@ BEGIN_PROVIDER [ double precision, expected_s2] END_PROVIDER -BEGIN_PROVIDER [ double precision, s2_values, (N_states) ] + BEGIN_PROVIDER [ double precision, s2_values, (N_states) ] +&BEGIN_PROVIDER [ double precision, s_values, (N_states) ] implicit none BEGIN_DOC ! array of the averaged values of the S^2 operator on the various states END_DOC integer :: i call u_0_S2_u_0(s2_values,psi_coef,n_det,psi_det,N_int,N_states,psi_det_size) + do i = 1, N_states + s_values(i) = 0.5d0 *(-1.d0 + dsqrt(1.d0 + 4 * s2_values(i))) + enddo END_PROVIDER diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index 04cf861f..3a33a37d 100644 --- a/src/determinants/slater_rules.irp.f +++ b/src/determinants/slater_rules.irp.f @@ -438,7 +438,7 @@ subroutine bitstring_to_list_ab( string, list, n_elements, Nint) use bitmasks implicit none BEGIN_DOC - ! Gives the inidices(+1) of the bits set to 1 in the bit string + ! Gives the indices(+1) of the bits set to 1 in the bit string ! For alpha/beta determinants. END_DOC integer, intent(in) :: Nint @@ -472,6 +472,35 @@ subroutine bitstring_to_list_ab( string, list, n_elements, Nint) end +!subroutine bitstring_to_list( string, list, n_elements, Nint) +! use bitmasks +! implicit none +! BEGIN_DOC +! ! Gives the indices(+1) of the bits set to 1 in the bit string +! END_DOC +! integer, intent(in) :: Nint +! integer(bit_kind), intent(in) :: string(Nint) +! integer, intent(out) :: list(Nint*bit_kind_size) +! integer, intent(out) :: n_elements +! +! integer :: i, j, ishift +! integer(bit_kind) :: l +! +! n_elements = 0 +! ishift = 1 +! do i=1,Nint +! l = string(i) +! do while (l /= 0_bit_kind) +! j = trailz(l) +! n_elements = n_elements + 1 +! l = ibclr(l,j) +! list(n_elements) = ishift+j +! enddo +! ishift = ishift + bit_kind_size +! enddo +! +!end + subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2) use bitmasks diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f index dea4a566..dd55e112 100644 --- a/src/determinants/spindeterminants.irp.f +++ b/src/determinants/spindeterminants.irp.f @@ -585,7 +585,7 @@ END_PROVIDER enddo !$OMP ENDDO !$OMP END PARALLEL - call i8radix_sort(to_sort, psi_bilinear_matrix_transp_order, N_det,-1) + call i8sort(to_sort, psi_bilinear_matrix_transp_order, N_det) call iset_order(psi_bilinear_matrix_transp_rows,psi_bilinear_matrix_transp_order,N_det) call iset_order(psi_bilinear_matrix_transp_columns,psi_bilinear_matrix_transp_order,N_det) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(l) diff --git a/src/determinants/utils.irp.f b/src/determinants/utils.irp.f index 3aec16f9..7b75d985 100644 --- a/src/determinants/utils.irp.f +++ b/src/determinants/utils.irp.f @@ -8,6 +8,7 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] double precision :: hij integer :: degree(N_det),idx(0:N_det) call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij) + print*,'Providing the H_matrix_all_dets ...' !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hij,degree,idx,k) & !$OMP SHARED (N_det, psi_det, N_int,H_matrix_all_dets) do i =1,N_det @@ -18,6 +19,26 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] enddo enddo !$OMP END PARALLEL DO + print*,'H_matrix_all_dets done ' +END_PROVIDER + +BEGIN_PROVIDER [ double precision, H_matrix_diag_all_dets,(N_det) ] + use bitmasks + implicit none + BEGIN_DOC + ! |H| matrix on the basis of the Slater determinants defined by psi_det + END_DOC + integer :: i + double precision :: hij + integer :: degree(N_det) + call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij) + !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,hij,degree) & + !$OMP SHARED (N_det, psi_det, N_int,H_matrix_diag_all_dets) + do i =1,N_det + call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hij) + H_matrix_diag_all_dets(i) = hij + enddo + !$OMP END PARALLEL DO END_PROVIDER diff --git a/src/dft_keywords/EZFIO.cfg b/src/dft_keywords/EZFIO.cfg index b452c863..a13cbac9 100644 --- a/src/dft_keywords/EZFIO.cfg +++ b/src/dft_keywords/EZFIO.cfg @@ -16,3 +16,8 @@ doc: Percentage of HF exchange in the DFT model interface: ezfio,provider,ocaml default: 0. +[mu_dft_type] +type: character*(32) +doc: type of mu(r) for rsdft [ cst ] +interface: ezfio, provider, ocaml +default: cst diff --git a/src/dft_one_e/NEED b/src/dft_one_e/NEED index 3a942f28..615ee97e 100644 --- a/src/dft_one_e/NEED +++ b/src/dft_one_e/NEED @@ -6,3 +6,4 @@ ao_one_e_ints ao_two_e_ints mo_two_e_erf_ints ao_two_e_erf_ints +mu_of_r diff --git a/src/dft_one_e/mu_erf_dft.irp.f b/src/dft_one_e/mu_erf_dft.irp.f index 53effcb6..0b870564 100644 --- a/src/dft_one_e/mu_erf_dft.irp.f +++ b/src/dft_one_e/mu_erf_dft.irp.f @@ -8,3 +8,73 @@ BEGIN_PROVIDER [double precision, mu_erf_dft] mu_erf_dft = mu_erf END_PROVIDER + +BEGIN_PROVIDER [double precision, mu_of_r_dft, (n_points_final_grid)] + implicit none + integer :: i + if(mu_dft_type == "Read")then + call ezfio_get_mu_of_r_mu_of_r_disk(mu_of_r_dft) + else + do i = 1, n_points_final_grid + if(mu_dft_type == "cst")then + mu_of_r_dft(i) = mu_erf_dft + else if(mu_dft_type == "hf")then + mu_of_r_dft(i) = mu_of_r_hf(i) + else if(mu_dft_type == "rsc")then + mu_of_r_dft(i) = mu_rsc_of_r(i) + else if(mu_dft_type == "grad_rho")then + mu_of_r_dft(i) = mu_grad_rho(i) + else + print*,'mu_dft_type is not of good type = ',mu_dft_type + print*,'it must be of type Read, cst, hf, rsc' + print*,'Stopping ...' + stop + endif + enddo + endif +END_PROVIDER + +BEGIN_PROVIDER [double precision, mu_rsc_of_r, (n_points_final_grid)] + implicit none + integer :: i + double precision :: mu_rs_c,rho,r(3), dm_a, dm_b + do i = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + call dm_dft_alpha_beta_at_r(r,dm_a,dm_b) + rho = dm_a + dm_b + mu_rsc_of_r(i) = mu_rs_c(rho) + enddo +END_PROVIDER + +BEGIN_PROVIDER [double precision, mu_grad_rho, (n_points_final_grid)] + implicit none + integer :: i + double precision :: mu_grad_rho_func, r(3) + do i = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + mu_grad_rho(i) = mu_grad_rho_func(r) + enddo +END_PROVIDER + + +BEGIN_PROVIDER [double precision, mu_of_r_dft_average] + implicit none + integer :: i + double precision :: mu_rs_c,rho,r(3), dm_a, dm_b + mu_of_r_dft_average = 0.d0 + do i = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + call dm_dft_alpha_beta_at_r(r,dm_a,dm_b) + rho = dm_a + dm_b + if(mu_of_r_dft(i).gt.1.d+3)cycle + mu_of_r_dft_average += rho * mu_of_r_dft(i) * final_weight_at_r_vector(i) + enddo + mu_of_r_dft_average = mu_of_r_dft_average / dble(elec_alpha_num + elec_beta_num) + print*,'mu_of_r_dft_average = ',mu_of_r_dft_average +END_PROVIDER diff --git a/src/dft_utils_func/mu_of_r_dft.irp.f b/src/dft_utils_func/mu_of_r_dft.irp.f new file mode 100644 index 00000000..0e9a0f1b --- /dev/null +++ b/src/dft_utils_func/mu_of_r_dft.irp.f @@ -0,0 +1,37 @@ +double precision function mu_rs_c(rho) + implicit none + double precision, intent(in) :: rho + include 'constants.include.F' + double precision :: cst_rs,alpha_rs,rs + cst_rs = (4.d0 * dacos(-1.d0)/3.d0)**(-1.d0/3.d0) + alpha_rs = 2.d0 * dsqrt((9.d0 * dacos(-1.d0)/4.d0)**(-1.d0/3.d0)) / sqpi + + rs = cst_rs * rho**(-1.d0/3.d0) + mu_rs_c = alpha_rs/dsqrt(rs) + +end + +double precision function mu_grad_rho_func(r) + implicit none + double precision , intent(in) :: r(3) + integer :: m + double precision :: rho, dm_a, dm_b, grad_dm_a(3), grad_dm_b(3) + double precision :: eta, grad_rho(3), grad_sqr + eta = mu_erf + call density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b) + rho = dm_a + dm_b + do m = 1,3 + grad_rho(m) = grad_dm_a(m) + grad_dm_b(m) + enddo + grad_sqr=0.d0 + do m = 1,3 + grad_sqr=grad_sqr+grad_rho(m)*grad_rho(m) + enddo + grad_sqr = dsqrt(grad_sqr) + if (rho<1.d-12) then + mu_grad_rho_func = 1.d-10 + else + mu_grad_rho_func = eta * grad_sqr / rho + endif + +end diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f index 6fa6a4c7..38478d21 100644 --- a/src/dft_utils_in_r/ao_in_r.irp.f +++ b/src/dft_utils_in_r/ao_in_r.irp.f @@ -91,7 +91,19 @@ enddo END_PROVIDER - BEGIN_PROVIDER[double precision, aos_lapl_in_r_array, (ao_num,n_points_final_grid,3)] + BEGIN_PROVIDER [double precision, aos_lapl_in_r_array_transp, (ao_num, n_points_final_grid,3)] + implicit none + integer :: i,j,m + do i = 1, n_points_final_grid + do j = 1, ao_num + do m = 1, 3 + aos_lapl_in_r_array_transp(j,i,m) = aos_lapl_in_r_array(m,j,i) + enddo + enddo + enddo + END_PROVIDER + + BEGIN_PROVIDER [double precision, aos_lapl_in_r_array, (3,ao_num,n_points_final_grid)] implicit none BEGIN_DOC ! aos_lapl_in_r_array(i,j,k) = value of the kth component of the laplacian of jth ao on the ith grid point @@ -100,20 +112,20 @@ END_DOC integer :: i,j,m double precision :: aos_array(ao_num), r(3) - double precision :: aos_grad_array(ao_num,3) - double precision :: aos_lapl_array(ao_num,3) + double precision :: aos_grad_array(3,ao_num) + double precision :: aos_lapl_array(3,ao_num) !$OMP PARALLEL DO & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,r,aos_array,aos_grad_array,aos_lapl_array,j,m) & !$OMP SHARED(aos_lapl_in_r_array,n_points_final_grid,ao_num,final_grid_points) - do m = 1, 3 - do i = 1, n_points_final_grid - r(1) = final_grid_points(1,i) - r(2) = final_grid_points(2,i) - r(3) = final_grid_points(3,i) - call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array) - do j = 1, ao_num - aos_lapl_in_r_array(j,i,m) = aos_lapl_array(j,m) + do i = 1, n_points_final_grid + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array) + do j = 1, ao_num + do m = 1, 3 + aos_lapl_in_r_array(m,j,i) = aos_lapl_array(m,j) enddo enddo enddo diff --git a/src/dft_utils_in_r/ints_grad.irp.f b/src/dft_utils_in_r/ints_grad.irp.f new file mode 100644 index 00000000..239fe554 --- /dev/null +++ b/src/dft_utils_in_r/ints_grad.irp.f @@ -0,0 +1,39 @@ + BEGIN_PROVIDER [ double precision, mo_grad_ints, (mo_num, mo_num,3)] + implicit none + BEGIN_DOC +! mo_grad_ints(i,j,m) = + END_DOC + integer :: i,j,ipoint,m + double precision :: weight + mo_grad_ints = 0.d0 + do m = 1, 3 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do j = 1, mo_num + do i = 1, mo_num + mo_grad_ints(i,j,m) += mos_grad_in_r_array(j,ipoint,m) * mos_in_r_array(i,ipoint) * weight + enddo + enddo + enddo + enddo + + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, mo_grad_ints_transp, (3,mo_num, mo_num)] + implicit none + BEGIN_DOC +! mo_grad_ints(i,j,m) = + END_DOC + integer :: i,j,ipoint,m + double precision :: weight + do m = 1, 3 + do j = 1, mo_num + do i = 1, mo_num + mo_grad_ints_transp(m,i,j) = mo_grad_ints(i,j,m) + enddo + enddo + enddo + + +END_PROVIDER diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f index 0a8b4d52..192cb25a 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -138,7 +138,7 @@ integer :: m mos_lapl_in_r_array = 0.d0 do m=1,3 - call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_lapl_in_r_array(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_num) + call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_lapl_in_r_array_transp(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_num) enddo END_PROVIDER diff --git a/src/dressing/alpha_factory.irp.f b/src/dressing/alpha_factory.irp.f index 5eeeb1a6..c7adffe3 100644 --- a/src/dressing/alpha_factory.irp.f +++ b/src/dressing/alpha_factory.irp.f @@ -1179,7 +1179,7 @@ subroutine bitstring_to_list_in_selection( string, list, n_elements, Nint) use bitmasks implicit none BEGIN_DOC - ! Gives the inidices(+1) of the bits set to 1 in the bit string + ! Gives the indices(+1) of the bits set to 1 in the bit string END_DOC integer, intent(in) :: Nint integer(bit_kind), intent(in) :: string(Nint) diff --git a/src/dressing/run_dress_slave.irp.f b/src/dressing/run_dress_slave.irp.f index a33fb1dd..08b654c9 100644 --- a/src/dressing/run_dress_slave.irp.f +++ b/src/dressing/run_dress_slave.irp.f @@ -72,7 +72,7 @@ subroutine run_dress_slave(thread,iproce,energy) provide psi_energy ending = dress_N_cp+1 ntask_tbd = 0 - call omp_set_max_active_levels(8) + call set_multiple_levels_omp(.True.) !$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(interesting, breve_delta_m, task_id) & @@ -84,7 +84,7 @@ subroutine run_dress_slave(thread,iproce,energy) zmq_socket_push = new_zmq_push_socket(thread) integer, external :: connect_to_taskserver !$OMP CRITICAL - call omp_set_max_active_levels(1) + call set_multiple_levels_omp(.False.) if (connect_to_taskserver(zmq_to_qp_run_socket,worker_id,thread) == -1) then print *, irp_here, ': Unable to connect to task server' stop -1 @@ -296,7 +296,7 @@ subroutine run_dress_slave(thread,iproce,energy) !$OMP END CRITICAL !$OMP END PARALLEL - call omp_set_max_active_levels(1) + call set_multiple_levels_omp(.False.) ! do i=0,dress_N_cp+1 ! call omp_destroy_lock(lck_sto(i)) ! end do diff --git a/src/ezfio_files/output.irp.f b/src/ezfio_files/output.irp.f index 48512f92..7b2663a0 100644 --- a/src/ezfio_files/output.irp.f +++ b/src/ezfio_files/output.irp.f @@ -25,7 +25,7 @@ subroutine write_time(iunit) ct = ct - output_cpu_time_0 call wall_time(wt) wt = wt - output_wall_time_0 - write(6,'(A,F14.6,A,F14.6,A)') & + write(6,'(A,F14.2,A,F14.2,A)') & '.. >>>>> [ WALL TIME: ', wt, ' s ] [ CPU TIME: ', ct, ' s ] <<<<< ..' write(6,*) end diff --git a/src/fci/40.fci.bats b/src/fci/40.fci.bats index f2c78067..1a0c5507 100644 --- a/src/fci/40.fci.bats +++ b/src/fci/40.fci.bats @@ -59,43 +59,43 @@ function run_stoch() { @test "HCO" { # 12.2868s qp set_file hco.ezfio - run -113.389297812482 6.e-4 100000 + run -113.393356604085 1.e-3 100000 } @test "H2O2" { # 12.9214s qp set_file h2o2.ezfio qp set_mo_class --core="[1-2]" --act="[3-24]" --del="[25-38]" - run -151.00467 1.e-4 100000 + run -151.005848404095 1.e-3 100000 } @test "HBO" { # 13.3144s [[ -n $TRAVIS ]] && skip qp set_file hbo.ezfio - run -100.212560384678 1.e-3 100000 + run -100.214099486337 1.e-3 100000 } @test "H2O" { # 11.3727s [[ -n $TRAVIS ]] && skip qp set_file h2o.ezfio - run -76.2361605151999 3.e-4 100000 + run -76.2361605151999 5.e-4 100000 } @test "ClO" { # 13.3755s [[ -n $TRAVIS ]] && skip qp set_file clo.ezfio - run -534.545616787223 3.e-4 100000 + run -534.546453546852 1.e-3 100000 } @test "SO" { # 13.4952s [[ -n $TRAVIS ]] && skip qp set_file so.ezfio - run -26.0096209515081 1.e-3 100000 + run -26.0176563764039 1.e-3 100000 } @test "H2S" { # 13.6745s [[ -n $TRAVIS ]] && skip qp set_file h2s.ezfio - run -398.859168655255 3.e-4 100000 + run -398.859577605891 5.e-4 100000 } @test "OH" { # 13.865s @@ -113,13 +113,13 @@ function run_stoch() { @test "H3COH" { # 14.7299s [[ -n $TRAVIS ]] && skip qp set_file h3coh.ezfio - run -115.205191406072 3.e-4 100000 + run -115.205632960026 1.e-3 100000 } @test "SiH3" { # 15.99s [[ -n $TRAVIS ]] && skip qp set_file sih3.ezfio - run -5.57241217753818 3.e-4 100000 + run -5.57241217753818 5.e-4 100000 } @test "CH4" { # 16.1612s @@ -132,28 +132,28 @@ function run_stoch() { @test "ClF" { # 16.8864s [[ -n $TRAVIS ]] && skip qp set_file clf.ezfio - run -559.169313755572 3.e-4 100000 + run -559.169748890031 1.5e-3 100000 } @test "SO2" { # 17.5645s [[ -n $TRAVIS ]] && skip qp set_file so2.ezfio qp set_mo_class --core="[1-8]" --act="[9-87]" - run -41.5746738713298 3.e-4 100000 + run -41.5746738713298 1.5e-3 100000 } @test "C2H2" { # 17.6827s [[ -n $TRAVIS ]] && skip qp set_file c2h2.ezfio qp set_mo_class --act="[1-30]" --del="[31-36]" - run -12.3685464085969 3.e-4 100000 + run -12.3685464085969 2.e-3 100000 } @test "N2" { # 18.0198s [[ -n $TRAVIS ]] && skip qp set_file n2.ezfio qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-60]" - run -109.28681540699360 3.e-4 100000 + run -109.287917088107 1.5e-3 100000 } @test "N2H4" { # 18.5006s @@ -167,7 +167,7 @@ function run_stoch() { [[ -n $TRAVIS ]] && skip qp set_file co2.ezfio qp set_mo_class --core="[1,2]" --act="[3-30]" --del="[31-42]" - run -187.968547952413 3.e-4 100000 + run -187.970184372047 1.5e-3 100000 } @@ -182,6 +182,6 @@ function run_stoch() { [[ -n $TRAVIS ]] && skip qp set_file hcn.ezfio qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-55]" - run -93.0771143355433 3.e-4 100000 + run -93.0777619629755 1.e-3 100000 } diff --git a/src/functionals/sr_lda.irp.f b/src/functionals/sr_lda.irp.f index 965a744c..bd062a02 100644 --- a/src/functionals/sr_lda.irp.f +++ b/src/functionals/sr_lda.irp.f @@ -21,7 +21,9 @@ weight = final_weight_at_r_vector(i) rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate) rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate) - call ex_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_x,vx_a,vx_b) + double precision :: mu_local + mu_local = mu_of_r_dft(i) + call ex_lda_sr(mu_local,rhoa(istate),rhob(istate),e_x,vx_a,vx_b) energy_x_sr_lda(istate) += weight * e_x enddo enddo @@ -48,7 +50,9 @@ weight = final_weight_at_r_vector(i) rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate) rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate) - call ec_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_c,vc_a,vc_b) + double precision :: mu_local + mu_local = mu_of_r_dft(i) + call ec_lda_sr(mu_local,rhoa(istate),rhob(istate),e_c,vc_a,vc_b) energy_c_sr_lda(istate) += weight * e_c enddo enddo @@ -122,8 +126,10 @@ END_PROVIDER weight = final_weight_at_r_vector(i) rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate) rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate) - call ec_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_c,sr_vc_a,sr_vc_b) - call ex_lda_sr(mu_erf_dft,rhoa(istate),rhob(istate),e_x,sr_vx_a,sr_vx_b) + double precision :: mu_local + mu_local = mu_of_r_dft(i) + call ec_lda_sr(mu_local,rhoa(istate),rhob(istate),e_c,sr_vc_a,sr_vc_b) + call ex_lda_sr(mu_local,rhoa(istate),rhob(istate),e_x,sr_vx_a,sr_vx_b) do j =1, ao_num aos_sr_vc_alpha_lda_w(j,i,istate) = sr_vc_a * aos_in_r_array(j,i)*weight aos_sr_vc_beta_lda_w(j,i,istate) = sr_vc_b * aos_in_r_array(j,i)*weight @@ -147,8 +153,6 @@ END_PROVIDER double precision :: mu,weight double precision :: e_c,sr_vc_a,sr_vc_b,e_x,sr_vx_a,sr_vx_b double precision, allocatable :: rhoa(:),rhob(:) - double precision :: mu_local - mu_local = mu_erf_dft allocate(rhoa(N_states), rhob(N_states)) do istate = 1, N_states do i = 1, n_points_final_grid @@ -158,6 +162,8 @@ END_PROVIDER weight = final_weight_at_r_vector(i) rhoa(istate) = one_e_dm_and_grad_alpha_in_r(4,i,istate) rhob(istate) = one_e_dm_and_grad_beta_in_r(4,i,istate) + double precision :: mu_local + mu_local = mu_of_r_dft(i) call ec_lda_sr(mu_local,rhoa(istate),rhob(istate),e_c,sr_vc_a,sr_vc_b) call ex_lda_sr(mu_local,rhoa(istate),rhob(istate),e_x,sr_vx_a,sr_vx_b) do j =1, ao_num diff --git a/src/functionals/sr_pbe.irp.f b/src/functionals/sr_pbe.irp.f index 93c51067..7053cfb6 100644 --- a/src/functionals/sr_pbe.irp.f +++ b/src/functionals/sr_pbe.irp.f @@ -35,9 +35,11 @@ grad_rho_b_2 += grad_rho_b(m) * grad_rho_b(m) grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m) enddo - + + double precision :: mu_local + mu_local = mu_of_r_dft(i) ! inputs - call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange + call GGA_sr_type_functionals(mu_local,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b ) energy_x_sr_pbe(istate) += ex * weight @@ -135,8 +137,10 @@ END_PROVIDER grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m) enddo + double precision :: mu_local + mu_local = mu_of_r_dft(i) ! inputs - call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange + call GGA_sr_type_functionals(mu_local,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b ) vx_rho_a *= weight @@ -292,8 +296,10 @@ END_PROVIDER grad_rho_a_b += grad_rho_a(m) * grad_rho_b(m) enddo + double precision :: mu_local + mu_local = mu_of_r_dft(i) ! inputs - call GGA_sr_type_functionals(mu_erf_dft,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange + call GGA_sr_type_functionals(mu_local,rho_a,rho_b,grad_rho_a_2,grad_rho_b_2,grad_rho_a_b, & ! outputs exchange ex,vx_rho_a,vx_rho_b,vx_grad_rho_a_2,vx_grad_rho_b_2,vx_grad_rho_a_b, & ! outputs correlation ec,vc_rho_a,vc_rho_b,vc_grad_rho_a_2,vc_grad_rho_b_2,vc_grad_rho_a_b ) vx_rho_a *= weight diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index 641ee209..8e6285e2 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -98,7 +98,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s enddo endif - call print_energy_components() +! call print_energy_components() end subroutine diff --git a/src/kohn_sham_rs/rs_ks_scf.irp.f b/src/kohn_sham_rs/rs_ks_scf.irp.f index 5d23544e..84b85136 100644 --- a/src/kohn_sham_rs/rs_ks_scf.irp.f +++ b/src/kohn_sham_rs/rs_ks_scf.irp.f @@ -17,7 +17,7 @@ program rs_ks_scf print*, '**************************' print*, 'mu_erf_dft = ',mu_erf_dft print*, '**************************' - call check_coherence_functional +! call check_coherence_functional call create_guess call orthonormalize_mos call run diff --git a/src/mo_basis/mos_in_r.irp.f b/src/mo_basis/mos_in_r.irp.f index ee2795d0..e5d3b243 100644 --- a/src/mo_basis/mos_in_r.irp.f +++ b/src/mo_basis/mos_in_r.irp.f @@ -1,6 +1,9 @@ subroutine give_all_mos_at_r(r,mos_array) implicit none + BEGIN_DOC +! mos_array(i) = ith MO function evaluated at "r" + END_DOC double precision, intent(in) :: r(3) double precision, intent(out) :: mos_array(mo_num) double precision :: aos_array(ao_num) diff --git a/src/mo_guess/h_core_guess_routine.irp.f b/src/mo_guess/h_core_guess_routine.irp.f index cbf23a9a..fcbdde49 100644 --- a/src/mo_guess/h_core_guess_routine.irp.f +++ b/src/mo_guess/h_core_guess_routine.irp.f @@ -7,7 +7,7 @@ subroutine hcore_guess label = 'Guess' call mo_as_eigvectors_of_mo_matrix(mo_one_e_integrals, & size(mo_one_e_integrals,1), & - size(mo_one_e_integrals,2),label,1,.false.) + size(mo_one_e_integrals,2),label,1,.true.) call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12 ) call save_mos TOUCH mo_coef mo_label diff --git a/src/mo_two_e_ints/map_integrals.irp.f b/src/mo_two_e_ints/map_integrals.irp.f index 9f73d518..272916e3 100644 --- a/src/mo_two_e_ints/map_integrals.irp.f +++ b/src/mo_two_e_ints/map_integrals.irp.f @@ -302,21 +302,21 @@ end integer(key_kind) :: idx double precision :: tmp - icount = 1 ! Avoid division by zero - do j=1,mo_num - do i=1,j-1 - call two_e_integrals_index(i,j,j,i,idx) - !DIR$ FORCEINLINE - call map_get(mo_integrals_map,idx,tmp) - banned_excitation(i,j) = dabs(tmp) < 1.d-14 - banned_excitation(j,i) = banned_excitation(i,j) - if (banned_excitation(i,j)) icount = icount+2 - enddo - enddo - use_banned_excitation = (mo_num*mo_num) / icount <= 100 !1% - if (use_banned_excitation) then - print *, 'Using sparsity of exchange integrals' - endif +!icount = 1 ! Avoid division by zero +!do j=1,mo_num +! do i=1,j-1 +! call two_e_integrals_index(i,j,j,i,idx) +! !DIR$ FORCEINLINE +! call map_get(mo_integrals_map,idx,tmp) +! banned_excitation(i,j) = dabs(tmp) < 1.d-14 +! banned_excitation(j,i) = banned_excitation(i,j) +! if (banned_excitation(i,j)) icount = icount+2 +! enddo +!enddo +!use_banned_excitation = (mo_num*mo_num) / icount <= 100 !1% +!if (use_banned_excitation) then +! print *, 'Using sparsity of exchange integrals' +!endif END_PROVIDER diff --git a/src/pseudo/EZFIO.cfg b/src/pseudo/EZFIO.cfg index 952ebe04..2b76a5c3 100644 --- a/src/pseudo/EZFIO.cfg +++ b/src/pseudo/EZFIO.cfg @@ -10,7 +10,7 @@ type:integer interface: ezfio,provider [pseudo_n_k] -doc: Number of gaussians in the local component +doc: Powers of r - 2 in the local component type: integer interface: ezfio,provider size: (nuclei.nucl_num,pseudo.pseudo_klocmax) @@ -38,7 +38,7 @@ type:integer interface: ezfio,provider [pseudo_n_kl] -doc: Number of functions in the non-local component +doc: Powers of r - 2 in the non-local component type: integer interface: ezfio,provider size: (nuclei.nucl_num,pseudo.pseudo_kmax,0:pseudo.pseudo_lmax) @@ -69,7 +69,7 @@ default: 1000 [pseudo_grid_rmax] type: double precision -doc: R_max of the QMC grid +doc: R_max of the QMC grid interface: ezfio,provider,ocaml default: 10.0 diff --git a/src/tools/NEED b/src/tools/NEED index c07c9109..0f4e17b0 100644 --- a/src/tools/NEED +++ b/src/tools/NEED @@ -2,3 +2,4 @@ fci mo_two_e_erf_ints aux_quantities hartree_fock +two_body_rdm diff --git a/src/tools/molden.irp.f b/src/tools/molden.irp.f index 417b25ad..830a141e 100644 --- a/src/tools/molden.irp.f +++ b/src/tools/molden.irp.f @@ -52,8 +52,8 @@ program molden l += 1 if (l > ao_num) exit enddo - write(i_unit_output,*)'' enddo + write(i_unit_output,*)'' enddo diff --git a/src/tools/print_dipole.irp.f b/src/tools/print_dipole.irp.f index 8351308e..8db9aa09 100644 --- a/src/tools/print_dipole.irp.f +++ b/src/tools/print_dipole.irp.f @@ -1,5 +1,7 @@ program print_dipole implicit none - call print_z_dipole_moment_only + read_wf = .True. + TOUCH read_wf + call print_dipole_moments end diff --git a/src/tools/print_wf.irp.f b/src/tools/print_wf.irp.f index 7e51caaf..64eb1a1f 100644 --- a/src/tools/print_wf.irp.f +++ b/src/tools/print_wf.irp.f @@ -32,8 +32,9 @@ subroutine routine double precision :: norm_mono_a,norm_mono_b double precision :: norm_mono_a_2,norm_mono_b_2 double precision :: norm_mono_a_pert_2,norm_mono_b_pert_2 - double precision :: norm_mono_a_pert,norm_mono_b_pert + double precision :: norm_mono_a_pert,norm_mono_b_pert,norm_double_1 double precision :: delta_e,coef_2_2 + norm_mono_a = 0.d0 norm_mono_b = 0.d0 norm_mono_a_2 = 0.d0 @@ -42,6 +43,7 @@ subroutine routine norm_mono_b_pert = 0.d0 norm_mono_a_pert_2 = 0.d0 norm_mono_b_pert_2 = 0.d0 + norm_double_1 = 0.d0 do i = 1, min(N_det_print_wf,N_det) print*,'' print*,'i = ',i @@ -93,6 +95,7 @@ subroutine routine print*,'h1,p1 = ',h1,p1 print*,'s2',s2 print*,'h2,p2 = ',h2,p2 + norm_double_1 += dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1)) endif print*,' = ',hij @@ -109,6 +112,7 @@ subroutine routine print*,'' print*,'L1 norm of mono alpha = ',norm_mono_a print*,'L1 norm of mono beta = ',norm_mono_b + print*,'L1 norm of double exc = ',norm_double_1 print*, '---' print*,'L2 norm of mono alpha = ',norm_mono_a_2 print*,'L2 norm of mono beta = ',norm_mono_b_2 diff --git a/src/tools/save_natorb_no_ref.irp.f b/src/tools/save_natorb_no_ref.irp.f new file mode 100644 index 00000000..9d253fa0 --- /dev/null +++ b/src/tools/save_natorb_no_ref.irp.f @@ -0,0 +1,24 @@ +program save_natorb + implicit none + BEGIN_DOC +! Save natural |MOs| into the |EZFIO|. +! +! This program reads the wave function stored in the |EZFIO| directory, +! extracts the corresponding natural orbitals and setd them as the new +! |MOs|. +! +! If this is a multi-state calculation, the density matrix that produces +! the natural orbitals is obtained from an average of the density +! matrices of each state with the corresponding +! :option:`determinants state_average_weight` + END_DOC + read_wf = .True. + touch read_wf + call save_natural_mos_canon_label + call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals('None') + call ezfio_set_mo_one_e_ints_io_mo_one_e_integrals('None') + call ezfio_set_mo_one_e_ints_io_mo_integrals_kinetic('None') + call ezfio_set_mo_one_e_ints_io_mo_integrals_n_e('None') + call ezfio_set_mo_one_e_ints_io_mo_integrals_pseudo('None') +end + diff --git a/src/tools/truncate_wf.irp.f b/src/tools/truncate_wf.irp.f new file mode 100644 index 00000000..6c66c8ec --- /dev/null +++ b/src/tools/truncate_wf.irp.f @@ -0,0 +1,98 @@ +program truncate_wf + implicit none + BEGIN_DOC +! Truncate the wave function + END_DOC + read_wf = .True. + if (s2_eig) then + call routine_s2 + else + call routine + endif +end + +subroutine routine + implicit none + integer :: ndet_max + print*, 'Max number of determinants ?' + read(5,*) ndet_max + integer(bit_kind), allocatable :: psi_det_tmp(:,:,:) + double precision, allocatable :: psi_coef_tmp(:,:) + allocate(psi_det_tmp(N_int,2,ndet_max),psi_coef_tmp(ndet_max, N_states)) + + integer :: i,j + double precision :: accu(N_states) + accu = 0.d0 + do i = 1, ndet_max + do j = 1, N_int + psi_det_tmp(j,1,i) = psi_det_sorted(j,1,i) + psi_det_tmp(j,2,i) = psi_det_sorted(j,2,i) + enddo + do j = 1, N_states + psi_coef_tmp(i,j) = psi_coef_sorted(i,j) + accu(j) += psi_coef_tmp(i,j) **2 + enddo + enddo + do j = 1, N_states + accu(j) = 1.d0/dsqrt(accu(j)) + enddo + do j = 1, N_states + do i = 1, ndet_max + psi_coef_tmp(i,j) = psi_coef_tmp(i,j) * accu(j) + enddo + enddo + + call save_wavefunction_general(ndet_max,N_states,psi_det_tmp,size(psi_coef_tmp,1),psi_coef_tmp) + +end + +subroutine routine_s2 + implicit none + integer :: ndet_max + double precision :: wmin + integer(bit_kind), allocatable :: psi_det_tmp(:,:,:) + double precision, allocatable :: psi_coef_tmp(:,:) + integer :: i,j,k + double precision :: accu(N_states) + + print *, 'Weights of the CFG' + do i=1,N_det + print *, i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:))) + enddo + print*, 'Min weight of the configuration?' + read(5,*) wmin + + ndet_max = 0 + do i=1,N_det + if (maxval(weight_configuration( det_to_configuration(i),:)) < wmin) cycle + ndet_max = ndet_max+1 + enddo + + allocate(psi_det_tmp(N_int,2,ndet_max),psi_coef_tmp(ndet_max, N_states)) + + accu = 0.d0 + k=0 + do i = 1, N_det + if (maxval(weight_configuration( det_to_configuration(i),:)) < wmin) cycle + k = k+1 + do j = 1, N_int + psi_det_tmp(j,1,k) = psi_det(j,1,i) + psi_det_tmp(j,2,k) = psi_det(j,2,i) + enddo + do j = 1, N_states + psi_coef_tmp(k,j) = psi_coef(i,j) + accu(j) += psi_coef_tmp(k,j)**2 + enddo + enddo + do j = 1, N_states + accu(j) = 1.d0/dsqrt(accu(j)) + enddo + do j = 1, N_states + do i = 1, ndet_max + psi_coef_tmp(i,j) = psi_coef_tmp(i,j) * accu(j) + enddo + enddo + + call save_wavefunction_general(ndet_max,N_states,psi_det_tmp,size(psi_coef_tmp,1),psi_coef_tmp) + +end diff --git a/src/davidson/print_e_components.irp.f b/src/two_body_rdm/print_e_components.irp.f similarity index 100% rename from src/davidson/print_e_components.irp.f rename to src/two_body_rdm/print_e_components.irp.f diff --git a/src/two_body_rdm/two_e_dm_mo.irp.f b/src/two_body_rdm/two_e_dm_mo.irp.f index 4dadd2e6..a4dea15f 100644 --- a/src/two_body_rdm/two_e_dm_mo.irp.f +++ b/src/two_body_rdm/two_e_dm_mo.irp.f @@ -1,9 +1,8 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)] implicit none BEGIN_DOC - ! two_e_dm_bb_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of beta/beta electrons - ! - ! + ! \sum_{\sigma \sigma'} + ! ! ! where the indices (i,j,k,l) belong to all MOs. ! @@ -12,7 +11,7 @@ BEGIN_PROVIDER [double precision, two_e_dm_mo, (mo_num,mo_num,mo_num,mo_num)] ! !!!!! WARNING !!!!! IF "no_core_density" then all elements involving at least one CORE MO are set to zero ! The state-averaged two-electron energy : ! - ! \sum_{i,j,k,l = 1, mo_num} two_e_dm_mo(i,j,k,l) * < ii jj | kk ll > + ! \sum_{i,j,k,l = 1, mo_num} two_e_dm_mo(i,j,k,l) * < kk ll | ii jj > END_DOC two_e_dm_mo = 0.d0 integer :: i,j,k,l,iorb,jorb,korb,lorb,istate diff --git a/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f b/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f index eb247dea..26ed5ae6 100644 --- a/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f +++ b/src/two_rdm_routines/davidson_like_state_av_2rdm.irp.f @@ -529,10 +529,14 @@ subroutine orb_range_2_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_ c_average += c_1(l) * c_1(l) * state_weights(l) enddo - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + if (nkeys > 0) then + call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + endif nkeys = 0 call orb_range_diag_to_all_2_rdm_dm_buffer(tmp_det,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) - call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + if (nkeys > 0) then + call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm) + endif nkeys = 0 end do diff --git a/src/utils/memory.irp.f b/src/utils/memory.irp.f index 3ea242b0..d5a066a1 100644 --- a/src/utils/memory.irp.f +++ b/src/utils/memory.irp.f @@ -114,7 +114,7 @@ subroutine print_memory_usage() call resident_memory(rss) call total_memory(mem) - write(*,'(A,F14.6,A,F14.6,A)') & + write(*,'(A,F14.3,A,F14.3,A)') & '.. >>>>> [ RES MEM : ', rss , & ' GB ] [ VIRT MEM : ', mem, ' GB ] <<<<< ..' end diff --git a/src/utils/set_multiple_levels_omp.irp.f b/src/utils/set_multiple_levels_omp.irp.f new file mode 100644 index 00000000..572a13f4 --- /dev/null +++ b/src/utils/set_multiple_levels_omp.irp.f @@ -0,0 +1,26 @@ +subroutine set_multiple_levels_omp(activate) + + BEGIN_DOC +! If true, activate OpenMP nested parallelism. If false, deactivate. + END_DOC + + implicit none + logical, intent(in) :: activate + + if (activate) then + call omp_set_max_active_levels(3) + + IRP_IF SET_NESTED + call omp_set_nested(.True.) + IRP_ENDIF + + else + + call omp_set_max_active_levels(1) + + IRP_IF SET_NESTED + call omp_set_nested(.False.) + IRP_ENDIF + end if + +end diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index a63eb4a3..ff40263c 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -356,7 +356,8 @@ BEGIN_TEMPLATE if ( isize < 32) then call insertion_$Xsort(x,iorder,isize) else - call $Xradix_sort(x,iorder,isize,-1) +! call $Xradix_sort(x,iorder,isize,-1) + call quick_$Xsort(x,iorder,isize) endif end subroutine $Xsort @@ -450,7 +451,8 @@ BEGIN_TEMPLATE if ( isize < 32) then call insertion_$Xsort(x,iorder,isize) else - call $Xradix_sort(x,iorder,isize,-1) +! call $Xradix_sort(x,iorder,isize,-1) + call quick_$Xsort(x,iorder,isize) endif end subroutine $Xsort