LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 11:03:29 +01:00
Code Releases Activity

Merge branch 'dev' of github.com:QuantumPackage/qp2 into dev

Browse Source
This commit is contained in:
Anthony Scemama 2021-07-28 17:51:40 +02:00
commit d90eb54b24
30 changed files with 1326 additions and 262 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitmodules vendored
View File

@ -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

4
AUTHORS
View File

@ -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

34
INSTALL.rst
View File

@ -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 <https://gmplib.org/>`_ : Gnu Multiple Precision Arithmetic Library
- |OCaml| compiler with |OPAM| package manager
- `Bubblewrap <https://github.com/projectatomic/bubblewrap>`_ : Sandboxing tool required by Opam
- `libcap <https://git.kernel.org/pub/scm/linux/kernel/git/morgan/libcap.git>`_ : 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

5
RELEASE_NOTES.org
View File

@ -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

63
bin/qp_convert_output_to_ezfio
View File

@ -120,6 +120,7 @@ def write_ezfio(res, filename):
exponent = []
res.convert_to_cartesian()
# ~#~#~#~#~#~#~ #
# P a r s i n g #
# ~#~#~#~#~#~#~ #
@ -177,6 +178,68 @@ def write_ezfio(res, filename):
print("OK")
# _
# |_) _. _ o _
# |_) (_| _> | _>
#
print("Basis\t\t...\t", end=' ')
# ~#~#~#~ #
# I n i t #
# ~#~#~#~ #
coefficient = []
exponent = []
# ~#~#~#~#~#~#~ #
# P a r s i n g #
# ~#~#~#~#~#~#~ #
nbasis = 0
nucl_center = []
curr_center = -1
nucl_shell_num = []
ang_mom = []
nshell = 0
shell_prim_index = [1]
shell_prim_num = []
for b in res.basis:
s = b.sym
if str.count(s, "y") + str.count(s, "x") == 0:
c = b.center
nshell += 1
if c != curr_center:
curr_center = c
nucl_center.append(nbasis+1)
nucl_shell_num.append(nshell)
nshell = 0
nbasis += 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:]
# ~#~#~#~#~ #
# 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_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)
print("OK")
# _
# |\/| _ _ |_) _. _ o _
# | | (_) _> |_) (_| _> | _>

184
configure vendored
View File

