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

Merge pull request #164 from QuantumPackage/cleaning_dft

Browse Source 
fixed 2 rdms
This commit is contained in:
Anthony Scemama 2021-06-11 16:31:05 +02:00 committed by GitHub
commit b636fbfee9
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
42 changed files with 1345 additions and 149 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
config/gfortran_debug.cfg
View File

@ -51,7 +51,8 @@ FCFLAGS : -Ofast
# -g : Extra debugging information
#
[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
#################

365
scripts/opam_installer_no_usr_bin.sh Executable file
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

7
src/ao_one_e_ints/ao_overlap.irp.f
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)
c = ao_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i)
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_y(i,j) += c * overlap_y
ao_overlap_z(i,j) += c * overlap_z

18
src/basis_correction/basis_correction.irp.f
View File

@ -7,22 +7,10 @@ program basis_correction
touch read_wf
no_core_density = .True.
touch no_core_density
provide ao_two_e_integrals_in_map
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_basis_correction
! call print_e_b
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

28
src/basis_correction/print_su_pbe_ot.irp.f Normal file
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

31
src/becke_numerical_grid/EZFIO.cfg
View File

@ -33,3 +33,34 @@ doc: Number of angular grid points given from input. Warning, this number cannot
interface: ezfio,provider,ocaml
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

104
src/becke_numerical_grid/angular_extra_grid.irp.f Normal file
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

178
src/becke_numerical_grid/extra_grid.irp.f Normal file
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

60
src/becke_numerical_grid/extra_grid_vector.irp.f Normal file
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

8
src/bitmask/modify_bitmasks.irp.f
View File

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

8
src/cas_based_on_top/example.irp.f
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
END_DOC
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
character*(128) :: output
@ -15,7 +16,7 @@ subroutine write_on_top_in_real_space
i_unit_output = getUnitAndOpen(output,'w')
zmax = 2.0d0
zmax = 5.0d0
print*,'nucl_coord(1,3) = ',nucl_coord(1,3)
print*,'nucl_coord(2,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
print*,'r(3) = ',r(3)
istate = 1
write(i_unit_output,*)" z, on-top(z), n(z) "
do i = 1, nz
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
r(3) += dz
enddo
end

20
src/density_for_dft/density_for_dft.irp.f
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)
enddo
elec_alpha_valence(istate) = elec_alpha_frozen_num/elec_alpha_valence(istate)
one_e_dm_mo_alpha_for_dft(:,:,istate) = one_e_dm_mo_alpha_for_dft(:,:,istate) * 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)
endif
enddo
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
END_DOC
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
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
@ -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(:,:)
endif
if(no_core_density)then
integer :: ii,i,j
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
enddo
enddo
double precision :: elec_beta_valence(N_states),elec_beta_frozen_num
integer :: istate
if(normalize_dm)then
double precision :: elec_beta_valence(N_states),elec_beta_frozen_num
elec_beta_frozen_num = elec_beta_num - n_core_orb
elec_beta_valence = 0.d0
integer :: istate
do istate = 1, N_states
do i = 1, mo_num
elec_beta_valence(istate) += one_e_dm_mo_beta_for_dft(i,i,istate)
enddo
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)
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)
one_e_dm_mo_beta_for_dft(:,:,istate) = one_e_dm_mo_beta_for_dft(:,:,istate) * elec_beta_valence(istate)
endif
enddo
endif
endif

65
src/determinants/dipole_moments.irp.f Normal file
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

30
src/determinants/psi_energy_mono_elec.irp.f
View File

@ -27,3 +27,33 @@
psi_energy_h_core(i) = psi_energy_h_core(i) * accu
enddo
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

60
src/dft_utils_in_r/ao_in_r.irp.f
View File

@ -120,43 +120,41 @@
!$OMP END PARALLEL DO
END_PROVIDER
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)]
BEGIN_PROVIDER[double precision, aos_grad_in_r_array_transp_bis, (n_points_final_grid,ao_num,3)]
implicit none
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
integer :: i,j,k
integer :: i,j,m
double precision :: aos_array(ao_num), r(3)
do k = 1, nucl_num
do i = 1, n_pts_per_atom(k)
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
aos_in_r_array_per_atom(j,i,k) = aos_array(j)
aos_in_r_array_per_atom_transp(i,j,k) = aos_array(j)
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_bis(i,j,m) = aos_grad_in_r_array(j,i,m)
enddo
enddo
enddo
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

56
src/dft_utils_in_r/dm_in_r_routines.irp.f
View File

@ -110,6 +110,62 @@ end
grad_dm_b *= 2.d0
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)

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

@ -79,7 +79,7 @@
END_DOC
integer :: m
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 j = 1, mo_num
do m = 1, 3
@ -89,6 +89,24 @@
enddo
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, beta_dens_kin_in_r, (n_points_final_grid)]
implicit none
@ -115,8 +133,6 @@
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_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
END_DOC
integer :: m
@ -127,3 +143,21 @@
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

