9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 03:23:29 +01:00

Merge pull request #178 from QuantumPackage/dev

Dev
This commit is contained in:
Anthony Scemama 2021-11-18 09:10:13 +01:00 committed by GitHub
commit 924dd3a65b
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
94 changed files with 2931 additions and 533 deletions

3
.gitmodules vendored
View File

@ -4,3 +4,6 @@
[submodule "external/irpf90"] [submodule "external/irpf90"]
path = external/irpf90 path = external/irpf90
url = https://gitlab.com/scemama/irpf90.git url = https://gitlab.com/scemama/irpf90.git
[submodule "external/qp2-dependencies"]
path = external/qp2-dependencies
url = https://github.com/QuantumPackage/qp2-dependencies.git

View File

@ -4,12 +4,14 @@
- Thomas Applencourt - Thomas Applencourt
- Anouar Benali - Anouar Benali
- Michel Caffarel - Michel Caffarel
- Vijay Gopal Chilkuri
- Yann Damour
- Grégoire David - Grégoire David
- Anthony Ferté - Anthony Ferté
- Yann Garniron - Yann Garniron
- Kevin Gasperich - Kevin Gasperich
- Vijay Gopal Chilkuri
- Emmanuel Giner - Emmanuel Giner
- Fabris Kossoski
- Pierre-François Loos - Pierre-François Loos
- Jean-Paul Malrieu - Jean-Paul Malrieu
- Julien Paquier - Julien Paquier

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 This script will create the :file:`quantum_package.rc` bash script, which
sets all the environment variables required for the normal operation of the 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 When all dependencies have been installed, (the :command:`configure` will
and need to be installed. inform 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 tell 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 Now all the requirements are met, you can compile the programs using
@ -51,8 +53,6 @@ Requirements
- |ZeroMQ| : networking library - |ZeroMQ| : networking library
- `GMP <https://gmplib.org/>`_ : Gnu Multiple Precision Arithmetic Library - `GMP <https://gmplib.org/>`_ : Gnu Multiple Precision Arithmetic Library
- |OCaml| compiler with |OPAM| package manager - |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 - |Ninja| : a parallel build system
- |pkg-config| : a tool which returns information about installed libraries - |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 * zeromq
* f77zmq * f77zmq
* gmp * gmp
* libcap * ocaml (:math:`\approx` 5 minutes)
* bwrap
* ocaml ( :math:`\approx` 10 minutes)
* ezfio * ezfio
* docopt * docopt
* resultsFile * resultsFile
@ -111,8 +109,9 @@ Example:
.. note:: .. note::
When installing the ocaml package, you will be asked the location of where it should be installed. When installing the ocaml package, you will be asked the location of where
A safe option is to enter the path proposed by the |QP|: 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
@ -122,8 +121,9 @@ Example:
If the :command:`configure` executable fails to install a specific dependency If the :command:`configure` executable fails to install a specific dependency
----------------------------------------------------------------------------- -----------------------------------------------------------------------------
If the :command:`configure` executable does not succeed to install a specific dependency, If the :command:`configure` executable does not succeed to install a specific
there are some proposition of how to download and install the minimal dependencies to compile and use the |QP|. 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 Before doing anything below, try to install the packages with your package manager

View File

@ -30,6 +30,7 @@
- Fixed bug in DIIS - Fixed bug in DIIS
- Fixed bug in molden (Au -> Angs) - Fixed bug in molden (Au -> Angs)
- Fixed bug with non-contiguous MOs in active space and deleter MOs - Fixed bug with non-contiguous MOs in active space and deleter MOs
- Complete network-free installation
*** User interface *** User interface
@ -83,9 +84,7 @@
- Added LIB file to add extra libs in plugin - Added LIB file to add extra libs in plugin
- Using Intel IPP for sorting when using Intel compiler - Using Intel IPP for sorting when using Intel compiler
- Removed parallelism in sorting - Removed parallelism in sorting
- Compute banned_excitations from exchange integrals to accelerate with local MOs
ao_one_e_integral_zero
banned_excitations

View File

@ -120,6 +120,7 @@ def write_ezfio(res, filename):
exponent = [] exponent = []
res.convert_to_cartesian() res.convert_to_cartesian()
# ~#~#~#~#~#~#~ # # ~#~#~#~#~#~#~ #
# P a r s i n g # # P a r s i n g #
# ~#~#~#~#~#~#~ # # ~#~#~#~#~#~#~ #
@ -177,6 +178,68 @@ def write_ezfio(res, filename):
print("OK") 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 _ # |\/| _ _ |_) _. _ o _
# | | (_) _> |_) (_| _> | _> # | | (_) _> |_) (_| _> | _>
@ -226,10 +289,15 @@ def write_ezfio(res, filename):
for i in range(mo_num): for i in range(mo_num):
energies.append(MOs[i].eigenvalue) energies.append(MOs[i].eigenvalue)
if res.occ_num is not None:
OccNum = [] OccNum = []
if res.occ_num is not None:
for i in MOindices: for i in MOindices:
OccNum.append(res.occ_num[MO_type][i]) OccNum.append(res.occ_num[MO_type][i])
else:
for i in range(res.num_beta):
OccNum.append(2.)
for i in range(res.num_beta,res.num_alpha):
OccNum.append(1.)
while len(OccNum) < mo_num: while len(OccNum) < mo_num:
OccNum.append(0.) OccNum.append(0.)
@ -254,8 +322,9 @@ def write_ezfio(res, filename):
# ~#~#~#~#~ # # ~#~#~#~#~ #
ezfio.set_mo_basis_mo_num(mo_num) ezfio.set_mo_basis_mo_num(mo_num)
ezfio.set_mo_basis_mo_occ(OccNum)
ezfio.set_mo_basis_mo_coef(MoMatrix) ezfio.set_mo_basis_mo_coef(MoMatrix)
ezfio.set_mo_basis_mo_occ(OccNum)
print("OK") print("OK")

View File

@ -51,7 +51,8 @@ FCFLAGS : -Ofast
# -g : Extra debugging information # -g : Extra debugging information
# #
[DEBUG] [DEBUG]
FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan #FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow -finit-real=nan
# OpenMP flags # OpenMP flags
################# #################

237
configure vendored
View File

@ -3,8 +3,6 @@
# Quantum Package configuration script # Quantum Package configuration script
# #
TEMP=$(getopt -o c:i:h -l config:,install:,help -n $0 -- "$@") || exit 1
eval set -- "$TEMP"
export QP_ROOT="$( cd "$(dirname "$0")" ; pwd -P )" export QP_ROOT="$( cd "$(dirname "$0")" ; pwd -P )"
echo "QP_ROOT="$QP_ROOT echo "QP_ROOT="$QP_ROOT
@ -18,35 +16,21 @@ export CC=gcc
git submodule init git submodule init
git submodule update 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() function help()
{ {
cat <<EOF cat <<EOF
Quantum Package configuration script. Quantum Package configuration script.
Usage: Usage:
$(basename $0) -c <file> | --config=<file> $(basename $0) -c <file>
$(basename $0) -h | --help $(basename $0) -h
$(basename $0) -i <package> | --install=<package> $(basename $0) -i <package>
Options: Options:
-c, --config=<file> Define a COMPILATION configuration file, -c <file> Define a COMPILATION configuration file,
in "${QP_ROOT}/config/". in "${QP_ROOT}/config/".
-h, --help Print the HELP message -h Print the HELP message
-i, --install=<package> INSTALL <package>. Use at your OWN RISK: -i <package> INSTALL <package>. Use at your OWN RISK:
no support will be provided for the installation of no support will be provided for the installation of
dependencies. dependencies.
@ -82,33 +66,31 @@ function execute () {
} }
PACKAGES="" PACKAGES=""
OCAML_PACKAGES="ocamlbuild cryptokit zmq sexplib ppx_sexp_conv ppx_deriving getopt" echo $@
while true ; do
case "$1" in while getopts "d:c:i:h" c ; do
-c|--config) case "$c" in
case "$2" in c)
case "$OPTARG" in
"") help ; break;; "") help ; break;;
*) if [[ -f $2 ]] ; then *) if [[ -f $OPTARG ]] ; then
CONFIG="$2" CONFIG="$OPTARG"
else else
error "error: configuration file $2 not found." error "error: configuration file $OPTARG not found."
exit 1 exit 1
fi fi
esac esac;;
shift 2;; i)
-i|--install) case "$OPTARG" in
case "$2" in
"") help ; break;; "") help ; break;;
*) PACKAGES="${PACKAGE} $2" *) PACKAGES="${PACKAGE} $OPTARG"
esac esac;;
shift 2;; h)
-h|-help|--help)
help help
exit 0;; exit 0;;
--) shift ; break ;;
*) *)
error $(basename $0)": unknown option $1, try --help" error $(basename $0)": unknown option $c, try --help"
exit 2;; exit 2;;
esac esac
done done
@ -134,16 +116,6 @@ function success() {
exit 0 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() { function not_found() {
echo 'not_found' echo 'not_found'
} }
@ -176,6 +148,10 @@ function find_dir() {
fi fi
} }
# Make program believe stdin is a tty
function faketty() {
script -qfc "$(printf "%q " "$@")" /dev/null
}
# Install IRPF90 if needed # Install IRPF90 if needed
IRPF90=$(find_exe irpf90) IRPF90=$(find_exe irpf90)
@ -205,7 +181,7 @@ if [[ "${PACKAGES}.x" != ".x" ]] ; then
fi fi
if [[ ${PACKAGES} = all ]] ; then 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 fi
@ -213,10 +189,9 @@ for PACKAGE in ${PACKAGES} ; do
if [[ ${PACKAGE} = ninja ]] ; then if [[ ${PACKAGE} = ninja ]] ; then
download ${NINJA_URL} "${QP_ROOT}"/external/ninja.zip
execute << EOF execute << EOF
rm -f "\${QP_ROOT}"/bin/ninja 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 EOF
@ -224,146 +199,62 @@ EOF
execute << EOF execute << EOF
cd "\${QP_ROOT}"/external 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 cd gmp-6.1.2
./configure --prefix=$QP_ROOT && make -j 8 ./configure --prefix=$QP_ROOT && make -j 8
make install make -j 8 install
EOF 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 elif [[ ${PACKAGE} = zeromq ]] ; then
download ${ZEROMQ_URL} "${QP_ROOT}"/external/zeromq.tar.gz
execute << EOF execute << EOF
export CC=gcc export CC=gcc
export CXX=g++ export CXX=g++
cd "\${QP_ROOT}"/external cd "\${QP_ROOT}"/external
tar --gunzip --extract --file zeromq.tar.gz tar --gunzip --extract --file qp2-dependencies/zeromq-4.2.5.tar.gz
rm zeromq.tar.gz
cd zeromq-* cd zeromq-*
./configure --prefix="\$QP_ROOT" --without-libsodium --enable-libunwind=no ./configure --prefix="\$QP_ROOT" --without-libsodium --enable-libunwind=no
make make -j 8
make install make install
EOF EOF
elif [[ ${PACKAGE} = f77zmq ]] ; then elif [[ ${PACKAGE} = f77zmq ]] ; then
download ${F77ZMQ_URL} "${QP_ROOT}"/external/f77_zmq.tar.gz
execute << EOF execute << EOF
cd "\${QP_ROOT}"/external cd "\${QP_ROOT}"/external
tar --gunzip --extract --file f77_zmq.tar.gz tar --gunzip --extract --file qp2-dependencies/f77-zmq-4.3.2.tar.gz
rm f77_zmq.tar.gz cd f77-zmq-*
cd f77_zmq-* ./configure --prefix=\$QP_ROOT
export ZMQ_H="\$QP_ROOT"/include/zmq.h export ZMQ_H="\$QP_ROOT"/include/zmq.h
make make && make check && make install
cp libf77zmq.a "\${QP_ROOT}"/lib
cp libf77zmq.so "\${QP_ROOT}"/lib
cp f77_zmq_free.h "\${QP_ROOT}"/include
EOF EOF
elif [[ ${PACKAGE} = ocaml ]] ; then elif [[ ${PACKAGE} = ocaml ]] ; then
download ${OCAML_URL} "${QP_ROOT}"/external/opam_installer.sh execute <<EOF
source "${QP_ROOT}"/quantum_package.rc
if [[ -n ${TRAVIS} ]] ; then cd "${QP_ROOT}"/external/
# Special commands for Travis CI tar --gunzip --extract --file qp2-dependencies/ocaml-bundle_x86.tar.gz
chmod +x "${QP_ROOT}"/external/opam_installer.sh echo "" | ./ocaml-bundle/bootstrap.sh "${QP_ROOT}"
rm --force ${QP_ROOT}/bin/opam ./ocaml-bundle/configure.sh "${QP_ROOT}"
if [[ -n ${NO_CACHE} ]] ; then echo "" | ./ocaml-bundle/compile.sh "${QP_ROOT}"
rm -rf ${HOME}/.opam
fi
export OPAMROOT=${HOME}/.opam
cat << EOF | bash ${QP_ROOT}/external/opam_installer.sh --no-backup
${QP_ROOT}/bin
EOF 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 elif [[ ${PACKAGE} = bse ]] ; then
download ${BSE_URL} "${QP_ROOT}"/external/bse.tar.gz
execute << EOF execute << EOF
cd "\${QP_ROOT}"/external 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-* pip install -e basis_set_exchange-*
EOF EOF
elif [[ ${PACKAGE} = zlib ]] ; then elif [[ ${PACKAGE} = zlib ]] ; then
download ${ZLIB_URL} "${QP_ROOT}"/external/zlib.tar.gz
execute << EOF execute << EOF
cd "\${QP_ROOT}"/external cd "\${QP_ROOT}"/external
tar --gunzip --extract --file zlib.tar.gz tar --gunzip --extract --file qp2-dependencies/zlib-1.2.11.tar.gz
rm zlib.tar.gz && \
cd zlib-*/ cd zlib-*/
./configure --prefix=${QP_ROOT} && \ ./configure --prefix=${QP_ROOT} && \
make && make install make && make install
@ -372,33 +263,27 @@ EOF
elif [[ ${PACKAGE} = docopt ]] ; then elif [[ ${PACKAGE} = docopt ]] ; then
download ${DOCOPT_URL} "${QP_ROOT}"/external/docopt.tar.gz
execute << EOF execute << EOF
cd "\${QP_ROOT}"/external 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" mv docopt-*/docopt.py "\${QP_ROOT}/external/Python"
rm --recursive --force -- docopt-*/ docopt.tar.gz
EOF EOF
elif [[ ${PACKAGE} = resultsFile ]] ; then elif [[ ${PACKAGE} = resultsFile ]] ; then
download ${RESULTS_URL} "${QP_ROOT}"/external/resultsFile.tar.gz
execute << EOF execute << EOF
cd "\${QP_ROOT}"/external 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/" mv resultsFile-*/resultsFile "\${QP_ROOT}/external/Python/"
rm --recursive --force resultsFile-* resultsFile.tar.gz
EOF EOF
elif [[ ${PACKAGE} = bats ]] ; then elif [[ ${PACKAGE} = bats ]] ; then
download ${BATS_URL} "${QP_ROOT}"/external/bats.tar.gz
execute << EOF execute << EOF
cd "\${QP_ROOT}"/external 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}) ( 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 EOF
else else
@ -417,12 +302,6 @@ if [[ ${NINJA} = $(not_found) ]] ; then
fail fail
fi fi
IRPF90=$(find_exe irpf90)
if [[ ${IRPF90} = $(not_found) ]] ; then
error "IRPF90 (irpf90) is not installed."
fail
fi
ZEROMQ=$(find_lib -lzmq) ZEROMQ=$(find_lib -lzmq)
if [[ ${ZEROMQ} = $(not_found) ]] ; then if [[ ${ZEROMQ} = $(not_found) ]] ; then
error "ZeroMQ (zeromq) is not installed." error "ZeroMQ (zeromq) is not installed."
@ -441,24 +320,6 @@ if [[ ${ZLIB} = $(not_found) ]] ; then
fail fail
fi 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) OCAML=$(find_exe ocaml)
if [[ ${OCAML} = $(not_found) ]] ; then if [[ ${OCAML} = $(not_found) ]] ; then
error "OCaml (ocaml) compiler is not installed." error "OCaml (ocaml) compiler is not installed."

1
etc/cflags.rc Normal file
View File

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

View File

@ -4,8 +4,10 @@ if [[ -z $OPAMROOT ]]
then then
# Comment these lines if you have a system-wide OCaml installation # 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 fi
source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true

1
etc/openmp.rc Normal file
View File

@ -0,0 +1 @@
export OMP_NESTED=True

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 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 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 C_INCLUDE_PATH=$(qp_prepend_export "C_INCLUDE_PATH" "${QP_ROOT}"/include)
export CPATH=$(qp_prepend_export "CPATH" "${QP_ROOT}"/include) export CPATH=$(qp_prepend_export "CPATH" "${QP_ROOT}"/include)

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/

Binary file not shown.

1
external/qp2-dependencies vendored Submodule

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

View File

@ -709,6 +709,11 @@ def save_subninja_file(path_module):
" description = Cleaning module {0}".format(path_module.rel), " description = Cleaning module {0}".format(path_module.rel),
""] ""]
l_string += ["rule make_tidy",
" command = module_handler.py tidy {0}".format(path_module.rel),
" description = Cleaning module {0}".format(path_module.rel),
""]
l_string += ["rule executables", l_string += ["rule executables",
" command = make -C {0} executables .gitignore qp_edit.native qp_run.native".format(join("$QP_ROOT","ocaml")), " command = make -C {0} executables .gitignore qp_edit.native qp_run.native".format(join("$QP_ROOT","ocaml")),
" description = Updating OCaml executables", " description = Updating OCaml executables",
@ -719,6 +724,7 @@ def save_subninja_file(path_module):
"build local: make_local_binaries dummy_target", "", "build local: make_local_binaries dummy_target", "",
"build executables: executables local dummy_target", "", "build executables: executables local dummy_target", "",
"default executables", "", "build clean: make_clean dummy_target", "default executables", "", "build clean: make_clean dummy_target",
"", "build tidy: make_tidy dummy_target",
""] ""]
path_ninja_cur = join(path_module.abs, "build.ninja") path_ninja_cur = join(path_module.abs, "build.ninja")
@ -745,6 +751,10 @@ def create_build_ninja_global():
" command = module_handler.py clean --all", " command = module_handler.py clean --all",
" description = Cleaning all modules", ""] " description = Cleaning all modules", ""]
l_string += ["rule make_tidy",
" command = module_handler.py tidy --all",
" description = Cleaning all modules", ""]
l_string += ["rule make_ocaml", l_string += ["rule make_ocaml",
" command = make -C {0}/ocaml".format("$QP_ROOT"), " command = make -C {0}/ocaml".format("$QP_ROOT"),
" pool = console", " pool = console",
@ -759,6 +769,8 @@ def create_build_ninja_global():
"default ocaml_target", "default ocaml_target",
"", "",
"build clean: make_clean dummy_target", "build clean: make_clean dummy_target",
"",
"build tidy: make_tidy dummy_target",
"", ] "", ]
path_ninja_cur = join(QP_ROOT, "build.ninja") path_ninja_cur = join(QP_ROOT, "build.ninja")