@ -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 <<EOF
@ -82,7 +68,6 @@ function execute () {
}
PACKAGES=""
OCAML_PACKAGES="ocamlbuild cryptokit zmq sexplib ppx_sexp_conv ppx_deriving getopt"
while true ; do
case "$1" in
@ -134,16 +119,6 @@ function success() {
exit 0
}
function download() {
echo "Downloading $1"
echo ""
printf "\e[0;34m"
wget --no-check-certificate $1 --output-document=$2 || error "Unable to download $1"
printf "\e[m"
echo "Saved dowloaded file as $2"
echo ""
}
function not_found() {
echo 'not_found'
}
@ -176,6 +151,10 @@ function find_dir() {
fi
}
# Make program believe stdin is a tty
function faketty() {
script -qfc "$(printf "%q " "$@")" /dev/null
}
# Install IRPF90 if needed
IRPF90=$(find_exe irpf90)
@ -205,7 +184,7 @@ if [[ "${PACKAGES}.x" != ".x" ]] ; then
fi
if [[ ${PACKAGES} = all ]] ; then
PACKAGES="zlib ninja irpf90 zeromq f77zmq gmp libcap bwrap ocaml docopt resultsFile bats"
PACKAGES="zlib ninja zeromq f77zmq gmp ocaml docopt resultsFile bats"
fi
@ -213,10 +192,9 @@ for PACKAGE in ${PACKAGES} ; do
if [[ ${PACKAGE} = ninja ]] ; then
download ${NINJA_URL} "${QP_ROOT}"/external/ninja.zip
execute << EOF
rm -f "\${QP_ROOT}"/bin/ninja
unzip "\${QP_ROOT}"/external/ninja.zip -d "\${QP_ROOT}"/bin
unzip "\${QP_ROOT}"/external/qp2-dependencies/ninja-linux.zip -d "\${QP_ROOT}"/bin
EOF
@ -224,70 +202,31 @@ EOF
execute << EOF
cd "\${QP_ROOT}"/external
tar --bzip2 --extract --file gmp-6.1.2.tar.bz2
tar --bzip2 --extract --file qp2-dependencies/gmp-6.1.2.tar.bz2
cd gmp-6.1.2
./configure --prefix=$QP_ROOT && make -j 8
make install
make -j 8 install
EOF
elif [[ ${PACKAGE} = libcap ]] ; then
download ${LIBCAP_URL} "${QP_ROOT}"/external/libcap.tar.gz
execute << EOF
cd "\${QP_ROOT}"/external
tar --gunzip --extract --file libcap.tar.gz
rm libcap.tar.gz
cd libcap-*/libcap
prefix=$QP_ROOT make BUILD_GPERF=no install
EOF
elif [[ ${PACKAGE} = bwrap ]] ; then
download ${BUBBLE_URL} "${QP_ROOT}"/external/bwrap.tar.xz
execute << EOF
cd "\${QP_ROOT}"/external
tar --xz --extract --file bwrap.tar.xz
rm bwrap.tar.xz
cd bubblewrap*
./configure --prefix=$QP_ROOT && make -j 8
make install-exec-am
EOF
elif [[ ${PACKAGE} = irpf90 ]] ; then
execute << EOF
cd "\${QP_ROOT}"/external
tar --gunzip --extract --file irpf90.tar.gz
rm irpf90.tar.gz
mv irpf90-* irpf90
cd irpf90
make
EOF
elif [[ ${PACKAGE} = zeromq ]] ; then
download ${ZEROMQ_URL} "${QP_ROOT}"/external/zeromq.tar.gz
execute << EOF
export CC=gcc
export CXX=g++
cd "\${QP_ROOT}"/external
tar --gunzip --extract --file zeromq.tar.gz
rm zeromq.tar.gz
tar --gunzip --extract --file qp2-dependencies/zeromq-4.2.5.tar.gz
cd zeromq-*
./configure --prefix="\$QP_ROOT" --without-libsodium --enable-libunwind=no
make
make -j 8
make install
EOF
elif [[ ${PACKAGE} = f77zmq ]] ; then
download ${F77ZMQ_URL} "${QP_ROOT}"/external/f77_zmq.tar.gz
execute << EOF
cd "\${QP_ROOT}"/external
tar --gunzip --extract --file f77_zmq.tar.gz
rm f77_zmq.tar.gz
tar --gunzip --extract --file qp2-dependencies/f77_zmq-4.2.5.tar.gz
cd f77_zmq-*
export ZMQ_H="\$QP_ROOT"/include/zmq.h
make
@ -299,71 +238,28 @@ EOF
elif [[ ${PACKAGE} = ocaml ]] ; then
download ${OCAML_URL} "${QP_ROOT}"/external/opam_installer.sh
if [[ -n ${TRAVIS} ]] ; then
# Special commands for Travis CI
chmod +x "${QP_ROOT}"/external/opam_installer.sh
rm --force ${QP_ROOT}/bin/opam
if [[ -n ${NO_CACHE} ]] ; then
rm -rf ${HOME}/.opam
fi
export OPAMROOT=${HOME}/.opam
cat << EOF | bash ${QP_ROOT}/external/opam_installer.sh --no-backup
${QP_ROOT}/bin
execute <<EOF
source "${QP_ROOT}"/quantum_package.rc
cd "${QP_ROOT}"/external/
tar --gunzip --extract --file qp2-dependencies/ocaml-bundle_x86.tar.gz
echo "" | ./ocaml-bundle/bootstrap.sh "${QP_ROOT}"
./ocaml-bundle/configure.sh "${QP_ROOT}"
echo "" | ./ocaml-bundle/compile.sh "${QP_ROOT}"
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
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."

1
etc/cflags.rc Normal file
View File

@ -0,0 +1 @@
export CFLAGS="$CFLAGS --std=gnu99"

6
etc/ocaml.rc
View File

@ -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

2
etc/paths.rc
View File

@ -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)

46
external/Python/.gitignore vendored
View File

@ -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/

BIN
external/gmp-6.1.2.tar.bz2 vendored
View File

Binary file not shown.

1
external/qp2-dependencies vendored Submodule