26
src/mo_basis/mos_in_r.irp.f
View File

@ -12,18 +12,18 @@ subroutine give_all_mos_and_grad_at_r(r,mos_array,mos_grad_array)
implicit none
double precision, intent(in) :: r(3)
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
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)
mos_array=0d0
mos_grad_array=0d0
do j = 1, mo_num
do k=1, ao_num
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(j,2) += mo_coef(k,j)*aos_grad_array(k,2)
mos_grad_array(j,3) += mo_coef(k,j)*aos_grad_array(k,3)
mos_grad_array(1,j) += mo_coef(k,j)*aos_grad_array(1,k)
mos_grad_array(2,j) += mo_coef(k,j)*aos_grad_array(2,k)
mos_grad_array(3,j) += mo_coef(k,j)*aos_grad_array(3,k)
enddo
enddo
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
double precision, intent(in) :: r(3)
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
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)
mos_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 k=1, ao_num
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(j,2) += mo_coef(k,j) * aos_grad_array(k,2)
mos_grad_array(j,3) += mo_coef(k,j) * aos_grad_array(k,3)
mos_lapl_array(j,1) += mo_coef(k,j) * aos_lapl_array(k,1)
mos_lapl_array(j,2) += mo_coef(k,j) * aos_lapl_array(k,2)
mos_lapl_array(j,3) += mo_coef(k,j) * aos_lapl_array(k,3)
mos_grad_array(1,j) += mo_coef(k,j) * aos_grad_array(1,k)
mos_grad_array(2,j) += mo_coef(k,j) * aos_grad_array(2,k)
mos_grad_array(3,j) += mo_coef(k,j) * aos_grad_array(3,k)
mos_lapl_array(1,j) += mo_coef(k,j) * aos_lapl_array(1,k)
mos_lapl_array(2,j) += mo_coef(k,j) * aos_lapl_array(2,k)
mos_lapl_array(3,j) += mo_coef(k,j) * aos_lapl_array(3,k)
enddo
enddo
end

2
src/mo_guess/h_core_guess_routine.irp.f
View File

@ -10,5 +10,5 @@ subroutine hcore_guess
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 save_mos
SOFT_TOUCH mo_coef mo_label
TOUCH mo_coef mo_label
end

2
src/mo_two_e_erf_ints/core_quantities_erf.irp.f
View File

@ -7,7 +7,7 @@ BEGIN_PROVIDER [double precision, core_energy_erf]
core_energy_erf = 0.d0
do i = 1, n_core_orb
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
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))

1
src/mu_of_r/f_hf_utils.irp.f
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
integer :: orb_i,orb_j,i,j,orb_m,orb_n,m,n
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
m = list_valence_orb_for_hf(orb_m,1)
do orb_n = 1, n_max_occ_val_orb_for_hf! electron 2

3
src/mu_of_r/f_psi_i_a_v_utils.irp.f
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
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
PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
do orb_m = 1, n_act_orb ! electron 1
m = list_act(orb_m)
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
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
PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
do orb_m = 1, n_act_orb ! electron 1
m = list_act(orb_m)
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
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)
PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals
do orb_m = 1, n_inact_orb ! electron 1
m = list_inact(orb_m)
do orb_n = 1, n_inact_orb ! electron 2

1
src/mu_of_r/mu_of_r_conditions.irp.f
View File

@ -148,3 +148,4 @@
mu_average_prov(istate) = mu_average_prov(istate) / elec_num_grid_becke(istate)
enddo
END_PROVIDER

7
src/scf_utils/EZFIO.cfg
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
interface: ezfio,provider,ocaml
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

21
src/scf_utils/diagonalize_fock.irp.f
View File

@ -31,6 +31,27 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num)
enddo
enddo
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
do i = elec_beta_num+1, elec_alpha_num

21
src/scf_utils/fock_matrix.irp.f
View File

@ -92,6 +92,27 @@
enddo
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

3
src/tools/hcore_guess.irp.f Normal file
View File

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

5
src/tools/huckel_guess.irp.f Normal file
View File

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

5
src/tools/print_dipole.irp.f Normal file
View File

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

11
src/tools/print_var_energy.irp.f Normal file
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

24
src/two_body_rdm/EZFIO.cfg
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]
type: Disk_access
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
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]
type: Disk_access
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
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]
type: Disk_access
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
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]
type: Disk_access
doc: Read/Write the active part of the two-body rdm for spin trace electrons from/to disk [ Write | Read | None ]

24
src/two_body_rdm/act_2_rdm.irp.f
View File