View File

@ -6,12 +6,18 @@ Module utilitary
Usage: Usage:
module_handler.py print_descendant [<module_name>...] module_handler.py print_descendant [<module_name>...]
module_handler.py clean [ --all | <module_name>...] module_handler.py clean [ --all | <module_name>...]
module_handler.py create_git_ignore [<module_name>...] module_handler.py tidy [ --all | <module_name>...]
module_handler.py create_git_ignore [ --all | <module_name>...]
Options: Options:
print_descendant Print the genealogy of the needed modules print_descendant Print the genealogy of the needed modules
clean Used for ninja clean
tidy A light version of clean, where only the intermediate
files are removed
create_git_ignore deprecated
NEED The path of NEED file. NEED The path of NEED file.
by default try to open the file in the current path by default try to open the file in the current path
""" """
import os import os
import sys import sys
@ -25,7 +31,7 @@ try:
from docopt import docopt from docopt import docopt
from qp_path import QP_SRC, QP_ROOT, QP_PLUGINS, QP_EZFIO from qp_path import QP_SRC, QP_ROOT, QP_PLUGINS, QP_EZFIO
except ImportError: except ImportError:
print("source .quantum_package.rc") print("source quantum_package.rc")
raise raise
@ -230,7 +236,7 @@ if __name__ == '__main__':
for module in l_module: for module in l_module:
print(" ".join(sorted(m.l_descendant_unique([module])))) print(" ".join(sorted(m.l_descendant_unique([module]))))
if arguments["clean"]: if arguments["clean"] or arguments["tidy"]:
l_dir = ['IRPF90_temp', 'IRPF90_man'] l_dir = ['IRPF90_temp', 'IRPF90_man']
l_file = ["irpf90_entities", "tags", "irpf90.make", "Makefile", l_file = ["irpf90_entities", "tags", "irpf90.make", "Makefile",
@ -242,7 +248,6 @@ if __name__ == '__main__':
l_symlink = m.l_descendant_unique([module]) l_symlink = m.l_descendant_unique([module])
l_exe = get_binaries(module_abs) l_exe = get_binaries(module_abs)
if arguments["clean"]:
for f in l_dir: for f in l_dir:
try: try:
shutil.rmtree(os.path.join(module_abs, f)) shutil.rmtree(os.path.join(module_abs, f))
@ -261,6 +266,7 @@ if __name__ == '__main__':
except: except:
pass pass
if arguments["clean"]:
for f in l_exe: for f in l_exe:
try: try:
@ -268,6 +274,4 @@ if __name__ == '__main__':
except: except:
pass pass
if arguments["create_git_ignore"]:
pass

View File

@ -0,0 +1,365 @@
#!/bin/sh
set -ue
# (c) Copyright Fabrice Le Fessant INRIA/OCamlPro 2013
# (c) Copyright Louis Gesbert OCamlPro 2014-2017
VERSION='2.0.7'
DEV_VERSION='2.1.0~beta2'
DEFAULT_BINDIR=/usr/local/bin
bin_sha512() {
case "$OPAM_BIN" in
opam-2.0.6-arm64-linux) echo "d2b3d92fd5fae7f053702b53ddbc7c224fcfbfc9b232247ba4e40cbf1cda28f160d8c14fde87aebeebfd2545e13265c0ee4a47e292f035767fb944b1b8ff5c90";;
opam-2.0.6-armhf-linux) echo "a42a7ad8c1afdb20ac5746934306576e6364f5453b176ccd42a3e5a116a5db05c2758cec31800ffab11411290cf671f9eee3f299df48c7ceca8e4d7e33dfedc8";;
opam-2.0.6-i686-linux) echo "6c0d965f89a2026ead3120e217d12b2df7426740d54bc94e2c46faaeff5893081e68aac162621bfa694ab597a18be28165f10cdda1217a4d73653789a9928b64";;
opam-2.0.6-x86_64-linux) echo "2b9d4a99aa28a193c88c7c6f6265203bd3cfeef98929d6f5cfce4b52cd9ddbd7be7eddc1d3d9c440f81d65074dd7851b8d29cd397fb06d2cfccffb54d3cdcc6a";;
opam-2.0.6-x86_64-macos) echo "cf02546b22ca91b1d97a3657b970b34d4acf4dc745696b7200ff185d25ebb5914ea8b6a94b503eb8c999634de6fdb944998a970105cd6a4c6df538c262b48b7f";;
opam-2.0.6-x86_64-openbsd) echo "2f58b3d4902d4c3fb823d251a50e034f9101b0c5a3827725876bb3bcb6c013c4f54138054d82abba0a9e917675275e26f05b98630cf7116c465d2110756f1309";;
opam-2.0.7-arm64-linux) echo "0dd4d80496545f684af39dc5b4b28867bc19a74186577c38bd2a8934d871c2cbcdb9891bfd41c080b5f12d5a3c8801e203df8a76d55e1e22fe80d31447402e46";;
opam-2.0.7-armhf-linux) echo "ea691bc9565acc1207dea3dfb89192b1865b5b5809efe804a329f39878640fb19771edcb05c5699f8e914e88e3155f31132b845c54b0095bedd3952d336bae0b";;
opam-2.0.7-i686-linux) echo "5fa8fb9664d36ead5760e7e1c337f6ae7b0fd4be5089ddfb50ae74028deec30893b1f4dee040402bc3f15da197ba89a45c7d626ecf6e5be80d176f43526c4bad";;
opam-2.0.7-x86_64-linux) echo "da75b0cb5ad50f95d31857a7d72f1836132a1fa1cdbfdedf684342b798e7107b4add4c74c05d5ce44881309fa1e57707538dbcda874e7f74b269b1bb204f3ae3";;
opam-2.0.7-x86_64-macos) echo "de1194c8e97e53956e5e47502c28881bbf26d1beaac4f33a43a922b8ca7ce97725533cfaf65a33fc0e183eab5a95e9ecd2e20f72faeaec333dc3850b79b5fe8a";;
opam-2.0.7-x86_64-openbsd) echo "b253809c4388847e1a33b5c4f1f5d72bef79a2f0c43b19ef65b40d0c10341aa0bee4a4b1f3a9ab70eb026e4cc220a63cfc56a18c035b6b0297c92f2bdb7f9a78";;
opam-2.1.0-alpha-arm64-linux) echo "1bf0acfa64aa01c3244e65eed60eef1caaa6de53aa8b32dd0d2446f91905a1e41591f53cd350e85b2b9f5edba9b137d723c32949115623e9753e77b707bb25b0";;
opam-2.1.0-alpha-armhf-linux) echo "87c12a422bd14a0d10a94ddaaa46de23700e3b89810a0c06232eff8d96b37c2fd43dcb5a8da5a2004aa8040d1b93293209f1ff1aab865ffd150364e24c87c716";;
opam-2.1.0-alpha-i686-linux) echo "b8369da6d4795a461ff1b49e687b027325d4e90bc8f19517e52a94ee3be167c4faaaf33bd0b3536be552d2add54865d0e33933acaa674f2e1a17249b022738af";;
opam-2.1.0-alpha-x86_64-linux) echo "2e22747829fb0bada3a74a23f5e0ff2228520d647fc4fe08a1ce76f3cb357cc7240f7b45e422c5f4b8eafe832ae3a8973ecbd4814ae0e8ce1096bcff39482020";;
opam-2.1.0-alpha-x86_64-macos) echo "c440e8ae1970fa7533e6e1b96ba3e3dd65b04432d41bc57ce4c768ed9b4229954546d59ec06f3d4ee49cbe00bb4bfd0b3f509d6d9a27de2db17725e097a61c86";;
opam-2.1.0-alpha-x86_64-openbsd) echo "d87afe99fee541a1c6fae30b72653db7a5ea2abdec3fa3b2b480daddf3fcd8d4096e2a40458310755faec3722119f29ed981ffbfa65142e618f99b70572f892f";;
opam-2.1.0-alpha2-arm64-linux) echo "b67520bb2a6c59f800da100278d74e58f2bbf66924f94643023dc46b97b16f17a30de95e439c6f9b032bd555c062ddba325f3e5169cac186615b959a8c434788";;
opam-2.1.0-alpha2-armhf-linux) echo "9a6312eb54d6c9c2036ca90f7816789c27c23f1b1d325cd69d27a910cdd8760b82f19c9e9b61b5b6214818f1f40f8b4d2ef081acb43f0dad68c976986a7c6a45";;
opam-2.1.0-alpha2-i686-linux) echo "0dc07f236405777ad74d58fcc6cb6c3247e7dfc31408df4a199599077d5cb41ec86895f1d0c5eaa2a9c70842a2a998226674f986ba0044c82896c073ac90b209";;
opam-2.1.0-alpha2-x86_64-linux) echo "21509e8abd8463f4e18a55398f690700772e25f0ddb9f3fd7644e2f9a9a89ebbf5c09efbeceafe4a0ab5015d0d03b2f29506be514aae813a2f3dac7dd01261f3";;
opam-2.1.0-alpha2-x86_64-macos) echo "1c1bd26621eebb5bf3783dec80d5555aa5ff02dcbf43eb44398798e6162c1964bc1964e3980391ea115e5c068c1bb66960f8ebdd91bc4f0bac844f3a61433f1e";;
opam-2.1.0-alpha2-x86_64-openbsd) echo "941f3e306bc36e8e44e4245ca5e635b04e0a54f33439d55d41875ced47384cad8c222b649027d3c4eacc3c2c569cf5006c872763b19c490d9b289c9cfe4f491a";;
opam-2.1.0-alpha3-arm64-linux) echo "ad906bb2ab764a92fabdf0b906310c5034bf5daf0ebfb2529e9b87661ddbf8fd14f51dee5ce75b4fd4bb5789e29c7be71063f1ebcc92e92333be12aa62efdff9";;
opam-2.1.0-alpha3-armhf-linux) echo "2a7022c1f5dbc855a0d067f29677b13253dccbc9792b8170fa72a743802bbcd6e41ce7512c4845091af0f73b8ba7573038ec53ea9aaf74be04367ac1767e7220";;
opam-2.1.0-alpha3-i686-linux) echo "6f2fce0c45ae700e7a1b32d0a24988645c9aed3afc45998c8fbe70e97a65e3ba5d824069914a892bb3f9b1336383cfd492c28678ff16db5cada863da924b07d8";;
opam-2.1.0-alpha3-x86_64-linux) echo "1d219dbf670e1550bf71c28e586d14f1d8af2605f0e13bea2f11ad52a7f176bd9a89637e44a91a024f0088db1b2aba8dc3207bc81fa930580e54f4031255c178";;
opam-2.1.0-alpha3-x86_64-macos) echo "93edb6c1151f8f5bd017f230ffd9277f6ad943e3f5032ea000c37f012738fb3ab4b4add172e1f624c37e6564963fef0716b876b0113c8e43f5943d77bbbc173c";;
opam-2.1.0-alpha3-x86_64-openbsd) echo "0e3b3761e877c57f5b333aacb70c86bf60f50eecdca6e9e1a552e3d666cea034d8873f3a87e585a5970b1aef7e540adb18c71e0e8fd8794843dd5d1d421a87ec";;
opam-2.1.0-beta-arm64-linux) echo "954670c74ea8244b440756e4f7755bd2b5548ab67428ce577c4c507fc33c8d00eb73c4d7b59ccb0ef800f4465b5c704573c63486b78a23e9568f3751bf9aef78";;
opam-2.1.0-beta-armhf-linux) echo "cc666f2c6b1ac07d1bc8a035c6b3a9455794b51a827c54bb92786ae1a75c6c55839d3f48b378508f42a66ac887fdc68f7628a67e2826813cb6df048c906755ca";;
opam-2.1.0-beta-i686-linux) echo "66ac48b298741f753ca868be362851ccd9bf84fd8772d18f3307e99cf72c8c68ac9fa17bf2d610d7f3b5dc6209eb8371bf0e10b363e963fc6c31d70e5938017f";;
opam-2.1.0-beta-x86_64-linux) echo "e316f1b5f1c668affba6c2819f692c28776e131a17fb64b2c0e23f8a3b7d456575a8109fcdcb9babfad13bc33c17fa619cbb4a48ca6198765f86296b7e611f24";;
opam-2.1.0-beta-x86_64-macos) echo "acb29b7c64df314c6629e14f6d8f079504d39b7fd3104867fd22df3395ccfea9f1014a3a87dff9c12bf03ca451e9ee2918b9d9d8f17ce1a6d7de0c0649452fa9";;
opam-2.1.0-beta-x86_64-openbsd) echo "ff9fa1ee0ae7e54b4e18999cf5ea9b899c0b4039b744a950e96221e3e86c21eaa50904bdbc836ff8103f7713506d0de3d32ec77b169561e0cd694bfeea812cae";;
opam-2.1.0-beta2-arm64-linux) echo "a58ba3ebb4431d3cabfe96b806c9897205153e8a546ebe74f0229982758d140b4fcbcea421db70589b1eb3080dc86534522a3cba0330ce82e0898a60048d51ba";;
opam-2.1.0-beta2-armhf-linux) echo "fc4e6b753ce6368f75a0d3005f4b21ce9606599d21607a67015db55a38b6ef473b4205f5b128c5808189feed8ae58f93bd79348988be7c5007ae1b39307a5cd0";;
opam-2.1.0-beta2-i686-linux) echo "a376a6e0e1e2b08ea4d0a5c1c38487e67984bef2e89f978536dd08283f945f74dd31ee287bc68d91690603ba0fa657e91ff0d30bea217743f79ed99d2390eba5";;
opam-2.1.0-beta2-x86_64-linux) echo "12c5e2b0087ed389fa12fdb0e1f6f7dc0b3df3f95c59e8bc576279b7780921d47bbc4ebcba6caddde30f4fb1cc9e4a873cc8a6aef80fcc48a878aba69be7af44";;
opam-2.1.0-beta2-x86_64-macos) echo "4acc12672a2e3ad7e78540634edcae2e7e84860057b86a56b1cdf7eacf8d97957aaa864f571d6fb8f61ee8280f8a4ed73df7881d91a22c9d8c2d73e8a558f61d";;
opam-2.1.0-beta2-x86_64-openbsd) echo "84d7d409220c72e3ed7e6acdd7cce3b5a208f2966d232648a57a48641ab8ce4fa58e94e40b7176201455d82260e6c501a6ba4a30b1426a552f8d09cfd027ddde";;
*) echo "no sha";;
esac
}
usage() {
echo "opam binary installer v.$VERSION"
echo "Downloads and installs a pre-compiled binary of opam $VERSION to the system."
echo "This can also be used to switch between opam versions"
echo
echo "Options:"
echo " --dev Install the latest alpha or beta instead: $DEV_VERSION"
echo " --no-backup Don't attempt to backup the current opam root"
echo " --backup Force the backup the current opam root (even if it"
echo " is from the 2.0 branch already)"
echo " --fresh Create the opam $VERSION root from scratch"
echo " --restore VERSION Restore a backed up opam binary and root"
echo
echo "The default is to backup if the current version of opam is 1.*, or when"
echo "using '--fresh' or '--dev'"
}
RESTORE=
NOBACKUP=
FRESH=
DOWNLOAD_ONLY=
while [ $# -gt 0 ]; do
case "$1" in
--dev)
VERSION=$DEV_VERSION
if [ -z "$NOBACKUP" ]; then NOBACKUP=0; fi;;
--restore)
if [ $# -lt 2 ]; then echo "Option $1 requires an argument"; exit 2; fi
shift;
RESTORE=$1;;
--no-backup)
NOBACKUP=1;;
--backup)
NOBACKUP=0;;
--fresh)
FRESH=1;;
--download-only)
DOWNLOAD_ONLY=1;;
--help|-h)
usage; exit 0;;
*)
usage; exit 2;;
esac
shift
done
TMP=${TMPDIR:-/tmp}
ARCH=$(uname -m || echo unknown)
case "$ARCH" in
x86|i?86) ARCH="i686";;
x86_64|amd64) ARCH="x86_64";;
ppc|powerpc|ppcle) ARCH="ppc";;
aarch64_be|aarch64|armv8b|armv8l) ARCH="arm64";;
armv5*|armv6*|earmv6*|armv7*|earmv7*) ARCH="armhf";;
*) ARCH=$(echo "$ARCH" | awk '{print tolower($0)}')
esac
OS=$( (uname -s || echo unknown) | awk '{print tolower($0)}')
if [ "$OS" = "darwin" ] ; then
OS=macos
fi
TAG=$(echo "$VERSION" | tr '~' '-')
OPAM_BIN_URL_BASE='https://github.com/ocaml/opam/releases/download/'
OPAM_BIN="opam-${TAG}-${ARCH}-${OS}"
OPAM_BIN_URL="${OPAM_BIN_URL_BASE}${TAG}/${OPAM_BIN}"
download() {
if command -v wget >/dev/null; then wget -q -O "$@"
else curl -s -L -o "$@"
fi
}
check_sha512() {
OPAM_BIN_LOC="$1"
if command -v openssl > /dev/null; then
sha512_devnull="cf83e1357eefb8bdf1542850d66d8007d620e4050b5715dc83f4a921d36ce9ce47d0d13c5d85f2b0ff8318d2877eec2f63b931bd47417a81a538327af927da3e"
sha512_check=`openssl sha512 2>&1 < /dev/null | cut -f 2 -d ' '`
if [ "x$sha512_devnull" = "x$sha512_check" ]; then
sha512=`openssl sha512 "$OPAM_BIN_LOC" 2> /dev/null | cut -f 2 -d ' '`
check=`bin_sha512`
test "x$sha512" = "x$check"
else
echo "openssl 512 option not handled, binary integrity check can't be performed."
return 0
fi
else
echo "openssl not found, binary integrity check can't be performed."
return 0
fi
}
download_and_check() {
OPAM_BIN_LOC="$1"
echo "## Downloading opam $VERSION for $OS on $ARCH..."
if ! download "$OPAM_BIN_LOC" "$OPAM_BIN_URL"; then
echo "There may not yet be a binary release for your architecture or OS, sorry."
echo "See https://github.com/ocaml/opam/releases/tag/$TAG for pre-compiled binaries,"
echo "or run 'make cold' from https://github.com/ocaml/opam/archive/$TAG.tar.gz"
echo "to build from scratch"
exit 10
else
if check_sha512 "$OPAM_BIN_LOC"; then
echo "## Downloaded."
else
echo "Checksum mismatch, a problem occurred during download."
exit 10
fi
fi
}
DOWNLOAD_ONLY=${DOWNLOAD_ONLY:-0}
if [ $DOWNLOAD_ONLY -eq 1 ]; then
OPAM_BIN_LOC="$PWD/$OPAM_BIN"
if [ -e "$OPAM_BIN_LOC" ]; then
echo "Found opam binary in $OPAM_BIN_LOC ..."
if check_sha512 "$OPAM_BIN_LOC" ; then
echo "... with matching sha512"
exit 0;
else
echo "... with mismatching sha512, download the good one."
fi
fi
download_and_check "$OPAM_BIN_LOC"
exit 0;
fi
EXISTING_OPAM=$(command -v opam || echo)
EXISTING_OPAMV=
if [ -n "$EXISTING_OPAM" ]; then
EXISTING_OPAMV=$("$EXISTING_OPAM" --version || echo "unknown")
fi
FRESH=${FRESH:-0}
OPAMROOT=${OPAMROOT:-$HOME/.opam}
if [ ! -d "$OPAMROOT" ]; then FRESH=1; fi
if [ -z "$NOBACKUP" ] && [ ! "$FRESH" = 1 ] && [ -z "$RESTORE" ]; then
case "$EXISTING_OPAMV" in
2.*) NOBACKUP=1;;
*) NOBACKUP=0;;
esac
fi
xsudo() {
local CMD=$1; shift
local DST
for DST in "$@"; do : ; done
local DSTDIR=$(dirname "$DST")
if [ ! -w "$DSTDIR" ]; then
echo "Write access to $DSTDIR required, using 'sudo'."
echo "Command: $CMD $@"
if [ "$CMD" = "install" ]; then
sudo "$CMD" -g 0 -o root "$@"
else
sudo "$CMD" "$@"
fi
else
"$CMD" "$@"
fi
}
if [ -n "$RESTORE" ]; then
OPAM=$(command -v opam)
OPAMV=$("$OPAM" --version)
OPAM_BAK="$OPAM.$RESTORE"
OPAMROOT_BAK="$OPAMROOT.$RESTORE"
if [ ! -e "$OPAM_BAK" ] || [ ! -d "$OPAMROOT_BAK" ]; then
echo "No backup of opam $RESTORE was found"
exit 1
fi
if [ "$NOBACKUP" = 1 ]; then
printf "## This will clear $OPAM and $OPAMROOT. Continue ? [Y/n] "
read R
case "$R" in
""|"y"|"Y"|"yes")
xsudo rm -f "$OPAM"
rm -rf "$OPAMROOT";;
*) exit 1
esac
else
xsudo mv "$OPAM" "$OPAM.$OPAMV"
mv "$OPAMROOT" "$OPAMROOT.$OPAMV"
fi
xsudo mv "$OPAM_BAK" "$OPAM"
mv "$OPAMROOT_BAK" "$OPAMROOT"
printf "## Opam $RESTORE and its root were restored."
if [ "$NOBACKUP" = 1 ]; then echo
else echo " Opam $OPAMV was backed up."
fi
exit 0
fi
#if [ -e "$TMP/$OPAM_BIN" ] && ! check_sha512 "$TMP/$OPAM_BIN" || [ ! -e "$TMP/$OPAM_BIN" ]; then
download_and_check "$TMP/$OPAM_BIN"
#else
# echo "## Using already downloaded \"$TMP/$OPAM_BIN\""
#fi
if [ -n "$EXISTING_OPAM" ]; then
DEFAULT_BINDIR=$(dirname "$EXISTING_OPAM")
fi
while true; do
printf "## Where should it be installed ? [$DEFAULT_BINDIR] "
read BINDIR
if [ -z "$BINDIR" ]; then BINDIR="$DEFAULT_BINDIR"; fi
if [ -d "$BINDIR" ]; then break
else
printf "## $BINDIR does not exist. Create ? [Y/n] "
read R
case "$R" in
""|"y"|"Y"|"yes")
xsudo mkdir -p $BINDIR
break;;
esac
fi
done
#if [ -e "$EXISTING_OPAM" ]; then
# if [ "$NOBACKUP" = 1 ]; then
# xsudo rm -f "$EXISTING_OPAM"
# else
# xsudo mv "$EXISTING_OPAM" "$EXISTING_OPAM.$EXISTING_OPAMV"
# echo "## $EXISTING_OPAM backed up as $(basename $EXISTING_OPAM).$EXISTING_OPAMV"
# fi
#fi
if [ -d "$OPAMROOT" ]; then
if [ "$FRESH" = 1 ]; then
if [ "$NOBACKUP" = 1 ]; then
printf "## This will clear $OPAMROOT. Continue ? [Y/n] "
read R
case "$R" in
""|"y"|"Y"|"yes")
rm -rf "$OPAMROOT";;
*) exit 1
esac
else
mv "$OPAMROOT" "$OPAMROOT.$EXISTING_OPAMV"
echo "## $OPAMROOT backed up as $(basename $OPAMROOT).$EXISTING_OPAMV"
fi
echo "## opam $VERSION installed. Please run 'opam init' to get started"
elif [ ! "$NOBACKUP" = 1 ]; then
echo "## Backing up $OPAMROOT to $(basename $OPAMROOT).$EXISTING_OPAMV (this may take a while)"
if [ -e "$OPAMROOT.$EXISTING_OPAMV" ]; then
echo "ERROR: there is already a backup at $OPAMROOT.$EXISTING_OPAMV"
echo "Please move it away or run with --no-backup"
fi
FREE=$(df -k "$OPAMROOT" | awk 'NR>1 {print $4}')
NEEDED=$(du -sk "$OPAMROOT" | awk '{print $1}')
if ! [ $NEEDED -lt $FREE ]; then
echo "Error: not enough free space to backup. You can retry with --no-backup,"
echo "--fresh, or remove '$OPAMROOT'"
exit 1
fi
cp -a "$OPAMROOT" "$OPAMROOT.$EXISTING_OPAMV"
echo "## $OPAMROOT backed up as $(basename $OPAMROOT).$EXISTING_OPAMV"
fi
rm -f "$OPAMROOT"/repo/*/*.tar.gz*
fi
xsudo install -m 755 "$TMP/$OPAM_BIN" "$BINDIR/opam"
echo "## opam $VERSION installed to $BINDIR"
if [ ! "$FRESH" = 1 ]; then
echo "## Converting the opam root format & updating"
"$BINDIR/opam" init --reinit -ni
fi
WHICH=$(command -v opam || echo notfound)
case "$WHICH" in
"$BINDIR/opam") ;;
notfound) echo "## Remember to add $BINDIR to your PATH";;
*)
echo "## WARNING: 'opam' command found in PATH does not match the installed one:"
echo " - Installed: '$BINDIR/opam'"
echo " - Found: '$WHICH'"
echo "Make sure to remove the second or fix your PATH to use the new opam"
echo
esac
if [ ! "$NOBACKUP" = 1 ]; then
echo "## Run this script again with '--restore $EXISTING_OPAMV' to revert."
fi
rm -f $TMP/$OPAM_BIN

View File

@ -21,6 +21,21 @@ BEGIN_PROVIDER [ integer, ao_shell, (ao_num) ]
enddo enddo
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 END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_coef_normalized, (ao_num,ao_prim_num_max) ] BEGIN_PROVIDER [ double precision, ao_coef_normalized, (ao_num,ao_prim_num_max) ]

View File

@ -54,6 +54,13 @@
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)
ao_overlap(i,j) += c * overlap ao_overlap(i,j) += c * overlap
if(isnan(ao_overlap(i,j)))then
print*,'i,j',i,j
print*,'l,n',l,n
print*,'c,overlap',c,overlap
print*,overlap_x,overlap_y,overlap_z
stop
endif
ao_overlap_x(i,j) += c * overlap_x ao_overlap_x(i,j) += c * overlap_x
ao_overlap_y(i,j) += c * overlap_y ao_overlap_y(i,j) += c * overlap_y
ao_overlap_z(i,j) += c * overlap_z ao_overlap_z(i,j) += c * overlap_z

View File

@ -28,6 +28,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals, (ao_num,ao_num)]
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)] BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
use omp_lib
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Local pseudo-potential ! Local pseudo-potential
@ -37,26 +38,57 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
integer :: num_A,num_B integer :: num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3) double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(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 :: 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 :: 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 ao_pseudo_integrals_local = 0.d0
print*, 'Providing the nuclear electron pseudo integrals (local)' print*, 'Providing the nuclear electron pseudo integrals (local)'
call wall_time(wall_1) ! Dummy iteration for OpenMP
call cpu_time(cpu_1) 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 thread_num = 0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& !$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 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 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, & !$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() !$ thread_num = omp_get_thread_num()
wall_0 = wall_1 wall_0 = wall_1
!$OMP DO SCHEDULE (guided) !$OMP DO
do j = 1, ao_num 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) do m=1,ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i) beta = ao_expo_ordered_transp(m,i)
double precision :: c
c = 0.d0 c = 0.d0
if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& 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 cycle
endif endif
do k = 1, nucl_num do k = 1, nucl_num
double precision :: Z
Z = nucl_charge(k) Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3) C_center(1:3) = nucl_coord(k,1:3)
@ -128,6 +158,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_non_local, (ao_num,ao_num)] BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_non_local, (ao_num,ao_num)]
use omp_lib
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Non-local pseudo-potential ! Non-local pseudo-potential
@ -137,25 +168,27 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integrals_local, (ao_num,ao_num)]
integer :: num_A,num_B integer :: num_A,num_B
double precision :: A_center(3),B_center(3),C_center(3) double precision :: A_center(3),B_center(3),C_center(3)
integer :: power_A(3),power_B(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 :: 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 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 ao_pseudo_integrals_non_local = 0.d0
print*, 'Providing the nuclear electron pseudo integrals (non-local)' print*, 'Providing the nuclear electron pseudo integrals (non-local)'
call wall_time(wall_1) call wall_time(wall_1)
call cpu_time(cpu_1)
thread_num = 0 thread_num = 0
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& !$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 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 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,& !$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) do m=1,ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i) beta = ao_expo_ordered_transp(m,i)
double precision :: c
c = 0.d0 c = 0.d0
if (dabs(ao_coef_normalized_ordered_transp(l,j)*ao_coef_normalized_ordered_transp(m,i))& 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 endif
do k = 1, nucl_num do k = 1, nucl_num
double precision :: Z
Z = nucl_charge(k) Z = nucl_charge(k)
C_center(1:3) = nucl_coord(k,1:3) C_center(1:3) = nucl_coord(k,1:3)

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) 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) 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 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 Vloc=0.d0
return return
endif endif
@ -1839,7 +1839,7 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
m_1 = m+m+1 m_1 = m+m+1
nlm = n+m+l nlm = n+m+l
pi=dacos(-1.d0) pi=dacos(-1.d0)
a_over_b_square = (a/b)**2 a_over_b_square = (a*a)/(b*b)
! First term of the sequence ! 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) qk = dble(q)
two_qkmp1 = 2.d0*(qk+mk)+1.d0 two_qkmp1 = 2.d0*(qk+mk)+1.d0
do k=0,q-1 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 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 sum=sum+s_q_k
two_qkmp1 = two_qkmp1-2.d0 two_qkmp1 = two_qkmp1-2.d0
qk = qk-1.d0 qk = qk-1.d0
enddo enddo
inverses(q) = a_over_b_square/(dble(q+n+q+n+3) * dble(q+1)) 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 int=int+sum
@ -1892,7 +1887,6 @@ double precision function int_prod_bessel(l,gam,n,m,a,b,arg)
else else
!Compute the s_q+1_0 !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) 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) if(mod(n+m+l,2).eq.1)s_q_0=s_q_0*dsqrt(pi*.5d0)

View File

@ -327,6 +327,8 @@ double precision function get_ao_two_e_integral(i,j,k,l,map) result(result)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Gets one AO bi-electronic integral from the AO map ! Gets one AO bi-electronic integral from the AO map
!
! i,j,k,l in physicist notation <ij|kl>
END_DOC END_DOC
integer, intent(in) :: i,j,k,l integer, intent(in) :: i,j,k,l
integer(key_kind) :: idx integer(key_kind) :: idx

View File

@ -39,7 +39,7 @@ interface: ezfio, provider
[shell_prim_index] [shell_prim_index]
type: integer type: integer
doc: Max number of primitives in a shell doc: Index of the first primitive of the shell
size: (basis.shell_num) size: (basis.shell_num)
interface: ezfio, provider interface: ezfio, provider

View File

@ -7,22 +7,10 @@ program basis_correction
touch read_wf touch read_wf
no_core_density = .True. no_core_density = .True.
touch no_core_density touch no_core_density
if(io_mo_two_e_integrals .ne. "Read")then
provide ao_two_e_integrals_in_map provide ao_two_e_integrals_in_map
endif
provide mo_two_e_integrals_in_map provide mo_two_e_integrals_in_map
call print_basis_correction call print_basis_correction
! call print_e_b
end end
subroutine print_e_b
implicit none
print *, 'Hello world'
print*,'ecmd_lda_mu_of_r = ',ecmd_lda_mu_of_r
print*,'ecmd_pbe_ueg_mu_of_r = ',ecmd_pbe_ueg_mu_of_r
print*,'ecmd_pbe_ueg_eff_xi_mu_of_r = ',ecmd_pbe_ueg_eff_xi_mu_of_r
print*,''
print*,'psi_energy + E^B_LDA = ',psi_energy + ecmd_lda_mu_of_r
print*,'psi_energy + E^B_PBE_UEG = ',psi_energy + ecmd_pbe_ueg_mu_of_r
print*,'psi_energy + E^B_PBE_UEG_Xi = ',psi_energy + ecmd_pbe_ueg_eff_xi_mu_of_r
print*,''
print*,'mu_average_prov = ',mu_average_prov
end

View File

@ -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) write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
enddo 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*, ''
print*,'Using a CAS-like two-body density to define mu(r)' 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 ' print*,'This assumes that the CAS is a qualitative representation of the wave function '

View File

@ -0,0 +1,28 @@
program basis_corr_su_pbe_ot
implicit none
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
read_wf = .True.
touch read_wf
no_core_density = .True.
touch no_core_density
if(io_mo_two_e_integrals .ne. "Read")then
provide ao_two_e_integrals_in_map
endif
provide mo_two_e_integrals_in_map
call print_su_pbe_ot
end
subroutine print_su_pbe_ot
implicit none
integer :: istate
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
enddo
do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD SU-PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate)
enddo
end

View File

@ -33,3 +33,34 @@ doc: Number of angular grid points given from input. Warning, this number cannot
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: 1202 default: 1202
[extra_grid_type_sgn]
type: integer
doc: Type of extra_grid used for the Becke's numerical extra_grid. Can be, by increasing accuracy: [ 0 | 1 | 2 | 3 ]
interface: ezfio,provider,ocaml
default: 0
[thresh_extra_grid]
type: double precision
doc: threshold on the weight of a given extra_grid point
interface: ezfio,provider,ocaml
default: 1.e-20
[my_extra_grid_becke]
type: logical
doc: if True, the number of angular and radial extra_grid points are read from EZFIO
interface: ezfio,provider,ocaml
default: False
[my_n_pt_r_extra_grid]
type: integer
doc: Number of radial extra_grid points given from input
interface: ezfio,provider,ocaml
default: 300
[my_n_pt_a_extra_grid]
type: integer
doc: Number of angular extra_grid points given from input. Warning, this number cannot be any integer. See file list_angular_extra_grid
interface: ezfio,provider,ocaml
default: 1202

View File

@ -0,0 +1,104 @@
BEGIN_PROVIDER [double precision, angular_quadrature_points_extra, (n_points_extra_integration_angular,3) ]
&BEGIN_PROVIDER [double precision, weights_angular_points_extra, (n_points_extra_integration_angular)]
implicit none
BEGIN_DOC
! weights and grid points_extra for the integration on the angular variables on
! the unit sphere centered on (0,0,0)
! According to the LEBEDEV scheme
END_DOC
include 'constants.include.F'
integer :: i
double precision :: accu
double precision :: degre_rad
double precision :: x(n_points_extra_integration_angular)
double precision :: y(n_points_extra_integration_angular)
double precision :: z(n_points_extra_integration_angular)
double precision :: w(n_points_extra_integration_angular)
degre_rad = pi/180.d0
accu = 0.d0
select case (n_points_extra_integration_angular)
case (0006)
call LD0006(X,Y,Z,W,n_points_extra_integration_angular)
case (0014)
call LD0014(X,Y,Z,W,n_points_extra_integration_angular)
case (0026)
call LD0026(X,Y,Z,W,n_points_extra_integration_angular)
case (0038)
call LD0038(X,Y,Z,W,n_points_extra_integration_angular)
case (0050)
call LD0050(X,Y,Z,W,n_points_extra_integration_angular)
case (0074)
call LD0074(X,Y,Z,W,n_points_extra_integration_angular)
case (0086)
call LD0086(X,Y,Z,W,n_points_extra_integration_angular)
case (0110)
call LD0110(X,Y,Z,W,n_points_extra_integration_angular)
case (0146)
call LD0146(X,Y,Z,W,n_points_extra_integration_angular)
case (0170)
call LD0170(X,Y,Z,W,n_points_extra_integration_angular)
case (0194)
call LD0194(X,Y,Z,W,n_points_extra_integration_angular)
case (0230)
call LD0230(X,Y,Z,W,n_points_extra_integration_angular)
case (0266)
call LD0266(X,Y,Z,W,n_points_extra_integration_angular)
case (0302)
call LD0302(X,Y,Z,W,n_points_extra_integration_angular)
case (0350)
call LD0350(X,Y,Z,W,n_points_extra_integration_angular)
case (0434)
call LD0434(X,Y,Z,W,n_points_extra_integration_angular)
case (0590)
call LD0590(X,Y,Z,W,n_points_extra_integration_angular)
case (0770)
call LD0770(X,Y,Z,W,n_points_extra_integration_angular)
case (0974)
call LD0974(X,Y,Z,W,n_points_extra_integration_angular)
case (1202)
call LD1202(X,Y,Z,W,n_points_extra_integration_angular)
case (1454)
call LD1454(X,Y,Z,W,n_points_extra_integration_angular)
case (1730)
call LD1730(X,Y,Z,W,n_points_extra_integration_angular)
case (2030)
call LD2030(X,Y,Z,W,n_points_extra_integration_angular)
case (2354)
call LD2354(X,Y,Z,W,n_points_extra_integration_angular)
case (2702)
call LD2702(X,Y,Z,W,n_points_extra_integration_angular)
case (3074)
call LD3074(X,Y,Z,W,n_points_extra_integration_angular)
case (3470)
call LD3470(X,Y,Z,W,n_points_extra_integration_angular)
case (3890)
call LD3890(X,Y,Z,W,n_points_extra_integration_angular)
case (4334)
call LD4334(X,Y,Z,W,n_points_extra_integration_angular)
case (4802)
call LD4802(X,Y,Z,W,n_points_extra_integration_angular)
case (5294)
call LD5294(X,Y,Z,W,n_points_extra_integration_angular)
case (5810)
call LD5810(X,Y,Z,W,n_points_extra_integration_angular)
case default
print *, irp_here//': wrong n_points_extra_integration_angular. See in ${QP_ROOT}/src/becke_numerical_grid/list_angular_grid to see the possible angular grid points_extra. Ex: '
print *, '[ 50 | 74 | 170 | 194 | 266 | 302 | 590 | 1202 | 2030 | 5810 ]'
stop -1
end select
do i = 1, n_points_extra_integration_angular
angular_quadrature_points_extra(i,1) = x(i)
angular_quadrature_points_extra(i,2) = y(i)
angular_quadrature_points_extra(i,3) = z(i)
weights_angular_points_extra(i) = w(i) * 4.d0 * pi
accu += w(i)
enddo
END_PROVIDER

View File

@ -0,0 +1,178 @@
BEGIN_PROVIDER [integer, n_points_extra_radial_grid]
&BEGIN_PROVIDER [integer, n_points_extra_integration_angular]
implicit none
BEGIN_DOC
! n_points_extra_radial_grid = number of radial grid points_extra per atom
!
! n_points_extra_integration_angular = number of angular grid points_extra per atom
!
! These numbers are automatically set by setting the grid_type_sgn parameter
END_DOC
if(.not.my_extra_grid_becke)then
select case (extra_grid_type_sgn)
case(0)
n_points_extra_radial_grid = 23
n_points_extra_integration_angular = 170
case(1)
n_points_extra_radial_grid = 50
n_points_extra_integration_angular = 194
case(2)
n_points_extra_radial_grid = 75
n_points_extra_integration_angular = 302
case(3)
n_points_extra_radial_grid = 99
n_points_extra_integration_angular = 590
case default
write(*,*) '!!! Quadrature grid not available !!!'
stop
end select
else
n_points_extra_radial_grid = my_n_pt_r_extra_grid
n_points_extra_integration_angular = my_n_pt_a_extra_grid
endif
END_PROVIDER
BEGIN_PROVIDER [integer, n_points_extra_grid_per_atom]
implicit none
BEGIN_DOC
! Number of grid points_extra per atom
END_DOC
n_points_extra_grid_per_atom = n_points_extra_integration_angular * n_points_extra_radial_grid
END_PROVIDER
BEGIN_PROVIDER [double precision, grid_points_extra_radial, (n_points_extra_radial_grid)]
&BEGIN_PROVIDER [double precision, dr_radial_extra_integral]
implicit none
BEGIN_DOC
! points_extra in [0,1] to map the radial integral [0,\infty]
END_DOC
dr_radial_extra_integral = 1.d0/dble(n_points_extra_radial_grid-1)
integer :: i
do i = 1, n_points_extra_radial_grid
grid_points_extra_radial(i) = dble(i-1) * dr_radial_extra_integral
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, grid_points_extra_per_atom, (3,n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num)]
BEGIN_DOC
! x,y,z coordinates of grid points_extra used for integration in 3d space
END_DOC
implicit none
integer :: i,j,k
double precision :: dr,x_ref,y_ref,z_ref
double precision :: knowles_function
do i = 1, nucl_num
x_ref = nucl_coord(i,1)
y_ref = nucl_coord(i,2)
z_ref = nucl_coord(i,3)
do j = 1, n_points_extra_radial_grid-1
double precision :: x,r
! x value for the mapping of the [0, +\infty] to [0,1]
x = grid_points_extra_radial(j)
! value of the radial coordinate for the integration
r = knowles_function(alpha_knowles(grid_atomic_number(i)),m_knowles,x)
! explicit values of the grid points_extra centered around each atom
do k = 1, n_points_extra_integration_angular
grid_points_extra_per_atom(1,k,j,i) = &
x_ref + angular_quadrature_points_extra(k,1) * r
grid_points_extra_per_atom(2,k,j,i) = &
y_ref + angular_quadrature_points_extra(k,2) * r
grid_points_extra_per_atom(3,k,j,i) = &
z_ref + angular_quadrature_points_extra(k,3) * r
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, weight_at_r_extra, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num) ]
BEGIN_DOC
! Weight function at grid points_extra : w_n(r) according to the equation (22)
! of Becke original paper (JCP, 88, 1988)
!
! The "n" discrete variable represents the nucleis which in this array is
! represented by the last dimension and the points_extra are labelled by the
! other dimensions.
END_DOC
implicit none
integer :: i,j,k,l,m
double precision :: r(3)
double precision :: accu,cell_function_becke
double precision :: tmp_array(nucl_num)
! run over all points_extra in space
! that are referred to each atom
do j = 1, nucl_num
!for each radial grid attached to the "jth" atom
do k = 1, n_points_extra_radial_grid -1
! for each angular point attached to the "jth" atom
do l = 1, n_points_extra_integration_angular
r(1) = grid_points_extra_per_atom(1,l,k,j)
r(2) = grid_points_extra_per_atom(2,l,k,j)
r(3) = grid_points_extra_per_atom(3,l,k,j)
accu = 0.d0
! For each of these points_extra in space, ou need to evaluate the P_n(r)
do i = 1, nucl_num
! function defined for each atom "i" by equation (13) and (21) with k == 3
tmp_array(i) = cell_function_becke(r,i) ! P_n(r)
! Then you compute the summ the P_n(r) function for each of the "r" points_extra
accu += tmp_array(i)
enddo
accu = 1.d0/accu
weight_at_r_extra(l,k,j) = tmp_array(j) * accu
if(isnan(weight_at_r_extra(l,k,j)))then
print*,'isnan(weight_at_r_extra(l,k,j))'
print*,l,k,j
accu = 0.d0
do i = 1, nucl_num
! function defined for each atom "i" by equation (13) and (21) with k == 3
tmp_array(i) = cell_function_becke(r,i) ! P_n(r)
print*,i,tmp_array(i)
! Then you compute the summ the P_n(r) function for each of the "r" points_extra
accu += tmp_array(i)
enddo
write(*,'(100(F16.10,X))')tmp_array(j) , accu
stop
endif
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, final_weight_at_r_extra, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num) ]
BEGIN_DOC
! Total weight on each grid point which takes into account all Lebedev, Voronoi and radial weights.
END_DOC
implicit none
integer :: i,j,k,l,m
double precision :: r(3)
double precision :: accu,cell_function_becke
double precision :: tmp_array(nucl_num)
double precision :: contrib_integration,x
double precision :: derivative_knowles_function,knowles_function
! run over all points_extra in space
do j = 1, nucl_num ! that are referred to each atom
do i = 1, n_points_extra_radial_grid -1 !for each radial grid attached to the "jth" atom
x = grid_points_extra_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1]
do k = 1, n_points_extra_integration_angular ! for each angular point attached to the "jth" atom
contrib_integration = derivative_knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)&
*knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)**2
final_weight_at_r_extra(k,i,j) = weights_angular_points_extra(k) * weight_at_r_extra(k,i,j) * contrib_integration * dr_radial_extra_integral
if(isnan(final_weight_at_r_extra(k,i,j)))then
print*,'isnan(final_weight_at_r_extra(k,i,j))'
print*,k,i,j
write(*,'(100(F16.10,X))')weights_angular_points_extra(k) , weight_at_r_extra(k,i,j) , contrib_integration , dr_radial_extra_integral
stop
endif
enddo
enddo
enddo
END_PROVIDER

View File

@ -0,0 +1,60 @@
BEGIN_PROVIDER [integer, n_points_extra_final_grid]
implicit none
BEGIN_DOC
! Number of points_extra which are non zero
END_DOC
integer :: i,j,k,l
n_points_extra_final_grid = 0
do j = 1, nucl_num
do i = 1, n_points_extra_radial_grid -1
do k = 1, n_points_extra_integration_angular
if(dabs(final_weight_at_r_extra(k,i,j)) < thresh_extra_grid)then
cycle
endif
n_points_extra_final_grid += 1
enddo
enddo
enddo
print*,'n_points_extra_final_grid = ',n_points_extra_final_grid
print*,'n max point = ',n_points_extra_integration_angular*(n_points_extra_radial_grid*nucl_num - 1)
! call ezfio_set_becke_numerical_grid_n_points_extra_final_grid(n_points_extra_final_grid)
END_PROVIDER
BEGIN_PROVIDER [double precision, final_grid_points_extra, (3,n_points_extra_final_grid)]
&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_extra, (n_points_extra_final_grid) ]
&BEGIN_PROVIDER [integer, index_final_points_extra, (3,n_points_extra_final_grid) ]
&BEGIN_PROVIDER [integer, index_final_points_extra_reverse, (n_points_extra_integration_angular,n_points_extra_radial_grid,nucl_num) ]
implicit none
BEGIN_DOC
! final_grid_points_extra(1:3,j) = (/ x, y, z /) of the jth grid point
!
! final_weight_at_r_vector_extra(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
!
! index_final_points_extra(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point
!
! index_final_points_extra_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices
END_DOC
integer :: i,j,k,l,i_count
double precision :: r(3)
i_count = 0
do j = 1, nucl_num
do i = 1, n_points_extra_radial_grid -1
do k = 1, n_points_extra_integration_angular
if(dabs(final_weight_at_r_extra(k,i,j)) < thresh_extra_grid)then
cycle
endif
i_count += 1
final_grid_points_extra(1,i_count) = grid_points_extra_per_atom(1,k,i,j)
final_grid_points_extra(2,i_count) = grid_points_extra_per_atom(2,k,i,j)
final_grid_points_extra(3,i_count) = grid_points_extra_per_atom(3,k,i,j)
final_weight_at_r_vector_extra(i_count) = final_weight_at_r_extra(k,i,j)
index_final_points_extra(1,i_count) = k
index_final_points_extra(2,i_count) = i
index_final_points_extra(3,i_count) = j
index_final_points_extra_reverse(k,i,j) = i_count
enddo
enddo
enddo
END_PROVIDER

View File

@ -143,10 +143,10 @@ subroutine print_generators_bitmasks_holes
key_tmp(j,1) = generators_bitmask(j,1,i) key_tmp(j,1) = generators_bitmask(j,1,i)
key_tmp(j,2) = generators_bitmask(j,2,i) key_tmp(j,2) = generators_bitmask(j,2,i)
enddo enddo
print*,'' ! print*,''
print*,'index hole = ',i ! print*,'index hole = ',i
call print_det(key_tmp,N_int) ! call print_det(key_tmp,N_int)
print*,'' ! print*,''
enddo enddo
deallocate(key_tmp) deallocate(key_tmp)

View File

@ -5,6 +5,7 @@ subroutine write_on_top_in_real_space
! This routines is a simple example of how to plot the on-top pair density on a simple 1D grid ! This routines is a simple example of how to plot the on-top pair density on a simple 1D grid
END_DOC END_DOC
double precision :: zmax,dz,r(3),on_top_in_r,total_density,zcenter,dist double precision :: zmax,dz,r(3),on_top_in_r,total_density,zcenter,dist
double precision :: core_dens, inact_dens,act_dens(2,1)
integer :: nz,i,istate integer :: nz,i,istate
character*(128) :: output character*(128) :: output
@ -15,7 +16,7 @@ subroutine write_on_top_in_real_space
i_unit_output = getUnitAndOpen(output,'w') i_unit_output = getUnitAndOpen(output,'w')
zmax = 2.0d0 zmax = 5.0d0
print*,'nucl_coord(1,3) = ',nucl_coord(1,3) print*,'nucl_coord(1,3) = ',nucl_coord(1,3)
print*,'nucl_coord(2,3) = ',nucl_coord(2,3) print*,'nucl_coord(2,3) = ',nucl_coord(2,3)
dist = dabs(nucl_coord(1,3) - nucl_coord(2,3)) dist = dabs(nucl_coord(1,3) - nucl_coord(2,3))
@ -29,13 +30,14 @@ subroutine write_on_top_in_real_space
r(3) = zcenter -zmax * 0.5d0 r(3) = zcenter -zmax * 0.5d0
print*,'r(3) = ',r(3) print*,'r(3) = ',r(3)
istate = 1 istate = 1
write(i_unit_output,*)" z, on-top(z), n(z) "
do i = 1, nz do i = 1, nz
call give_on_top_in_r_one_state(r,istate,on_top_in_r) call give_on_top_in_r_one_state(r,istate,on_top_in_r)
call give_cas_density_in_r(r,total_density) call give_cas_density_in_r(core_dens,inact_dens,act_dens,total_density,r)
write(i_unit_output,*)r(3),on_top_in_r,total_density write(i_unit_output,*)r(3),on_top_in_r,total_density
r(3) += dz r(3) += dz
enddo enddo
end end

View File

@ -1,7 +1,6 @@
perturbation perturbation
zmq zmq
mpi mpi
davidson_undressed
iterations iterations
two_body_rdm two_body_rdm
csf csf

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_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted
PROVIDE psi_det_hii selection_weight pseudo_sym 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 if (h0_type == 'CFG') then
PROVIDE psi_configuration_hii det_to_configuration PROVIDE psi_configuration_hii det_to_configuration

View File

@ -13,7 +13,7 @@ subroutine run_stochastic_cipsi
double precision :: rss double precision :: rss
double precision, external :: memory_of_double double precision, external :: memory_of_double
PROVIDE H_apply_buffer_allocated PROVIDE H_apply_buffer_allocated distributed_davidson mo_two_e_integrals_in_map
N_iter = 1 N_iter = 1
threshold_generators = 1.d0 threshold_generators = 1.d0

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_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_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym 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') call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection')

View File

@ -62,7 +62,8 @@ void getIthBF(Node *inode, int isomo, bool foundBF, int NSOMOMax, int getaddr, i
if(isomo == NSOMOMax){ if(isomo == NSOMOMax){
if(inode->addr == getaddr){ if(inode->addr == getaddr){
for(int i = NSOMOMax-1; i > -1; i--){ int i;
for(i = NSOMOMax-1; i > -1; i--){
vecBF[i] = inode->cpl; vecBF[i] = inode->cpl;
inode = inode->PREV; inode = inode->PREV;
} }
@ -150,7 +151,8 @@ void getIthDet(Node *inode, int isomo, bool foundBF, int NSOMOMax, int getaddr,
if(isomo == NSOMOMax){ if(isomo == NSOMOMax){
if(inode->addr == getaddr){ if(inode->addr == getaddr){
for(int i = NSOMOMax-1; i > -1; i--){ int i;
for(i = NSOMOMax-1; i > -1; i--){
vecBF[i] = inode->cpl; vecBF[i] = inode->cpl;
inode = inode->PREV; inode = inode->PREV;
} }
@ -224,7 +226,8 @@ void getDetlist(Node *inode, int isomo, int NSOMOMax, int *vecBF, int *detlist){
if(isomo == NSOMOMax){ if(isomo == NSOMOMax){
int idet=0; int idet=0;
for(int k=0;k<NSOMOMax;k++){ int k;
for(k=0;k<NSOMOMax;k++){
if(vecBF[k] == 1) idet = idet | (1<<(NSOMOMax-1-k)); if(vecBF[k] == 1) idet = idet | (1<<(NSOMOMax-1-k));
} }
detlist[inode->addr]=idet; detlist[inode->addr]=idet;

1
src/dav_general_mat/NEED Normal file
View File

@ -0,0 +1 @@
davidson_undressed

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.

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

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

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

View File

@ -34,6 +34,12 @@ doc: If |true|, a memory-mapped file may be used to store the W and S2 vectors i
default: True default: True
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
[csf_based]
type: logical
doc: If |true|, use the CSF-based algorithm
default: False
interface: ezfio,provider,ocaml
[distributed_davidson] [distributed_davidson]
type: logical type: logical
doc: If |true|, use the distributed algorithm doc: If |true|, use the distributed algorithm
@ -52,3 +58,8 @@ doc: Maximum number of determinants where |H| is fully diagonalized
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: 1000 default: 1000
[without_diagonal]
type: logical
doc: If |true|, don't use denominator
default: False
interface: ezfio,provider,ocaml

View File

@ -447,6 +447,7 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
! Compute residual vector and davidson step ! Compute residual vector and davidson step
! ----------------------------------------- ! -----------------------------------------
if (without_diagonal) then
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k) !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
do k=1,N_st_diag do k=1,N_st_diag
do i=1,sze do i=1,sze
@ -455,6 +456,15 @@ subroutine davidson_diag_csf_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
enddo enddo
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
else
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,k)
do k=1,N_st_diag
do i=1,sze
U(i,k) = (lambda(k) * U(i,k) - W(i,k) )
enddo
enddo
!$OMP END PARALLEL DO
endif
do k=1,N_st do k=1,N_st
residual_norm(k) = u_dot_u(U(1,k),sze) residual_norm(k) = u_dot_u(U(1,k),sze)

View File

@ -42,7 +42,7 @@ END_PROVIDER
logical :: converged logical :: converged
logical :: do_csf 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 ! Guess values for the "N_states" states of the |CI| eigenvectors
do j=1,min(N_states,N_det) do j=1,min(N_states,N_det)
do i=1,N_det do i=1,N_det
@ -56,9 +56,7 @@ END_PROVIDER
enddo enddo
enddo enddo
! Deactivated temporarily: bug in N_csf do_csf = s2_eig .and. only_expected_s2 .and. csf_based
! do_csf = s2_eig .and. only_expected_s2 .and. (expected_s2 == 0.d0)
do_csf = .False.
if (diag_algorithm == "Davidson") then if (diag_algorithm == "Davidson") then

View File

@ -46,7 +46,7 @@ subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint,N_st,sze)
do i=1,N_st do i=1,N_st
norm = u_dot_u(u_0(1,i),n) norm = u_dot_u(u_0(1,i),n)
if (norm /= 0.d0) then if (norm /= 0.d0) then
e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n) e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n) / dsqrt(norm)
else else
e_0(i) = 0.d0 e_0(i) = 0.d0
endif endif

View File

@ -75,8 +75,8 @@ subroutine u_0_HS2_u_0(e_0,s_0,u_0,n,keys_tmp,Nint,N_st,sze)
do i=1,N_st do i=1,N_st
norm = u_dot_u(u_0(1,i),n) norm = u_dot_u(u_0(1,i),n)
if (norm /= 0.d0) then if (norm /= 0.d0) then
e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n) e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/norm
s_0(i) = u_dot_v(s_vec(1,i),u_0(1,i),n) s_0(i) = u_dot_v(s_vec(1,i),u_0(1,i),n)/norm
else else
e_0(i) = 0.d0 e_0(i) = 0.d0
s_0(i) = 0.d0 s_0(i) = 0.d0

View File

@ -60,9 +60,11 @@ END_PROVIDER
CI_eigenvectors_dressed(i,j) = psi_coef(i,j) CI_eigenvectors_dressed(i,j) = psi_coef(i,j)
enddo enddo
enddo enddo
logical :: converged
converged = .False.
call davidson_diag_HS2(psi_det,CI_eigenvectors_dressed, CI_eigenvectors_s2_dressed,& call davidson_diag_HS2(psi_det,CI_eigenvectors_dressed, CI_eigenvectors_s2_dressed,&
size(CI_eigenvectors_dressed,1), CI_electronic_energy_dressed,& size(CI_eigenvectors_dressed,1), CI_electronic_energy_dressed,&
N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1) N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,1,converged)
else if (diag_algorithm == "Lapack") then else if (diag_algorithm == "Lapack") then
@ -156,7 +158,8 @@ subroutine diagonalize_CI_dressed
! eigenstates of the CI matrix ! eigenstates of the CI matrix
END_DOC END_DOC
integer :: i,j integer :: i,j
PROVIDE delta_ij ! PROVIDE delta_ij
PROVIDE dressing_column_h
do j=1,N_states do j=1,N_states
do i=1,N_det do i=1,N_det
psi_coef(i,j) = CI_eigenvectors_dressed(i,j) psi_coef(i,j) = CI_eigenvectors_dressed(i,j)

View File

@ -41,7 +41,11 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_alpha_for_dft, (mo_num,mo_num, N_s
elec_alpha_valence(istate) += one_e_dm_mo_alpha_for_dft(i,i,istate) elec_alpha_valence(istate) += one_e_dm_mo_alpha_for_dft(i,i,istate)
enddo enddo
elec_alpha_valence(istate) = elec_alpha_frozen_num/elec_alpha_valence(istate) elec_alpha_valence(istate) = elec_alpha_frozen_num/elec_alpha_valence(istate)
if( dabs(elec_alpha_valence(istate)) .lt.1.d-12)then
one_e_dm_mo_alpha_for_dft = 0.d0
else
one_e_dm_mo_alpha_for_dft(:,:,istate) = one_e_dm_mo_alpha_for_dft(:,:,istate) * elec_alpha_valence(istate) one_e_dm_mo_alpha_for_dft(:,:,istate) = one_e_dm_mo_alpha_for_dft(:,:,istate) * elec_alpha_valence(istate)
endif
enddo enddo
endif endif
@ -55,6 +59,7 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st
! density matrix for beta electrons in the MO basis used for all DFT calculations based on the density ! density matrix for beta electrons in the MO basis used for all DFT calculations based on the density
END_DOC END_DOC
double precision :: delta_beta(mo_num,mo_num,N_states) double precision :: delta_beta(mo_num,mo_num,N_states)
one_e_dm_mo_beta_for_dft = 0.d0
if(density_for_dft .EQ. "damping_rs_dft")then if(density_for_dft .EQ. "damping_rs_dft")then
delta_beta = one_e_dm_mo_beta - data_one_e_dm_beta_mo delta_beta = one_e_dm_mo_beta - data_one_e_dm_beta_mo
one_e_dm_mo_beta_for_dft = data_one_e_dm_beta_mo + damping_for_rs_dft * delta_beta one_e_dm_mo_beta_for_dft = data_one_e_dm_beta_mo + damping_for_rs_dft * delta_beta
@ -73,6 +78,7 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st
one_e_dm_mo_beta_for_dft(:,:,1) = one_e_dm_mo_beta_average(:,:) one_e_dm_mo_beta_for_dft(:,:,1) = one_e_dm_mo_beta_average(:,:)
endif endif
if(no_core_density)then if(no_core_density)then
integer :: ii,i,j integer :: ii,i,j
do ii = 1, n_core_orb do ii = 1, n_core_orb
@ -82,17 +88,21 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st
one_e_dm_mo_beta_for_dft(i,j,:) = 0.d0 one_e_dm_mo_beta_for_dft(i,j,:) = 0.d0
enddo enddo
enddo enddo
if(normalize_dm)then
double precision :: elec_beta_valence(N_states),elec_beta_frozen_num double precision :: elec_beta_valence(N_states),elec_beta_frozen_num
integer :: istate
if(normalize_dm)then
elec_beta_frozen_num = elec_beta_num - n_core_orb elec_beta_frozen_num = elec_beta_num - n_core_orb
elec_beta_valence = 0.d0 elec_beta_valence = 0.d0
integer :: istate
do istate = 1, N_states do istate = 1, N_states
do i = 1, mo_num do i = 1, mo_num
elec_beta_valence(istate) += one_e_dm_mo_beta_for_dft(i,i,istate) elec_beta_valence(istate) += one_e_dm_mo_beta_for_dft(i,i,istate)
enddo enddo
if(dabs(elec_beta_valence(istate)).lt.1.d-12)then
one_e_dm_mo_beta_for_dft = 0.d0
else
elec_beta_valence(istate) = elec_beta_frozen_num/elec_beta_valence(istate) elec_beta_valence(istate) = elec_beta_frozen_num/elec_beta_valence(istate)
one_e_dm_mo_beta_for_dft(:,:,istate) = one_e_dm_mo_beta_for_dft(:,:,istate) * elec_beta_valence(istate) one_e_dm_mo_beta_for_dft(:,:,istate) = one_e_dm_mo_beta_for_dft(:,:,istate) * elec_beta_valence(istate)
endif
enddo enddo
endif endif
endif endif

View File

@ -262,17 +262,48 @@ subroutine set_natural_mos
iorb = list_virt(i) iorb = list_virt(i)
do j = 1, n_core_inact_act_orb do j = 1, n_core_inact_act_orb
jorb = list_core_inact_act(j) jorb = list_core_inact_act(j)
if(one_e_dm_mo(iorb,jorb).ne. 0.d0)then
print*,'AHAHAH'
print*,iorb,jorb,one_e_dm_mo(iorb,jorb)
stop
endif
enddo enddo
enddo enddo
call mo_as_svd_vectors_of_mo_matrix_eig(one_e_dm_mo,size(one_e_dm_mo,1),mo_num,mo_num,mo_occ,label) call mo_as_svd_vectors_of_mo_matrix_eig(one_e_dm_mo,size(one_e_dm_mo,1),mo_num,mo_num,mo_occ,label)
soft_touch mo_occ soft_touch mo_occ
end end
subroutine set_natorb_no_ov_rot
implicit none
BEGIN_DOC
! Set natural orbitals, obtained by diagonalization of the one-body density matrix
! in the |MO| basis
END_DOC
character*(64) :: label
double precision, allocatable :: tmp(:,:)
allocate(tmp(mo_num, mo_num))
label = "Natural"
tmp = one_e_dm_mo
integer :: i,j,iorb,jorb
do i = 1, n_virt_orb
iorb = list_virt(i)
do j = 1, n_core_inact_act_orb
jorb = list_core_inact_act(j)
tmp(iorb, jorb) = 0.d0
tmp(jorb, iorb) = 0.d0
enddo
enddo
call mo_as_svd_vectors_of_mo_matrix_eig(tmp,size(tmp,1),mo_num,mo_num,mo_occ,label)
soft_touch mo_occ
end
subroutine save_natural_mos_no_ov_rot
implicit none
BEGIN_DOC
! Save natural orbitals, obtained by diagonalization of the one-body density matrix in
! the |MO| basis
END_DOC
call set_natorb_no_ov_rot
call nullify_small_elements(ao_num,mo_num,mo_coef,size(mo_coef,1),1.d-10)
call orthonormalize_mos
call save_mos
end
subroutine save_natural_mos subroutine save_natural_mos
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -384,6 +415,14 @@ END_PROVIDER
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, one_e_dm_ao, (ao_num, ao_num)]
implicit none
BEGIN_DOC
! one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta
END_DOC
one_e_dm_ao = one_e_dm_ao_alpha + one_e_dm_ao_beta
END_PROVIDER
subroutine get_occupation_from_dets(istate,occupation) subroutine get_occupation_from_dets(istate,occupation)
implicit none implicit none

View File

@ -0,0 +1,65 @@
BEGIN_PROVIDER [double precision, z_dipole_moment, (N_states)]
&BEGIN_PROVIDER [double precision, y_dipole_moment, (N_states)]
&BEGIN_PROVIDER [double precision, x_dipole_moment, (N_states)]
implicit none
BEGIN_DOC
! blablabla
END_DOC
integer :: ipoint,istate,i,j
double precision :: weight, r(3)
double precision :: cpu0,cpu1,nuclei_part_z,nuclei_part_y,nuclei_part_x
call cpu_time(cpu0)
z_dipole_moment = 0.d0
y_dipole_moment = 0.d0
x_dipole_moment = 0.d0
do istate = 1, N_states
do i = 1, mo_num
z_dipole_moment(istate) += -(one_e_dm_mo_alpha(i,i,istate)+one_e_dm_mo_beta(i,i,istate)) * mo_dipole_z(i,i)
y_dipole_moment(istate) += -(one_e_dm_mo_alpha(i,i,istate)+one_e_dm_mo_beta(i,i,istate)) * mo_dipole_y(i,i)
x_dipole_moment(istate) += -(one_e_dm_mo_alpha(i,i,istate)+one_e_dm_mo_beta(i,i,istate)) * mo_dipole_x(i,i)
do j = i+1, mo_num
z_dipole_moment(istate) += - 2.d0 * (one_e_dm_mo_alpha(j,i,istate)+one_e_dm_mo_beta(j,i,istate)) * mo_dipole_z(j,i)
y_dipole_moment(istate) += - 2.d0 * (one_e_dm_mo_alpha(j,i,istate)+one_e_dm_mo_beta(j,i,istate)) * mo_dipole_y(j,i)
x_dipole_moment(istate) += - 2.d0 * (one_e_dm_mo_alpha(j,i,istate)+one_e_dm_mo_beta(j,i,istate)) * mo_dipole_x(j,i)
enddo
enddo
enddo
print*,'electron part for z_dipole = ',z_dipole_moment
print*,'electron part for y_dipole = ',y_dipole_moment
print*,'electron part for x_dipole = ',x_dipole_moment
nuclei_part_z = 0.d0
nuclei_part_y = 0.d0
nuclei_part_x = 0.d0
do i = 1,nucl_num
nuclei_part_z += nucl_charge(i) * nucl_coord(i,3)
nuclei_part_y += nucl_charge(i) * nucl_coord(i,2)
nuclei_part_x += nucl_charge(i) * nucl_coord(i,1)
enddo
print*,'nuclei part for z_dipole = ',nuclei_part_z
print*,'nuclei part for y_dipole = ',nuclei_part_y
print*,'nuclei part for x_dipole = ',nuclei_part_x
do istate = 1, N_states
z_dipole_moment(istate) += nuclei_part_z
y_dipole_moment(istate) += nuclei_part_y
x_dipole_moment(istate) += nuclei_part_x
enddo
call cpu_time(cpu1)
print*,'Time to provide the dipole moment :',cpu1-cpu0
END_PROVIDER
subroutine print_z_dipole_moment_only
implicit none
print*, ''
print*, ''
print*, '****************************************'
print*, 'z_dipole_moment = ',z_dipole_moment
print*, '****************************************'
end

View File

@ -27,3 +27,33 @@
psi_energy_h_core(i) = psi_energy_h_core(i) * accu psi_energy_h_core(i) = psi_energy_h_core(i) * accu
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, v_ne_psi_energy, (N_states) ]
implicit none
integer :: i
integer :: j,k
double precision :: tmp(mo_num,mo_num),mono_ints(mo_num,mo_num)
BEGIN_DOC
! v_ne_psi_energy = $\langle \Psi | v_ne |\Psi \rangle$
!
! computed using the :c:data:`one_e_dm_mo_alpha` +
! :c:data:`one_e_dm_mo_beta` and :c:data:`mo_one_e_integrals`
END_DOC
v_ne_psi_energy = 0.d0
do i = 1, N_states
do j = 1, mo_num
do k = 1, mo_num
v_ne_psi_energy(i) += mo_integrals_n_e(k,j) * (one_e_dm_mo_alpha(k,j,i) + one_e_dm_mo_beta(k,j,i))
enddo
enddo
enddo
double precision :: accu
do i = 1, N_states
accu = 0.d0
do j = 1, mo_num
accu += one_e_dm_mo_alpha(j,j,i) + one_e_dm_mo_beta(j,j,i)
enddo
accu = (elec_alpha_num + elec_beta_num ) / accu
v_ne_psi_energy(i) = v_ne_psi_energy(i) * accu
enddo
END_PROVIDER

View File

@ -37,13 +37,15 @@ double precision function g0_UEG_mu_inf(rho_a,rho_b)
rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19 rs = (3d0 / (4d0*pi*rho))**(1d0/3d0) ! JT: serious bug fixed 20/03/19
x = -d2*rs x = -d2*rs
if(dabs(x).lt.50.d0)then if(dabs(x).lt.50.d0)then
g0_UEG_mu_inf= 0.5d0 * (1d0- B*rs + C*rs**2 + D*rs**3 + E*rs**4)*dexp(x) ! g0_UEG_mu_inf= 0.5d0 * (1d0- B*rs + C*rs**2 + D*rs**3 + E*rs**4)*dexp(x)
g0_UEG_mu_inf= 0.5d0 * (1d0+ rs* (-B + rs*(C + rs*(D + rs*E))))*dexp(x)
else else
g0_UEG_mu_inf= 0.d0 g0_UEG_mu_inf= 0.d0
endif endif
else else
g0_UEG_mu_inf= 0.d0 g0_UEG_mu_inf= 0.d0
endif endif
g0_UEG_mu_inf = max(g0_UEG_mu_inf,1.d-14)
end end

View File

@ -120,43 +120,41 @@
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_bis, (n_points_final_grid,ao_num,3)]
BEGIN_PROVIDER[double precision, aos_lapl_in_r_array_transp, (n_points_final_grid,ao_num,3)]
implicit none
!
! aos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplacian of jth ao on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z
integer :: i,j,m
do m = 1, 3
do i = 1, n_points_final_grid
do j = 1, ao_num
aos_lapl_in_r_array_transp(i,j,m) = aos_lapl_in_r_array(j,i,m)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom, (ao_num,n_pts_max_per_atom,nucl_num)]
&BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom_transp, (n_pts_max_per_atom,ao_num,nucl_num)]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! aos_in_r_array_per_atom(i,j,k) = value of the ith ao on the jth grid point attached on the kth atom ! Transposed gradients
!
END_DOC END_DOC
integer :: i,j,k integer :: i,j,m
double precision :: aos_array(ao_num), r(3) double precision :: aos_array(ao_num), r(3)
do k = 1, nucl_num double precision :: aos_grad_array(3,ao_num)
do i = 1, n_pts_per_atom(k) do m = 1, 3
r(1) = final_grid_points_per_atom(1,i,k)
r(2) = final_grid_points_per_atom(2,i,k)
r(3) = final_grid_points_per_atom(3,i,k)
call give_all_aos_at_r(r,aos_array)
do j = 1, ao_num do j = 1, ao_num
aos_in_r_array_per_atom(j,i,k) = aos_array(j) do i = 1, n_points_final_grid
aos_in_r_array_per_atom_transp(i,j,k) = aos_array(j) aos_grad_in_r_array_transp_bis(i,j,m) = aos_grad_in_r_array(j,i,m)
enddo enddo
enddo enddo
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_3, (3,n_points_final_grid,ao_num)]
implicit none
BEGIN_DOC
! Transposed gradients
!
END_DOC
integer :: i,j,m
double precision :: aos_array(ao_num), r(3)
double precision :: aos_grad_array(3,ao_num)
do m = 1, 3
do j = 1, ao_num
do i = 1, n_points_final_grid
aos_grad_in_r_array_transp_3(m,i,j) = aos_grad_in_r_array(j,i,m)
enddo
enddo
enddo
END_PROVIDER

View File

@ -110,6 +110,62 @@ end
grad_dm_b *= 2.d0 grad_dm_b *= 2.d0
end end
subroutine density_and_grad_alpha_beta(r,dm_a,dm_b, grad_dm_a, grad_dm_b)
implicit none
BEGIN_DOC
! input:
!
! * r(1) ==> r(1) = x, r(2) = y, r(3) = z
!
! output:
!
! * dm_a = alpha density evaluated at r
! * dm_b = beta density evaluated at r
! * grad_dm_a(1) = X gradient of the alpha density evaluated in r
! * grad_dm_a(1) = X gradient of the beta density evaluated in r
!
END_DOC
double precision, intent(in) :: r(3)
double precision, intent(out) :: dm_a(N_states),dm_b(N_states)
double precision, intent(out) :: grad_dm_a(3,N_states),grad_dm_b(3,N_states)
double precision :: grad_aos_array(3,ao_num)
integer :: i,j,istate
double precision :: aos_array(ao_num),aos_array_bis(ao_num),u_dot_v
double precision :: aos_grad_array(ao_num,3), aos_grad_array_bis(ao_num,3)
call give_all_aos_and_grad_at_r(r,aos_array,grad_aos_array)
do i = 1, ao_num
do j = 1, 3
aos_grad_array(i,j) = grad_aos_array(j,i)
enddo
enddo
do istate = 1, N_states
! alpha density
! aos_array_bis = \rho_ao * aos_array
call dsymv('U',ao_num,1.d0,one_e_dm_alpha_ao_for_dft(1,1,istate),size(one_e_dm_alpha_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
dm_a(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i)
grad_dm_a(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
grad_dm_a(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
grad_dm_a(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
! aos_grad_array_bis = \rho_ao * aos_grad_array
! beta density
call dsymv('U',ao_num,1.d0,one_e_dm_beta_ao_for_dft(1,1,istate),size(one_e_dm_beta_ao_for_dft,1),aos_array,1,0.d0,aos_array_bis,1)
dm_b(istate) = u_dot_v(aos_array,aos_array_bis,ao_num)
! grad_dm(1) = \sum_i aos_grad_array(i,1) * aos_array_bis(i)
grad_dm_b(1,istate) = u_dot_v(aos_grad_array(1,1),aos_array_bis,ao_num)
grad_dm_b(2,istate) = u_dot_v(aos_grad_array(1,2),aos_array_bis,ao_num)
grad_dm_b(3,istate) = u_dot_v(aos_grad_array(1,3),aos_array_bis,ao_num)
! aos_grad_array_bis = \rho_ao * aos_grad_array
enddo
grad_dm_a *= 2.d0
grad_dm_b *= 2.d0
end
subroutine density_and_grad_lapl_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, grad_dm_a, grad_dm_b, lapl_dm_a, lapl_dm_b, aos_array, grad_aos_array, lapl_aos_array) subroutine density_and_grad_lapl_alpha_beta_and_all_aos_and_grad_aos_at_r(r,dm_a,dm_b, grad_dm_a, grad_dm_b, lapl_dm_a, lapl_dm_b, aos_array, grad_aos_array, lapl_aos_array)

View File

@ -79,7 +79,7 @@
END_DOC END_DOC
integer :: m integer :: m
integer :: i,j integer :: i,j
mos_grad_in_r_array = 0.d0 mos_grad_in_r_array_tranp = 0.d0
do i = 1, n_points_final_grid do i = 1, n_points_final_grid
do j = 1, mo_num do j = 1, mo_num
do m = 1, 3 do m = 1, 3
@ -89,6 +89,24 @@
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER[double precision, mos_grad_in_r_array_transp_bis, (n_points_final_grid,mo_num,3)]
implicit none
BEGIN_DOC
! Transposed gradients
!
END_DOC
integer :: i,j,m
do m = 1, 3
do j = 1, mo_num
do i = 1, n_points_final_grid
mos_grad_in_r_array_transp_bis(i,j,m) = mos_grad_in_r_array(j,i,m)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, alpha_dens_kin_in_r, (n_points_final_grid)] BEGIN_PROVIDER [double precision, alpha_dens_kin_in_r, (n_points_final_grid)]
&BEGIN_PROVIDER [double precision, beta_dens_kin_in_r, (n_points_final_grid)] &BEGIN_PROVIDER [double precision, beta_dens_kin_in_r, (n_points_final_grid)]
implicit none implicit none
@ -115,8 +133,6 @@
BEGIN_DOC BEGIN_DOC
! mos_lapl_in_r_array(i,j,k) = value of the kth component of the laplacian of ith mo on the jth grid point ! mos_lapl_in_r_array(i,j,k) = value of the kth component of the laplacian of ith mo on the jth grid point
! !
! mos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplacian of jth mo on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z ! k = 1 : x, k= 2, y, k 3, z
END_DOC END_DOC
integer :: m integer :: m
@ -127,3 +143,41 @@
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER[double precision, mos_lapl_in_r_array_tranp,(3,mo_num,n_points_final_grid)]
implicit none
BEGIN_DOC
! mos_lapl_in_r_array_transp(i,j,k) = value of the kth component of the laplient of jth mo on the ith grid point
!
! k = 1 : x, k= 2, y, k 3, z
END_DOC
integer :: m
integer :: i,j
mos_lapl_in_r_array_tranp = 0.d0
do i = 1, n_points_final_grid
do j = 1, mo_num
do m = 1, 3
mos_lapl_in_r_array_tranp(m,j,i) = mos_lapl_in_r_array(j,i,m)
enddo
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

View File

@ -3,6 +3,7 @@
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! \Delta_{state-specific}. \Psi ! \Delta_{state-specific}. \Psi
! Diagonal element is divided by 2 because Delta = D + D^t
END_DOC END_DOC
integer :: i,ii,k,j, l integer :: i,ii,k,j, l

View File

@ -1,3 +1,4 @@
cipsi cipsi
davidson_undressed
selectors_full selectors_full
generators_full generators_full

View File

@ -12,18 +12,18 @@ subroutine give_all_mos_and_grad_at_r(r,mos_array,mos_grad_array)
implicit none implicit none
double precision, intent(in) :: r(3) double precision, intent(in) :: r(3)
double precision, intent(out) :: mos_array(mo_num) double precision, intent(out) :: mos_array(mo_num)
double precision, intent(out) :: mos_grad_array(mo_num,3) double precision, intent(out) :: mos_grad_array(3,mo_num)
integer :: i,j,k integer :: i,j,k
double precision :: aos_array(ao_num),aos_grad_array(ao_num,3) double precision :: aos_array(ao_num),aos_grad_array(3,ao_num)
call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array) call give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
mos_array=0d0 mos_array=0d0
mos_grad_array=0d0 mos_grad_array=0d0
do j = 1, mo_num do j = 1, mo_num
do k=1, ao_num do k=1, ao_num
mos_array(j) += mo_coef(k,j)*aos_array(k) mos_array(j) += mo_coef(k,j)*aos_array(k)
mos_grad_array(j,1) += mo_coef(k,j)*aos_grad_array(k,1) mos_grad_array(1,j) += mo_coef(k,j)*aos_grad_array(1,k)
mos_grad_array(j,2) += mo_coef(k,j)*aos_grad_array(k,2) mos_grad_array(2,j) += mo_coef(k,j)*aos_grad_array(2,k)
mos_grad_array(j,3) += mo_coef(k,j)*aos_grad_array(k,3) mos_grad_array(3,j) += mo_coef(k,j)*aos_grad_array(3,k)
enddo enddo
enddo enddo
end end
@ -33,9 +33,9 @@ subroutine give_all_mos_and_grad_and_lapl_at_r(r,mos_array,mos_grad_array,mos_la
implicit none implicit none
double precision, intent(in) :: r(3) double precision, intent(in) :: r(3)
double precision, intent(out) :: mos_array(mo_num) double precision, intent(out) :: mos_array(mo_num)
double precision, intent(out) :: mos_grad_array(mo_num,3),mos_lapl_array(mo_num,3) double precision, intent(out) :: mos_grad_array(3,mo_num),mos_lapl_array(3,mo_num)
integer :: i,j,k integer :: i,j,k
double precision :: aos_array(ao_num),aos_grad_array(ao_num,3),aos_lapl_array(ao_num,3) double precision :: aos_array(ao_num),aos_grad_array(3,ao_num),aos_lapl_array(3,ao_num)
call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array) call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
mos_array = 0.d0 mos_array = 0.d0
mos_grad_array = 0.d0 mos_grad_array = 0.d0
@ -43,12 +43,12 @@ subroutine give_all_mos_and_grad_and_lapl_at_r(r,mos_array,mos_grad_array,mos_la
do j = 1, mo_num do j = 1, mo_num
do k=1, ao_num do k=1, ao_num
mos_array(j) += mo_coef(k,j) * aos_array(k) mos_array(j) += mo_coef(k,j) * aos_array(k)
mos_grad_array(j,1) += mo_coef(k,j) * aos_grad_array(k,1) mos_grad_array(1,j) += mo_coef(k,j) * aos_grad_array(1,k)
mos_grad_array(j,2) += mo_coef(k,j) * aos_grad_array(k,2) mos_grad_array(2,j) += mo_coef(k,j) * aos_grad_array(2,k)
mos_grad_array(j,3) += mo_coef(k,j) * aos_grad_array(k,3) mos_grad_array(3,j) += mo_coef(k,j) * aos_grad_array(3,k)
mos_lapl_array(j,1) += mo_coef(k,j) * aos_lapl_array(k,1) mos_lapl_array(1,j) += mo_coef(k,j) * aos_lapl_array(1,k)
mos_lapl_array(j,2) += mo_coef(k,j) * aos_lapl_array(k,2) mos_lapl_array(2,j) += mo_coef(k,j) * aos_lapl_array(2,k)
mos_lapl_array(j,3) += mo_coef(k,j) * aos_lapl_array(k,3) mos_lapl_array(3,j) += mo_coef(k,j) * aos_lapl_array(3,k)
enddo enddo
enddo enddo
end end

View File

@ -10,5 +10,5 @@ subroutine hcore_guess
size(mo_one_e_integrals,2),label,1,.false.) size(mo_one_e_integrals,2),label,1,.false.)
call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12 ) call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12 )
call save_mos call save_mos
SOFT_TOUCH mo_coef mo_label TOUCH mo_coef mo_label
end end

View File

@ -1,4 +1,3 @@
BEGIN_PROVIDER [ double precision, mo_overlap,(mo_num,mo_num) ] BEGIN_PROVIDER [ double precision, mo_overlap,(mo_num,mo_num) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC

View File

@ -7,7 +7,7 @@ BEGIN_PROVIDER [double precision, core_energy_erf]
core_energy_erf = 0.d0 core_energy_erf = 0.d0
do i = 1, n_core_orb do i = 1, n_core_orb
j = list_core(i) j = list_core(i)
core_energy_erf += 2.d0 * mo_one_e_integrals(j,j) + mo_two_e_int_erf_jj(j,j) core_energy_erf += mo_two_e_int_erf_jj(j,j)
do k = i+1, n_core_orb do k = i+1, n_core_orb
l = list_core(k) l = list_core(k)
core_energy_erf += 2.d0 * (2.d0 * mo_two_e_int_erf_jj(j,l) - mo_two_e_int_erf_jj_exchange(j,l)) core_energy_erf += 2.d0 * (2.d0 * mo_two_e_int_erf_jj(j,l) - mo_two_e_int_erf_jj_exchange(j,l))

View File

@ -51,6 +51,13 @@ end
subroutine four_idx_novvvv 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 use map_module
implicit none implicit none
BEGIN_DOC BEGIN_DOC

View File

@ -296,6 +296,8 @@ end
! If true, the excitation is banned in the selection. Useful with local MOs. ! If true, the excitation is banned in the selection. Useful with local MOs.
END_DOC END_DOC
banned_excitation = .False. banned_excitation = .False.
use_banned_excitation = .False.
integer :: i,j, icount integer :: i,j, icount
integer(key_kind) :: idx integer(key_kind) :: idx
double precision :: tmp double precision :: tmp

View File

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

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

View File

@ -6,7 +6,7 @@ size: (becke_numerical_grid.n_points_final_grid,determinants.n_states)
[mu_of_r_potential] [mu_of_r_potential]
type: character*(32) 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 interface: ezfio, provider, ocaml
default: hf default: hf

View File

@ -76,7 +76,11 @@ BEGIN_PROVIDER [integer, n_basis_orb]
! !
! It corresponds to all MOs except those defined as "deleted" ! It corresponds to all MOs except those defined as "deleted"
END_DOC END_DOC
if(mu_of_r_potential == "pure_act")then
n_basis_orb = n_act_orb
else
n_basis_orb = n_all_but_del_orb n_basis_orb = n_all_but_del_orb
endif
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [integer, list_basis, (n_basis_orb)] 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" ! It corresponds to all MOs except those defined as "deleted"
END_DOC END_DOC
integer :: i integer :: i
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 do i = 1, n_all_but_del_orb
list_basis(i) = list_all_but_del_orb(i) list_basis(i) = list_all_but_del_orb(i)
enddo enddo
endif
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [double precision, basis_mos_in_r_array, (n_basis_orb,n_points_final_grid)] BEGIN_PROVIDER [double precision, basis_mos_in_r_array, (n_basis_orb,n_points_final_grid)]

View File

@ -9,6 +9,7 @@ BEGIN_PROVIDER [double precision, two_e_int_hf_f, (n_basis_orb,n_basis_orb,n_max
END_DOC END_DOC
integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
double precision :: get_two_e_integral double precision :: get_two_e_integral
PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
do orb_m = 1, n_max_occ_val_orb_for_hf! electron 1 do orb_m = 1, n_max_occ_val_orb_for_hf! electron 1
m = list_valence_orb_for_hf(orb_m,1) m = list_valence_orb_for_hf(orb_m,1)
do orb_n = 1, n_max_occ_val_orb_for_hf! electron 2 do orb_n = 1, n_max_occ_val_orb_for_hf! electron 2

View File

@ -235,6 +235,7 @@ BEGIN_PROVIDER [double precision, two_e_int_aa_f, (n_basis_orb,n_basis_orb,n_act
END_DOC END_DOC
integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
double precision :: integrals_array(mo_num,mo_num),get_two_e_integral double precision :: integrals_array(mo_num,mo_num),get_two_e_integral
PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
do orb_m = 1, n_act_orb ! electron 1 do orb_m = 1, n_act_orb ! electron 1
m = list_act(orb_m) m = list_act(orb_m)
do orb_n = 1, n_act_orb ! electron 2 do orb_n = 1, n_act_orb ! electron 2
@ -264,6 +265,7 @@ BEGIN_PROVIDER [double precision, two_e_int_ia_f, (n_basis_orb,n_basis_orb,n_ina
END_DOC END_DOC
integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
double precision :: integrals_array(mo_num,mo_num),get_two_e_integral double precision :: integrals_array(mo_num,mo_num),get_two_e_integral
PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
do orb_m = 1, n_act_orb ! electron 1 do orb_m = 1, n_act_orb ! electron 1
m = list_act(orb_m) m = list_act(orb_m)
do orb_n = 1, n_inact_orb ! electron 2 do orb_n = 1, n_inact_orb ! electron 2
@ -293,6 +295,7 @@ BEGIN_PROVIDER [double precision, two_e_int_ii_f, (n_basis_orb,n_basis_orb,n_ina
END_DOC END_DOC
integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
double precision :: get_two_e_integral,integrals_array(mo_num,mo_num) double precision :: get_two_e_integral,integrals_array(mo_num,mo_num)
PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
do orb_m = 1, n_inact_orb ! electron 1 do orb_m = 1, n_inact_orb ! electron 1
m = list_inact(orb_m) m = list_inact(orb_m)
do orb_n = 1, n_inact_orb ! electron 2 do orb_n = 1, n_inact_orb ! electron 2

View File

@ -26,7 +26,7 @@
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
if(mu_of_r_potential.EQ."hf")then if(mu_of_r_potential.EQ."hf")then
mu_of_r_prov(ipoint,istate) = mu_of_r_hf(ipoint) 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) mu_of_r_prov(ipoint,istate) = mu_of_r_psi_cas(ipoint,istate)
else else
print*,'you requested the following mu_of_r_potential' print*,'you requested the following mu_of_r_potential'
@ -148,3 +148,4 @@
mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate) mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate)
enddo enddo
END_PROVIDER END_PROVIDER

View File

@ -50,3 +50,10 @@ type: logical
doc: If true, leave untouched all the orbitals defined as core and optimize all the orbitals defined as active with qp_set_mo_class doc: If true, leave untouched all the orbitals defined as core and optimize all the orbitals defined as active with qp_set_mo_class
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
[no_oa_or_av_opt]
type: logical
doc: If true, you set to zero all Fock elements between the orbital set to active and all the other orbitals
interface: ezfio,provider,ocaml
default: False

View File

@ -31,6 +31,27 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num)
enddo enddo
enddo enddo
endif endif
if(no_oa_or_av_opt)then
do i = 1, n_act_orb
iorb = list_act(i)
do j = 1, n_inact_orb
jorb = list_inact(j)
F(iorb,jorb) = 0.d0
F(jorb,iorb) = 0.d0
enddo
do j = 1, n_virt_orb
jorb = list_virt(j)
F(iorb,jorb) = 0.d0
F(jorb,iorb) = 0.d0
enddo
do j = 1, n_core_orb
jorb = list_core(j)
F(iorb,jorb) = 0.d0
F(jorb,iorb) = 0.d0
enddo
enddo
endif
! Insert level shift here ! Insert level shift here
do i = elec_beta_num+1, elec_alpha_num do i = elec_beta_num+1, elec_alpha_num

View File

@ -92,6 +92,27 @@
enddo enddo
endif endif
if(no_oa_or_av_opt)then
do i = 1, n_act_orb
iorb = list_act(i)
do j = 1, n_inact_orb
jorb = list_inact(j)
Fock_matrix_mo(iorb,jorb) = 0.d0
Fock_matrix_mo(jorb,iorb) = 0.d0
enddo
do j = 1, n_virt_orb
jorb = list_virt(j)
Fock_matrix_mo(iorb,jorb) = 0.d0
Fock_matrix_mo(jorb,iorb) = 0.d0
enddo
do j = 1, n_core_orb
jorb = list_core(j)
Fock_matrix_mo(iorb,jorb) = 0.d0
Fock_matrix_mo(jorb,iorb) = 0.d0
enddo
enddo
endif
END_PROVIDER END_PROVIDER

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

View File

@ -0,0 +1,3 @@
program hcore_guess_prog
call hcore_guess
end

View File

@ -0,0 +1,5 @@
program pouet
implicit none
call huckel_guess
end

View File

@ -0,0 +1,5 @@
program print_dipole
implicit none
call print_z_dipole_moment_only
end

View File

@ -0,0 +1,11 @@
program print_var_energy
implicit none
read_wf = .True.
touch read_wf
call routine
end
subroutine routine
implicit none
print*,'psi_energy = ',psi_energy
end

View File

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

View File

@ -1,45 +1,21 @@
[two_rdm_ab_disk]
type: double precision
doc: active part of the two body rdm alpha/beta stored on disk
interface: ezfio
size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
[io_two_body_rdm_ab] [io_two_body_rdm_ab]
type: Disk_access type: Disk_access
doc: Read/Write the active part of the two-body rdm for alpha/beta electrons from/to disk [ Write | Read | None ] doc: Read/Write the active part of the two-body rdm for alpha/beta electrons from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: None default: None
[two_rdm_aa_disk]
type: double precision
doc: active part of the two body rdm alpha/alpha stored on disk
interface: ezfio
size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
[io_two_body_rdm_aa] [io_two_body_rdm_aa]
type: Disk_access type: Disk_access
doc: Read/Write the active part of the two-body rdm for alpha/alpha electrons from/to disk [ Write | Read | None ] doc: Read/Write the active part of the two-body rdm for alpha/alpha electrons from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: None default: None
[two_rdm_bb_disk]
type: double precision
doc: active part of the two body rdm beta/beta stored on disk
interface: ezfio
size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
[io_two_body_rdm_bb] [io_two_body_rdm_bb]
type: Disk_access type: Disk_access
doc: Read/Write the active part of the two-body rdm for beta/beta electrons from/to disk [ Write | Read | None ] doc: Read/Write the active part of the two-body rdm for beta/beta electrons from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: None default: None
[two_rdm_spin_trace_disk]
type: double precision
doc: active part of the two body rdm spin trace stored on disk
interface: ezfio
size: (bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,bitmask.n_act_orb,determinants.n_states)
[io_two_body_rdm_spin_trace] [io_two_body_rdm_spin_trace]
type: Disk_access type: Disk_access
doc: Read/Write the active part of the two-body rdm for spin trace electrons from/to disk [ Write | Read | None ] doc: Read/Write the active part of the two-body rdm for spin trace electrons from/to disk [ Write | Read | None ]

View File

@ -22,6 +22,8 @@
END_DOC END_DOC
integer :: ispin integer :: ispin
double precision :: wall_1, wall_2 double precision :: wall_1, wall_2
character*(128) :: name_file
name_file = 'act_2_rdm_ab_mo'
! condition for alpha/beta spin ! condition for alpha/beta spin
print*,'' print*,''
print*,'Providing act_2_rdm_ab_mo ' print*,'Providing act_2_rdm_ab_mo '
@ -31,13 +33,13 @@
call wall_time(wall_1) call wall_time(wall_1)
if(read_two_body_rdm_ab)then if(read_two_body_rdm_ab)then
print*,'Reading act_2_rdm_ab_mo from disk ...' print*,'Reading act_2_rdm_ab_mo from disk ...'
call ezfio_get_two_body_rdm_two_rdm_ab_disk(act_2_rdm_ab_mo) call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_ab_mo,name_file)
else else
call orb_range_2_rdm_openmp(act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call orb_range_2_rdm_openmp(act_2_rdm_ab_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
endif endif
if(write_two_body_rdm_ab)then if(write_two_body_rdm_ab)then
print*,'Writing act_2_rdm_ab_mo on disk ...' print*,'Writing act_2_rdm_ab_mo on disk ...'
call ezfio_set_two_body_rdm_two_rdm_ab_disk(act_2_rdm_ab_mo) call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_ab_mo,name_file)
call ezfio_set_two_body_rdm_io_two_body_rdm_ab("Read") call ezfio_set_two_body_rdm_io_two_body_rdm_ab("Read")
endif endif
call wall_time(wall_2) call wall_time(wall_2)
@ -63,18 +65,20 @@
! condition for alpha/beta spin ! condition for alpha/beta spin
print*,'' print*,''
print*,'Providing act_2_rdm_aa_mo ' print*,'Providing act_2_rdm_aa_mo '
character*(128) :: name_file
name_file = 'act_2_rdm_aa_mo'
ispin = 1 ispin = 1
act_2_rdm_aa_mo = 0.d0 act_2_rdm_aa_mo = 0.d0
call wall_time(wall_1) call wall_time(wall_1)
if(read_two_body_rdm_aa)then if(read_two_body_rdm_aa)then
print*,'Reading act_2_rdm_aa_mo from disk ...' print*,'Reading act_2_rdm_aa_mo from disk ...'
call ezfio_get_two_body_rdm_two_rdm_aa_disk(act_2_rdm_aa_mo) call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_aa_mo,name_file)
else else
call orb_range_2_rdm_openmp(act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call orb_range_2_rdm_openmp(act_2_rdm_aa_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
endif endif
if(write_two_body_rdm_aa)then if(write_two_body_rdm_aa)then
print*,'Writing act_2_rdm_aa_mo on disk ...' print*,'Writing act_2_rdm_aa_mo on disk ...'
call ezfio_set_two_body_rdm_two_rdm_aa_disk(act_2_rdm_aa_mo) call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_aa_mo,name_file)
call ezfio_set_two_body_rdm_io_two_body_rdm_aa("Read") call ezfio_set_two_body_rdm_io_two_body_rdm_aa("Read")
endif endif
@ -101,18 +105,20 @@
! condition for beta/beta spin ! condition for beta/beta spin
print*,'' print*,''
print*,'Providing act_2_rdm_bb_mo ' print*,'Providing act_2_rdm_bb_mo '
character*(128) :: name_file
name_file = 'act_2_rdm_bb_mo'
ispin = 2 ispin = 2
act_2_rdm_bb_mo = 0.d0 act_2_rdm_bb_mo = 0.d0
call wall_time(wall_1) call wall_time(wall_1)
if(read_two_body_rdm_bb)then if(read_two_body_rdm_bb)then
print*,'Reading act_2_rdm_bb_mo from disk ...' print*,'Reading act_2_rdm_bb_mo from disk ...'
call ezfio_get_two_body_rdm_two_rdm_bb_disk(act_2_rdm_bb_mo) call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_bb_mo,name_file)
else else
call orb_range_2_rdm_openmp(act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call orb_range_2_rdm_openmp(act_2_rdm_bb_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
endif endif
if(write_two_body_rdm_bb)then if(write_two_body_rdm_bb)then
print*,'Writing act_2_rdm_bb_mo on disk ...' print*,'Writing act_2_rdm_bb_mo on disk ...'
call ezfio_set_two_body_rdm_two_rdm_bb_disk(act_2_rdm_bb_mo) call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_bb_mo,name_file)
call ezfio_set_two_body_rdm_io_two_body_rdm_bb("Read") call ezfio_set_two_body_rdm_io_two_body_rdm_bb("Read")
endif endif
@ -138,18 +144,20 @@
! condition for beta/beta spin ! condition for beta/beta spin
print*,'' print*,''
print*,'Providing act_2_rdm_spin_trace_mo ' print*,'Providing act_2_rdm_spin_trace_mo '
character*(128) :: name_file
name_file = 'act_2_rdm_spin_trace_mo'
ispin = 4 ispin = 4
act_2_rdm_spin_trace_mo = 0.d0 act_2_rdm_spin_trace_mo = 0.d0
call wall_time(wall_1) call wall_time(wall_1)
if(read_two_body_rdm_spin_trace)then if(read_two_body_rdm_spin_trace)then
print*,'Reading act_2_rdm_spin_trace_mo from disk ...' print*,'Reading act_2_rdm_spin_trace_mo from disk ...'
call ezfio_get_two_body_rdm_two_rdm_spin_trace_disk(act_2_rdm_spin_trace_mo) call read_array_two_rdm(n_act_orb,N_states,act_2_rdm_spin_trace_mo,name_file)
else else
call orb_range_2_rdm_openmp(act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) call orb_range_2_rdm_openmp(act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
endif endif
if(write_two_body_rdm_spin_trace)then if(write_two_body_rdm_spin_trace)then
print*,'Writing act_2_rdm_spin_trace_mo on disk ...' print*,'Writing act_2_rdm_spin_trace_mo on disk ...'
call ezfio_set_two_body_rdm_two_rdm_spin_trace_disk(act_2_rdm_spin_trace_mo) call write_array_two_rdm(n_act_orb,n_states,act_2_rdm_spin_trace_mo,name_file)
call ezfio_set_two_body_rdm_io_two_body_rdm_spin_trace("Read") call ezfio_set_two_body_rdm_io_two_body_rdm_spin_trace("Read")
endif endif

View File

@ -15,7 +15,7 @@ subroutine routine_active_only
double precision :: wee_aa_st_av, rdm_aa_st_av double precision :: wee_aa_st_av, rdm_aa_st_av
double precision :: wee_bb_st_av, rdm_bb_st_av double precision :: wee_bb_st_av, rdm_bb_st_av
double precision :: wee_ab_st_av, rdm_ab_st_av double precision :: wee_ab_st_av, rdm_ab_st_av
double precision :: wee_tot_st_av, rdm_tot_st_av double precision :: wee_tot_st_av, rdm_tot_st_av,spin_trace
double precision :: wee_aa_st_av_2,wee_ab_st_av_2,wee_bb_st_av_2,wee_tot_st_av_2,wee_tot_st_av_3 double precision :: wee_aa_st_av_2,wee_ab_st_av_2,wee_bb_st_av_2,wee_tot_st_av_2,wee_tot_st_av_3
wee_ab = 0.d0 wee_ab = 0.d0
@ -38,6 +38,27 @@ subroutine routine_active_only
provide act_2_rdm_ab_mo act_2_rdm_aa_mo act_2_rdm_bb_mo act_2_rdm_spin_trace_mo provide act_2_rdm_ab_mo act_2_rdm_aa_mo act_2_rdm_bb_mo act_2_rdm_spin_trace_mo
provide state_av_act_2_rdm_ab_mo state_av_act_2_rdm_aa_mo provide state_av_act_2_rdm_ab_mo state_av_act_2_rdm_aa_mo
provide state_av_act_2_rdm_bb_mo state_av_act_2_rdm_spin_trace_mo provide state_av_act_2_rdm_bb_mo state_av_act_2_rdm_spin_trace_mo
i = 1
j = 2
! print*,'testing stuffs'
! istate = 1
! print*,'alpha/beta'
! print*,'',j,i,j,i
! print*,act_2_rdm_ab_mo(j,i,j,i,istate)
! print*,'',i,j,i,j
! print*,act_2_rdm_ab_mo(i,j,i,j,istate)
! print*,'alpha/alpha'
! print*,'',j,i,j,i
! print*,act_2_rdm_aa_mo(j,i,j,i,istate)
! print*,'',i,j,i,j
! print*,act_2_rdm_aa_mo(i,j,i,j,istate)
! print*,'spin_trace'
! print*,'',j,i,j,i
! print*,act_2_rdm_spin_trace_mo(j,i,j,i,istate)
! print*,'',i,j,i,j
! print*,act_2_rdm_spin_trace_mo(i,j,i,j,istate)
! stop
!
print*,'**************************' print*,'**************************'
print*,'**************************' print*,'**************************'
do istate = 1, N_states do istate = 1, N_states
@ -51,6 +72,19 @@ subroutine routine_active_only
korb = list_act(k) korb = list_act(k)
do l = 1, n_act_orb do l = 1, n_act_orb
lorb = list_act(l) lorb = list_act(l)
if(dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(j,i,l,k,istate)).gt.1.d-10)then
print*,'Error in act_2_rdm_spin_trace_mo'
print*,"dabs(act_2_rdm_spin_trace_mo(i,j,k,l) - act_2_rdm_spin_trace_mo(j,i,l,k)).gt.1.d-10"
print*,i,j,k,l
print*,act_2_rdm_spin_trace_mo(i,j,k,l,istate),act_2_rdm_spin_trace_mo(j,i,l,k,istate),dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(j,i,l,k,istate))
endif
if(dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(k,l,i,j,istate)).gt.1.d-10)then
print*,'Error in act_2_rdm_spin_trace_mo'
print*,"dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(k,l,i,j,istate),istate).gt.1.d-10"
print*,i,j,k,l
print*,act_2_rdm_spin_trace_mo(i,j,k,l,istate),act_2_rdm_spin_trace_mo(k,l,i,j,istate),dabs(act_2_rdm_spin_trace_mo(i,j,k,l,istate) - act_2_rdm_spin_trace_mo(k,l,i,j,istate))
endif
vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map)
@ -59,6 +93,17 @@ subroutine routine_active_only
rdmaa = act_2_rdm_aa_mo(l,k,j,i,istate) rdmaa = act_2_rdm_aa_mo(l,k,j,i,istate)
rdmbb = act_2_rdm_bb_mo(l,k,j,i,istate) rdmbb = act_2_rdm_bb_mo(l,k,j,i,istate)
rdmtot = act_2_rdm_spin_trace_mo(l,k,j,i,istate) rdmtot = act_2_rdm_spin_trace_mo(l,k,j,i,istate)
spin_trace = rdmaa + rdmbb + rdmab
if(dabs(rdmtot- spin_trace).gt.1.d-10)then
print*,'Error in non state average !!!!'
print*,l,k,j,i
print*,lorb,korb,jorb,iorb
print*,spin_trace,rdmtot,dabs(spin_trace - rdmtot)
print*,'rdmab,rdmaa,rdmbb'
print*, rdmab,rdmaa,rdmbb
endif
wee_ab(istate) += vijkl * rdmab wee_ab(istate) += vijkl * rdmab
@ -71,8 +116,8 @@ subroutine routine_active_only
enddo enddo
enddo enddo
wee_aa_st_av_2 += wee_aa(istate) * state_average_weight(istate) wee_aa_st_av_2 += wee_aa(istate) * state_average_weight(istate)
wee_bb_st_av_2 += wee_aa(istate) * state_average_weight(istate) wee_bb_st_av_2 += wee_bb(istate) * state_average_weight(istate)
wee_ab_st_av_2 += wee_aa(istate) * state_average_weight(istate) wee_ab_st_av_2 += wee_ab(istate) * state_average_weight(istate)
wee_tot_st_av_2 += wee_tot(istate) * state_average_weight(istate) wee_tot_st_av_2 += wee_tot(istate) * state_average_weight(istate)
wee_tot_st_av_3 += psi_energy_two_e(istate) * state_average_weight(istate) wee_tot_st_av_3 += psi_energy_two_e(istate) * state_average_weight(istate)
print*,'' print*,''
@ -87,7 +132,6 @@ subroutine routine_active_only
print*,'Full energy ' print*,'Full energy '
print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate) print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate)
enddo enddo
wee_aa_st_av = 0.d0 wee_aa_st_av = 0.d0
wee_bb_st_av = 0.d0 wee_bb_st_av = 0.d0
wee_ab_st_av = 0.d0 wee_ab_st_av = 0.d0
@ -103,10 +147,30 @@ subroutine routine_active_only
vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map)
if(dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(j,i,l,k)).gt.1.d-10)then
print*,'Error in state_av_act_2_rdm_spin_trace_mo'
print*,"dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(j,i,l,k)).gt.1.d-10"
print*,i,j,k,l
print*,state_av_act_2_rdm_spin_trace_mo(i,j,k,l),state_av_act_2_rdm_spin_trace_mo(j,i,l,k),dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(j,i,l,k))
endif
if(dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(k,l,i,j)).gt.1.d-10)then
print*,'Error in state_av_act_2_rdm_spin_trace_mo'
print*,"dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(k,l,i,j)).gt.1.d-10"
print*,i,j,k,l
print*,state_av_act_2_rdm_spin_trace_mo(i,j,k,l),state_av_act_2_rdm_spin_trace_mo(k,l,i,j),dabs(state_av_act_2_rdm_spin_trace_mo(i,j,k,l) - state_av_act_2_rdm_spin_trace_mo(k,l,i,j))
endif
rdm_aa_st_av = state_av_act_2_rdm_aa_mo(l,k,j,i) rdm_aa_st_av = state_av_act_2_rdm_aa_mo(l,k,j,i)
rdm_bb_st_av = state_av_act_2_rdm_bb_mo(l,k,j,i) rdm_bb_st_av = state_av_act_2_rdm_bb_mo(l,k,j,i)
rdm_ab_st_av = state_av_act_2_rdm_ab_mo(l,k,j,i) rdm_ab_st_av = state_av_act_2_rdm_ab_mo(l,k,j,i)
spin_trace = rdm_aa_st_av + rdm_bb_st_av + rdm_ab_st_av
rdm_tot_st_av = state_av_act_2_rdm_spin_trace_mo(l,k,j,i) rdm_tot_st_av = state_av_act_2_rdm_spin_trace_mo(l,k,j,i)
if(dabs(spin_trace - rdm_tot_st_av).gt.1.d-10)then
print*,'Error !!!!'
print*,l,k,j,i
print*,spin_trace,rdm_tot_st_av,dabs(spin_trace - rdm_tot_st_av)
endif
wee_aa_st_av += vijkl * rdm_aa_st_av wee_aa_st_av += vijkl * rdm_aa_st_av
wee_bb_st_av += vijkl * rdm_bb_st_av wee_bb_st_av += vijkl * rdm_bb_st_av

View File

@ -0,0 +1,29 @@
subroutine write_array_two_rdm(n_orb,nstates,array_tmp,name_file)
implicit none
integer, intent(in) :: n_orb,nstates
character*(128), intent(in) :: name_file
double precision, intent(in) :: array_tmp(n_orb,n_orb,n_orb,n_orb,nstates)
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'W')
write(i_unit_output)array_tmp
close(unit=i_unit_output)
end
subroutine read_array_two_rdm(n_orb,nstates,array_tmp,name_file)
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
integer, intent(in) :: n_orb,nstates
character*(128), intent(in) :: name_file
double precision, intent(out) :: array_tmp(n_orb,n_orb,n_orb,n_orb,N_states)
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'/work/'//trim(name_file)
i_unit_output = getUnitAndOpen(output,'R')
read(i_unit_output)array_tmp
close(unit=i_unit_output)
end

View File

@ -119,7 +119,11 @@
call wall_time(wall_1) call wall_time(wall_1)
double precision :: wall_1, wall_2 double precision :: wall_1, wall_2
print*,'providing state_av_act_2_rdm_spin_trace_mo ' print*,'providing state_av_act_2_rdm_spin_trace_mo '
call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1)) state_av_act_2_rdm_spin_trace_mo = state_av_act_2_rdm_ab_mo &
+ state_av_act_2_rdm_aa_mo &
+ state_av_act_2_rdm_bb_mo
! call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
call wall_time(wall_2) call wall_time(wall_2)
print*,'Time to provide state_av_act_2_rdm_spin_trace_mo',wall_2 - wall_1 print*,'Time to provide state_av_act_2_rdm_spin_trace_mo',wall_2 - wall_1

View File

@ -61,14 +61,24 @@
! Therefore you don't necessayr have symmetry between electron 1 and 2 ! Therefore you don't necessayr have symmetry between electron 1 and 2
nkeys += 1 nkeys += 1
do istate = 1, N_st do istate = 1, N_st
values(istate,nkeys) = c_1(istate) values(istate,nkeys) = 0.5d0 * c_1(istate)
enddo enddo
keys(1,nkeys) = h1 keys(1,nkeys) = h1
keys(2,nkeys) = h2 keys(2,nkeys) = h2
keys(3,nkeys) = h1 keys(3,nkeys) = h1
keys(4,nkeys) = h2 keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate)
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = h1
enddo enddo
enddo enddo
else if (alpha_alpha)then else if (alpha_alpha)then
do i = 1, n_occ_ab(1) do i = 1, n_occ_ab(1)
i1 = occ(i,1) i1 = occ(i,1)
@ -259,12 +269,20 @@
if(alpha_beta)then if(alpha_beta)then
nkeys += 1 nkeys += 1
do istate = 1, N_st do istate = 1, N_st
values(istate,nkeys) = c_1(istate) * phase values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo enddo
keys(1,nkeys) = h1 keys(1,nkeys) = h1
keys(2,nkeys) = h2 keys(2,nkeys) = h2
keys(3,nkeys) = p1 keys(3,nkeys) = p1
keys(4,nkeys) = p2 keys(4,nkeys) = p2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p2
keys(4,nkeys) = p1
else if(spin_trace)then else if(spin_trace)then
nkeys += 1 nkeys += 1
do istate = 1, N_st do istate = 1, N_st
@ -278,10 +296,10 @@
do istate = 1, N_st do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo enddo
keys(1,nkeys) = p1 keys(1,nkeys) = h2
keys(2,nkeys) = p2 keys(2,nkeys) = h1
keys(3,nkeys) = h1 keys(3,nkeys) = p2
keys(4,nkeys) = h2 keys(4,nkeys) = p1
endif endif
end end
@ -356,12 +374,20 @@
h2 = list_orb_reverse(h2) h2 = list_orb_reverse(h2)
nkeys += 1 nkeys += 1
do istate = 1, N_st do istate = 1, N_st
values(istate,nkeys) = c_1(istate) * phase values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo enddo
keys(1,nkeys) = h1 keys(1,nkeys) = h1
keys(2,nkeys) = h2 keys(2,nkeys) = h2
keys(3,nkeys) = p1 keys(3,nkeys) = p1
keys(4,nkeys) = h2 keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo enddo
else else
! Mono beta ! Mono beta
@ -377,12 +403,20 @@
h2 = list_orb_reverse(h2) h2 = list_orb_reverse(h2)
nkeys += 1 nkeys += 1
do istate = 1, N_st do istate = 1, N_st
values(istate,nkeys) = c_1(istate) * phase values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo enddo
keys(1,nkeys) = h1 keys(1,nkeys) = h1
keys(2,nkeys) = h2 keys(2,nkeys) = h2
keys(3,nkeys) = p1 keys(3,nkeys) = p1
keys(4,nkeys) = h2 keys(4,nkeys) = h2
nkeys += 1
do istate = 1, N_st
values(istate,nkeys) = 0.5d0 * c_1(istate) * phase
enddo
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo enddo
endif endif
else if(spin_trace)then else if(spin_trace)then

View File

@ -60,11 +60,18 @@
! If alpha/beta, electron 1 is alpha, electron 2 is beta ! If alpha/beta, electron 1 is alpha, electron 2 is beta
! Therefore you don't necessayr have symmetry between electron 1 and 2 ! Therefore you don't necessayr have symmetry between electron 1 and 2
nkeys += 1 nkeys += 1
values(nkeys) = c_1 values(nkeys) = 0.5d0 * c_1
keys(1,nkeys) = h1 keys(1,nkeys) = h1
keys(2,nkeys) = h2 keys(2,nkeys) = h2
keys(3,nkeys) = h1 keys(3,nkeys) = h1
keys(4,nkeys) = h2 keys(4,nkeys) = h2
nkeys += 1
values(nkeys) = 0.5d0 * c_1
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = h1
enddo enddo
enddo enddo
else if (alpha_alpha)then else if (alpha_alpha)then
@ -236,11 +243,17 @@
p2 = list_orb_reverse(p2) p2 = list_orb_reverse(p2)
if(alpha_beta)then if(alpha_beta)then
nkeys += 1 nkeys += 1
values(nkeys) = c_1 * phase values(nkeys) = 0.5d0 * c_1 * phase
keys(1,nkeys) = h1 keys(1,nkeys) = h1
keys(2,nkeys) = h2 keys(2,nkeys) = h2
keys(3,nkeys) = p1 keys(3,nkeys) = p1
keys(4,nkeys) = p2 keys(4,nkeys) = p2
nkeys += 1
values(nkeys) = 0.5d0 * c_1 * phase
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = p2
keys(4,nkeys) = p1
else if(spin_trace)then else if(spin_trace)then
nkeys += 1 nkeys += 1
values(nkeys) = 0.5d0 * c_1 * phase values(nkeys) = 0.5d0 * c_1 * phase
@ -250,10 +263,10 @@
keys(4,nkeys) = p2 keys(4,nkeys) = p2
nkeys += 1 nkeys += 1
values(nkeys) = 0.5d0 * c_1 * phase values(nkeys) = 0.5d0 * c_1 * phase
keys(1,nkeys) = p1 keys(1,nkeys) = h2
keys(2,nkeys) = p2 keys(2,nkeys) = h1
keys(3,nkeys) = h1 keys(3,nkeys) = p2
keys(4,nkeys) = h2 keys(4,nkeys) = p1
endif endif
end end
@ -327,11 +340,17 @@
if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2) h2 = list_orb_reverse(h2)
nkeys += 1 nkeys += 1
values(nkeys) = c_1 * phase values(nkeys) = 0.5d0 * c_1 * phase
keys(1,nkeys) = h1 keys(1,nkeys) = h1
keys(2,nkeys) = h2 keys(2,nkeys) = h2
keys(3,nkeys) = p1 keys(3,nkeys) = p1
keys(4,nkeys) = h2 keys(4,nkeys) = h2
nkeys += 1
values(nkeys) = 0.5d0 * c_1 * phase
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo enddo
else else
! Mono beta ! Mono beta
@ -346,11 +365,17 @@
if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
h2 = list_orb_reverse(h2) h2 = list_orb_reverse(h2)
nkeys += 1 nkeys += 1
values(nkeys) = c_1 * phase values(nkeys) = 0.5d0 * c_1 * phase
keys(1,nkeys) = h1 keys(1,nkeys) = h1
keys(2,nkeys) = h2 keys(2,nkeys) = h2
keys(3,nkeys) = p1 keys(3,nkeys) = p1
keys(4,nkeys) = h2 keys(4,nkeys) = h2
nkeys += 1
values(nkeys) = 0.5d0 * c_1 * phase
keys(1,nkeys) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo enddo
endif endif
else if(spin_trace)then else if(spin_trace)then

View File

@ -3,6 +3,7 @@ integer, parameter :: SIMD_vector = 32
integer, parameter :: N_int_max = 32 integer, parameter :: N_int_max = 32
double precision, parameter :: pi = dacos(-1.d0) double precision, parameter :: pi = dacos(-1.d0)
double precision, parameter :: inv_pi = 1.d0/dacos(-1.d0)
double precision, parameter :: sqpi = dsqrt(dacos(-1.d0)) double precision, parameter :: sqpi = dsqrt(dacos(-1.d0))
double precision, parameter :: pi_5_2 = 34.9868366552d0 double precision, parameter :: pi_5_2 = 34.9868366552d0
double precision, parameter :: dfour_pi = 4.d0*dacos(-1.d0) double precision, parameter :: dfour_pi = 4.d0*dacos(-1.d0)

View File

@ -30,7 +30,11 @@ subroutine give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_k,iorder,alph
ab = alpha * beta ab = alpha * beta
d_AB = (A_center - B_center) * (A_center - B_center) d_AB = (A_center - B_center) * (A_center - B_center)
P_center = (alpha * A_center + beta * B_center) * p_inv P_center = (alpha * A_center + beta * B_center) * p_inv
if(dabs(ab*p_inv * d_AB).lt.50.d0)then
fact_k = exp(-ab*p_inv * d_AB) fact_k = exp(-ab*p_inv * d_AB)
else
fact_k = 0.d0
endif
! Recenter the polynomials P_a and P_b on x ! Recenter the polynomials P_a and P_b on x
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
@ -51,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 ) ! 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_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 ) ! * [ 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 END_DOC
implicit none implicit none
include 'constants.include.F' include 'constants.include.F'
@ -78,6 +86,13 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center) call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center)
if (fact_k < thresh) then 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+15
P_new = 0.d0
iorder = 0
fact_k = 0.d0 fact_k = 0.d0
return return
endif endif

View File

@ -1589,9 +1589,9 @@ subroutine restore_symmetry(m,n,A,LDA,thresh)
thresh2 = dsqrt(thresh) thresh2 = dsqrt(thresh)
call nullify_small_elements(m,n,A,LDA,thresh) call nullify_small_elements(m,n,A,LDA,thresh)
if (.not.restore_symm) then ! if (.not.restore_symm) then
return ! return
endif ! endif
! TODO: Costs O(n^4), but can be improved to (2 n^2 * log(n)): ! TODO: Costs O(n^4), but can be improved to (2 n^2 * log(n)):
! - copy all values in a 1D array ! - copy all values in a 1D array

View File

@ -15,10 +15,10 @@ double precision function overlap_gaussian_x(A_center,B_center,alpha,beta,power_
call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,& call give_explicit_poly_and_gaussian_x(P_new,P_center,p,fact_p,iorder_p,alpha,&
beta,power_A,power_B,A_center,B_center,dim) beta,power_A,power_B,A_center,B_center,dim)
! if(fact_p.lt.0.000001d0)then if(fact_p.lt.1.d-20)then
! overlap_gaussian_x = 0.d0 overlap_gaussian_x = 0.d0
! return return
! endif endif
overlap_gaussian_x = 0.d0 overlap_gaussian_x = 0.d0
integer :: i integer :: i
@ -53,13 +53,13 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
integer :: iorder_p(3) integer :: iorder_p(3)
call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim) call give_explicit_poly_and_gaussian(P_new,P_center,p,fact_p,iorder_p,alpha,beta,power_A,power_B,A_center,B_center,dim)
! if(fact_p.lt.1d-20)then if(fact_p.lt.1d-20)then
! overlap_x = 0.d0 overlap_x = 1.d-10
! overlap_y = 0.d0 overlap_y = 1.d-10
! overlap_z = 0.d0 overlap_z = 1.d-10
! overlap = 0.d0 overlap = 1.d-10
! return return
! endif endif
integer :: nmax integer :: nmax
double precision :: F_integral double precision :: F_integral
nmax = maxval(iorder_p) nmax = maxval(iorder_p)

View File

@ -1,3 +1,15 @@
double precision function derf_mu_x(mu,x)
implicit none
include 'utils/constants.include.F'
double precision, intent(in) :: mu,x
if(dabs(x).gt.1.d-6)then
derf_mu_x = derf(mu * x)/x
else
derf_mu_x = inv_sq_pi * 2.d0 * mu * (1.d0 - mu*mu*x*x/3.d0)
endif
end
double precision function binom_func(i,j) double precision function binom_func(i,j)
implicit none implicit none
BEGIN_DOC BEGIN_DOC
@ -288,12 +300,12 @@ subroutine wall_time(t)
end end
BEGIN_PROVIDER [ integer, nproc ] BEGIN_PROVIDER [ integer, nproc ]
use omp_lib
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Number of current OpenMP threads ! Number of current OpenMP threads
END_DOC END_DOC
integer :: omp_get_num_threads
nproc = 1 nproc = 1
!$OMP PARALLEL !$OMP PARALLEL
!$OMP MASTER !$OMP MASTER

View File

@ -127,9 +127,9 @@ function zmq_port(ishift)
END_DOC END_DOC
integer, intent(in) :: ishift integer, intent(in) :: ishift
character*(8) :: zmq_port character*(8) :: zmq_port
!$OMP CRITICAL(write) !$OMP CRITICAL
write(zmq_port,'(I8)') zmq_port_start+ishift write(zmq_port,'(I8)') zmq_port_start+ishift
!$OMP END CRITICAL(write) !$OMP END CRITICAL
zmq_port = adjustl(trim(zmq_port)) zmq_port = adjustl(trim(zmq_port))
end end
@ -520,9 +520,9 @@ subroutine new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in)
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
zmq_socket_pull = new_zmq_pull_socket () zmq_socket_pull = new_zmq_pull_socket ()
!$OMP CRITICAL(write) !$OMP CRITICAL
write(name,'(A,I8.8)') trim(name_in)//'.', icount write(name,'(A,I8.8)') trim(name_in)//'.', icount
!$OMP END CRITICAL(write) !$OMP END CRITICAL
sze = len(trim(name)) sze = len(trim(name))
zmq_state = trim(name) zmq_state = trim(name)
call lowercase(name,sze) call lowercase(name,sze)
@ -586,9 +586,9 @@ subroutine end_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,name_in)
integer, save :: icount=0 integer, save :: icount=0
icount = icount+1 icount = icount+1
!$OMP CRITICAL(write) !$OMP CRITICAL
write(name,'(A,I8.8)') trim(name_in)//'.', icount write(name,'(A,I8.8)') trim(name_in)//'.', icount
!$OMP END CRITICAL(write) !$OMP END CRITICAL
sze = len(trim(name)) sze = len(trim(name))
call lowercase(name,sze) call lowercase(name,sze)
if (name /= zmq_state) then if (name /= zmq_state) then
@ -710,9 +710,9 @@ integer function disconnect_from_taskserver_state(zmq_to_qp_run_socket, worker_i
disconnect_from_taskserver_state = -1 disconnect_from_taskserver_state = -1
!$OMP CRITICAL(write) !$OMP CRITICAL
write(message,*) 'disconnect '//trim(state), worker_id write(message,*) 'disconnect '//trim(state), worker_id
!$OMP END CRITICAL(write) !$OMP END CRITICAL
sze = min(510,len(trim(message))) sze = min(510,len(trim(message)))
rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0) rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
@ -789,9 +789,9 @@ integer function zmq_abort(zmq_to_qp_run_socket)
character*(512) :: message character*(512) :: message
zmq_abort = 0 zmq_abort = 0
!$OMP CRITICAL(write) !$OMP CRITICAL
write(message,*) 'abort ' write(message,*) 'abort '
!$OMP END CRITICAL(write) !$OMP END CRITICAL
sze = len(trim(message)) sze = len(trim(message))
@ -833,9 +833,9 @@ integer function task_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_i
task_done_to_taskserver = 0 task_done_to_taskserver = 0
!$OMP CRITICAL(write) !$OMP CRITICAL
write(message,*) 'task_done '//trim(zmq_state), worker_id, task_id write(message,*) 'task_done '//trim(zmq_state), worker_id, task_id
!$OMP END CRITICAL(write) !$OMP END CRITICAL
sze = len(trim(message)) sze = len(trim(message))
rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0) rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
@ -868,11 +868,11 @@ integer function tasks_done_to_taskserver(zmq_to_qp_run_socket, worker_id, task_
tasks_done_to_taskserver = 0 tasks_done_to_taskserver = 0
!$OMP CRITICAL(write) !$OMP CRITICAL
allocate(character(LEN=64+n_tasks*12) :: message) allocate(character(LEN=64+n_tasks*12) :: message)
write(fmt,*) '(A,X,A,I10,X,', n_tasks, '(I11,1X))' write(fmt,*) '(A,X,A,I10,X,', n_tasks, '(I11,1X))'
write(message,*) 'task_done '//trim(zmq_state), worker_id, (task_id(k), k=1,n_tasks) write(message,*) 'task_done '//trim(zmq_state), worker_id, (task_id(k), k=1,n_tasks)
!$OMP END CRITICAL(write) !$OMP END CRITICAL
sze = len(trim(message)) sze = len(trim(message))
rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0) rc = f77_zmq_send(zmq_to_qp_run_socket, trim(message), sze, 0)
@ -914,9 +914,9 @@ integer function get_task_from_taskserver(zmq_to_qp_run_socket,worker_id,task_id
get_task_from_taskserver = 0 get_task_from_taskserver = 0
!$OMP CRITICAL(write) !$OMP CRITICAL
write(message,*) 'get_task '//trim(zmq_state), worker_id write(message,*) 'get_task '//trim(zmq_state), worker_id
!$OMP END CRITICAL(write) !$OMP END CRITICAL
sze = len(trim(message)) sze = len(trim(message))
rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0) rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0)
@ -977,9 +977,9 @@ integer function get_tasks_from_taskserver(zmq_to_qp_run_socket,worker_id,task_i
get_tasks_from_taskserver = 0 get_tasks_from_taskserver = 0
!$OMP CRITICAL(write) !$OMP CRITICAL
write(message,'(A,A,X,I10,I10)') 'get_tasks ', trim(zmq_state), worker_id, n_tasks write(message,'(A,A,X,I10,I10)') 'get_tasks ', trim(zmq_state), worker_id, n_tasks
!$OMP END CRITICAL(write) !$OMP END CRITICAL
sze = len(trim(message)) sze = len(trim(message))
rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0) rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0)
@ -1079,9 +1079,9 @@ integer function zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,mo
zmq_delete_task = 0 zmq_delete_task = 0
!$OMP CRITICAL(write) !$OMP CRITICAL
write(message,*) 'del_task ', zmq_state, task_id write(message,*) 'del_task ', zmq_state, task_id
!$OMP END CRITICAL(write) !$OMP END CRITICAL
rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0) rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)
if (rc /= len(trim(message))) then if (rc /= len(trim(message))) then
zmq_delete_task = -1 zmq_delete_task = -1
@ -1121,9 +1121,9 @@ integer function zmq_delete_task_async_send(zmq_to_qp_run_socket,task_id,sending
endif endif
zmq_delete_task_async_send = 0 zmq_delete_task_async_send = 0
!$OMP CRITICAL(write) !$OMP CRITICAL
write(message,*) 'del_task ', zmq_state, task_id write(message,*) 'del_task ', zmq_state, task_id
!$OMP END CRITICAL(write) !$OMP END CRITICAL
rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0) rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)
if (rc /= len(trim(message))) then if (rc /= len(trim(message))) then
zmq_delete_task_async_send = -1 zmq_delete_task_async_send = -1
@ -1181,10 +1181,10 @@ integer function zmq_delete_tasks(zmq_to_qp_run_socket,zmq_socket_pull,task_id,n
allocate(character(LEN=64+n_tasks*12) :: message) allocate(character(LEN=64+n_tasks*12) :: message)
!$OMP CRITICAL(write) !$OMP CRITICAL
write(fmt,*) '(A,1X,A,1X,', n_tasks, '(I11,1X))' write(fmt,*) '(A,1X,A,1X,', n_tasks, '(I11,1X))'
write(message,*) 'del_task '//trim(zmq_state), (task_id(k), k=1,n_tasks) write(message,*) 'del_task '//trim(zmq_state), (task_id(k), k=1,n_tasks)
!$OMP END CRITICAL(write) !$OMP END CRITICAL
rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0) rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)
@ -1230,10 +1230,10 @@ integer function zmq_delete_tasks_async_send(zmq_to_qp_run_socket,task_id,n_task
allocate(character(LEN=64+n_tasks*12) :: message) allocate(character(LEN=64+n_tasks*12) :: message)
!$OMP CRITICAL(write) !$OMP CRITICAL
write(fmt,*) '(A,1X,A,1X,', n_tasks, '(I11,1X))' write(fmt,*) '(A,1X,A,1X,', n_tasks, '(I11,1X))'
write(message,*) 'del_task '//trim(zmq_state), (task_id(k), k=1,n_tasks) write(message,*) 'del_task '//trim(zmq_state), (task_id(k), k=1,n_tasks)
!$OMP END CRITICAL(write) !$OMP END CRITICAL
rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0) rc = f77_zmq_send(zmq_to_qp_run_socket,trim(message),len(trim(message)),0)

View File

@ -2,7 +2,7 @@
# Stage 1 # Stage 1
# Configure QP2 # Configure QP2
./configure --install all --config ./config/travis.cfg || exit -1 ./configure --download all --install all --config ./config/travis.cfg || exit -1
# Create cache # Create cache
cd ../ cd ../