@ -0,0 +1 @@
Subproject commit ae9397c3b4e689a487fdd4b1425af5f519d3ea82

2
scripts/module/module_handler.py
View File

@ -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

59
src/ao_one_e_ints/pot_ao_pseudo_ints.irp.f
View File

@ -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)

18
src/ao_one_e_ints/pseudopot.f90
View File

@ -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)

2
src/cipsi/pt2_stoch_routines.irp.f
View File

@ -132,6 +132,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
if (h0_type == 'CFG') then
PROVIDE psi_configuration_hii det_to_configuration

3
src/cipsi/zmq_selection.irp.f
View File

@ -21,7 +21,8 @@ subroutine ZMQ_selection(N_in, pt2_data)
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
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
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection')

1
src/dav_general_mat/NEED Normal file
View File

@ -0,0 +1 @@
davidson_undressed

13
src/dav_general_mat/README.rst Normal file
View File

@ -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.

447
src/dav_general_mat/dav_ext_rout.irp.f Normal file
View File

@ -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><i|H|u_k>
! -----------------------------------
! 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 = <u_k | W_l> = <u_k| H |u_l>
! -------------------------------------------
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).and.(U(k,j) == 0.d0))
k = k+1
enddo
if (U(k,j) * u_in(k,j) < 0.d0) then
do i=1,sze
W(i,j) = -W(i,j)
enddo
endif
enddo
enddo
do k=1,N_st
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, h, &
y, &
lambda &
)
deallocate(overlap)
FREE nthreads_davidson
end
subroutine hcalc_template(v,u,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! Template of routine for the application of H
!
! Here, it is done with the Hamiltonian matrix
!
! on the set of determinants of psi_det
!
! Computes $v = H | u \rangle$
!
END_DOC
integer, intent(in) :: N_st,sze
double precision, intent(in) :: u(sze,N_st)
double precision, intent(inout) :: v(sze,N_st)
integer :: i,j,istate
v = 0.d0
do istate = 1, N_st
do i = 1, sze
do j = 1, sze
v(i,istate) += H_matrix_all_dets(j,i) * u(j,istate)
enddo
enddo
do i = 1, sze
v(i,istate) += u(i,istate) * nuclear_repulsion
enddo
enddo
end

435
src/dav_general_mat/dav_general.irp.f Normal file
View File

@ -0,0 +1,435 @@
subroutine davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,h_mat)
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
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><i|H|u_k>
! -----------------------------------
! 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 = <u_k | W_l> = <u_k| H |u_l>
! -------------------------------------------
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).and.(U(k,j) == 0.d0))
k = k+1
enddo
if (U(k,j) * u_in(k,j) < 0.d0) then
do i=1,sze
W(i,j) = -W(i,j)
enddo
endif
enddo
enddo
do k=1,N_st
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, h, &
y, &
lambda &
)
deallocate(overlap)
FREE nthreads_davidson
end
subroutine hpsi(v,u,N_st,sze,h_mat)
use bitmasks
implicit none
BEGIN_DOC
! Computes $v = H | u \rangle$ and
END_DOC
integer, intent(in) :: N_st,sze
double precision, intent(in) :: u(sze,N_st),h_mat(sze,sze)
double precision, intent(inout) :: v(sze,N_st)
integer :: i,j,istate
v = 0.d0
do istate = 1, N_st
do i = 1, sze
do j = 1, sze
v(i,istate) += h_mat(j,i) * u(j,istate)
enddo
enddo
enddo
end

48
src/dav_general_mat/test_dav.irp.f Normal file
View File