@ -22,6 +22,8 @@
END_DOC
integer :: ispin
double precision :: wall_1, wall_2
character*(128) :: name_file
name_file = 'act_2_rdm_ab_mo'
! condition for alpha/beta spin
print*,''
print*,'Providing act_2_rdm_ab_mo '
@ -31,13 +33,13 @@
call wall_time(wall_1)
if(read_two_body_rdm_ab)then
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
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
if(write_two_body_rdm_ab)then
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")
endif
call wall_time(wall_2)
@ -63,18 +65,20 @@
! condition for alpha/beta spin
print*,''
print*,'Providing act_2_rdm_aa_mo '
character*(128) :: name_file
name_file = 'act_2_rdm_aa_mo'
ispin = 1
act_2_rdm_aa_mo = 0.d0
call wall_time(wall_1)
if(read_two_body_rdm_aa)then
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
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
if(write_two_body_rdm_aa)then
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")
endif
@ -101,18 +105,20 @@
! condition for beta/beta spin
print*,''
print*,'Providing act_2_rdm_bb_mo '
character*(128) :: name_file
name_file = 'act_2_rdm_bb_mo'
ispin = 2
act_2_rdm_bb_mo = 0.d0
call wall_time(wall_1)
if(read_two_body_rdm_bb)then
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
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
if(write_two_body_rdm_bb)then
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")
endif
@ -138,18 +144,20 @@
! condition for beta/beta spin
print*,''
print*,'Providing act_2_rdm_spin_trace_mo '
character*(128) :: name_file
name_file = 'act_2_rdm_spin_trace_mo'
ispin = 4
act_2_rdm_spin_trace_mo = 0.d0
call wall_time(wall_1)
if(read_two_body_rdm_spin_trace)then
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
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
if(write_two_body_rdm_spin_trace)then
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")
endif

72
src/two_body_rdm/example.irp.f
View File

@ -15,7 +15,7 @@ subroutine routine_active_only
double precision :: wee_aa_st_av, rdm_aa_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_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
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 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
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*,'**************************'
do istate = 1, N_states
@ -51,6 +72,19 @@ subroutine routine_active_only
korb = list_act(k)
do l = 1, n_act_orb
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)
@ -59,6 +93,17 @@ subroutine routine_active_only
rdmaa = act_2_rdm_aa_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)
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
@ -71,8 +116,8 @@ subroutine routine_active_only
enddo
enddo
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_ab_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_ab(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)
print*,''
@ -87,7 +132,6 @@ subroutine routine_active_only
print*,'Full energy '
print*,'psi_energy_two_e(istate)= ',psi_energy_two_e(istate)
enddo
wee_aa_st_av = 0.d0
wee_bb_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)
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_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)
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)
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_bb_st_av += vijkl * rdm_bb_st_av

29
src/two_body_rdm/io_two_rdm.irp.f Normal file
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

6
src/two_body_rdm/state_av_act_2rdm.irp.f
View File

@ -119,7 +119,11 @@
call wall_time(wall_1)
double precision :: wall_1, wall_2
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)
print*,'Time to provide state_av_act_2_rdm_spin_trace_mo',wall_2 - wall_1

50
src/two_rdm_routines/update_rdm.irp.f
View File

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

41
src/two_rdm_routines/update_state_av_rdm.irp.f
View File

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

1
src/utils/constants.include.F
View File

@ -3,6 +3,7 @@ integer, parameter :: SIMD_vector = 32
integer, parameter :: N_int_max = 32
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 :: pi_5_2 = 34.9868366552d0
double precision, parameter :: dfour_pi = 4.d0*dacos(-1.d0)

10
src/utils/integration.irp.f
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
d_AB = (A_center - B_center) * (A_center - B_center)
P_center = (alpha * A_center + beta * B_center) * p_inv
fact_k = exp(-ab*p_inv * d_AB)
if(dabs(ab*p_inv * d_AB).lt.50.d0)then
fact_k = exp(-ab*p_inv * d_AB)
else
fact_k = 0.d0
endif
! Recenter the polynomials P_a and P_b on x
!DIR$ FORCEINLINE
@ -78,6 +82,10 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
!DIR$ FORCEINLINE
call gaussian_product(alpha,A_center,beta,B_center,fact_k,p,P_center)
if (fact_k < thresh) then
P_center = 0.d0
p = 1.d-10
P_new = 0.d0
iorder = -1
fact_k = 0.d0
return
endif

6
src/utils/linear_algebra.irp.f
View File

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

22
src/utils/one_e_integration.irp.f
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,&
beta,power_A,power_B,A_center,B_center,dim)
! if(fact_p.lt.0.000001d0)then
! overlap_gaussian_x = 0.d0
! return
! endif
if(fact_p.lt.1.d-20)then
overlap_gaussian_x = 0.d0
return
endif
overlap_gaussian_x = 0.d0
integer :: i
@ -53,13 +53,13 @@ subroutine overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,&
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)
! if(fact_p.lt.1d-20)then
! overlap_x = 0.d0
! overlap_y = 0.d0
! overlap_z = 0.d0
! overlap = 0.d0
! return
! endif
if(fact_p.lt.1d-20)then
overlap_x = 1.d-10
overlap_y = 1.d-10
overlap_z = 1.d-10
overlap = 1.d-10
return
endif
integer :: nmax
double precision :: F_integral
nmax = maxval(iorder_p)

12
src/utils/util.irp.f
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)
implicit none
BEGIN_DOC