diff --git a/.gitmodules b/.gitmodules index 73299cba..d250eddf 100644 --- a/.gitmodules +++ b/.gitmodules @@ -4,3 +4,6 @@ [submodule "external/irpf90"] path = external/irpf90 url = https://gitlab.com/scemama/irpf90.git +[submodule "external/qp2-dependencies"] + path = external/qp2-dependencies + url = https://github.com/QuantumPackage/qp2-dependencies.git diff --git a/AUTHORS b/AUTHORS index 4e7cce4b..18d99e25 100644 --- a/AUTHORS +++ b/AUTHORS @@ -4,12 +4,14 @@ - Thomas Applencourt - Anouar Benali - Michel Caffarel +- Vijay Gopal Chilkuri +- Yann Damour - Grégoire David - Anthony Ferté - Yann Garniron - Kevin Gasperich -- Vijay Gopal Chilkuri - Emmanuel Giner +- Fabris Kossoski - Pierre-François Loos - Jean-Paul Malrieu - Julien Paquier diff --git a/INSTALL.rst b/INSTALL.rst index e12142db..229bf40a 100644 --- a/INSTALL.rst +++ b/INSTALL.rst @@ -20,13 +20,15 @@ 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*. +*Quantum Package*. 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 +internet connection is not needed any more. -Running this script will also tell you which external dependencies are missing -and need to be installed. - -When all dependencies have been installed, ( the :command:`configure` will tell you) -source the :file:`quantum_package.rc` in order to load all environment variables and compile the |QP|. +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|. Now all the requirements are met, you can compile the programs using @@ -51,8 +53,6 @@ Requirements - |ZeroMQ| : networking library - `GMP `_ : Gnu Multiple Precision Arithmetic Library - |OCaml| compiler with |OPAM| package manager -- `Bubblewrap `_ : Sandboxing tool required by Opam -- `libcap `_ : POSIX capabilities required by Bubblewrap - |Ninja| : a parallel build system - |pkg-config| : a tool which returns information about installed libraries @@ -95,9 +95,7 @@ The following packages are supported by the :command:`configure` installer: * zeromq * f77zmq * gmp -* libcap -* bwrap -* ocaml ( :math:`\approx` 10 minutes) +* ocaml (:math:`\approx` 5 minutes) * ezfio * docopt * resultsFile @@ -111,19 +109,21 @@ Example: .. 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|: + 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 + QP>> Please install it here: /your_quantum_package_directory/bin - So just enter the proposition of the |QP| and press enter. + 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 to install a specific +dependency, there are some proposition of how to download and install the +minimal dependencies to compile and use the |QP|. Before doing anything below, try to install the packages with your package manager diff --git a/RELEASE_NOTES.org b/RELEASE_NOTES.org index 758618a9..98830f3f 100644 --- a/RELEASE_NOTES.org +++ b/RELEASE_NOTES.org @@ -30,6 +30,7 @@ - Fixed bug in DIIS - Fixed bug in molden (Au -> Angs) - Fixed bug with non-contiguous MOs in active space and deleter MOs + - Complete network-free installation *** User interface @@ -83,9 +84,7 @@ - Added LIB file to add extra libs in plugin - Using Intel IPP for sorting when using Intel compiler - Removed parallelism in sorting - - ao_one_e_integral_zero - banned_excitations + - Compute banned_excitations from exchange integrals to accelerate with local MOs diff --git a/configure b/configure index b45bfd27..0debde04 100755 --- a/configure +++ b/configure @@ -3,7 +3,7 @@ # Quantum Package configuration script # -TEMP=$(getopt -o c:i:h -l config:,install:,help -n $0 -- "$@") || exit 1 +TEMP=$(getopt -o d:c:i:h -l download:,config:,install:,help -n $0 -- "$@") || exit 1 eval set -- "$TEMP" export QP_ROOT="$( cd "$(dirname "$0")" ; pwd -P )" @@ -18,20 +18,6 @@ export CC=gcc git submodule init git submodule update -# /!\ When updating version, update also etc files -BATS_URL="https://github.com/bats-core/bats-core/archive/v1.1.0.tar.gz" -BUBBLE_URL="https://github.com/projectatomic/bubblewrap/releases/download/v0.3.3/bubblewrap-0.3.3.tar.xz" -DOCOPT_URL="https://github.com/docopt/docopt/archive/0.6.2.tar.gz" -BSE_URL="https://github.com/MolSSI-BSE/basis_set_exchange/archive/v0.8.11.tar.gz" -F77ZMQ_URL="https://github.com/scemama/f77_zmq/archive/v4.2.5.tar.gz" -LIBCAP_URL="https://git.kernel.org/pub/scm/linux/kernel/git/morgan/libcap.git/snapshot/libcap-2.25.tar.gz" -NINJA_URL="https://github.com/ninja-build/ninja/releases/download/v1.8.2/ninja-linux.zip" -OCAML_URL="https://raw.githubusercontent.com/ocaml/opam/master/shell/install.sh" -RESULTS_URL="https://gitlab.com/scemama/resultsFile/-/archive/v2.0/resultsFile-v2.0.tar.gz" -ZEROMQ_URL="https://github.com/zeromq/libzmq/releases/download/v4.2.5/zeromq-4.2.5.tar.gz" -ZLIB_URL="https://www.zlib.net/zlib-1.2.11.tar.gz" - - function help() { cat < /dev/null 2> /dev/null || true -# -# opam switch create ocaml-base-compiler.4.11.1 - opam init --verbose --yes --compiler=4.11.1 --disable-sandboxing - - eval $(opam env) - opam install -y ${OCAML_PACKAGES} || exit 1 - - else - # Conventional commands - execute << EOF - chmod +x "${QP_ROOT}"/external/opam_installer.sh - "${QP_ROOT}"/external/opam_installer.sh --no-backup -EOF - execute << EOF - rm --force ${QP_ROOT}/bin/opam - export OPAMROOT=${OPAMROOT:-${QP_ROOT}/external/opam} - echo ${QP_ROOT}/bin \ - | sh ${QP_ROOT}/external/opam_installer.sh -EOF - rm ${QP_ROOT}/external/opam_installer.sh -# source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true -# opam switch create ocaml-base-compiler.4.11.1 || exit 1 - - opam init --verbose --yes --compiler=4.11.1 --disable-sandboxing - eval $(opam env) - execute << EOF - opam install -y \${OCAML_PACKAGES} || exit 1 -EOF - fi - - elif [[ ${PACKAGE} = bse ]] ; then - download ${BSE_URL} "${QP_ROOT}"/external/bse.tar.gz execute << EOF cd "\${QP_ROOT}"/external - tar --gunzip --extract --file bse.tar.gz + tar --gunzip --extract --file qp2-dependencies/bse-v0.8.11.tar.gz pip install -e basis_set_exchange-* EOF + elif [[ ${PACKAGE} = zlib ]] ; then - download ${ZLIB_URL} "${QP_ROOT}"/external/zlib.tar.gz execute << EOF cd "\${QP_ROOT}"/external - tar --gunzip --extract --file zlib.tar.gz - rm zlib.tar.gz && \ + tar --gunzip --extract --file qp2-dependencies/zlib-1.2.11.tar.gz cd zlib-*/ ./configure --prefix=${QP_ROOT} && \ make && make install @@ -372,33 +268,27 @@ EOF elif [[ ${PACKAGE} = docopt ]] ; then - download ${DOCOPT_URL} "${QP_ROOT}"/external/docopt.tar.gz execute << EOF cd "\${QP_ROOT}"/external - tar --gunzip --extract --file docopt.tar.gz + tar --gunzip --extract --file qp2-dependencies/docopt-0.6.2.tar.gz mv docopt-*/docopt.py "\${QP_ROOT}/external/Python" - rm --recursive --force -- docopt-*/ docopt.tar.gz EOF elif [[ ${PACKAGE} = resultsFile ]] ; then - download ${RESULTS_URL} "${QP_ROOT}"/external/resultsFile.tar.gz execute << EOF cd "\${QP_ROOT}"/external - tar --gunzip --extract --file resultsFile.tar.gz + tar --gunzip --extract --file qp2-dependencies/resultsFile-v2.0.tar.gz mv resultsFile-*/resultsFile "\${QP_ROOT}/external/Python/" - rm --recursive --force resultsFile-* resultsFile.tar.gz EOF elif [[ ${PACKAGE} = bats ]] ; then - download ${BATS_URL} "${QP_ROOT}"/external/bats.tar.gz execute << EOF cd "\${QP_ROOT}"/external - tar -zxf bats.tar.gz + tar -zxf qp2-dependencies/bats-v1.1.0.tar.gz ( cd bats-core-1.1.0/ ; ./install.sh \${QP_ROOT}) - rm --recursive --force -- bats-core-1.1.0 \ "\${QP_ROOT}"/external/bats.tar.gz EOF else @@ -417,12 +307,6 @@ if [[ ${NINJA} = $(not_found) ]] ; then fail fi -IRPF90=$(find_exe irpf90) -if [[ ${IRPF90} = $(not_found) ]] ; then - error "IRPF90 (irpf90) is not installed." - fail -fi - ZEROMQ=$(find_lib -lzmq) if [[ ${ZEROMQ} = $(not_found) ]] ; then error "ZeroMQ (zeromq) is not installed." @@ -441,24 +325,6 @@ if [[ ${ZLIB} = $(not_found) ]] ; then fail fi -LIBCAP=$(find_lib -lcap) -if [[ ${LIBCAP} = $(not_found) ]] ; then - error "Libcap (libcap) is not installed." - fail -fi - -BWRAP=$(find_exe bwrap) -if [[ ${BWRAP} = $(not_found) ]] ; then - error "Bubblewrap (bwrap) is not installed." - fail -fi - -OPAM=$(find_exe opam) -if [[ ${OPAM} = $(not_found) ]] ; then - error "OPAM (ocaml) package manager is not installed." - fail -fi - OCAML=$(find_exe ocaml) if [[ ${OCAML} = $(not_found) ]] ; then error "OCaml (ocaml) compiler is not installed." diff --git a/etc/cflags.rc b/etc/cflags.rc new file mode 100644 index 00000000..89ab1d7f --- /dev/null +++ b/etc/cflags.rc @@ -0,0 +1 @@ +export CFLAGS="$CFLAGS --std=gnu99" diff --git a/etc/ocaml.rc b/etc/ocaml.rc index 1a1e5612..da6de03f 100644 --- a/etc/ocaml.rc +++ b/etc/ocaml.rc @@ -4,8 +4,10 @@ if [[ -z $OPAMROOT ]] then # Comment these lines if you have a system-wide OCaml installation - export OPAMROOT=${QP_ROOT}/external/opam - + export PATH="${QP_ROOT}/external/ocaml-bundle/bootstrap/bin:$PATH" + if [[ -f "${QP_ROOT}/external/ocaml-bundle/bootstrap/bin/opam" ]] ; then + eval $(opam env --root "${QP_ROOT}/external/ocaml-bundle/opam" --set-root) + fi fi source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true diff --git a/etc/openmp.rc b/etc/openmp.rc new file mode 100644 index 00000000..7f71c9b8 --- /dev/null +++ b/etc/openmp.rc @@ -0,0 +1 @@ +export OMP_NESTED=True diff --git a/etc/paths.rc b/etc/paths.rc index aff62f6e..84c2d12f 100644 --- a/etc/paths.rc +++ b/etc/paths.rc @@ -34,9 +34,9 @@ export PATH=$(qp_prepend_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROO export LD_LIBRARY_PATH=$(qp_prepend_export "LD_LIBRARY_PATH" "${QP_ROOT}"/lib) - export LIBRARY_PATH=$(qp_prepend_export "LIBRARY_PATH" "${QP_ROOT}"/lib:"${QP_ROOT}"/lib64) +export PKG_CONFIG_PATH=$(qp_prepend_export "PKG_CONFIG_PATH" "${QP_ROOT}"/lib/pkgconfig) export C_INCLUDE_PATH=$(qp_prepend_export "C_INCLUDE_PATH" "${QP_ROOT}"/include) export CPATH=$(qp_prepend_export "CPATH" "${QP_ROOT}"/include) diff --git a/external/Python/.gitignore b/external/Python/.gitignore index 1cd79fd6..e69de29b 100644 --- a/external/Python/.gitignore +++ b/external/Python/.gitignore @@ -1,46 +0,0 @@ -# Byte-compiled / optimized / DLL files -__pycache__/ -*.py[cod] - -# C extensions -*.so - -# Distribution / packaging -.Python -env/ -build/ -develop-eggs/ -dist/ -downloads/ -eggs/ - -# PyInstaller -# Usually these files are written by a python script from a template -# before PyInstaller builds the exe, so as to inject date/other infos into it. -*.manifest -*.spec - -# Installer logs -pip-log.txt -pip-delete-this-directory.txt - -# Unit test / coverage reports -htmlcov/ -.tox/ -.coverage -.cache -nosetests.xml -coverage.xml - -# Translations -*.mo -*.pot - -# Django stuff: -*.log - -# Sphinx documentation -docs/_build/ - -# PyBuilder -target/ diff --git a/external/gmp-6.1.2.tar.bz2 b/external/gmp-6.1.2.tar.bz2 deleted file mode 100644 index 3b9b275c..00000000 Binary files a/external/gmp-6.1.2.tar.bz2 and /dev/null differ diff --git a/external/qp2-dependencies b/external/qp2-dependencies new file mode 160000 index 00000000..0e1ca913 --- /dev/null +++ b/external/qp2-dependencies @@ -0,0 +1 @@ +Subproject commit 0e1ca91313e4b6ba3ea042b6378c3ff483781fb1 diff --git a/scripts/module/module_handler.py b/scripts/module/module_handler.py index 707a6734..d66918e2 100755 --- a/scripts/module/module_handler.py +++ b/scripts/module/module_handler.py @@ -31,7 +31,7 @@ try: from docopt import docopt from qp_path import QP_SRC, QP_ROOT, QP_PLUGINS, QP_EZFIO except ImportError: - print("source .quantum_package.rc") + print("source quantum_package.rc") raise diff --git a/src/ao_basis/aos.irp.f b/src/ao_basis/aos.irp.f index d2ce1ab2..1cbd3976 100644 --- a/src/ao_basis/aos.irp.f +++ b/src/ao_basis/aos.irp.f @@ -21,6 +21,21 @@ BEGIN_PROVIDER [ integer, ao_shell, (ao_num) ] enddo enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer, ao_first_of_shell, (shell_num) ] + implicit none + BEGIN_DOC + ! Index of the shell to which the AO corresponds + END_DOC + integer :: i, j, k, n + k=1 + do i=1,shell_num + ao_first_of_shell(i) = k + n = shell_ang_mom(i)+1 + k = k+(n*(n+1))/2 + enddo + END_PROVIDER BEGIN_PROVIDER [ double precision, ao_coef_normalized, (ao_num,ao_prim_num_max) ] diff --git a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f index 988bbe0a..24f43311 100644 --- a/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f @@ -37,26 +37,58 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] integer :: num_A,num_B double precision :: A_center(3),B_center(3),C_center(3) integer :: power_A(3),power_B(3) - integer :: i,j,k,l,n_pt_in,m + integer :: i,j,k,l,m double precision :: Vloc, Vpseudo - double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 + double precision :: wall_1, wall_2, wall_0 integer :: thread_num integer :: omp_get_thread_num + double precision :: c + double precision :: Z + + PROVIDE ao_coef_normalized_ordered_transp + PROVIDE pseudo_v_k_transp pseudo_n_k_transp pseudo_klocmax pseudo_dz_k_transp ao_pseudo_integrals_local = 0.d0 print*, 'Providing the nuclear electron pseudo integrals (local)' - call wall_time(wall_1) - call cpu_time(cpu_1) + ! Dummy iteration for OpenMP + j=1 + i=1 + l=1 + m=1 + num_A = ao_nucl(j) + power_A(1:3)= ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + num_B = ao_nucl(i) + power_B(1:3)= ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + alpha = ao_expo_ordered_transp(l,j) + beta = ao_expo_ordered_transp(m,i) + c = 0.d0 + do k = 1, nucl_num + Z = nucl_charge(k) + C_center(1:3) = nucl_coord(k,1:3) + + c = c + Vloc(pseudo_klocmax, & + pseudo_v_k_transp (1,k), & + pseudo_n_k_transp (1,k), & + pseudo_dz_k_transp(1,k), & + A_center,power_A,alpha,B_center,power_B,beta,C_center) + + enddo + + + ao_pseudo_integrals_local = 0.d0 + call wall_time(wall_1) thread_num = 0 !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& - !$OMP num_A,num_B,Z,c,n_pt_in, & + !$OMP num_A,num_B,Z,c, & !$OMP wall_0,wall_2,thread_num) & !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,& !$OMP ao_pseudo_integrals_local,nucl_num,nucl_charge, & @@ -66,7 +98,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] !$ thread_num = omp_get_thread_num() wall_0 = wall_1 - !$OMP DO SCHEDULE (guided) + !$OMP DO do j = 1, ao_num @@ -85,7 +117,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] do m=1,ao_prim_num(i) beta = ao_expo_ordered_transp(m,i) - double precision :: c c = 0.d0 if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& @@ -93,7 +124,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] cycle endif do k = 1, nucl_num - double precision :: Z Z = nucl_charge(k) C_center(1:3) = nucl_coord(k,1:3) @@ -137,25 +167,28 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] integer :: num_A,num_B double precision :: A_center(3),B_center(3),C_center(3) integer :: power_A(3),power_B(3) - integer :: i,j,k,l,n_pt_in,m + integer :: i,j,k,l,m double precision :: Vloc, Vpseudo integer :: omp_get_thread_num - double precision :: cpu_1, cpu_2, wall_1, wall_2, wall_0 + double precision :: wall_1, wall_2, wall_0 integer :: thread_num + double precision :: c + double precision :: Z + PROVIDE ao_coef_normalized_ordered_transp + PROVIDE pseudo_lmax pseudo_kmax pseudo_v_kl_transp pseudo_n_kl_transp pseudo_dz_kl_transp ao_pseudo_integrals_non_local = 0.d0 print*, 'Providing the nuclear electron pseudo integrals (non-local)' call wall_time(wall_1) - call cpu_time(cpu_1) thread_num = 0 !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& - !$OMP num_A,num_B,Z,c,n_pt_in, & + !$OMP num_A,num_B,Z,c, & !$OMP wall_0,wall_2,thread_num) & !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,nucl_coord,ao_coef_normalized_ordered_transp,& !$OMP ao_pseudo_integrals_non_local,nucl_num,nucl_charge,& @@ -184,7 +217,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] do m=1,ao_prim_num(i) beta = ao_expo_ordered_transp(m,i) - double precision :: c c = 0.d0 if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& @@ -193,7 +225,6 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] endif do k = 1, nucl_num - double precision :: Z Z = nucl_charge(k) C_center(1:3) = nucl_coord(k,1:3) diff --git a/src/ao_one_e_ints/pseudopot.f90 b/src/ao_one_e_ints/pseudopot.f90 index 48e3803e..141292d8 100644 --- a/src/ao_one_e_ints/pseudopot.f90 +++ b/src/ao_one_e_ints/pseudopot.f90 @@ -666,7 +666,7 @@ double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) arg=g_a*ac**2+g_b*bc**2 - if(arg.gt.-dlog(10.d-20))then + if(arg.gt.-dlog(1.d-20))then Vloc=0.d0 return endif @@ -1839,7 +1839,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg) m_1 = m+m+1 nlm = n+m+l pi=dacos(-1.d0) - a_over_b_square = (a/b)**2 + a_over_b_square = (a*a)/(b*b) ! First term of the sequence @@ -1869,21 +1869,16 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg) qk = dble(q) two_qkmp1 = 2.d0*(qk+mk)+1.d0 do k=0,q-1 - ! possible FPE here. To be checked + if (s_q_k < 1.d-32) then + s_q_k = 0.d0 + exit + endif s_q_k = two_qkmp1*qk*inverses(k)*s_q_k -! if (s_q_k < 1.d-32) then -! s_q_k = 0.d0 -! exit -! endif sum=sum+s_q_k two_qkmp1 = two_qkmp1-2.d0 qk = qk-1.d0 enddo inverses(q) = a_over_b_square/(dble(q+n+q+n+3) * dble(q+1)) -! do k=0,q -! sum=sum+s_q_k -! s_q_k = a_over_b_square * ( dble(2*(q-k+m)+1)*dble(q-k)/(dble(2*(k+n)+3) * dble(k+1)) ) * s_q_k -! enddo int=int+sum @@ -1892,7 +1887,6 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg) else !Compute the s_q+1_0 -! s_q_0=s_q_0*(2.d0*q+nlm+1)*b**2/((2.d0*(m+q)+3)*4.d0*(q+1)*gam) s_q_0=s_q_0*(q+q+nlm+1)*b*b/(dble(8*(m+q)+12)*(q+1)*gam) if(mod(n+m+l,2).eq.1)s_q_0=s_q_0*dsqrt(pi*.5d0) diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index 55b2d5e2..c3b206e1 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -327,6 +327,8 @@ double precision function get_ao_two_e_integral(i,j,k,l,map) result(result) implicit none BEGIN_DOC ! Gets one AO bi-electronic integral from the AO map + ! + ! i,j,k,l in physicist notation END_DOC integer, intent(in) :: i,j,k,l integer(key_kind) :: idx diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg index 7f2ede4c..a6df8e7a 100644 --- a/src/basis/EZFIO.cfg +++ b/src/basis/EZFIO.cfg @@ -39,7 +39,7 @@ interface: ezfio, provider [shell_prim_index] type: integer -doc: Max number of primitives in a shell +doc: Index of the first primitive of the shell size: (basis.shell_num) interface: ezfio, provider diff --git a/src/basis_correction/print_routine.irp.f b/src/basis_correction/print_routine.irp.f index 05fbbf60..67c5c6c2 100644 --- a/src/basis_correction/print_routine.irp.f +++ b/src/basis_correction/print_routine.irp.f @@ -38,7 +38,7 @@ subroutine print_basis_correction write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate) enddo - else if(mu_of_r_potential.EQ."cas_ful")then + else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then print*, '' print*,'Using a CAS-like two-body density to define mu(r)' print*,'This assumes that the CAS is a qualitative representation of the wave function ' diff --git a/src/cipsi/NEED b/src/cipsi/NEED index 6c14b4b6..bfbc559a 100644 --- a/src/cipsi/NEED +++ b/src/cipsi/NEED @@ -1,7 +1,6 @@ perturbation zmq mpi -davidson_undressed iterations two_body_rdm csf diff --git a/src/dav_general_mat/NEED b/src/dav_general_mat/NEED new file mode 100644 index 00000000..711fbf96 --- /dev/null +++ b/src/dav_general_mat/NEED @@ -0,0 +1 @@ +davidson_undressed diff --git a/src/dav_general_mat/README.rst b/src/dav_general_mat/README.rst new file mode 100644 index 00000000..e528109f --- /dev/null +++ b/src/dav_general_mat/README.rst @@ -0,0 +1,13 @@ +=============== +dav_general_mat +=============== + +This modules allows to use the Davidson Algorithm for general squared symmetric matrices +You have two options : + a) the routine "davidson_general" to whom you pass the matrix you want to diagonalize + b) the routine "davidson_general_ext_rout" to whom you pass the subroutine that realizes v = H u + +See the routines in "test_dav.irp.f" for a clear example. + + + diff --git a/src/dav_general_mat/dav_ext_rout.irp.f b/src/dav_general_mat/dav_ext_rout.irp.f new file mode 100644 index 00000000..a2aad413 --- /dev/null +++ b/src/dav_general_mat/dav_ext_rout.irp.f @@ -0,0 +1,447 @@ + +subroutine davidson_general_ext_rout(u_in,H_jj,energies,dim_in,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 + ! + ! 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 : 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) :: dim_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(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 + + !--------------- + + 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) + + +! 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( & + ! 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) + 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(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 + END_DOC + integer, intent(in) :: dim_in, sze, N_st, N_st_diag_in + double precision, intent(in) :: H_jj(sze),h_mat(sze,sze) + double precision, intent(inout) :: u_in(dim_in,N_st_diag_in) + double precision, intent(out) :: energies(N_st) + + 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 :: 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 + 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 + + !--------------- + + 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) + + +! 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( & + ! 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-Smitt 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) + +! call H_S2_u_0_nstates_openmp(W(1,shift+1),U(1,shift+1),N_st_diag,sze) + call hpsi(W(1,shift+1),U(1,shift+1),N_st_diag,sze,h_mat) + 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(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 + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,4) = core_inact_act_bitmask_4(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I V V !!!!!!!!!!!!!!!!!!!! + ! (core+inact+act) ^ 2 (virt) ^2 + ! = J_iv + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = virt_bitmask(i,1) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + + ! (core+inact+act) ^ 2 (virt) ^2 + ! = (iv|iv) + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,3) = virt_bitmask(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! V V V !!!!!!!!!!!!!!!!!!!!!!! +! if(.not.no_vvv_integrals)then + print*, '' + print*, ' and ' + do i = 1,N_int + mask_ijk(i,1) = virt_bitmask(i,1) + mask_ijk(i,2) = virt_bitmask(i,1) + mask_ijk(i,3) = virt_bitmask(i,1) + enddo + call add_integrals_to_map_three_indices(mask_ijk) +! endif + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I V !!!!!!!!!!!!!!!!!!!! + ! (core+inact+act) ^ 3 (virt) ^1 + ! + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map(mask_ijkl) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I V V V !!!!!!!!!!!!!!!!!!!! + ! (core+inact+act) ^ 1 (virt) ^3 + ! +! if(.not.no_ivvv_integrals)then + print*, '' + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = virt_bitmask(i,1) + mask_ijkl(i,3) = virt_bitmask(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map_no_exit_34(mask_ijkl) +end diff --git a/src/mu_of_r/EZFIO.cfg b/src/mu_of_r/EZFIO.cfg index 5677b3ab..c774ec82 100644 --- a/src/mu_of_r/EZFIO.cfg +++ b/src/mu_of_r/EZFIO.cfg @@ -6,7 +6,7 @@ size: (becke_numerical_grid.n_points_final_grid,determinants.n_states) [mu_of_r_potential] type: character*(32) -doc: type of potential for the mu(r) interaction: can be [ hf| cas_ful | cas_truncated] +doc: type of potential for the mu(r) interaction: can be [ hf| cas_ful | cas_truncated | pure_act] interface: ezfio, provider, ocaml default: hf diff --git a/src/mu_of_r/basis_def.irp.f b/src/mu_of_r/basis_def.irp.f index 4da27cb0..fff9f581 100644 --- a/src/mu_of_r/basis_def.irp.f +++ b/src/mu_of_r/basis_def.irp.f @@ -76,7 +76,11 @@ BEGIN_PROVIDER [integer, n_basis_orb] ! ! It corresponds to all MOs except those defined as "deleted" END_DOC - n_basis_orb = n_all_but_del_orb + if(mu_of_r_potential == "pure_act")then + n_basis_orb = n_act_orb + else + n_basis_orb = n_all_but_del_orb + endif END_PROVIDER BEGIN_PROVIDER [integer, list_basis, (n_basis_orb)] @@ -89,9 +93,15 @@ BEGIN_PROVIDER [integer, list_basis, (n_basis_orb)] ! It corresponds to all MOs except those defined as "deleted" END_DOC integer :: i - do i = 1, n_all_but_del_orb - list_basis(i) = list_all_but_del_orb(i) - enddo + if(mu_of_r_potential == "pure_act")then + do i = 1, n_act_orb + list_basis(i) = list_act(i) + enddo + else + do i = 1, n_all_but_del_orb + list_basis(i) = list_all_but_del_orb(i) + enddo + endif END_PROVIDER BEGIN_PROVIDER [double precision, basis_mos_in_r_array, (n_basis_orb,n_points_final_grid)] diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f index 148c65b3..5c41acdc 100644 --- a/src/mu_of_r/mu_of_r_conditions.irp.f +++ b/src/mu_of_r/mu_of_r_conditions.irp.f @@ -26,7 +26,7 @@ do ipoint = 1, n_points_final_grid if(mu_of_r_potential.EQ."hf")then mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint) - else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then + else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate) else print*,'you requested the following mu_of_r_potential' diff --git a/src/tools/fcidump_pyscf.irp.f b/src/tools/fcidump_pyscf.irp.f new file mode 100644 index 00000000..9cbf733a --- /dev/null +++ b/src/tools/fcidump_pyscf.irp.f @@ -0,0 +1,78 @@ +program fcidump + implicit none + BEGIN_DOC +! Produce a regular `FCIDUMP` file from the |MOs| stored in the |EZFIO| +! directory. +! +! To specify an active space, the class of the |MOs| have to set in the +! |EZFIO| directory (see :ref:`qp_set_mo_class`). +! +! The :ref:`fcidump` program supports 3 types of |MO| classes : +! +! * the *core* orbitals which are always doubly occupied in the +! calculation +! +! * the *deleted* orbitals that are never occupied in the calculation +! +! * the *active* orbitals that are occupied with a varying number of +! electrons +! + END_DOC + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + output=trim(ezfio_filename)//'.FCIDUMP' + i_unit_output = getUnitAndOpen(output,'w') + + integer :: i,j,k,l + integer :: i1,j1,k1,l1 + integer :: i2,j2,k2,l2 + integer*8 :: m + character*(2), allocatable :: A(:) + + write(i_unit_output,*) '&FCI NORB=', n_act_orb, ', NELEC=', elec_num-n_core_orb*2, & + ', MS2=', (elec_alpha_num-elec_beta_num), ',' + allocate (A(n_act_orb)) + A = '1,' + write(i_unit_output,*) 'ORBSYM=', (A(i), i=1,n_act_orb) + write(i_unit_output,*) 'ISYM=0,' + write(i_unit_output,*) '&end' + deallocate(A) + + integer(key_kind), allocatable :: keys(:) + double precision, allocatable :: values(:) + integer(cache_map_size_kind) :: n_elements, n_elements_max + PROVIDE mo_two_e_integrals_in_map + + double precision :: get_two_e_integral, integral + + do l=1,n_act_orb + l1 = list_act(l) + do k=1,n_act_orb + k1 = list_act(k) + do j=l,n_act_orb + j1 = list_act(j) + do i=k,n_act_orb + i1 = list_act(i) + if (i1>=j1) then + integral = get_two_e_integral(i1,j1,k1,l1,mo_integrals_map) + if (dabs(integral) > mo_integrals_threshold) then + write(i_unit_output,*) integral, i,k,j,l + endif + end if + enddo + enddo + enddo + enddo + + do j=1,n_act_orb + j1 = list_act(j) + do i=j,n_act_orb + i1 = list_act(i) + integral = mo_one_e_integrals(i1,j1) + core_fock_operator(i1,j1) + if (dabs(integral) > mo_integrals_threshold) then + write(i_unit_output,*) integral, i,j,0,0 + endif + enddo + enddo + write(i_unit_output,*) core_energy, 0, 0, 0, 0 +end diff --git a/src/tools/save_natorb_no_ov_rot.irp.f b/src/tools/save_natorb_no_ov_rot.irp.f new file mode 100644 index 00000000..e5b69fbf --- /dev/null +++ b/src/tools/save_natorb_no_ov_rot.irp.f @@ -0,0 +1,25 @@ +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_no_ov_rot + call save_ref_determinant + 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/utils/integration.irp.f b/src/utils/integration.irp.f index 5fe7438b..38e198dc 100644 --- a/src/utils/integration.irp.f +++ b/src/utils/integration.irp.f @@ -55,6 +55,10 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha, ! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 ) ! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 ) ! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 ) + ! + ! WARNING ::: IF fact_k is too smal then: + ! returns a "s" function centered in zero + ! with an inifinite exponent and a zero polynom coef END_DOC implicit none include 'constants.include.F' @@ -82,10 +86,13 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha, !DIR$ FORCEINLINE call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center) if (fact_k < thresh) then + ! IF fact_k is too smal then: + ! returns a "s" function centered in zero + ! with an inifinite exponent and a zero polynom coef P_center = 0.d0 - p = 1.d-10 + p = 1.d+15 P_new = 0.d0 - iorder = -1 + iorder = 0 fact_k = 0.d0 return endif diff --git a/travis/configuration.sh b/travis/configuration.sh index 7b3f5423..f925107d 100755 --- a/travis/configuration.sh +++ b/travis/configuration.sh @@ -2,7 +2,7 @@ # Stage 1 # Configure QP2 -./configure --install all --config ./config/travis.cfg || exit -1 +./configure --download all --install all --config ./config/travis.cfg || exit -1 # Create cache cd ../