@ -0,0 +1,48 @@
program test_dav
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
print *, 'Hello world'
read_wf = .True.
touch read_wf
PROVIDE threshold_davidson nthreads_davidson
call routine
end
subroutine routine
implicit none
double precision, allocatable :: u_in(:,:), H_jj(:), energies(:),h_mat(:,:)
integer :: dim_in,sze,N_st,N_st_diag_in,dressing_state
logical :: converged
integer :: i,j
external hcalc_template
N_st = N_states
N_st_diag_in = N_states_diag
sze = N_det
dim_in = sze
dressing_state = 0
!!!! MARK THAT u_in mut dimensioned with "N_st_diag_in" as a second dimension
allocate(u_in(dim_in,N_st_diag_in),H_jj(sze),h_mat(sze,sze),energies(N_st))
u_in = 0.d0
do i = 1, N_st
u_in(1,i) = 1.d0
enddo
!!! Matrix "h_mat" is the matrix we want to diagonalize with the first routine
!!! "davidson_general"
do i = 1, sze
do j = 1, sze
h_mat(j,i) = H_matrix_all_dets(j,i)
enddo
H_jj(i) = H_mat(i,i) + nuclear_repulsion
h_mat(i,i) = H_mat(i,i) + nuclear_repulsion
enddo
provide nthreads_davidson
call davidson_general(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,h_mat)
print*,'energies = ',energies
!!! hcalc_template is the routine that computes v = H u
!!! and you can use the routine "davidson_general_ext_rout"
call davidson_general_ext_rout(u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag_in,converged,hcalc_template)
print*,'energies = ',energies
end

2
src/davidson/diagonalize_ci.irp.f
View File

@ -42,7 +42,7 @@ END_PROVIDER
logical :: converged
logical :: do_csf
PROVIDE threshold_davidson nthreads_davidson
PROVIDE threshold_davidson nthreads_davidson distributed_davidson
! Guess values for the "N_states" states of the |CI| eigenvectors
do j=1,min(N_states,N_det)
do i=1,N_det

20
src/dft_utils_in_r/mo_in_r.irp.f
View File

@ -161,3 +161,23 @@
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, mos_grad_in_r_array_transp_3, (3,n_points_final_grid,mo_num)]
implicit none
BEGIN_DOC
! Transposed gradients
!
END_DOC
integer :: i,j,m
double precision :: mos_array(mo_num), r(3)
double precision :: mos_grad_array(3,mo_num)
do m = 1, 3
do j = 1, mo_num
do i = 1, n_points_final_grid
mos_grad_in_r_array_transp_3(m,i,j) = mos_grad_in_r_array(j,i,m)
enddo
enddo
enddo
END_PROVIDER

7
src/mo_two_e_ints/four_idx_novvvv.irp.f
View File

@ -51,6 +51,13 @@ end
subroutine four_idx_novvvv
print*,'********'
print*,'********'
print*,'********'
print*,'WARNING :: Using four_idx_novvvv, and we are not sure that this routine is not bugged ...'
print*,'********'
print*,'********'
print*,'********'
use map_module
implicit none
BEGIN_DOC

3
src/mo_two_e_ints/mo_bi_integrals.irp.f
View File

@ -50,7 +50,8 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
call cpu_time(cpu_1)
if(no_vvvv_integrals)then
call four_idx_novvvv
! call four_idx_novvvv
call four_idx_novvvv_old
else
call add_integrals_to_map(full_ijkl_bitmask_4)
endif

88
src/mo_two_e_ints/no_vvvv.irp.f Normal file
View File

@ -0,0 +1,88 @@
subroutine four_idx_novvvv_old
use map_module
use bitmasks
implicit none
BEGIN_DOC
! Retransform MO integrals for next CAS-SCF step
END_DOC
integer(bit_kind) :: mask_ijkl(N_int,4)
integer(bit_kind) :: mask_ijk(N_int,3)
print*,'Using partial transformation'
print*,'It will not transform all integrals with at least 3 indices within the virtuals'
integer :: i,j,k,l
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! I I I I !!!!!!!!!!!!!!!!!!!!
! (core+inact+act) ^ 4
! <ii|ii>
print*, ''
print*, '<ii|ii>'
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
! <iv|iv> = J_iv
print*, ''
print*, '<iv|iv>'
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
! <ii|vv> = (iv|iv)
print*, ''
print*, '<ii|vv>'
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*, '<rv|sv> and <rv|vs>'
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
! <iv|ii>
print*, ''
print*, '<iv|ii>'
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
! <iv|vv>
! if(.not.no_ivvv_integrals)then
print*, ''
print*, '<iv|vv>'
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

78
src/tools/fcidump_pyscf.irp.f Normal file
View File

@ -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

11
src/utils/integration.irp.f
View File

@ -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

2
travis/configuration.sh
View File

@ -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 ../