PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-21 20:52:28 +02:00
Code Releases Activity

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

Browse Source
This commit is contained in:
Emmanuel Giner LCT 2024-01-25 14:43:29 +01:00
commit d35131967c
376 changed files with 8877 additions and 6858 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
bin/python Executable file
View File

@ -0,0 +1,4 @@
#!/bin/bash
exec python3 $@

1
bin/qp_convert_output_to_ezfio
View File

@ -256,6 +256,7 @@ def write_ezfio(res, filename):
MoTag = res.determinants_mo_type MoTag = res.determinants_mo_type
ezfio.set_mo_basis_mo_label('Orthonormalized') ezfio.set_mo_basis_mo_label('Orthonormalized')
ezfio.set_determinants_mo_label('Orthonormalized')
MO_type = MoTag MO_type = MoTag
allMOs = res.mo_sets[MO_type] allMOs = res.mo_sets[MO_type]

1
bin/qp_plugins
View File

@ -127,6 +127,7 @@ def main(arguments):
l_repository = list(d_tmp.keys()) l_repository = list(d_tmp.keys())
if l_repository == []: if l_repository == []:
l_result = [] l_result = []
l_plugins = []
else: else:
m_instance = ModuleHandler(l_repository) m_instance = ModuleHandler(l_repository)
l_plugins = [module for module in m_instance.l_module] l_plugins = [module for module in m_instance.l_module]

2
bin/qp_reset
View File

@ -97,6 +97,8 @@ if [[ $dets -eq 1 ]] ; then
rm --force -- ${ezfio}/determinants/psi_{det,coef}.gz rm --force -- ${ezfio}/determinants/psi_{det,coef}.gz
rm --force -- ${ezfio}/determinants/n_det_qp_edit rm --force -- ${ezfio}/determinants/n_det_qp_edit
rm --force -- ${ezfio}/determinants/psi_{det,coef}_qp_edit.gz rm --force -- ${ezfio}/determinants/psi_{det,coef}_qp_edit.gz
rm --force -- ${ezfio}/tc_bi_ortho/psi_{l,r}_coef_bi_ortho.gz
fi fi
if [[ $mos -eq 1 ]] ; then if [[ $mos -eq 1 ]] ; then

3
bin/qpsh
View File

@ -1,6 +1,7 @@
#!/bin/bash #!/bin/bash
export QP_ROOT=$(dirname "$(readlink -f "$0")")/.. REALPATH=$( cd "$(dirname "$0")" ; pwd -P )
export QP_ROOT=${REALPATH}/..
bash --init-file <(cat << EOF bash --init-file <(cat << EOF
[[ -f /etc/bashrc ]] && source /etc/bashrc [[ -f /etc/bashrc ]] && source /etc/bashrc

23
bin/zcat Executable file
View File

@ -0,0 +1,23 @@
#!/bin/bash
# On Darwin: try gzcat if available, otherwise use Python
if [[ $(uname -s) = Darwin ]] ; then
which gzcat &> /dev/null
if [[ $? -eq 0 ]] ; then
exec gzcat $@
else
exec python3 << EOF
import sys
import gzip
with gzip.open("$1", "rt") as f:
print(f.read())
EOF
fi
else
SCRIPTPATH="$( cd -- "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"
command=$(which -a zcat | grep -v "$SCRIPTPATH/" | head -1)
exec $command $@
fi

2
config/gfortran.cfg
View File

@ -10,7 +10,7 @@
# #
# #
[COMMON] [COMMON]
FC : gfortran -g -ffree-line-length-none -I . -fPIC -march=native FC : gfortran -g -ffree-line-length-none -I . -fPIC -march=native -std=legacy
LAPACK_LIB : -lblas -llapack LAPACK_LIB : -lblas -llapack
IRPF90 : irpf90 IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED

2
config/gfortran_armpl.cfg
View File

@ -13,7 +13,7 @@
# #
# #
[COMMON] [COMMON]
FC : gfortran -g -ffree-line-length-none -I . -fPIC -march=native FC : gfortran -g -ffree-line-length-none -I . -fPIC -march=native -std=legacy
LAPACK_LIB : -larmpl_lp64_mp LAPACK_LIB : -larmpl_lp64_mp
IRPF90 : irpf90 IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED

2
config/gfortran_avx.cfg
View File

@ -10,7 +10,7 @@
# #
# #
[COMMON] [COMMON]
FC : gfortran -ffree-line-length-none -I . -mavx -g -fPIC FC : gfortran -ffree-line-length-none -I . -mavx -g -fPIC -std=legacy
LAPACK_LIB : -llapack -lblas LAPACK_LIB : -llapack -lblas
IRPF90 : irpf90 IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DSET_NESTED IRPF90_FLAGS : --ninja --align=32 -DSET_NESTED

2
config/gfortran_debug.cfg
View File

@ -10,7 +10,7 @@
# #
# #
[COMMON] [COMMON]
FC : gfortran -g -ffree-line-length-none -I . -fPIC FC : gfortran -g -ffree-line-length-none -I . -fPIC -std=legacy
LAPACK_LIB : -lblas -llapack LAPACK_LIB : -lblas -llapack
IRPF90 : irpf90 IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED

62
config/gfortran_macos.cfg Normal file
View File

@ -0,0 +1,62 @@
# Common flags
##############
#
# -ffree-line-length-none : Needed for IRPF90 which produces long lines
# -lblas -llapack : Link with libblas and liblapack libraries provided by the system
# -I . : Include the curent directory (Mandatory)
#
# --ninja : Allow the utilisation of ninja. (Mandatory)
# --align=32 : Align all provided arrays on a 32-byte boundary
#
#
[COMMON]
FC : gfortran -ffree-line-length-none -I . -g -fPIC -std=legacy
LAPACK_LIB : -llapack -lblas
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DSET_NESTED -DMACOS
# Global options
################
#
# 1 : Activate
# 0 : Deactivate
#
[OPTION]
MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below
CACHE : 0 ; Enable cache_compile.py
OPENMP : 1 ; Append OpenMP flags
# Optimization flags
####################
#
# -Ofast : Disregard strict standards compliance. Enables all -O3 optimizations.
# It also enables optimizations that are not valid
# for all standard-compliant programs. It turns on
# -ffast-math and the Fortran-specific
# -fno-protect-parens and -fstack-arrays.
[OPT]
FCFLAGS : -Ofast -march=native
# Profiling flags
#################
#
[PROFILE]
FC : -p -g
FCFLAGS : -Ofast
# Debugging flags
#################
#
# -fcheck=all : Checks uninitialized variables, array subscripts, etc...
# -g : Extra debugging information
#
[DEBUG]
FCFLAGS : -fcheck=all -g
# OpenMP flags
#################
#
[OPENMP]
FC : -fopenmp
IRPF90_FLAGS : --openmp

2
config/gfortran_mpi.cfg
View File

@ -10,7 +10,7 @@
# #
# #
[COMMON] [COMMON]
FC : mpif90 -ffree-line-length-none -I . -g -fPIC FC : mpif90 -ffree-line-length-none -I . -g -fPIC -std=legacy
LAPACK_LIB : -lblas -llapack LAPACK_LIB : -lblas -llapack
IRPF90 : irpf90 IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DMPI -DSET_NESTED IRPF90_FLAGS : --ninja --align=32 -DMPI -DSET_NESTED

2
config/gfortran_openblas.cfg
View File

@ -10,7 +10,7 @@
# #
# #
[COMMON] [COMMON]
FC : gfortran -g -ffree-line-length-none -I . -fPIC -march=native FC : gfortran -g -ffree-line-length-none -I . -fPIC -march=native -std=legacy
LAPACK_LIB : -lopenblas LAPACK_LIB : -lopenblas
IRPF90 : irpf90 IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED

38
configure vendored
View File

@ -19,7 +19,11 @@ git submodule init
git submodule update git submodule update
# Update ARM or x86 dependencies # Update ARM or x86 dependencies
ARCHITECTURE=$(uname -m) SYSTEM=$(uname -s)
if [[ $SYSTEM = "Linux" ]] ; then
SYSTEM=""
fi
ARCHITECTURE=$(uname -m)$SYSTEM
cd ${QP_ROOT}/external/qp2-dependencies cd ${QP_ROOT}/external/qp2-dependencies
git checkout master git checkout master
git pull git pull
@ -191,7 +195,7 @@ if [[ "${PACKAGES}.x" != ".x" ]] ; then
fi fi
if [[ ${PACKAGES} = all ]] ; then if [[ ${PACKAGES} = all ]] ; then
PACKAGES="zlib ninja zeromq f77zmq gmp ocaml docopt resultsFile bats trexio qmckl" PACKAGES="zlib ninja zeromq f77zmq gmp ocaml docopt resultsFile bats trexio"
fi fi
@ -211,6 +215,7 @@ EOF
execute << EOF execute << EOF
cd "\${QP_ROOT}"/external cd "\${QP_ROOT}"/external
wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz
rm -rf trexio-${VERSION}
tar -zxf trexio-${VERSION}.tar.gz && rm trexio-${VERSION}.tar.gz tar -zxf trexio-${VERSION}.tar.gz && rm trexio-${VERSION}.tar.gz
cd trexio-${VERSION} cd trexio-${VERSION}
./configure --prefix=\${QP_ROOT} --without-hdf5 CFLAGS='-g' ./configure --prefix=\${QP_ROOT} --without-hdf5 CFLAGS='-g'
@ -224,6 +229,7 @@ EOF
execute << EOF execute << EOF
cd "\${QP_ROOT}"/external cd "\${QP_ROOT}"/external
wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz wget https://github.com/TREX-CoE/trexio/releases/download/v${VERSION}/trexio-${VERSION}.tar.gz
rm -rf trexio-${VERSION}
tar -zxf trexio-${VERSION}.tar.gz && rm trexio-${VERSION}.tar.gz tar -zxf trexio-${VERSION}.tar.gz && rm trexio-${VERSION}.tar.gz
cd trexio-${VERSION} cd trexio-${VERSION}
./configure --prefix=\${QP_ROOT} CFLAGS="-g" ./configure --prefix=\${QP_ROOT} CFLAGS="-g"
@ -231,15 +237,28 @@ EOF
EOF EOF
elif [[ ${PACKAGE} = qmckl ]] ; then elif [[ ${PACKAGE} = qmckl ]] ; then
VERSION=0.5.3 VERSION=0.5.4
execute << EOF execute << EOF
cd "\${QP_ROOT}"/external cd "\${QP_ROOT}"/external
wget https://github.com/TREX-CoE/qmckl/releases/download/v${VERSION}/qmckl-${VERSION}.tar.gz wget https://github.com/TREX-CoE/qmckl/releases/download/v${VERSION}/qmckl-${VERSION}.tar.gz
rm -rf qmckl-${VERSION}
tar -zxf qmckl-${VERSION}.tar.gz && rm qmckl-${VERSION}.tar.gz tar -zxf qmckl-${VERSION}.tar.gz && rm qmckl-${VERSION}.tar.gz
cd qmckl-${VERSION} cd qmckl-${VERSION}
./configure --prefix=\${QP_ROOT} --enable-hpc --disable-doc CFLAGS='-g' ./configure --prefix=\${QP_ROOT} --enable-hpc --disable-doc CFLAGS='-g'
make && make -j 4 check && make install make && make -j 4 check && make install
EOF EOF
elif [[ ${PACKAGE} = qmckl-intel ]] ; then
VERSION=0.5.4
execute << EOF
cd "\${QP_ROOT}"/external
wget https://github.com/TREX-CoE/qmckl/releases/download/v${VERSION}/qmckl-${VERSION}.tar.gz
rm -rf qmckl-${VERSION}
tar -zxf qmckl-${VERSION}.tar.gz && rm qmckl-${VERSION}.tar.gz
cd qmckl-${VERSION}
./configure --prefix=\${QP_ROOT} --enable-hpc --disable-doc --with-icc --with-ifort CFLAGS='-g'
make && make -j 4 check && make install
EOF
elif [[ ${PACKAGE} = gmp ]] ; then elif [[ ${PACKAGE} = gmp ]] ; then
@ -260,6 +279,7 @@ EOF
cd "\${QP_ROOT}"/external cd "\${QP_ROOT}"/external
tar --gunzip --extract --file qp2-dependencies/zeromq-4.2.5.tar.gz tar --gunzip --extract --file qp2-dependencies/zeromq-4.2.5.tar.gz
cd zeromq-* cd zeromq-*
[[ "${SYSTEM}" = "Darwin" ]] && ./autogen.sh
./configure --prefix="\$QP_ROOT" --without-libsodium --enable-libunwind=no ./configure --prefix="\$QP_ROOT" --without-libsodium --enable-libunwind=no
make -j 8 make -j 8
make install make install
@ -378,15 +398,15 @@ fi
TREXIO=$(find_lib -ltrexio) TREXIO=$(find_lib -ltrexio)
if [[ ${TREXIO} = $(not_found) ]] ; then if [[ ${TREXIO} = $(not_found) ]] ; then
error "TREXIO (trexio,trexio-nohdf5) is not installed. If you don't have HDF5, use trexio-nohdf5" error "TREXIO (trexio | trexio-nohdf5) is not installed. If you don't have HDF5, use trexio-nohdf5"
fail fail
fi fi
QMCKL=$(find_lib -lqmckl) #QMCKL=$(find_lib -lqmckl)
if [[ ${QMCKL} = $(not_found) ]] ; then #if [[ ${QMCKL} = $(not_found) ]] ; then
error "QMCkl (qmckl) is not installed." # error "QMCkl (qmckl | qmckl-intel) is not installed."
fail # fail
fi #fi
F77ZMQ=$(find_lib -lzmq -lf77zmq -lpthread) F77ZMQ=$(find_lib -lzmq -lf77zmq -lpthread)
if [[ ${F77ZMQ} = $(not_found) ]] ; then if [[ ${F77ZMQ} = $(not_found) ]] ; then

2
etc/paths.rc
View File

@ -32,7 +32,7 @@ export PYTHONPATH=$(qp_prepend_export "PYTHONPATH" "${QP_EZFIO}/Python":"${QP_PY
export PATH=$(qp_prepend_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml) export PATH=$(qp_prepend_export "PATH" "${QP_PYTHON}":"${QP_ROOT}"/bin:"${QP_ROOT}"/ocaml)
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:"${QP_ROOT}"/lib64)
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)

2
external/irpf90 vendored

@ -1 +1 @@
Subproject commit 0007f72f677fe7d61c5e1ed461882cb239517102 Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6

33
ocaml/Input_determinants_by_hand.ml
View File

@ -13,6 +13,7 @@ module Determinants_by_hand : sig
psi_coef : Det_coef.t array; psi_coef : Det_coef.t array;
psi_det : Determinant.t array; psi_det : Determinant.t array;
state_average_weight : Positive_float.t array; state_average_weight : Positive_float.t array;
mo_label : MO_label.t;
} [@@deriving sexp] } [@@deriving sexp]
val read : ?full:bool -> unit -> t option val read : ?full:bool -> unit -> t option
val write : ?force:bool -> t -> unit val write : ?force:bool -> t -> unit
@ -34,11 +35,21 @@ end = struct
psi_coef : Det_coef.t array; psi_coef : Det_coef.t array;
psi_det : Determinant.t array; psi_det : Determinant.t array;
state_average_weight : Positive_float.t array; state_average_weight : Positive_float.t array;
mo_label : MO_label.t;
} [@@deriving sexp] } [@@deriving sexp]
;; ;;
let get_default = Qpackage.get_ezfio_default "determinants";; let get_default = Qpackage.get_ezfio_default "determinants";;
let read_mo_label () =
if not (Ezfio.has_determinants_mo_label ()) then
if Ezfio.has_mo_basis_mo_label () then (
let label = Ezfio.get_mo_basis_mo_label () in
Ezfio.set_determinants_mo_label label) ;
Ezfio.get_determinants_mo_label ()
|> MO_label.of_string
;;
let read_n_int () = let read_n_int () =
if not (Ezfio.has_determinants_n_int()) then if not (Ezfio.has_determinants_n_int()) then
Ezfio.get_mo_basis_mo_num () Ezfio.get_mo_basis_mo_num ()
@ -222,7 +233,7 @@ end = struct
and n_states = and n_states =
States_number.to_int n_states States_number.to_int n_states
in in
let r = let r =
Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c Ezfio.ezfio_array_of_list ~rank:2 ~dim:[| n_det ; n_states |] ~data:c
in in
Ezfio.set_determinants_psi_coef r; Ezfio.set_determinants_psi_coef r;
@ -283,19 +294,23 @@ end = struct
|> Array.concat |> Array.concat
|> Array.to_list |> Array.to_list
in in
let r = let r =
Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; Det_number.to_int n_det |] ~data:data Ezfio.ezfio_array_of_list ~rank:3 ~dim:[| N_int_number.to_int n_int ; 2 ; Det_number.to_int n_det |] ~data:data
in in
Ezfio.set_determinants_psi_det r; Ezfio.set_determinants_psi_det r;
Ezfio.set_determinants_psi_det_qp_edit r Ezfio.set_determinants_psi_det_qp_edit r
;; ;;
let write_mo_label a =
MO_label.to_string a
|> Ezfio.set_determinants_mo_label
let read ?(full=true) () = let read ?(full=true) () =
let n_det_qp_edit = read_n_det_qp_edit () in let n_det_qp_edit = read_n_det_qp_edit () in
let n_det = read_n_det () in let n_det = read_n_det () in
let read_only = let read_only =
if full then false else n_det_qp_edit <> n_det if full then false else n_det_qp_edit <> n_det
in in
@ -311,6 +326,7 @@ end = struct
psi_det = read_psi_det ~read_only () ; psi_det = read_psi_det ~read_only () ;
n_states = read_n_states () ; n_states = read_n_states () ;
state_average_weight = read_state_average_weight () ; state_average_weight = read_state_average_weight () ;
mo_label = read_mo_label () ;
} }
with _ -> None with _ -> None
else else
@ -328,6 +344,7 @@ end = struct
psi_det ; psi_det ;
n_states ; n_states ;
state_average_weight ; state_average_weight ;
mo_label ;
} = } =
write_n_int n_int ; write_n_int n_int ;
write_bit_kind bit_kind; write_bit_kind bit_kind;
@ -340,7 +357,9 @@ end = struct
write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ; write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ;
write_psi_det ~n_int:n_int ~n_det:n_det psi_det write_psi_det ~n_int:n_int ~n_det:n_det psi_det
end; end;
write_state_average_weight state_average_weight write_state_average_weight state_average_weight ;
write_mo_label mo_label ;
()
;; ;;
@ -439,7 +458,7 @@ psi_det = %s
in in
(* Split into header and determinants data *) (* Split into header and determinants data *)
let idx = let idx =
match String_ext.substr_index r ~pos:0 ~pattern:"\nDeterminants" with match String_ext.substr_index r ~pos:0 ~pattern:"\nDeterminants" with
| Some x -> x | Some x -> x
| None -> assert false | None -> assert false
@ -545,6 +564,8 @@ psi_det = %s
let bitkind = let bitkind =
Printf.sprintf "(bit_kind %d)" (Lazy.force Qpackage.bit_kind Printf.sprintf "(bit_kind %d)" (Lazy.force Qpackage.bit_kind
|> Bit_kind.to_int) |> Bit_kind.to_int)
and mo_label =
Printf.sprintf "(mo_label %s)" (MO_label.to_string @@ read_mo_label ())
and n_int = and n_int =
Printf.sprintf "(n_int %d)" (N_int_number.get_max ()) Printf.sprintf "(n_int %d)" (N_int_number.get_max ())
and n_states = and n_states =
@ -553,7 +574,7 @@ psi_det = %s
Printf.sprintf "(n_det_qp_edit %d)" (Det_number.to_int @@ read_n_det_qp_edit ()) Printf.sprintf "(n_det_qp_edit %d)" (Det_number.to_int @@ read_n_det_qp_edit ())
in in
let s = let s =
String.concat "" [ header ; bitkind ; n_int ; n_states ; psi_coef ; psi_det ; n_det_qp_edit ] String.concat "" [ header ; mo_label ; bitkind ; n_int ; n_states ; psi_coef ; psi_det ; n_det_qp_edit ]
in in

21
ocaml/Zmatrix.ml
View File

@ -58,17 +58,32 @@ let int_of_atom_id : atom_id -> int = fun x -> x
let float_of_distance : float StringMap.t -> distance -> float = let float_of_distance : float StringMap.t -> distance -> float =
fun map -> function fun map -> function
| Value x -> x | Value x -> x
| Label s -> StringMap.find s map | Label s -> begin
try StringMap.find s map with
| Not_found ->
Printf.sprintf "Zmatrix error: distance %s undefined" s
|> failwith
end
let float_of_angle : float StringMap.t -> angle -> float = let float_of_angle : float StringMap.t -> angle -> float =
fun map -> function fun map -> function
| Value x -> x | Value x -> x
| Label s -> StringMap.find s map | Label s -> begin
try StringMap.find s map with
| Not_found ->
Printf.sprintf "Zmatrix error: angle %s undefined" s
|> failwith
end
let float_of_dihedral : float StringMap.t -> dihedral -> float = let float_of_dihedral : float StringMap.t -> dihedral -> float =
fun map -> function fun map -> function
| Value x -> x | Value x -> x
| Label s -> StringMap.find s map | Label s -> begin
try StringMap.find s map with
| Not_found ->
Printf.sprintf "Zmatrix error: dihedral %s undefined" s
|> failwith
end
type line = type line =

4
ocaml/qptypes_generator.ml
View File

@ -154,8 +154,8 @@ let input_ezfio = "
* N_int_number : int * N_int_number : int
determinants_n_int determinants_n_int
1 : 30 1 : 128
N_int > 30 N_int > 128
* Det_number : int * Det_number : int
determinants_n_det determinants_n_det

1
plugins/.gitignore vendored
View File

@ -1,2 +1 @@
*

1
src/ao_many_one_e_ints/NEED → plugins/local/ao_many_one_e_ints/NEED
View File

@ -4,3 +4,4 @@ becke_numerical_grid
mo_one_e_ints mo_one_e_ints
dft_utils_in_r dft_utils_in_r
tc_keywords tc_keywords
hamiltonian

0
src/ao_many_one_e_ints/README.rst → plugins/local/ao_many_one_e_ints/README.rst
View File

38
src/ao_many_one_e_ints/ao_erf_gauss.irp.f → plugins/local/ao_many_one_e_ints/ao_erf_gauss.irp.f
View File

@ -98,7 +98,7 @@ double precision function phi_j_erf_mu_r_phi(i, j, mu_in, C_center)
enddo enddo
enddo enddo
end function phi_j_erf_mu_r_phi end
! --- ! ---
@ -201,7 +201,7 @@ subroutine erf_mu_gauss_ij_ao(i, j, mu, C_center, delta, gauss_ints)
enddo enddo
enddo enddo
end subroutine erf_mu_gauss_ij_ao end
! --- ! ---
@ -266,7 +266,7 @@ subroutine NAI_pol_x_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
enddo enddo
enddo enddo
end subroutine NAI_pol_x_mult_erf_ao end
! --- ! ---
@ -340,7 +340,7 @@ subroutine NAI_pol_x_mult_erf_ao_v0(i_ao, j_ao, mu_in, C_center, LD_C, ints, LD_
deallocate(integral) deallocate(integral)
end subroutine NAI_pol_x_mult_erf_ao_v0 end
! --- ! ---
@ -420,7 +420,7 @@ subroutine NAI_pol_x_mult_erf_ao_v(i_ao, j_ao, mu_in, C_center, LD_C, ints, LD_i
deallocate(integral) deallocate(integral)
end subroutine NAI_pol_x_mult_erf_ao_v end
! --- ! ---
@ -479,7 +479,7 @@ double precision function NAI_pol_x_mult_erf_ao_x(i_ao, j_ao, mu_in, C_center)
enddo enddo
enddo enddo
end function NAI_pol_x_mult_erf_ao_x end
! --- ! ---
@ -538,7 +538,7 @@ double precision function NAI_pol_x_mult_erf_ao_y(i_ao, j_ao, mu_in, C_center)
enddo enddo
enddo enddo
end function NAI_pol_x_mult_erf_ao_y end
! --- ! ---
@ -597,7 +597,7 @@ double precision function NAI_pol_x_mult_erf_ao_z(i_ao, j_ao, mu_in, C_center)
enddo enddo
enddo enddo
end function NAI_pol_x_mult_erf_ao_z end
! --- ! ---
@ -667,7 +667,7 @@ double precision function NAI_pol_x_mult_erf_ao_with1s_x(i_ao, j_ao, beta, B_cen
enddo enddo
enddo enddo
end function NAI_pol_x_mult_erf_ao_with1s_x end
! --- ! ---
@ -737,7 +737,7 @@ double precision function NAI_pol_x_mult_erf_ao_with1s_y(i_ao, j_ao, beta, B_cen
enddo enddo
enddo enddo
end function NAI_pol_x_mult_erf_ao_with1s_y end
! --- ! ---
@ -807,7 +807,7 @@ double precision function NAI_pol_x_mult_erf_ao_with1s_z(i_ao, j_ao, beta, B_cen
enddo enddo
enddo enddo
end function NAI_pol_x_mult_erf_ao_with1s_z end
! --- ! ---
@ -880,7 +880,7 @@ subroutine NAI_pol_x_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_cen
enddo enddo
enddo enddo
end subroutine NAI_pol_x_mult_erf_ao_with1s end
! --- ! ---
@ -967,7 +967,7 @@ subroutine NAI_pol_x_mult_erf_ao_with1s_v0(i_ao, j_ao, beta, B_center, LD_B, mu_
deallocate(integral) deallocate(integral)
end subroutine NAI_pol_x_mult_erf_ao_with1s_v0 end
! --- ! ---
@ -1057,7 +1057,7 @@ subroutine NAI_pol_x_mult_erf_ao_with1s_v(i_ao, j_ao, beta, B_center, LD_B, mu_i
deallocate(integral) deallocate(integral)
end subroutine NAI_pol_x_mult_erf_ao_with1s_v end
! --- ! ---
@ -1175,7 +1175,7 @@ subroutine NAI_pol_x2_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_ce
enddo enddo
enddo enddo
end subroutine NAI_pol_x2_mult_erf_ao_with1s end
! --- ! ---
@ -1241,7 +1241,7 @@ subroutine NAI_pol_x2_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
enddo enddo
enddo enddo
end subroutine NAI_pol_x2_mult_erf_ao end
! --- ! ---
@ -1320,7 +1320,7 @@ subroutine NAI_pol_012_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_c
enddo enddo
enddo enddo
end subroutine NAI_pol_012_mult_erf_ao_with1s end
! --- ! ---
@ -1328,7 +1328,7 @@ subroutine NAI_pol_012_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
BEGIN_DOC BEGIN_DOC
! !
! Computes the following integral : ! Computes the following integrals :
! !
! int(1) = $\int_{-\infty}^{infty} dr x^0 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$. ! int(1) = $\int_{-\infty}^{infty} dr x^0 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
! !
@ -1395,7 +1395,7 @@ subroutine NAI_pol_012_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
enddo enddo
enddo enddo
end subroutine NAI_pol_012_mult_erf_ao end
! --- ! ---

0
src/ao_many_one_e_ints/ao_erf_gauss_grad.irp.f → plugins/local/ao_many_one_e_ints/ao_erf_gauss_grad.irp.f
View File

89
src/ao_many_one_e_ints/ao_gaus_gauss.irp.f → plugins/local/ao_many_one_e_ints/ao_gaus_gauss.irp.f
View File

@ -152,7 +152,7 @@ double precision function overlap_gauss_r12_ao(D_center, delta, i, j)
enddo enddo
enddo enddo
end function overlap_gauss_r12_ao end
! -- ! --
@ -199,7 +199,7 @@ double precision function overlap_abs_gauss_r12_ao(D_center, delta, i, j)
enddo enddo
enddo enddo
end function overlap_gauss_r12_ao end
! -- ! --
@ -257,7 +257,7 @@ subroutine overlap_gauss_r12_ao_v(D_center, LD_D, delta, i, j, resv, LD_resv, n_
deallocate(analytical_j) deallocate(analytical_j)
end subroutine overlap_gauss_r12_ao_v end
! --- ! ---
@ -327,7 +327,7 @@ double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center,
enddo enddo
enddo enddo
end function overlap_gauss_r12_ao_with1s end
! --- ! ---
@ -420,7 +420,86 @@ subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, LD_D, delta,
deallocate(fact_g, G_center, analytical_j) deallocate(fact_g, G_center, analytical_j)
end subroutine overlap_gauss_r12_ao_with1s_v end
! ---
subroutine overlap_gauss_r12_ao_012(D_center, delta, i, j, ints)
BEGIN_DOC
!
! Computes the following integrals :
!
! ints(1) = $\int_{-\infty}^{infty} dr x^0 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
!
! ints(2) = $\int_{-\infty}^{infty} dr x^1 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(3) = $\int_{-\infty}^{infty} dr y^1 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(4) = $\int_{-\infty}^{infty} dr z^1 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
!
! ints(5) = $\int_{-\infty}^{infty} dr x^2 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(6) = $\int_{-\infty}^{infty} dr y^2 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(7) = $\int_{-\infty}^{infty} dr z^2 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: i, j
double precision, intent(in) :: delta, D_center(3)
double precision, intent(out) :: ints(7)
integer :: k, l, m
integer :: power_A(3), power_B(3), power_A1(3), power_A2(3)
double precision :: A_center(3), B_center(3), alpha, beta, coef1, coef
double precision :: integral0, integral1, integral2
double precision, external :: overlap_gauss_r12
ints = 0.d0
if(ao_overlap_abs(j,i).lt.1.d-12) then
return
endif
power_A(1:3) = ao_power(i,1:3)
power_B(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(ao_nucl(i),1:3)
B_center(1:3) = nucl_coord(ao_nucl(j),1:3)
do l = 1, ao_prim_num(i)
alpha = ao_expo_ordered_transp (l,i)
coef1 = ao_coef_normalized_ordered_transp(l,i)
do k = 1, ao_prim_num(j)
beta = ao_expo_ordered_transp(k,j)
coef = coef1 * ao_coef_normalized_ordered_transp(k,j)
if(dabs(coef) .lt. 1d-12) cycle
integral0 = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta)
ints(1) += coef * integral0
do m = 1, 3
power_A1 = power_A
power_A1(m) += 1
integral1 = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A1, power_B, alpha, beta)
ints(1+m) += coef * (integral1 + A_center(m)*integral0)
power_A2 = power_A
power_A2(m) += 2
integral2 = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A2, power_B, alpha, beta)
ints(4+m) += coef * (integral2 + A_center(m) * (2.d0*integral1 + A_center(m)*integral0))
enddo
enddo ! k
enddo ! l
return
end
! --- ! ---

198
src/ao_many_one_e_ints/grad2_jmu_manu.irp.f → plugins/local/ao_many_one_e_ints/grad2_jmu_manu.irp.f
View File

@ -1,11 +1,11 @@
! --- ! ---
BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, int2_grad1u2_grad2u2_env2_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2 ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 [1 - erf(mu r12)]^2
! !
END_DOC END_DOC
@ -15,30 +15,30 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
double precision :: coef, beta, B_center(3) double precision :: coef, beta, B_center(3)
double precision :: tmp double precision :: tmp
double precision :: wall0, wall1 double precision :: wall0, wall1
double precision :: int_gauss, dsqpi_3_2, int_j1b double precision :: int_gauss, dsqpi_3_2, int_env
double precision :: factor_ij_1s, beta_ij, center_ij_1s(3), sq_pi_3_2 double precision :: factor_ij_1s, beta_ij, center_ij_1s(3), sq_pi_3_2
double precision, allocatable :: int_fit_v(:) double precision, allocatable :: int_fit_v(:)
double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing int2_grad1u2_grad2u2_j1b2_test ...' print*, ' providing int2_grad1u2_grad2u2_env2_test ...'
sq_pi_3_2 = (dacos(-1.d0))**(1.5d0) sq_pi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points_transp j1b_pen List_comb_thr_b3_coef provide mu_erf final_grid_points_transp List_comb_thr_b3_coef
call wall_time(wall0) call wall_time(wall0)
int2_grad1u2_grad2u2_j1b2_test(:,:,:) = 0.d0 int2_grad1u2_grad2u2_env2_test(:,:,:) = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit_v, tmp,int_gauss,int_j1b,factor_ij_1s,beta_ij,center_ij_1s) & !$OMP coef_fit, expo_fit, int_fit_v, tmp,int_gauss,int_env,factor_ij_1s,beta_ij,center_ij_1s) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points,List_comb_thr_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, final_grid_points,List_comb_thr_b3_size, &
!$OMP final_grid_points_transp, ng_fit_jast, & !$OMP final_grid_points_transp, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b, & !$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_env2_test, ao_abs_comb_b3_env, &
!$OMP ao_overlap_abs,sq_pi_3_2,thrsh_cycle_tc) !$OMP ao_overlap_abs,sq_pi_3_2,thrsh_cycle_tc)
!$OMP DO SCHEDULE(dynamic) !$OMP DO SCHEDULE(dynamic)
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint) r(1) = final_grid_points(1,ipoint)
@ -54,13 +54,13 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
! i_1s = 1 ! i_1s = 1
! --- --- --- ! --- --- ---
int_j1b = ao_abs_comb_b3_j1b(1,j,i) int_env = ao_abs_comb_b3_env(1,j,i)
do i_fit = 1, ng_fit_jast do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_1_erf_x_2(i_fit) expo_fit = expo_gauss_1_erf_x_2(i_fit)
coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit)
! if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle ! if(dabs(coef_fit*int_env*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
int_gauss = overlap_gauss_r12_ao(r, expo_fit, i, j) int_gauss = overlap_gauss_r12_ao(r, expo_fit, i, j)
int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss int2_grad1u2_grad2u2_env2_test(j,i,ipoint) += coef_fit * int_gauss
enddo enddo
! --- --- --- ! --- --- ---
@ -71,7 +71,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
coef = List_comb_thr_b3_coef (i_1s,j,i) coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) int_env = ao_abs_comb_b3_env(i_1s,j,i)
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -81,11 +81,11 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef
! if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle ! if(dabs(coef_fit*factor_ij_1s*int_env*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, & ! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, &
! expo_fit, i, j, int_fit_v, n_points_final_grid) ! expo_fit, i, j, int_fit_v, n_points_final_grid)
int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss int2_grad1u2_grad2u2_env2_test(j,i,ipoint) += coef_fit * int_gauss
enddo enddo
enddo enddo
@ -98,26 +98,26 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 1, ao_num do i = 1, ao_num
do j = 1, i-1 do j = 1, i-1
int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) int2_grad1u2_grad2u2_env2_test(j,i,ipoint) = int2_grad1u2_grad2u2_env2_test(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test', wall1 - wall0 print*, ' wall time for int2_grad1u2_grad2u2_env2_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, int2_grad1u2_grad2u2_env2_test_v, (ao_num, ao_num, n_points_final_grid)]
!
! BEGIN_DOC BEGIN_DOC
! ! !
! ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2 ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 [1 - erf(mu r12)]^2
! ! !
! END_DOC END_DOC
!
implicit none implicit none
integer :: i, j, ipoint, i_1s, i_fit integer :: i, j, ipoint, i_1s, i_fit
double precision :: r(3), expo_fit, coef_fit double precision :: r(3), expo_fit, coef_fit
@ -128,24 +128,24 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
double precision, allocatable :: int_fit_v(:),big_array(:,:,:) double precision, allocatable :: int_fit_v(:),big_array(:,:,:)
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing int2_grad1u2_grad2u2_j1b2_test_v ...' print*, ' providing int2_grad1u2_grad2u2_env2_test_v ...'
provide mu_erf final_grid_points_transp j1b_pen provide mu_erf final_grid_points_transp
call wall_time(wall0) call wall_time(wall0)
double precision :: int_j1b double precision :: int_env
big_array(:,:,:) = 0.d0 big_array(:,:,:) = 0.d0
allocate(big_array(n_points_final_grid,ao_num, ao_num)) allocate(big_array(n_points_final_grid,ao_num, ao_num))
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,& !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,&
!$OMP coef_fit, expo_fit, int_fit_v, tmp,int_j1b) & !$OMP coef_fit, expo_fit, int_fit_v, tmp,int_env) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size,& !$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size,&
!$OMP final_grid_points_transp, ng_fit_jast, & !$OMP final_grid_points_transp, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, big_array,& !$OMP List_comb_thr_b3_cent, big_array,&
!$OMP ao_abs_comb_b3_j1b,ao_overlap_abs,thrsh_cycle_tc) !$OMP ao_abs_comb_b3_env,ao_overlap_abs,thrsh_cycle_tc)
! !
allocate(int_fit_v(n_points_final_grid)) allocate(int_fit_v(n_points_final_grid))
!$OMP DO SCHEDULE(dynamic) !$OMP DO SCHEDULE(dynamic)
do i = 1, ao_num do i = 1, ao_num
@ -159,7 +159,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
coef = List_comb_thr_b3_coef (i_1s,j,i) coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) int_env = ao_abs_comb_b3_env(i_1s,j,i)
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -187,7 +187,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
do i = 1, ao_num do i = 1, ao_num
do j = i, ao_num do j = i, ao_num
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
int2_grad1u2_grad2u2_j1b2_test_v(j,i,ipoint) = big_array(ipoint,j,i) int2_grad1u2_grad2u2_env2_test_v(j,i,ipoint) = big_array(ipoint,j,i)
enddo enddo
enddo enddo
enddo enddo
@ -195,23 +195,23 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
int2_grad1u2_grad2u2_j1b2_test_v(j,i,ipoint) = big_array(ipoint,i,j) int2_grad1u2_grad2u2_env2_test_v(j,i,ipoint) = big_array(ipoint,i,j)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test_v', wall1 - wall0 print*, ' wall time for int2_grad1u2_grad2u2_env2_test_v (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, int2_u2_env2_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [u_12^mu]^2 ! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 [u_12^mu]^2
! !
END_DOC END_DOC
@ -219,29 +219,29 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
integer :: i, j, ipoint, i_1s, i_fit integer :: i, j, ipoint, i_1s, i_fit
double precision :: r(3), int_fit, expo_fit, coef_fit double precision :: r(3), int_fit, expo_fit, coef_fit
double precision :: coef, beta, B_center(3), tmp double precision :: coef, beta, B_center(3), tmp
double precision :: wall0, wall1,int_j1b double precision :: wall0, wall1,int_env
double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
double precision :: factor_ij_1s,beta_ij,center_ij_1s(3),sq_pi_3_2 double precision :: factor_ij_1s,beta_ij,center_ij_1s(3),sq_pi_3_2
print*, ' providing int2_u2_j1b2_test ...' print*, ' providing int2_u2_env2_test ...'
sq_pi_3_2 = (dacos(-1.d0))**(1.5d0) sq_pi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points j1b_pen provide mu_erf final_grid_points
call wall_time(wall0) call wall_time(wall0)
int2_u2_j1b2_test = 0.d0 int2_u2_env2_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp, int_j1b,factor_ij_1s,beta_ij,center_ij_1s) & !$OMP coef_fit, expo_fit, int_fit, tmp, int_env,factor_ij_1s,beta_ij,center_ij_1s) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, &
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, & !$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo,sq_pi_3_2, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo,sq_pi_3_2, &
!$OMP List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b,thrsh_cycle_tc) !$OMP List_comb_thr_b3_cent, int2_u2_env2_test,ao_abs_comb_b3_env,thrsh_cycle_tc)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint) r(1) = final_grid_points(1,ipoint)
@ -257,12 +257,12 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
! i_1s = 1 ! i_1s = 1
! --- --- --- ! --- --- ---
int_j1b = ao_abs_comb_b3_j1b(1,j,i) int_env = ao_abs_comb_b3_env(1,j,i)
if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle if(dabs(int_env).lt.thrsh_cycle_tc) cycle
do i_fit = 1, ng_fit_jast do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_x_2(i_fit) expo_fit = expo_gauss_j_mu_x_2(i_fit)
coef_fit = coef_gauss_j_mu_x_2(i_fit) coef_fit = coef_gauss_j_mu_x_2(i_fit)
! if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle ! if(dabs(coef_fit*int_env*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
tmp += coef_fit * int_fit tmp += coef_fit * int_fit
enddo enddo
@ -275,8 +275,8 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
coef = List_comb_thr_b3_coef (i_1s,j,i) coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) int_env = ao_abs_comb_b3_env(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle ! if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -286,13 +286,13 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
coef_fit = coef_gauss_j_mu_x_2(i_fit) coef_fit = coef_gauss_j_mu_x_2(i_fit)
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
! if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle ! if(dabs(coef_fit*coef*factor_ij_1s*int_env*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
tmp += coef * coef_fit * int_fit tmp += coef * coef_fit * int_fit
enddo enddo
enddo enddo
int2_u2_j1b2_test(j,i,ipoint) = tmp int2_u2_env2_test(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -302,23 +302,23 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
int2_u2_j1b2_test(j,i,ipoint) = int2_u2_j1b2_test(i,j,ipoint) int2_u2_env2_test(j,i,ipoint) = int2_u2_env2_test(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_u2_j1b2_test', wall1 - wall0 print*, ' wall time for int2_u2_env2_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n_points_final_grid, 3)] BEGIN_PROVIDER [double precision, int2_u_grad1u_x_env2_test, (ao_num,ao_num,n_points_final_grid,3)]
BEGIN_DOC BEGIN_DOC
! !
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu] r2 ! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 u_12^mu [\grad_1 u_12^mu] r2
! !
END_DOC END_DOC
@ -327,27 +327,27 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
double precision :: r(3), int_fit(3), expo_fit, coef_fit double precision :: r(3), int_fit(3), expo_fit, coef_fit
double precision :: coef, beta, B_center(3), dist double precision :: coef, beta, B_center(3), dist
double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coef_tmp double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coef_tmp
double precision :: tmp_x, tmp_y, tmp_z, int_j1b double precision :: tmp_x, tmp_y, tmp_z, int_env
double precision :: wall0, wall1, sq_pi_3_2,sq_alpha double precision :: wall0, wall1, sq_pi_3_2,sq_alpha
print*, ' providing int2_u_grad1u_x_j1b2_test ...' print*, ' providing int2_u_grad1u_x_env2_test ...'
sq_pi_3_2 = dacos(-1.D0)**(1.d0) sq_pi_3_2 = dacos(-1.D0)**(1.d0)
provide mu_erf final_grid_points j1b_pen provide mu_erf final_grid_points
call wall_time(wall0) call wall_time(wall0)
int2_u_grad1u_x_j1b2_test = 0.d0 int2_u_grad1u_x_env2_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, alpha_1s, dist, & !$OMP coef_fit, expo_fit, int_fit, alpha_1s, dist, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp, & !$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp, &
!$OMP tmp_x, tmp_y, tmp_z,int_j1b,sq_alpha) & !$OMP tmp_x, tmp_y, tmp_z,int_env,sq_alpha) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, &
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, & !$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b,sq_pi_3_2,thrsh_cycle_tc) !$OMP List_comb_thr_b3_cent, int2_u_grad1u_x_env2_test,ao_abs_comb_b3_env,sq_pi_3_2,thrsh_cycle_tc)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -365,8 +365,8 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
coef = List_comb_thr_b3_coef (i_1s,j,i) coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) int_env = ao_abs_comb_b3_env(i_1s,j,i)
if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -389,7 +389,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
coef_tmp = coef * coef_fit * dexp(-expo_coef_1s) coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
sq_alpha = alpha_1s_inv * dsqrt(alpha_1s_inv) sq_alpha = alpha_1s_inv * dsqrt(alpha_1s_inv)
! if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .lt. thrsh_cycle_tc) cycle ! if(dabs(coef_tmp*int_env*sq_pi_3_2*sq_alpha) .lt. thrsh_cycle_tc) cycle
call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit) call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)
@ -402,9 +402,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
enddo enddo
int2_u_grad1u_x_j1b2_test(j,i,ipoint,1) = tmp_x int2_u_grad1u_x_env2_test(j,i,ipoint,1) = tmp_x
int2_u_grad1u_x_j1b2_test(j,i,ipoint,2) = tmp_y int2_u_grad1u_x_env2_test(j,i,ipoint,2) = tmp_y
int2_u_grad1u_x_j1b2_test(j,i,ipoint,3) = tmp_z int2_u_grad1u_x_env2_test(j,i,ipoint,3) = tmp_z
enddo enddo
enddo enddo
enddo enddo
@ -414,24 +414,25 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
int2_u_grad1u_x_j1b2_test(j,i,ipoint,1) = int2_u_grad1u_x_j1b2_test(i,j,ipoint,1) int2_u_grad1u_x_env2_test(j,i,ipoint,1) = int2_u_grad1u_x_env2_test(i,j,ipoint,1)
int2_u_grad1u_x_j1b2_test(j,i,ipoint,2) = int2_u_grad1u_x_j1b2_test(i,j,ipoint,2) int2_u_grad1u_x_env2_test(j,i,ipoint,2) = int2_u_grad1u_x_env2_test(i,j,ipoint,2)
int2_u_grad1u_x_j1b2_test(j,i,ipoint,3) = int2_u_grad1u_x_j1b2_test(i,j,ipoint,3) int2_u_grad1u_x_env2_test(j,i,ipoint,3) = int2_u_grad1u_x_env2_test(i,j,ipoint,3)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_u_grad1u_x_j1b2_test', wall1 - wall0 print*, ' wall time for int2_u_grad1u_x_env2_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, int2_u_grad1u_env2_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu] ! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 u_12^mu [\grad_1 u_12^mu]
! !
END_DOC END_DOC
@ -442,31 +443,31 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, tmp double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, tmp
double precision :: wall0, wall1 double precision :: wall0, wall1
double precision, external :: NAI_pol_mult_erf_ao_with1s double precision, external :: NAI_pol_mult_erf_ao_with1s
double precision :: j12_mu_r12,int_j1b double precision :: j12_mu_r12,int_env
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2 double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2
double precision :: beta_ij,center_ij_1s(3),factor_ij_1s double precision :: beta_ij,center_ij_1s(3),factor_ij_1s
print*, ' providing int2_u_grad1u_j1b2_test ...' print*, ' providing int2_u_grad1u_env2_test ...'
dsqpi_3_2 = (dacos(-1.d0))**(1.5d0) dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points j1b_pen ao_overlap_abs List_comb_thr_b3_cent provide mu_erf final_grid_points ao_overlap_abs List_comb_thr_b3_cent
call wall_time(wall0) call wall_time(wall0)
int2_u_grad1u_j1b2_test = 0.d0 int2_u_grad1u_env2_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, & !$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, &
!$OMP beta_ij,center_ij_1s,factor_ij_1s, & !$OMP beta_ij,center_ij_1s,factor_ij_1s, &
!$OMP int_j1b,alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) & !$OMP int_env,alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, &
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, & !$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2, & !$OMP ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, ao_abs_comb_b3_j1b, & !$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, ao_abs_comb_b3_env, &
!$OMP List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test,thrsh_cycle_tc) !$OMP List_comb_thr_b3_cent, int2_u_grad1u_env2_test,thrsh_cycle_tc)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 1, ao_num do i = 1, ao_num
@ -484,11 +485,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
! i_1s = 1 ! i_1s = 1
! --- --- --- ! --- --- ---
int_j1b = ao_abs_comb_b3_j1b(1,j,i) int_env = ao_abs_comb_b3_env(1,j,i)
! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
do i_fit = 1, ng_fit_jast do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_1_erf(i_fit) expo_fit = expo_gauss_j_mu_1_erf(i_fit)
! if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.thrsh_cycle_tc) cycle
coef_fit = coef_gauss_j_mu_1_erf(i_fit) coef_fit = coef_gauss_j_mu_1_erf(i_fit)
int_fit = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r) int_fit = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r)
tmp += coef_fit * int_fit tmp += coef_fit * int_fit
@ -502,8 +501,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
coef = List_comb_thr_b3_coef (i_1s,j,i) coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i) beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i) int_env = ao_abs_comb_b3_env(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -513,7 +511,6 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
do i_fit = 1, ng_fit_jast do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_1_erf(i_fit) expo_fit = expo_gauss_j_mu_1_erf(i_fit)
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s) call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
! if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.thrsh_cycle_tc)cycle
coef_fit = coef_gauss_j_mu_1_erf(i_fit) coef_fit = coef_gauss_j_mu_1_erf(i_fit)
alpha_1s = beta + expo_fit alpha_1s = beta + expo_fit
@ -533,7 +530,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
enddo enddo
enddo enddo
int2_u_grad1u_j1b2_test(j,i,ipoint) = tmp int2_u_grad1u_env2_test(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -543,14 +540,15 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
int2_u_grad1u_j1b2_test(j,i,ipoint) = int2_u_grad1u_j1b2_test(i,j,ipoint) int2_u_grad1u_env2_test(j,i,ipoint) = int2_u_grad1u_env2_test(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_u_grad1u_j1b2_test', wall1 - wall0 print*, ' wall time for int2_u_grad1u_env2_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---

203
src/ao_many_one_e_ints/grad2_jmu_modif.irp.f → plugins/local/ao_many_one_e_ints/grad2_jmu_modif.irp.f
View File

@ -6,7 +6,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_poin
BEGIN_DOC BEGIN_DOC
! !
! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) [1 - erf(mu r12)]^2 ! \frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) [1 - erf(mu r12)]^2
! !
END_DOC END_DOC
@ -21,7 +21,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_poin
print*, ' providing int2_grad1u2_grad2u2 ...' print*, ' providing int2_grad1u2_grad2u2 ...'
call wall_time(wall0) call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen provide mu_erf
provide final_grid_points
int2_grad1u2_grad2u2 = 0.d0 int2_grad1u2_grad2u2 = 0.d0
@ -44,7 +45,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_poin
expo_fit = expo_gauss_1_erf_x_2(i_fit) expo_fit = expo_gauss_1_erf_x_2(i_fit)
coef_fit = coef_gauss_1_erf_x_2(i_fit) coef_fit = coef_gauss_1_erf_x_2(i_fit)
tmp += -0.25d0 * coef_fit * overlap_gauss_r12_ao(r, expo_fit, i, j) tmp += 0.25d0 * coef_fit * overlap_gauss_r12_ao(r, expo_fit, i, j)
enddo enddo
int2_grad1u2_grad2u2(j,i,ipoint) = tmp int2_grad1u2_grad2u2(j,i,ipoint) = tmp
@ -63,17 +64,17 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_poin
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_grad1u2_grad2u2 =', wall1 - wall0 print*, ' wall time for int2_grad1u2_grad2u2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, int2_grad1u2_grad2u2_env2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2 ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 [1 - erf(mu r12)]^2
! !
END_DOC END_DOC
@ -87,21 +88,22 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing int2_grad1u2_grad2u2_j1b2 ...' print*, ' providing int2_grad1u2_grad2u2_env2 ...'
call wall_time(wall0) call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen provide mu_erf
provide final_grid_points
int2_grad1u2_grad2u2_j1b2 = 0.d0 int2_grad1u2_grad2u2_env2 = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp) & !$OMP coef_fit, expo_fit, int_fit, tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_env1s_square_size, &
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, & !$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & !$OMP List_env1s_square_coef, List_env1s_square_expo, &
!$OMP List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2) !$OMP List_env1s_square_cent, int2_grad1u2_grad2u2_env2)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint) r(1) = final_grid_points(1,ipoint)
@ -125,14 +127,14 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
! --- ! ---
do i_1s = 2, List_all_comb_b3_size do i_1s = 2, List_env1s_square_size
coef = List_all_comb_b3_coef (i_1s) coef = List_env1s_square_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0 if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b3_expo (i_1s) beta = List_env1s_square_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s) B_center(1) = List_env1s_square_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s) B_center(2) = List_env1s_square_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s) B_center(3) = List_env1s_square_cent(3,i_1s)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
@ -143,7 +145,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
enddo enddo
int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = tmp int2_grad1u2_grad2u2_env2(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -153,23 +155,23 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2(i,j,ipoint) int2_grad1u2_grad2u2_env2(j,i,ipoint) = int2_grad1u2_grad2u2_env2(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_grad1u2_grad2u2_j1b2 =', wall1 - wall0 print*, ' wall time for int2_grad1u2_grad2u2_env2 (min) =', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, int2_u2_env2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [u_12^mu]^2 ! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 [u_12^mu]^2
! !
END_DOC END_DOC
@ -182,21 +184,22 @@ BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_
double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing int2_u2_j1b2 ...' print*, ' providing int2_u2_env2 ...'
call wall_time(wall0) call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen provide mu_erf
provide final_grid_points
int2_u2_j1b2 = 0.d0 int2_u2_env2 = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp) & !$OMP coef_fit, expo_fit, int_fit, tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_env1s_square_size, &
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, & !$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & !$OMP List_env1s_square_coef, List_env1s_square_expo, &
!$OMP List_all_comb_b3_cent, int2_u2_j1b2) !$OMP List_env1s_square_cent, int2_u2_env2)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint) r(1) = final_grid_points(1,ipoint)
@ -220,14 +223,14 @@ BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_
! --- ! ---
do i_1s = 2, List_all_comb_b3_size do i_1s = 2, List_env1s_square_size
coef = List_all_comb_b3_coef (i_1s) coef = List_env1s_square_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0 if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b3_expo (i_1s) beta = List_env1s_square_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s) B_center(1) = List_env1s_square_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s) B_center(2) = List_env1s_square_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s) B_center(3) = List_env1s_square_cent(3,i_1s)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
@ -238,7 +241,7 @@ BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_
enddo enddo
int2_u2_j1b2(j,i,ipoint) = tmp int2_u2_env2(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -248,23 +251,23 @@ BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
int2_u2_j1b2(j,i,ipoint) = int2_u2_j1b2(i,j,ipoint) int2_u2_env2(j,i,ipoint) = int2_u2_env2(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_u2_j1b2', wall1 - wall0 print*, ' wall time for int2_u2_env2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_points_final_grid, 3)] BEGIN_PROVIDER [double precision, int2_u_grad1u_x_env2, (ao_num, ao_num, n_points_final_grid, 3)]
BEGIN_DOC BEGIN_DOC
! !
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu] r2 ! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 u_12^mu [\grad_1 u_12^mu] r2
! !
END_DOC END_DOC
@ -276,23 +279,24 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_poin
double precision :: tmp_x, tmp_y, tmp_z double precision :: tmp_x, tmp_y, tmp_z
double precision :: wall0, wall1 double precision :: wall0, wall1
print*, ' providing int2_u_grad1u_x_j1b2 ...' print*, ' providing int2_u_grad1u_x_env2 ...'
call wall_time(wall0) call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen provide mu_erf
provide final_grid_points
int2_u_grad1u_x_j1b2 = 0.d0 int2_u_grad1u_x_env2 = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, alpha_1s, dist, & !$OMP coef_fit, expo_fit, int_fit, alpha_1s, dist, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp, & !$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp, &
!$OMP tmp_x, tmp_y, tmp_z) & !$OMP tmp_x, tmp_y, tmp_z) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_env1s_square_size, &
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, & !$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & !$OMP List_env1s_square_coef, List_env1s_square_expo, &
!$OMP List_all_comb_b3_cent, int2_u_grad1u_x_j1b2) !$OMP List_env1s_square_cent, int2_u_grad1u_x_env2)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -321,14 +325,14 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_poin
! --- ! ---
do i_1s = 2, List_all_comb_b3_size do i_1s = 2, List_env1s_square_size
coef = List_all_comb_b3_coef (i_1s) coef = List_env1s_square_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0 if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b3_expo (i_1s) beta = List_env1s_square_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s) B_center(1) = List_env1s_square_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s) B_center(2) = List_env1s_square_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s) B_center(3) = List_env1s_square_cent(3,i_1s)
dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) & dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) &
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) & + (B_center(2) - r(2)) * (B_center(2) - r(2)) &
+ (B_center(3) - r(3)) * (B_center(3) - r(3)) + (B_center(3) - r(3)) * (B_center(3) - r(3))
@ -355,9 +359,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_poin
enddo enddo
int2_u_grad1u_x_j1b2(j,i,ipoint,1) = tmp_x int2_u_grad1u_x_env2(j,i,ipoint,1) = tmp_x
int2_u_grad1u_x_j1b2(j,i,ipoint,2) = tmp_y int2_u_grad1u_x_env2(j,i,ipoint,2) = tmp_y
int2_u_grad1u_x_j1b2(j,i,ipoint,3) = tmp_z int2_u_grad1u_x_env2(j,i,ipoint,3) = tmp_z
enddo enddo
enddo enddo
enddo enddo
@ -367,25 +371,25 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_poin
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
int2_u_grad1u_x_j1b2(j,i,ipoint,1) = int2_u_grad1u_x_j1b2(i,j,ipoint,1) int2_u_grad1u_x_env2(j,i,ipoint,1) = int2_u_grad1u_x_env2(i,j,ipoint,1)
int2_u_grad1u_x_j1b2(j,i,ipoint,2) = int2_u_grad1u_x_j1b2(i,j,ipoint,2) int2_u_grad1u_x_env2(j,i,ipoint,2) = int2_u_grad1u_x_env2(i,j,ipoint,2)
int2_u_grad1u_x_j1b2(j,i,ipoint,3) = int2_u_grad1u_x_j1b2(i,j,ipoint,3) int2_u_grad1u_x_env2(j,i,ipoint,3) = int2_u_grad1u_x_env2(i,j,ipoint,3)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_u_grad1u_x_j1b2 = ', wall1 - wall0 print*, ' wall time for int2_u_grad1u_x_env2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [ double precision, int2_u_grad1u_env2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu] ! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 u_12^mu [\grad_1 u_12^mu]
! !
END_DOC END_DOC
@ -397,22 +401,23 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
double precision :: wall0, wall1 double precision :: wall0, wall1
double precision, external :: NAI_pol_mult_erf_ao_with1s double precision, external :: NAI_pol_mult_erf_ao_with1s
print*, ' providing int2_u_grad1u_j1b2 ...' print*, ' providing int2_u_grad1u_env2 ...'
call wall_time(wall0) call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen provide mu_erf
provide final_grid_points
int2_u_grad1u_j1b2 = 0.d0 int2_u_grad1u_env2 = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, & !$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) & !$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_env1s_square_size, &
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, & !$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, & !$OMP List_env1s_square_coef, List_env1s_square_expo, &
!$OMP List_all_comb_b3_cent, int2_u_grad1u_j1b2) !$OMP List_env1s_square_cent, int2_u_grad1u_env2)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 1, ao_num do i = 1, ao_num
@ -436,14 +441,14 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
! --- ! ---
do i_1s = 2, List_all_comb_b3_size do i_1s = 2, List_env1s_square_size
coef = List_all_comb_b3_coef (i_1s) coef = List_env1s_square_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0 if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b3_expo (i_1s) beta = List_env1s_square_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s) B_center(1) = List_env1s_square_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s) B_center(2) = List_env1s_square_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s) B_center(3) = List_env1s_square_cent(3,i_1s)
dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) & dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) &
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) & + (B_center(2) - r(2)) * (B_center(2) - r(2)) &
+ (B_center(3) - r(3)) * (B_center(3) - r(3)) + (B_center(3) - r(3)) * (B_center(3) - r(3))
@ -468,7 +473,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
enddo enddo
int2_u_grad1u_j1b2(j,i,ipoint) = tmp int2_u_grad1u_env2(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -478,13 +483,13 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
int2_u_grad1u_j1b2(j,i,ipoint) = int2_u_grad1u_j1b2(i,j,ipoint) int2_u_grad1u_env2(j,i,ipoint) = int2_u_grad1u_env2(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for int2_u_grad1u_j1b2', wall1 - wall0 print*, ' wall time for int2_u_grad1u_env2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER

115
src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f → plugins/local/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f
View File

@ -1,11 +1,11 @@
! --- ! ---
BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, v_ij_erf_rk_cst_mu_env_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! int dr phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R| - 1) / |r - R| ! int dr phi_i(r) phi_j(r) 1s_env(r) (erf(mu(R) |r - R| - 1) / |r - R|
! !
END_DOC END_DOC
@ -13,24 +13,23 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
integer :: i, j, ipoint, i_1s integer :: i, j, ipoint, i_1s
double precision :: r(3), int_mu, int_coulomb double precision :: r(3), int_mu, int_coulomb
double precision :: coef, beta, B_center(3) double precision :: coef, beta, B_center(3)
double precision :: tmp,int_j1b double precision :: tmp,int_env
double precision :: wall0, wall1 double precision :: wall0, wall1
double precision, external :: NAI_pol_mult_erf_ao_with1s double precision, external :: NAI_pol_mult_erf_ao_with1s
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2 double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2
print*, ' providing v_ij_erf_rk_cst_mu_j1b_test ...' print*, ' providing v_ij_erf_rk_cst_mu_env_test ...'
dsqpi_3_2 = (dacos(-1.d0))**(1.5d0) dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0) call wall_time(wall0)
v_ij_erf_rk_cst_mu_j1b_test = 0.d0 v_ij_erf_rk_cst_mu_env_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp, int_j1b)& !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp, int_env)&
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points, & !$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_j1b, & !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_env, &
!$OMP v_ij_erf_rk_cst_mu_j1b_test, mu_erf, & !$OMP v_ij_erf_rk_cst_mu_env_test, mu_erf, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc) !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
!$OMP DO !$OMP DO
!do ipoint = 1, 10 !do ipoint = 1, 10
@ -48,8 +47,8 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
coef = List_comb_thr_b2_coef (i_1s,j,i) coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) int_env = ao_abs_comb_b2_env(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle ! if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@ -60,7 +59,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
tmp += coef * (int_mu - int_coulomb) tmp += coef * (int_mu - int_coulomb)
enddo enddo
v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) = tmp v_ij_erf_rk_cst_mu_env_test(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -70,22 +69,22 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) = v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint) v_ij_erf_rk_cst_mu_env_test(j,i,ipoint) = v_ij_erf_rk_cst_mu_env_test(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for v_ij_erf_rk_cst_mu_j1b_test', wall1 - wall0 print*, ' wall time for v_ij_erf_rk_cst_mu_env_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid, 3)] BEGIN_PROVIDER [double precision, x_v_ij_erf_rk_cst_mu_env_test, (ao_num, ao_num, n_points_final_grid, 3)]
BEGIN_DOC BEGIN_DOC
! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R| ! int dr x phi_i(r) phi_j(r) 1s_env(r) (erf(mu(R) |r - R|) - 1)/|r - R|
END_DOC END_DOC
implicit none implicit none
@ -93,23 +92,23 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3) double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3)
double precision :: tmp_x, tmp_y, tmp_z double precision :: tmp_x, tmp_y, tmp_z
double precision :: wall0, wall1 double precision :: wall0, wall1
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b,factor_ij_1s,beta_ij,center_ij_1s double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_env,factor_ij_1s,beta_ij,center_ij_1s
print*, ' providing x_v_ij_erf_rk_cst_mu_j1b_test ...' print*, ' providing x_v_ij_erf_rk_cst_mu_env_test ...'
dsqpi_3_2 = (dacos(-1.d0))**(1.5d0) dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
provide expo_erfc_mu_gauss ao_prod_sigma ao_prod_center provide expo_erfc_mu_gauss ao_prod_sigma ao_prod_center
call wall_time(wall0) call wall_time(wall0)
x_v_ij_erf_rk_cst_mu_j1b_test = 0.d0 x_v_ij_erf_rk_cst_mu_env_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, & !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, &
!$OMP int_j1b, tmp_x, tmp_y, tmp_z,factor_ij_1s,beta_ij,center_ij_1s) & !$OMP int_env, tmp_x, tmp_y, tmp_z,factor_ij_1s,beta_ij,center_ij_1s) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points,& !$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points,&
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent, & !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent, &
!$OMP x_v_ij_erf_rk_cst_mu_j1b_test, mu_erf,ao_abs_comb_b2_j1b, & !$OMP x_v_ij_erf_rk_cst_mu_env_test, mu_erf,ao_abs_comb_b2_env, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,thrsh_cycle_tc) !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,thrsh_cycle_tc)
! !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,expo_erfc_mu_gauss) ! !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,expo_erfc_mu_gauss)
!$OMP DO !$OMP DO
@ -129,8 +128,8 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
coef = List_comb_thr_b2_coef (i_1s,j,i) coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) int_env = ao_abs_comb_b2_env(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle ! if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@ -143,9 +142,9 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
tmp_z += coef * (ints(3) - ints_coulomb(3)) tmp_z += coef * (ints(3) - ints_coulomb(3))
enddo enddo
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,1) = tmp_x x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,1) = tmp_x
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,2) = tmp_y x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,2) = tmp_y
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,3) = tmp_z x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,3) = tmp_z
enddo enddo
enddo enddo
enddo enddo
@ -155,26 +154,26 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,1) x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_env_test(i,j,ipoint,1)
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,2) x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_env_test(i,j,ipoint,2)
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,3) x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_env_test(i,j,ipoint,3)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for x_v_ij_erf_rk_cst_mu_j1b_test', wall1 - wall0 print*, ' wall time for x_v_ij_erf_rk_cst_mu_env_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
! TODO analytically ! TODO analytically
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_env_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) ! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2) u(mu, r12)
! !
END_DOC END_DOC
@ -185,29 +184,28 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
double precision :: tmp double precision :: tmp
double precision :: wall0, wall1 double precision :: wall0, wall1
double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot
double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_j1b double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_env
double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing v_ij_u_cst_mu_j1b_test ...' print*, ' providing v_ij_u_cst_mu_env_test ...'
dsqpi_3_2 = (dacos(-1.d0))**(1.5d0) dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0) call wall_time(wall0)
v_ij_u_cst_mu_j1b_test = 0.d0 v_ij_u_cst_mu_env_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, & !$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, &
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) & !$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_env) &
!$OMP SHARED (n_points_final_grid, ao_num, & !$OMP SHARED (n_points_final_grid, ao_num, &
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, &
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, & !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_env_test,ao_abs_comb_b2_env, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc) !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -225,8 +223,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
! i_1s = 1 ! i_1s = 1
! --- --- --- ! --- --- ---
int_j1b = ao_abs_comb_b2_j1b(1,j,i) int_env = ao_abs_comb_b2_env(1,j,i)
! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle ! if(dabs(int_env).lt.thrsh_cycle_tc) cycle
do i_fit = 1, ng_fit_jast do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_x(i_fit) expo_fit = expo_gauss_j_mu_x(i_fit)
coef_fit = coef_gauss_j_mu_x(i_fit) coef_fit = coef_gauss_j_mu_x(i_fit)
@ -242,8 +240,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
do i_1s = 2, List_comb_thr_b2_size(j,i) do i_1s = 2, List_comb_thr_b2_size(j,i)
coef = List_comb_thr_b2_coef (i_1s,j,i) coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) int_env = ao_abs_comb_b2_env(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle ! if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@ -259,7 +257,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
enddo enddo
enddo enddo
v_ij_u_cst_mu_j1b_test(j,i,ipoint) = tmp v_ij_u_cst_mu_env_test(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -269,23 +267,23 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
v_ij_u_cst_mu_j1b_test(j,i,ipoint) = v_ij_u_cst_mu_j1b_test(i,j,ipoint) v_ij_u_cst_mu_env_test(j,i,ipoint) = v_ij_u_cst_mu_env_test(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_test', wall1 - wall0 print*, ' wall time for v_ij_u_cst_mu_env_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_env_ng_1_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) with u(mu,r12) \approx 1/2 mu e^{-2.5 * mu (r12)^2} ! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2) u(mu, r12) with u(mu,r12) \approx 1/2 mu e^{-2.5 * mu (r12)^2}
! !
END_DOC END_DOC
@ -296,27 +294,26 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
double precision :: tmp double precision :: tmp
double precision :: wall0, wall1 double precision :: wall0, wall1
double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot
double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_j1b double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_env
double precision, external :: overlap_gauss_r12_ao double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
dsqpi_3_2 = (dacos(-1.d0))**(1.5d0) dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0) call wall_time(wall0)
v_ij_u_cst_mu_j1b_ng_1_test = 0.d0 v_ij_u_cst_mu_env_ng_1_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
!$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, & !$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, &
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) & !$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_env) &
!$OMP SHARED (n_points_final_grid, ao_num, & !$OMP SHARED (n_points_final_grid, ao_num, &
!$OMP final_grid_points, expo_good_j_mu_1gauss,coef_good_j_mu_1gauss, & !$OMP final_grid_points, expo_good_j_mu_1gauss,coef_good_j_mu_1gauss, &
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & !$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, &
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, & !$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_env_ng_1_test,ao_abs_comb_b2_env, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc) !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -334,8 +331,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
! i_1s = 1 ! i_1s = 1
! --- --- --- ! --- --- ---
int_j1b = ao_abs_comb_b2_j1b(1,j,i) int_env = ao_abs_comb_b2_env(1,j,i)
! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle ! if(dabs(int_env).lt.thrsh_cycle_tc) cycle
expo_fit = expo_good_j_mu_1gauss expo_fit = expo_good_j_mu_1gauss
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
tmp += int_fit tmp += int_fit
@ -347,8 +344,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
do i_1s = 2, List_comb_thr_b2_size(j,i) do i_1s = 2, List_comb_thr_b2_size(j,i)
coef = List_comb_thr_b2_coef (i_1s,j,i) coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i) beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i) int_env = ao_abs_comb_b2_env(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle ! if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i) B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i) B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i) B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@ -364,7 +361,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
! enddo ! enddo
enddo enddo
v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint) = tmp v_ij_u_cst_mu_env_ng_1_test(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -374,13 +371,13 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint) = v_ij_u_cst_mu_j1b_ng_1_test(i,j,ipoint) v_ij_u_cst_mu_env_ng_1_test(j,i,ipoint) = v_ij_u_cst_mu_env_ng_1_test(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_ng_1_test', wall1 - wall0 print*, ' wall time for v_ij_u_cst_mu_env_ng_1_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER

237
src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f → plugins/local/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f
View File

@ -1,11 +1,11 @@
! --- ! ---
BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, v_ij_erf_rk_cst_mu_env, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! int dr phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R| - 1) / |r - R| ! int dr phi_i(r) phi_j(r) 1s_env(r) (erf(mu(R) |r - R| - 1) / |r - R|
! !
END_DOC END_DOC
@ -17,18 +17,20 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
double precision :: wall0, wall1 double precision :: wall0, wall1
double precision, external :: NAI_pol_mult_erf_ao_with1s double precision, external :: NAI_pol_mult_erf_ao_with1s
print *, ' providing v_ij_erf_rk_cst_mu_j1b ...' PROVIDE mu_erf
PROVIDE final_grid_points
PROVIDE env_expo
print *, ' providing v_ij_erf_rk_cst_mu_env ...'
call wall_time(wall0) call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen v_ij_erf_rk_cst_mu_env = 0.d0
v_ij_erf_rk_cst_mu_j1b = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp) & !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points, & !$OMP SHARED (n_points_final_grid, ao_num, List_env1s_size, final_grid_points, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & !$OMP List_env1s_coef, List_env1s_expo, List_env1s_cent, &
!$OMP v_ij_erf_rk_cst_mu_j1b, mu_erf) !$OMP v_ij_erf_rk_cst_mu_env, mu_erf)
!$OMP DO !$OMP DO
!do ipoint = 1, 10 !do ipoint = 1, 10
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -43,28 +45,27 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
! --- ! ---
coef = List_all_comb_b2_coef (1) coef = List_env1s_coef (1)
beta = List_all_comb_b2_expo (1) beta = List_env1s_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1) B_center(1) = List_env1s_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1) B_center(2) = List_env1s_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1) B_center(3) = List_env1s_cent(3,1)
int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r) int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r) int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
! if(dabs(coef)*dabs(int_mu - int_coulomb) .lt. 1d-12) cycle
tmp += coef * (int_mu - int_coulomb) tmp += coef * (int_mu - int_coulomb)
! --- ! ---
do i_1s = 2, List_all_comb_b2_size do i_1s = 2, List_env1s_size
coef = List_all_comb_b2_coef (i_1s) coef = List_env1s_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0 if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b2_expo (i_1s) beta = List_env1s_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s) B_center(1) = List_env1s_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s) B_center(2) = List_env1s_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s) B_center(3) = List_env1s_cent(3,i_1s)
int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r) int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r) int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
@ -74,7 +75,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
! --- ! ---
v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) = tmp v_ij_erf_rk_cst_mu_env(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -84,22 +85,22 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) = v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) v_ij_erf_rk_cst_mu_env(j,i,ipoint) = v_ij_erf_rk_cst_mu_env(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for v_ij_erf_rk_cst_mu_j1b', wall1 - wall0 print*, ' wall time for v_ij_erf_rk_cst_mu_env (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid, 3)] BEGIN_PROVIDER [double precision, x_v_ij_erf_rk_cst_mu_env, (ao_num, ao_num, n_points_final_grid, 3)]
BEGIN_DOC BEGIN_DOC
! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R| ! int dr x phi_i(r) phi_j(r) 1s_env(r) (erf(mu(R) |r - R|) - 1)/|r - R|
END_DOC END_DOC
implicit none implicit none
@ -108,17 +109,17 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_
double precision :: tmp_x, tmp_y, tmp_z double precision :: tmp_x, tmp_y, tmp_z
double precision :: wall0, wall1 double precision :: wall0, wall1
print*, ' providing x_v_ij_erf_rk_cst_mu_j1b ...' print*, ' providing x_v_ij_erf_rk_cst_mu_env ...'
call wall_time(wall0) call wall_time(wall0)
x_v_ij_erf_rk_cst_mu_j1b = 0.d0 x_v_ij_erf_rk_cst_mu_env = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, & !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, &
!$OMP tmp_x, tmp_y, tmp_z) & !$OMP tmp_x, tmp_y, tmp_z) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points,& !$OMP SHARED (n_points_final_grid, ao_num, List_env1s_size, final_grid_points,&
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, & !$OMP List_env1s_coef, List_env1s_expo, List_env1s_cent, &
!$OMP x_v_ij_erf_rk_cst_mu_j1b, mu_erf) !$OMP x_v_ij_erf_rk_cst_mu_env, mu_erf)
!$OMP DO !$OMP DO
!do ipoint = 1, 10 !do ipoint = 1, 10
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -135,11 +136,11 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_
! --- ! ---
coef = List_all_comb_b2_coef (1) coef = List_env1s_coef (1)
beta = List_all_comb_b2_expo (1) beta = List_env1s_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1) B_center(1) = List_env1s_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1) B_center(2) = List_env1s_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1) B_center(3) = List_env1s_cent(3,1)
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints ) call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints )
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb) call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb)
@ -152,14 +153,14 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_
! --- ! ---
do i_1s = 2, List_all_comb_b2_size do i_1s = 2, List_env1s_size
coef = List_all_comb_b2_coef (i_1s) coef = List_env1s_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0 if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b2_expo (i_1s) beta = List_env1s_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s) B_center(1) = List_env1s_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s) B_center(2) = List_env1s_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s) B_center(3) = List_env1s_cent(3,i_1s)
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints ) call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints )
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb) call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb)
@ -171,9 +172,9 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_
! --- ! ---
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,1) = tmp_x x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,1) = tmp_x
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,2) = tmp_y x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,2) = tmp_y
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,3) = tmp_z x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,3) = tmp_z
enddo enddo
enddo enddo
enddo enddo
@ -183,25 +184,25 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,1)
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,2)
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,3)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for x_v_ij_erf_rk_cst_mu_j1b =', wall1 - wall0 print*, ' wall time for x_v_ij_erf_rk_cst_mu_env (min) =', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_env_fit, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) ! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2) u(mu, r12)
! !
END_DOC END_DOC
@ -214,23 +215,23 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_poi
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing v_ij_u_cst_mu_j1b_fit ...' print*, ' providing v_ij_u_cst_mu_env_fit ...'
call wall_time(wall0) call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen provide mu_erf final_grid_points env_expo
PROVIDE ng_fit_jast expo_gauss_j_mu_x coef_gauss_j_mu_x PROVIDE ng_fit_jast expo_gauss_j_mu_x coef_gauss_j_mu_x
PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent PROVIDE List_env1s_size List_env1s_coef List_env1s_expo List_env1s_cent
v_ij_u_cst_mu_j1b_fit = 0.d0 v_ij_u_cst_mu_env_fit = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp) & !$OMP coef_fit, expo_fit, int_fit, tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_env1s_size, &
!$OMP final_grid_points, ng_fit_jast, & !$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, & !$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & !$OMP List_env1s_coef, List_env1s_expo, &
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_fit) !$OMP List_env1s_cent, v_ij_u_cst_mu_env_fit)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint) r(1) = final_grid_points(1,ipoint)
@ -247,11 +248,11 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_poi
! --- ! ---
coef = List_all_comb_b2_coef (1) coef = List_env1s_coef (1)
beta = List_all_comb_b2_expo (1) beta = List_env1s_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1) B_center(1) = List_env1s_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1) B_center(2) = List_env1s_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1) B_center(3) = List_env1s_cent(3,1)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
@ -259,14 +260,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_poi
! --- ! ---
do i_1s = 2, List_all_comb_b2_size do i_1s = 2, List_env1s_size
coef = List_all_comb_b2_coef (i_1s) coef = List_env1s_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0 if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b2_expo (i_1s) beta = List_env1s_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s) B_center(1) = List_env1s_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s) B_center(2) = List_env1s_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s) B_center(3) = List_env1s_cent(3,i_1s)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j) int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
@ -277,7 +278,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_poi
enddo enddo
v_ij_u_cst_mu_j1b_fit(j,i,ipoint) = tmp v_ij_u_cst_mu_env_fit(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -287,23 +288,23 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_poi
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
v_ij_u_cst_mu_j1b_fit(j,i,ipoint) = v_ij_u_cst_mu_j1b_fit(i,j,ipoint) v_ij_u_cst_mu_env_fit(j,i,ipoint) = v_ij_u_cst_mu_env_fit(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_fit', wall1 - wall0 print*, ' wall time for v_ij_u_cst_mu_env_fit (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_env_an_old, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) ! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2) u(mu, r12)
! !
END_DOC END_DOC
@ -322,24 +323,24 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_p
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
double precision, external :: NAI_pol_mult_erf_ao_with1s double precision, external :: NAI_pol_mult_erf_ao_with1s
print*, ' providing v_ij_u_cst_mu_j1b_an_old ...' print*, ' providing v_ij_u_cst_mu_env_an_old ...'
call wall_time(wall0) call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen provide mu_erf final_grid_points env_expo
PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent PROVIDE List_env1s_size List_env1s_coef List_env1s_expo List_env1s_cent
ct = inv_sq_pi_2 / mu_erf ct = inv_sq_pi_2 / mu_erf
v_ij_u_cst_mu_j1b_an_old = 0.d0 v_ij_u_cst_mu_env_an_old = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
!$OMP r1_2, tmp, int_c1, int_e1, int_o, int_c2, & !$OMP r1_2, tmp, int_c1, int_e1, int_o, int_c2, &
!$OMP int_e2, int_c3, int_e3) & !$OMP int_e2, int_c3, int_e3) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_env1s_size, &
!$OMP final_grid_points, mu_erf, ct, & !$OMP final_grid_points, mu_erf, ct, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & !$OMP List_env1s_coef, List_env1s_expo, &
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an_old) !$OMP List_env1s_cent, v_ij_u_cst_mu_env_an_old)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -353,11 +354,11 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_p
! --- ! ---
coef = List_all_comb_b2_coef (1) coef = List_env1s_coef (1)
beta = List_all_comb_b2_expo (1) beta = List_env1s_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1) B_center(1) = List_env1s_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1) B_center(2) = List_env1s_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1) B_center(3) = List_env1s_cent(3,1)
int_c1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r) int_c1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
int_e1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r) int_e1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
@ -379,14 +380,14 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_p
! --- ! ---
do i_1s = 2, List_all_comb_b2_size do i_1s = 2, List_env1s_size
coef = List_all_comb_b2_coef (i_1s) coef = List_env1s_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0 if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b2_expo (i_1s) beta = List_env1s_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s) B_center(1) = List_env1s_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s) B_center(2) = List_env1s_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s) B_center(3) = List_env1s_cent(3,i_1s)
int_c1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r) int_c1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
int_e1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r) int_e1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
@ -410,7 +411,7 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_p
! --- ! ---
v_ij_u_cst_mu_j1b_an_old(j,i,ipoint) = tmp v_ij_u_cst_mu_env_an_old(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -420,23 +421,23 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_p
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
v_ij_u_cst_mu_j1b_an_old(j,i,ipoint) = v_ij_u_cst_mu_j1b_an_old(i,j,ipoint) v_ij_u_cst_mu_env_an_old(j,i,ipoint) = v_ij_u_cst_mu_env_an_old(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_an_old', wall1 - wall0 print*, ' wall time for v_ij_u_cst_mu_env_an_old (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_points_final_grid)] BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_env_an, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC BEGIN_DOC
! !
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) ! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2) u(mu, r12)
! !
END_DOC END_DOC
@ -454,23 +455,23 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_point
double precision, external :: overlap_gauss_r12_ao_with1s double precision, external :: overlap_gauss_r12_ao_with1s
double precision, external :: NAI_pol_mult_erf_ao_with1s double precision, external :: NAI_pol_mult_erf_ao_with1s
print*, ' providing v_ij_u_cst_mu_j1b_an ...' print*, ' providing v_ij_u_cst_mu_env_an ...'
call wall_time(wall0) call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen provide mu_erf final_grid_points env_expo
PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent PROVIDE List_env1s_size List_env1s_coef List_env1s_expo List_env1s_cent
ct = inv_sq_pi_2 / mu_erf ct = inv_sq_pi_2 / mu_erf
v_ij_u_cst_mu_j1b_an = 0.d0 v_ij_u_cst_mu_env_an = 0.d0
!$OMP PARALLEL DEFAULT (NONE) & !$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, & !$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
!$OMP r1_2, tmp, int_c, int_e, int_o) & !$OMP r1_2, tmp, int_c, int_e, int_o) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, & !$OMP SHARED (n_points_final_grid, ao_num, List_env1s_size, &
!$OMP final_grid_points, mu_erf, ct, & !$OMP final_grid_points, mu_erf, ct, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, & !$OMP List_env1s_coef, List_env1s_expo, &
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an) !$OMP List_env1s_cent, v_ij_u_cst_mu_env_an)
!$OMP DO !$OMP DO
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -484,11 +485,11 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_point
! --- ! ---
coef = List_all_comb_b2_coef (1) coef = List_env1s_coef (1)
beta = List_all_comb_b2_expo (1) beta = List_env1s_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1) B_center(1) = List_env1s_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1) B_center(2) = List_env1s_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1) B_center(3) = List_env1s_cent(3,1)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c) call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e) call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e)
@ -504,14 +505,14 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_point
! --- ! ---
do i_1s = 2, List_all_comb_b2_size do i_1s = 2, List_env1s_size
coef = List_all_comb_b2_coef (i_1s) coef = List_env1s_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0 if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b2_expo (i_1s) beta = List_env1s_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s) B_center(1) = List_env1s_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s) B_center(2) = List_env1s_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s) B_center(3) = List_env1s_cent(3,i_1s)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c) call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e) call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e)
@ -529,7 +530,7 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_point
! --- ! ---
v_ij_u_cst_mu_j1b_an(j,i,ipoint) = tmp v_ij_u_cst_mu_env_an(j,i,ipoint) = tmp
enddo enddo
enddo enddo
enddo enddo
@ -539,13 +540,13 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_point
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
do i = 2, ao_num do i = 2, ao_num
do j = 1, i-1 do j = 1, i-1
v_ij_u_cst_mu_j1b_an(j,i,ipoint) = v_ij_u_cst_mu_j1b_an(i,j,ipoint) v_ij_u_cst_mu_env_an(j,i,ipoint) = v_ij_u_cst_mu_env_an(i,j,ipoint)
enddo enddo
enddo enddo
enddo enddo
call wall_time(wall1) call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_an', wall1 - wall0 print*, ' wall time for v_ij_u_cst_mu_env_an (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER END_PROVIDER

0
src/ao_many_one_e_ints/grad_related_ints.irp.f → plugins/local/ao_many_one_e_ints/grad_related_ints.irp.f
View File

574
plugins/local/ao_many_one_e_ints/lin_fc_rsdft.irp.f Normal file
View File

@ -0,0 +1,574 @@
! ---
BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du_0, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du_x, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du_y, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du_z, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du_2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! Ir2_Mu_long_Du_0 = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) [(1 - erf(mu r_12) / r_12]
!
! Ir2_Mu_long_Du_x = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) [(1 - erf(mu r_12) / r_12] * x2
! Ir2_Mu_long_Du_y = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) [(1 - erf(mu r_12) / r_12] * y2
! Ir2_Mu_long_Du_z = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) [(1 - erf(mu r_12) / r_12] * z2
!
! Ir2_Mu_long_Du_2 = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) [(1 - erf(mu r_12) / r_12] * r2^2
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s
double precision :: r(3), int_clb(7), int_erf(7)
double precision :: c_1s, e_1s, R_1s(3)
double precision :: tmp_Du_0, tmp_Du_x, tmp_Du_y, tmp_Du_z, tmp_Du_2
double precision :: wall0, wall1
PROVIDE mu_erf
PROVIDE final_grid_points
PROVIDE List_env1s_size List_env1s_expo List_env1s_coef List_env1s_cent
print *, ' providing Ir2_Mu_long_Du ...'
call wall_time(wall0)
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, c_1s, e_1s, R_1s, int_erf, int_clb, &
!$OMP tmp_Du_0, tmp_Du_x, tmp_Du_y, tmp_Du_z, tmp_Du_2) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, mu_erf, &
!$OMP List_env1s_size, List_env1s_expo, &
!$OMP List_env1s_coef, List_env1s_cent, &
!$OMP Ir2_Mu_long_Du_0, Ir2_Mu_long_Du_x, &
!$OMP Ir2_Mu_long_Du_y, Ir2_Mu_long_Du_z, &
!$OMP Ir2_Mu_long_Du_2)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
do i = 1, ao_num
do j = i, ao_num
call NAI_pol_012_mult_erf_ao(i, j, 1.d+9, r, int_clb)
call NAI_pol_012_mult_erf_ao(i, j, mu_erf, r, int_erf)
tmp_Du_0 = int_clb(1) - int_erf(1)
tmp_Du_x = int_clb(2) - int_erf(2)
tmp_Du_y = int_clb(3) - int_erf(3)
tmp_Du_z = int_clb(4) - int_erf(4)
tmp_Du_2 = int_clb(5) + int_clb(6) + int_clb(7) - int_erf(5) - int_erf(6) - int_erf(7)
do i_1s = 2, List_env1s_size
e_1s = List_env1s_expo(i_1s)
c_1s = List_env1s_coef(i_1s)
R_1s(1) = List_env1s_cent(1,i_1s)
R_1s(2) = List_env1s_cent(2,i_1s)
R_1s(3) = List_env1s_cent(3,i_1s)
call NAI_pol_012_mult_erf_ao_with1s(i, j, e_1s, R_1s, 1.d+9, r, int_clb)
call NAI_pol_012_mult_erf_ao_with1s(i, j, e_1s, R_1s, mu_erf, r, int_erf)
tmp_Du_0 = tmp_Du_0 + c_1s * (int_clb(1) - int_erf(1))
tmp_Du_x = tmp_Du_x + c_1s * (int_clb(2) - int_erf(2))
tmp_Du_y = tmp_Du_y + c_1s * (int_clb(3) - int_erf(3))
tmp_Du_z = tmp_Du_z + c_1s * (int_clb(4) - int_erf(4))
tmp_Du_2 = tmp_Du_2 + c_1s * (int_clb(5) + int_clb(6) + int_clb(7) - int_erf(5) - int_erf(6) - int_erf(7))
enddo
Ir2_Mu_long_Du_0(j,i,ipoint) = tmp_Du_0
Ir2_Mu_long_Du_x(j,i,ipoint) = tmp_Du_x
Ir2_Mu_long_Du_y(j,i,ipoint) = tmp_Du_y
Ir2_Mu_long_Du_z(j,i,ipoint) = tmp_Du_z
Ir2_Mu_long_Du_2(j,i,ipoint) = tmp_Du_2
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
Ir2_Mu_long_Du_0(j,i,ipoint) = Ir2_Mu_long_Du_0(i,j,ipoint)
Ir2_Mu_long_Du_x(j,i,ipoint) = Ir2_Mu_long_Du_x(i,j,ipoint)
Ir2_Mu_long_Du_y(j,i,ipoint) = Ir2_Mu_long_Du_y(i,j,ipoint)
Ir2_Mu_long_Du_z(j,i,ipoint) = Ir2_Mu_long_Du_z(i,j,ipoint)
Ir2_Mu_long_Du_2(j,i,ipoint) = Ir2_Mu_long_Du_2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for Ir2_Mu_long_Du (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, Ir2_Mu_gauss_Du, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! Ir2_Mu_gauss_Du = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) e^{-(mu r_12)^2}
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s
double precision :: r(3)
double precision :: coef, beta, B_center(3)
double precision :: tmp_Du
double precision :: mu_sq, dx, dy, dz, tmp_arg, rmu_sq(3)
double precision :: e_1s, c_1s, R_1s(3)
double precision :: wall0, wall1
double precision, external :: overlap_gauss_r12_ao
PROVIDE mu_erf
PROVIDE final_grid_points
PROVIDE List_env1s_size List_env1s_expo List_env1s_coef List_env1s_cent
print *, ' providing Ir2_Mu_gauss_Du ...'
call wall_time(wall0)
mu_sq = mu_erf * mu_erf
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, dx, dy, dz, r, tmp_arg, coef, &
!$OMP rmu_sq, e_1s, c_1s, R_1s, beta, B_center, tmp_Du) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, mu_sq, &
!$OMP List_env1s_size, List_env1s_expo, &
!$OMP List_env1s_coef, List_env1s_cent, &
!$OMP Ir2_Mu_gauss_Du)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
rmu_sq(1) = mu_sq * r(1)
rmu_sq(2) = mu_sq * r(2)
rmu_sq(3) = mu_sq * r(3)
do i = 1, ao_num
do j = i, ao_num
tmp_Du = overlap_gauss_r12_ao(r, mu_sq, j, i)
do i_1s = 2, List_env1s_size
e_1s = List_env1s_expo(i_1s)
c_1s = List_env1s_coef(i_1s)
R_1s(1) = List_env1s_cent(1,i_1s)
R_1s(2) = List_env1s_cent(2,i_1s)
R_1s(3) = List_env1s_cent(3,i_1s)
dx = r(1) - R_1s(1)
dy = r(2) - R_1s(2)
dz = r(3) - R_1s(3)
beta = mu_sq + e_1s
tmp_arg = mu_sq * e_1s * (dx*dx + dy*dy + dz*dz) / beta
coef = c_1s * dexp(-tmp_arg)
B_center(1) = (rmu_sq(1) + e_1s * R_1s(1)) / beta
B_center(2) = (rmu_sq(2) + e_1s * R_1s(2)) / beta
B_center(3) = (rmu_sq(3) + e_1s * R_1s(3)) / beta
tmp_Du += coef * overlap_gauss_r12_ao(B_center, beta, j, i)
enddo
Ir2_Mu_gauss_Du(j,i,ipoint) = tmp_Du
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
Ir2_Mu_gauss_Du(j,i,ipoint) = Ir2_Mu_gauss_Du(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for Ir2_Mu_gauss_Du (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du2_0, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du2_x, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du2_y, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du2_z, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du2_2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! Ir2_Mu_long_Du2_0 = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12) / r_12]
!
! Ir2_Mu_long_Du2_x = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12) / r_12] * x2
! Ir2_Mu_long_Du2_y = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12) / r_12] * y2
! Ir2_Mu_long_Du2_z = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12) / r_12] * z2
!
! Ir2_Mu_long_Du2_2 = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12) / r_12] * r2^2
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s
double precision :: r(3), int_clb(7), int_erf(7)
double precision :: coef, beta, B_center(3)
double precision :: tmp_Du2_0, tmp_Du2_x, tmp_Du2_y, tmp_Du2_z, tmp_Du2_2
double precision :: mu_sq, tmp_arg, dx, dy, dz, rmu_sq(3)
double precision :: e_1s, c_1s, R_1s(3)
double precision :: wall0, wall1
PROVIDE mu_erf
PROVIDE final_grid_points
PROVIDE List_env1s_square_size List_env1s_square_expo List_env1s_square_coef List_env1s_square_cent
print *, ' providing Ir2_Mu_long_Du2 ...'
call wall_time(wall0)
mu_sq = mu_erf * mu_erf
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, rmu_sq, dx, dy, dz, &
!$OMP e_1s, c_1s, R_1s, tmp_arg, coef, beta, B_center, &
!$OMP int_erf, int_clb, &
!$OMP tmp_Du2_0, tmp_Du2_x, tmp_Du2_y, tmp_Du2_z, tmp_Du2_2) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, mu_sq, &
!$OMP mu_erf, List_env1s_square_size, List_env1s_square_expo, &
!$OMP List_env1s_square_coef, List_env1s_square_cent, &
!$OMP Ir2_Mu_long_Du2_0, Ir2_Mu_long_Du2_x, &
!$OMP Ir2_Mu_long_Du2_y, Ir2_Mu_long_Du2_z, &
!$OMP Ir2_Mu_long_Du2_2)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
rmu_sq(1) = mu_sq * r(1)
rmu_sq(2) = mu_sq * r(2)
rmu_sq(3) = mu_sq * r(3)
do i = 1, ao_num
do j = i, ao_num
call NAI_pol_012_mult_erf_ao_with1s(i, j, mu_sq, r, 1.d+9, r, int_clb)
call NAI_pol_012_mult_erf_ao_with1s(i, j, mu_sq, r, mu_erf, r, int_erf)
tmp_Du2_0 = int_clb(1) - int_erf(1)
tmp_Du2_x = int_clb(2) - int_erf(2)
tmp_Du2_y = int_clb(3) - int_erf(3)
tmp_Du2_z = int_clb(4) - int_erf(4)
tmp_Du2_2 = int_clb(5) + int_clb(6) + int_clb(7) - int_erf(5) - int_erf(6) - int_erf(7)
do i_1s = 2, List_env1s_square_size
e_1s = List_env1s_square_expo(i_1s)
c_1s = List_env1s_square_coef(i_1s)
R_1s(1) = List_env1s_square_cent(1,i_1s)
R_1s(2) = List_env1s_square_cent(2,i_1s)
R_1s(3) = List_env1s_square_cent(3,i_1s)
dx = r(1) - R_1s(1)
dy = r(2) - R_1s(2)
dz = r(3) - R_1s(3)
beta = mu_sq + e_1s
tmp_arg = mu_sq * e_1s * (dx*dx + dy*dy + dz*dz) / beta
coef = c_1s * dexp(-tmp_arg)
B_center(1) = (rmu_sq(1) + e_1s * R_1s(1)) / beta
B_center(2) = (rmu_sq(2) + e_1s * R_1s(2)) / beta
B_center(3) = (rmu_sq(3) + e_1s * R_1s(3)) / beta
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_clb)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_erf)
tmp_Du2_0 = tmp_Du2_0 + coef * (int_clb(1) - int_erf(1))
tmp_Du2_x = tmp_Du2_x + coef * (int_clb(2) - int_erf(2))
tmp_Du2_y = tmp_Du2_y + coef * (int_clb(3) - int_erf(3))
tmp_Du2_z = tmp_Du2_z + coef * (int_clb(4) - int_erf(4))
tmp_Du2_2 = tmp_Du2_2 + coef * (int_clb(5) + int_clb(6) + int_clb(7) - int_erf(5) - int_erf(6) - int_erf(7))
enddo
Ir2_Mu_long_Du2_0(j,i,ipoint) = tmp_Du2_0
Ir2_Mu_long_Du2_x(j,i,ipoint) = tmp_Du2_x
Ir2_Mu_long_Du2_y(j,i,ipoint) = tmp_Du2_y
Ir2_Mu_long_Du2_z(j,i,ipoint) = tmp_Du2_z
Ir2_Mu_long_Du2_2(j,i,ipoint) = tmp_Du2_2
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
Ir2_Mu_long_Du2_0(j,i,ipoint) = Ir2_Mu_long_Du2_0(i,j,ipoint)
Ir2_Mu_long_Du2_x(j,i,ipoint) = Ir2_Mu_long_Du2_x(i,j,ipoint)
Ir2_Mu_long_Du2_y(j,i,ipoint) = Ir2_Mu_long_Du2_y(i,j,ipoint)
Ir2_Mu_long_Du2_z(j,i,ipoint) = Ir2_Mu_long_Du2_z(i,j,ipoint)
Ir2_Mu_long_Du2_2(j,i,ipoint) = Ir2_Mu_long_Du2_2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for Ir2_Mu_long_Du2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, Ir2_Mu_gauss_Du2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! Ir2_Mu_gauss_Du2 = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 e^{-(mu r_12)^2}
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s
double precision :: r(3)
double precision :: coef, beta, B_center(3)
double precision :: tmp_Du2
double precision :: mu_sq, dx, dy, dz, tmp_arg, rmu_sq(3)
double precision :: e_1s, c_1s, R_1s(3)
double precision :: wall0, wall1
double precision, external :: overlap_gauss_r12_ao
PROVIDE mu_erf
PROVIDE final_grid_points
PROVIDE List_env1s_square_size List_env1s_square_expo List_env1s_square_coef List_env1s_square_cent
print *, ' providing Ir2_Mu_gauss_Du2 ...'
call wall_time(wall0)
mu_sq = 2.d0 * mu_erf * mu_erf
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, dx, dy, dz, r, tmp_arg, coef, &
!$OMP rmu_sq, e_1s, c_1s, R_1s, beta, B_center, tmp_Du2) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, mu_sq, &
!$OMP List_env1s_square_size, List_env1s_square_expo, &
!$OMP List_env1s_square_coef, List_env1s_square_cent, &
!$OMP Ir2_Mu_gauss_Du2)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
rmu_sq(1) = mu_sq * r(1)
rmu_sq(2) = mu_sq * r(2)
rmu_sq(3) = mu_sq * r(3)
do i = 1, ao_num
do j = i, ao_num
tmp_Du2 = overlap_gauss_r12_ao(r, mu_sq, j, i)
do i_1s = 2, List_env1s_square_size
e_1s = List_env1s_square_expo(i_1s)
c_1s = List_env1s_square_coef(i_1s)
R_1s(1) = List_env1s_square_cent(1,i_1s)
R_1s(2) = List_env1s_square_cent(2,i_1s)
R_1s(3) = List_env1s_square_cent(3,i_1s)
dx = r(1) - R_1s(1)
dy = r(2) - R_1s(2)
dz = r(3) - R_1s(3)
beta = mu_sq + e_1s
tmp_arg = mu_sq * e_1s * (dx*dx + dy*dy + dz*dz) / beta
coef = c_1s * dexp(-tmp_arg)
B_center(1) = (rmu_sq(1) + e_1s * R_1s(1)) / beta
B_center(2) = (rmu_sq(2) + e_1s * R_1s(2)) / beta
B_center(3) = (rmu_sq(3) + e_1s * R_1s(3)) / beta
tmp_Du2 += coef * overlap_gauss_r12_ao(B_center, beta, j, i)
enddo
Ir2_Mu_gauss_Du2(j,i,ipoint) = tmp_Du2
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
Ir2_Mu_gauss_Du2(j,i,ipoint) = Ir2_Mu_gauss_Du2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for Ir2_Mu_gauss_Du2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, Ir2_Mu_short_Du2_0, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_short_Du2_x, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_short_Du2_y, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_short_Du2_z, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_short_Du2_2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! Ir2_Mu_short_Du2_0 = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12)]^2
!
! Ir2_Mu_short_Du2_x = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12)]^2 * x2
! Ir2_Mu_short_Du2_y = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12)]^2 * y2
! Ir2_Mu_short_Du2_z = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12)]^2 * z2
!
! Ir2_Mu_short_Du2_2 = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12)]^2 * r2^2
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s, i_fit
double precision :: r(3), ints(7)
double precision :: coef, beta, B_center(3)
double precision :: tmp_Du2_0, tmp_Du2_x, tmp_Du2_y, tmp_Du2_z, tmp_Du2_2
double precision :: tmp_arg, dx, dy, dz
double precision :: expo_fit, coef_fit, e_1s, c_1s, R_1s(3)
double precision :: wall0, wall1
PROVIDE final_grid_points
PROVIDE List_env1s_square_size List_env1s_square_expo List_env1s_square_coef List_env1s_square_cent
PROVIDE ng_fit_jast expo_gauss_1_erf_x_2 coef_gauss_1_erf_x_2
print *, ' providing Ir2_Mu_short_Du2 ...'
call wall_time(wall0)
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, dx, dy, dz, &
!$OMP expo_fit, coef_fit, e_1s, c_1s, R_1s, &
!$OMP tmp_arg, coef, beta, B_center, ints, &
!$OMP tmp_Du2_0, tmp_Du2_x, tmp_Du2_y, tmp_Du2_z, tmp_Du2_2) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, &
!$OMP ng_fit_jast, expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_env1s_square_size, List_env1s_square_expo, &
!$OMP List_env1s_square_coef, List_env1s_square_cent, &
!$OMP Ir2_Mu_short_Du2_0, Ir2_Mu_short_Du2_x, &
!$OMP Ir2_Mu_short_Du2_y, Ir2_Mu_short_Du2_z, &
!$OMP Ir2_Mu_short_Du2_2)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
do i = 1, ao_num
do j = i, ao_num
tmp_Du2_0 = 0.d0
tmp_Du2_x = 0.d0
tmp_Du2_y = 0.d0
tmp_Du2_z = 0.d0
tmp_Du2_2 = 0.d0
do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_1_erf_x_2(i_fit)
coef_fit = coef_gauss_1_erf_x_2(i_fit)
call overlap_gauss_r12_ao_012(r, expo_fit, i, j, ints)
tmp_Du2_0 += coef_fit * ints(1)
tmp_Du2_x += coef_fit * ints(2)
tmp_Du2_y += coef_fit * ints(3)
tmp_Du2_z += coef_fit * ints(4)
tmp_Du2_2 += coef_fit * (ints(5) + ints(6) + ints(7))
do i_1s = 2, List_env1s_square_size
e_1s = List_env1s_square_expo(i_1s)
c_1s = List_env1s_square_coef(i_1s)
R_1s(1) = List_env1s_square_cent(1,i_1s)
R_1s(2) = List_env1s_square_cent(2,i_1s)
R_1s(3) = List_env1s_square_cent(3,i_1s)
dx = r(1) - R_1s(1)
dy = r(2) - R_1s(2)
dz = r(3) - R_1s(3)
beta = expo_fit + e_1s
tmp_arg = expo_fit * e_1s * (dx*dx + dy*dy + dz*dz) / beta
coef = coef_fit * c_1s * dexp(-tmp_arg)
B_center(1) = (expo_fit * r(1) + e_1s * R_1s(1)) / beta
B_center(2) = (expo_fit * r(2) + e_1s * R_1s(2)) / beta
B_center(3) = (expo_fit * r(3) + e_1s * R_1s(3)) / beta
call overlap_gauss_r12_ao_012(B_center, beta, i, j, ints)
tmp_Du2_0 += coef * ints(1)
tmp_Du2_x += coef * ints(2)
tmp_Du2_y += coef * ints(3)
tmp_Du2_z += coef * ints(4)
tmp_Du2_2 += coef * (ints(5) + ints(6) + ints(7))
enddo ! i_1s
enddo ! i_fit
Ir2_Mu_short_Du2_0(j,i,ipoint) = tmp_Du2_0
Ir2_Mu_short_Du2_x(j,i,ipoint) = tmp_Du2_x
Ir2_Mu_short_Du2_y(j,i,ipoint) = tmp_Du2_y
Ir2_Mu_short_Du2_z(j,i,ipoint) = tmp_Du2_z
Ir2_Mu_short_Du2_2(j,i,ipoint) = tmp_Du2_2
enddo ! j
enddo ! i
enddo ! ipoint
!$OMP END DO
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
Ir2_Mu_short_Du2_0(j,i,ipoint) = Ir2_Mu_short_Du2_0(i,j,ipoint)
Ir2_Mu_short_Du2_x(j,i,ipoint) = Ir2_Mu_short_Du2_x(i,j,ipoint)
Ir2_Mu_short_Du2_y(j,i,ipoint) = Ir2_Mu_short_Du2_y(i,j,ipoint)
Ir2_Mu_short_Du2_z(j,i,ipoint) = Ir2_Mu_short_Du2_z(i,j,ipoint)
Ir2_Mu_short_Du2_2(j,i,ipoint) = Ir2_Mu_short_Du2_2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for Ir2_Mu_short_Du2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---

0
src/ao_many_one_e_ints/list_grid.irp.f → plugins/local/ao_many_one_e_ints/list_grid.irp.f
View File

197
plugins/local/ao_many_one_e_ints/listj1b_sorted.irp.f Normal file
View File

@ -0,0 +1,197 @@
! ---
BEGIN_PROVIDER [integer, List_comb_thr_b2_size, (ao_num, ao_num)]
&BEGIN_PROVIDER [integer, max_List_comb_thr_b2_size]
implicit none
integer :: i_1s, i, j, ipoint
integer :: list(ao_num)
double precision :: coef,beta,center(3),int_env
double precision :: r(3),weight,dist
List_comb_thr_b2_size = 0
print*,'List_env1s_size = ',List_env1s_size
do i = 1, ao_num
do j = i, ao_num
do i_1s = 1, List_env1s_size
coef = List_env1s_coef(i_1s)
if(dabs(coef).lt.thrsh_cycle_tc) cycle
beta = List_env1s_expo(i_1s)
beta = max(beta,1.d-12)
center(1:3) = List_env1s_cent(1:3,i_1s)
int_env = 0.d0
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(ipoint)
dist = ( center(1) - r(1) )*( center(1) - r(1) )
dist += ( center(2) - r(2) )*( center(2) - r(2) )
dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_env += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo
if(dabs(coef)*dabs(int_env).gt.thrsh_cycle_tc)then
List_comb_thr_b2_size(j,i) += 1
endif
enddo
enddo
enddo
do i = 1, ao_num
do j = 1, i-1
List_comb_thr_b2_size(j,i) = List_comb_thr_b2_size(i,j)
enddo
enddo
do i = 1, ao_num
list(i) = maxval(List_comb_thr_b2_size(:,i))
enddo
max_List_comb_thr_b2_size = maxval(list)
print*, ' max_List_comb_thr_b2_size = ',max_List_comb_thr_b2_size
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, List_comb_thr_b2_coef, ( max_List_comb_thr_b2_size,ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_expo, ( max_List_comb_thr_b2_size,ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_cent, (3,max_List_comb_thr_b2_size,ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_env , ( max_List_comb_thr_b2_size,ao_num,ao_num)]
implicit none
integer :: i_1s,i,j,ipoint,icount
double precision :: coef,beta,center(3),int_env
double precision :: r(3),weight,dist
ao_abs_comb_b2_env = 10000000.d0
do i = 1, ao_num
do j = i, ao_num
icount = 0
do i_1s = 1, List_env1s_size
coef = List_env1s_coef (i_1s)
if(dabs(coef).lt.thrsh_cycle_tc)cycle
beta = List_env1s_expo (i_1s)
center(1:3) = List_env1s_cent(1:3,i_1s)
int_env = 0.d0
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(ipoint)
dist = ( center(1) - r(1) )*( center(1) - r(1) )
dist += ( center(2) - r(2) )*( center(2) - r(2) )
dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_env += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo
if(dabs(coef)*dabs(int_env).gt.thrsh_cycle_tc)then
icount += 1
List_comb_thr_b2_coef(icount,j,i) = coef
List_comb_thr_b2_expo(icount,j,i) = beta
List_comb_thr_b2_cent(1:3,icount,j,i) = center(1:3)
ao_abs_comb_b2_env(icount,j,i) = int_env
endif
enddo
enddo
enddo
do i = 1, ao_num
do j = 1, i-1
do icount = 1, List_comb_thr_b2_size(j,i)
List_comb_thr_b2_coef(icount,j,i) = List_comb_thr_b2_coef(icount,i,j)
List_comb_thr_b2_expo(icount,j,i) = List_comb_thr_b2_expo(icount,i,j)
List_comb_thr_b2_cent(1:3,icount,j,i) = List_comb_thr_b2_cent(1:3,icount,i,j)
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [integer, List_comb_thr_b3_size, (ao_num,ao_num)]
&BEGIN_PROVIDER [integer, max_List_comb_thr_b3_size]
implicit none
integer :: i_1s,i,j,ipoint
integer :: list(ao_num)
double precision :: coef,beta,center(3),int_env
double precision :: r(3),weight,dist
List_comb_thr_b3_size = 0
print*,'List_env1s_square_size = ',List_env1s_square_size
do i = 1, ao_num
do j = 1, ao_num
do i_1s = 1, List_env1s_square_size
coef = List_env1s_square_coef (i_1s)
beta = List_env1s_square_expo (i_1s)
center(1:3) = List_env1s_square_cent(1:3,i_1s)
if(dabs(coef).lt.thrsh_cycle_tc)cycle
int_env = 0.d0
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(ipoint)
dist = ( center(1) - r(1) )*( center(1) - r(1) )
dist += ( center(2) - r(2) )*( center(2) - r(2) )
dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_env += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo
if(dabs(coef)*dabs(int_env).gt.thrsh_cycle_tc) then
List_comb_thr_b3_size(j,i) += 1
endif
enddo
enddo
enddo
do i = 1, ao_num
list(i) = maxval(List_comb_thr_b3_size(:,i))
enddo
max_List_comb_thr_b3_size = maxval(list)
print*, ' max_List_comb_thr_b3_size = ',max_List_comb_thr_b3_size
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, List_comb_thr_b3_coef, ( max_List_comb_thr_b3_size,ao_num,ao_num)]
&BEGIN_PROVIDER [double precision, List_comb_thr_b3_expo, ( max_List_comb_thr_b3_size,ao_num,ao_num)]
&BEGIN_PROVIDER [double precision, List_comb_thr_b3_cent, (3, max_List_comb_thr_b3_size,ao_num,ao_num)]
&BEGIN_PROVIDER [double precision, ao_abs_comb_b3_env , ( max_List_comb_thr_b3_size,ao_num,ao_num)]
implicit none
integer :: i_1s,i,j,ipoint,icount
double precision :: coef,beta,center(3),int_env
double precision :: r(3),weight,dist
ao_abs_comb_b3_env = 10000000.d0
do i = 1, ao_num
do j = 1, ao_num
icount = 0
do i_1s = 1, List_env1s_square_size
coef = List_env1s_square_coef (i_1s)
beta = List_env1s_square_expo (i_1s)
beta = max(beta,1.d-12)
center(1:3) = List_env1s_square_cent(1:3,i_1s)
if(dabs(coef).lt.thrsh_cycle_tc)cycle
int_env = 0.d0
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(ipoint)
dist = ( center(1) - r(1) )*( center(1) - r(1) )
dist += ( center(2) - r(2) )*( center(2) - r(2) )
dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_env += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo
if(dabs(coef)*dabs(int_env).gt.thrsh_cycle_tc)then
icount += 1
List_comb_thr_b3_coef(icount,j,i) = coef
List_comb_thr_b3_expo(icount,j,i) = beta
List_comb_thr_b3_cent(1:3,icount,j,i) = center(1:3)
ao_abs_comb_b3_env(icount,j,i) = int_env
endif
enddo
enddo
enddo
END_PROVIDER
! ---

0
src/ao_many_one_e_ints/prim_int_erf_gauss.irp.f → plugins/local/ao_many_one_e_ints/prim_int_erf_gauss.irp.f
View File

2
src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f → plugins/local/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f
View File

@ -200,7 +200,7 @@ subroutine overlap_gauss_r12_v(D_center, LD_D, delta, A_center, B_center, power_
deallocate(A_new, A_center_new, fact_a_new, iorder_a_new, overlap) deallocate(A_new, A_center_new, fact_a_new, iorder_a_new, overlap)
end subroutine overlap_gauss_r12_v end
!--- !---

0
src/ao_many_one_e_ints/stg_gauss_int.irp.f → plugins/local/ao_many_one_e_ints/stg_gauss_int.irp.f
View File

0
src/ao_many_one_e_ints/taylor_exp.irp.f → plugins/local/ao_many_one_e_ints/taylor_exp.irp.f
View File

0
src/ao_many_one_e_ints/xyz_grad_xyz_ao_pol.irp.f → plugins/local/ao_many_one_e_ints/xyz_grad_xyz_ao_pol.irp.f
View File

4
src/ao_tc_eff_map/NEED → plugins/local/ao_tc_eff_map/NEED
View File

@ -1,5 +1,7 @@
ao_two_e_erf_ints ao_two_e_ints
mo_one_e_ints mo_one_e_ints
ao_many_one_e_ints ao_many_one_e_ints
dft_utils_in_r dft_utils_in_r
tc_keywords tc_keywords
hamiltonian
jastrow

0
src/ao_tc_eff_map/README.rst → plugins/local/ao_tc_eff_map/README.rst
View File

11
src/ao_tc_eff_map/compute_ints_eff_pot.irp.f → plugins/local/ao_tc_eff_map/compute_ints_eff_pot.irp.f
View File

@ -23,10 +23,9 @@ subroutine compute_ao_tc_sym_two_e_pot_jl(j, l, n_integrals, buffer_i, buffer_va
logical, external :: ao_two_e_integral_zero logical, external :: ao_two_e_integral_zero
double precision :: ao_tc_sym_two_e_pot, ao_two_e_integral_erf double precision :: ao_tc_sym_two_e_pot, ao_two_e_integral_erf
double precision :: j1b_gauss_2e_j1, j1b_gauss_2e_j2 double precision :: env_gauss_2e_j1, env_gauss_2e_j2
PROVIDE j1b_type
thr = ao_integrals_threshold thr = ao_integrals_threshold
@ -53,14 +52,6 @@ subroutine compute_ao_tc_sym_two_e_pot_jl(j, l, n_integrals, buffer_i, buffer_va
integral_erf = ao_two_e_integral_erf(i, k, j, l) integral_erf = ao_two_e_integral_erf(i, k, j, l)
integral = integral_erf + integral_pot integral = integral_erf + integral_pot
!if( j1b_type .eq. 1 ) then
! !print *, ' j1b type 1 is added'
! integral = integral + j1b_gauss_2e_j1(i, k, j, l)
!elseif( j1b_type .eq. 2 ) then
! !print *, ' j1b type 2 is added'
! integral = integral + j1b_gauss_2e_j2(i, k, j, l)
!endif
if(abs(integral) < thr) then if(abs(integral) < thr) then
cycle cycle
endif endif

0
src/ao_tc_eff_map/integrals_eff_pot_in_map_slave.irp.f → plugins/local/ao_tc_eff_map/integrals_eff_pot_in_map_slave.irp.f
View File

0
src/ao_tc_eff_map/map_integrals_eff_pot.irp.f → plugins/local/ao_tc_eff_map/map_integrals_eff_pot.irp.f
View File

145
src/ao_tc_eff_map/one_e_1bgauss_grad2.irp.f → plugins/local/ao_tc_eff_map/one_e_1bgauss_grad2.irp.f
View File

@ -1,10 +1,10 @@
! --- ! ---
BEGIN_PROVIDER [ double precision, j1b_gauss_hermII, (ao_num,ao_num)] BEGIN_PROVIDER [double precision, env_gauss_hermII, (ao_num,ao_num)]
BEGIN_DOC BEGIN_DOC
! !
! :math:`\langle \chi_A | -0.5 \grad \tau_{1b} \cdot \grad \tau_{1b} | \chi_B \rangle` ! :math:`\langle \chi_A | -0.5 \grad \tau_{env} \cdot \grad \tau_{env} | \chi_B \rangle`
! !
END_DOC END_DOC
@ -22,8 +22,6 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermII, (ao_num,ao_num)]
double precision :: int_gauss_4G double precision :: int_gauss_4G
PROVIDE j1b_type j1b_pen j1b_coeff
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
! -- Dummy call to provide everything ! -- Dummy call to provide everything
dim1 = 100 dim1 = 100
@ -38,10 +36,7 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermII, (ao_num,ao_num)]
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
j1b_gauss_hermII(1:ao_num,1:ao_num) = 0.d0 env_gauss_hermII(1:ao_num,1:ao_num) = 0.d0
if(j1b_type .eq. 1) then
! \tau_1b = \sum_iA -[1 - exp(-alpha_A r_iA^2)]
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
@ -51,113 +46,51 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermII, (ao_num,ao_num)]
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, & !$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, & !$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, & !$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_hermII) !$OMP nucl_num, env_expo, env_gauss_hermII)
!$OMP DO SCHEDULE (dynamic) !$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num do j = 1, ao_num
num_A = ao_nucl(j) num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3) power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3) A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num do i = 1, ao_num
num_B = ao_nucl(i) num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3) power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3) B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j) do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j) alpha = ao_expo_ordered_transp(l,j)
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)
c = 0.d0 c = 0.d0
do k1 = 1, nucl_num do k1 = 1, nucl_num
gama1 = j1b_pen(k1) gama1 = env_expo(k1)
C_center1(1:3) = nucl_coord(k1,1:3) C_center1(1:3) = nucl_coord(k1,1:3)
do k2 = 1, nucl_num do k2 = 1, nucl_num
gama2 = j1b_pen(k2) gama2 = env_expo(k2)
C_center2(1:3) = nucl_coord(k2,1:3) C_center2(1:3) = nucl_coord(k2,1:3)
! < XA | exp[-gama1 r_C1^2 -gama2 r_C2^2] r_C1 \cdot r_C2 | XB > ! < XA | exp[-gama1 r_C1^2 -gama2 r_C2^2] r_C1 \cdot r_C2 | XB >
c1 = int_gauss_4G( A_center, B_center, C_center1, C_center2 & c1 = int_gauss_4G( A_center, B_center, C_center1, C_center2 &
, power_A, power_B, alpha, beta, gama1, gama2 ) , power_A, power_B, alpha, beta, gama1, gama2 )
c = c - 2.d0 * gama1 * gama2 * c1 c = c - 2.d0 * gama1 * gama2 * c1
enddo
enddo enddo
j1b_gauss_hermII(i,j) = j1b_gauss_hermII(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo enddo
env_gauss_hermII(i,j) = env_gauss_hermII(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo enddo
enddo enddo
enddo enddo
enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
elseif(j1b_type .eq. 2) then
! \tau_1b = \sum_iA [c_A exp(-alpha_A r_iA^2)]
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k1, k2, l, m, alpha, beta, gama1, gama2, &
!$OMP A_center, B_center, C_center1, C_center2, &
!$OMP power_A, power_B, num_A, num_B, c1, c, &
!$OMP coef1, coef2) &
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_hermII, &
!$OMP j1b_coeff)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k1 = 1, nucl_num
gama1 = j1b_pen (k1)
coef1 = j1b_coeff(k1)
C_center1(1:3) = nucl_coord(k1,1:3)
do k2 = 1, nucl_num
gama2 = j1b_pen (k2)
coef2 = j1b_coeff(k2)
C_center2(1:3) = nucl_coord(k2,1:3)
! < XA | exp[-gama1 r_C1^2 -gama2 r_C2^2] r_C1 \cdot r_C2 | XB >
c1 = int_gauss_4G( A_center, B_center, C_center1, C_center2 &
, power_A, power_B, alpha, beta, gama1, gama2 )
c = c - 2.d0 * gama1 * gama2 * coef1 * coef2 * c1
enddo
enddo
j1b_gauss_hermII(i,j) = j1b_gauss_hermII(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
END_PROVIDER END_PROVIDER

142
src/ao_tc_eff_map/one_e_1bgauss_lap.irp.f → plugins/local/ao_tc_eff_map/one_e_1bgauss_lap.irp.f
View File

@ -1,10 +1,10 @@
! --- ! ---
BEGIN_PROVIDER [ double precision, j1b_gauss_hermI, (ao_num,ao_num)] BEGIN_PROVIDER [double precision, env_gauss_hermI, (ao_num,ao_num)]
BEGIN_DOC BEGIN_DOC
! !
! :math:`\langle \chi_A | -0.5 \Delta \tau_{1b} | \chi_B \rangle` ! :math:`\langle \chi_A | -0.5 \Delta \tau_{env} | \chi_B \rangle`
! !
END_DOC END_DOC
@ -22,8 +22,6 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermI, (ao_num,ao_num)]
double precision :: int_gauss_r0, int_gauss_r2 double precision :: int_gauss_r0, int_gauss_r2
PROVIDE j1b_type j1b_pen j1b_coeff
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
! -- Dummy call to provide everything ! -- Dummy call to provide everything
dim1 = 100 dim1 = 100
@ -37,10 +35,7 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermI, (ao_num,ao_num)]
, overlap_y, d_a_2, overlap_z, overlap, dim1 ) , overlap_y, d_a_2, overlap_z, overlap, dim1 )
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
j1b_gauss_hermI(1:ao_num,1:ao_num) = 0.d0 env_gauss_hermI(1:ao_num,1:ao_num) = 0.d0
if(j1b_type .eq. 1) then
! \tau_1b = \sum_iA -[1 - exp(-alpha_A r_iA^2)]
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
@ -50,109 +45,50 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermI, (ao_num,ao_num)]
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, & !$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, & !$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, & !$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_hermI) !$OMP nucl_num, env_expo, env_gauss_hermI)
!$OMP DO SCHEDULE (dynamic) !$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num do j = 1, ao_num
num_A = ao_nucl(j) num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3) power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3) A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num do i = 1, ao_num
num_B = ao_nucl(i) num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3) power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3) B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j) do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j) alpha = ao_expo_ordered_transp(l,j)
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)
c = 0.d0 c = 0.d0
do k = 1, nucl_num do k = 1, nucl_num
gama = j1b_pen(k) gama = env_expo(k)
C_center(1:3) = nucl_coord(k,1:3) C_center(1:3) = nucl_coord(k,1:3)
! < XA | exp[-gama r_C^2] | XB > ! < XA | exp[-gama r_C^2] | XB >
c1 = int_gauss_r0( A_center, B_center, C_center & c1 = int_gauss_r0( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama ) , power_A, power_B, alpha, beta, gama )
! < XA | r_A^2 exp[-gama r_C^2] | XB > ! < XA | r_A^2 exp[-gama r_C^2] | XB >
c2 = int_gauss_r2( A_center, B_center, C_center & c2 = int_gauss_r2( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama ) , power_A, power_B, alpha, beta, gama )
c = c + 3.d0 * gama * c1 - 2.d0 * gama * gama * c2 c = c + 3.d0 * gama * c1 - 2.d0 * gama * gama * c2
enddo
j1b_gauss_hermI(i,j) = j1b_gauss_hermI(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo enddo
env_gauss_hermI(i,j) = env_gauss_hermI(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo enddo
enddo enddo
enddo enddo
enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
elseif(j1b_type .eq. 2) then
! \tau_1b = \sum_iA [c_A exp(-alpha_A r_iA^2)]
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k, l, m, alpha, beta, gama, coef, &
!$OMP A_center, B_center, C_center, power_A, power_B, &
!$OMP num_A, num_B, c1, c2, c) &
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_hermI, &
!$OMP j1b_coeff)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k = 1, nucl_num
gama = j1b_pen (k)
coef = j1b_coeff(k)
C_center(1:3) = nucl_coord(k,1:3)
! < XA | exp[-gama r_C^2] | XB >
c1 = int_gauss_r0( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
! < XA | r_A^2 exp[-gama r_C^2] | XB >
c2 = int_gauss_r2( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
c = c + 3.d0 * gama * coef * c1 - 2.d0 * gama * gama * coef * c2
enddo
j1b_gauss_hermI(i,j) = j1b_gauss_hermI(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
END_PROVIDER END_PROVIDER

130
src/ao_tc_eff_map/one_e_1bgauss_nonherm.irp.f → plugins/local/ao_tc_eff_map/one_e_1bgauss_nonherm.irp.f
View File

@ -1,10 +1,11 @@
! --- ! ---
BEGIN_PROVIDER [ double precision, j1b_gauss_nonherm, (ao_num,ao_num)] BEGIN_PROVIDER [double precision, env_gauss_nonherm, (ao_num,ao_num)]
BEGIN_DOC BEGIN_DOC
! !
! j1b_gauss_nonherm(i,j) = \langle \chi_j | - grad \tau_{1b} \cdot grad | \chi_i \rangle ! env_gauss_nonherm(i,j) = \langle \chi_j | - grad \tau_{env} \cdot grad | \chi_i \rangle
! !
END_DOC END_DOC
@ -22,8 +23,6 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_nonherm, (ao_num,ao_num)]
double precision :: int_gauss_deriv double precision :: int_gauss_deriv
PROVIDE j1b_type j1b_pen j1b_coeff
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
! -- Dummy call to provide everything ! -- Dummy call to provide everything
dim1 = 100 dim1 = 100
@ -38,10 +37,8 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_nonherm, (ao_num,ao_num)]
! -------------------------------------------------------------------------------- ! --------------------------------------------------------------------------------
j1b_gauss_nonherm(1:ao_num,1:ao_num) = 0.d0 env_gauss_nonherm(1:ao_num,1:ao_num) = 0.d0
if(j1b_type .eq. 1) then
! \tau_1b = \sum_iA -[1 - exp(-alpha_A r_iA^2)]
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
@ -51,101 +48,46 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_nonherm, (ao_num,ao_num)]
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, & !$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, & !$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, & !$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_nonherm) !$OMP nucl_num, env_expo, env_gauss_nonherm)
!$OMP DO SCHEDULE (dynamic) !$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num do j = 1, ao_num
num_A = ao_nucl(j) num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3) power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3) A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num do i = 1, ao_num
num_B = ao_nucl(i) num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3) power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3) B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j) do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j) alpha = ao_expo_ordered_transp(l,j)
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)
c = 0.d0 c = 0.d0
do k = 1, nucl_num do k = 1, nucl_num
gama = j1b_pen(k) gama = env_expo(k)
C_center(1:3) = nucl_coord(k,1:3) C_center(1:3) = nucl_coord(k,1:3)
! \langle \chi_A | exp[-gama r_C^2] r_C \cdot grad | \chi_B \rangle ! \langle \chi_A | exp[-gama r_C^2] r_C \cdot grad | \chi_B \rangle
c1 = int_gauss_deriv( A_center, B_center, C_center & c1 = int_gauss_deriv( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama ) , power_A, power_B, alpha, beta, gama )
c = c + 2.d0 * gama * c1 c = c + 2.d0 * gama * c1
enddo
j1b_gauss_nonherm(i,j) = j1b_gauss_nonherm(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo enddo
env_gauss_nonherm(i,j) = env_gauss_nonherm(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo enddo
enddo enddo
enddo enddo
enddo
!$OMP END DO !$OMP END DO
!$OMP END PARALLEL !$OMP END PARALLEL
elseif(j1b_type .eq. 2) then
! \tau_1b = \sum_iA [c_A exp(-alpha_A r_iA^2)]
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k, l, m, alpha, beta, gama, coef, &
!$OMP A_center, B_center, C_center, power_A, power_B, &
!$OMP num_A, num_B, c1, c) &
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_nonherm, &
!$OMP j1b_coeff)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k = 1, nucl_num
gama = j1b_pen (k)
coef = j1b_coeff(k)
C_center(1:3) = nucl_coord(k,1:3)
! \langle \chi_A | exp[-gama r_C^2] r_C \cdot grad | \chi_B \rangle
c1 = int_gauss_deriv( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
c = c + 2.d0 * gama * coef * c1
enddo
j1b_gauss_nonherm(i,j) = j1b_gauss_nonherm(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
END_PROVIDER END_PROVIDER

3
src/ao_tc_eff_map/providers_ao_eff_pot.irp.f → plugins/local/ao_tc_eff_map/providers_ao_eff_pot.irp.f
View File

@ -22,9 +22,6 @@ BEGIN_PROVIDER [ logical, ao_tc_sym_two_e_pot_in_map ]
integer :: kk, m, j1, i1, lmax integer :: kk, m, j1, i1, lmax
character*(64) :: fmt character*(64) :: fmt
!double precision :: j1b_gauss_coul_debug
!integral = j1b_gauss_coul_debug(1,1,1,1)
integral = ao_tc_sym_two_e_pot(1,1,1,1) integral = ao_tc_sym_two_e_pot(1,1,1,1)
double precision :: map_mb double precision :: map_mb

27
src/ao_tc_eff_map/two_e_1bgauss_j1.irp.f → plugins/local/ao_tc_eff_map/two_e_1bgauss_j1.irp.f
View File

@ -1,6 +1,6 @@
! --- ! ---
double precision function j1b_gauss_2e_j1(i, j, k, l) double precision function env_gauss_2e_j1(i, j, k, l)
BEGIN_DOC BEGIN_DOC
! !
@ -36,10 +36,10 @@ double precision function j1b_gauss_2e_j1(i, j, k, l)
double precision :: I_center(3), J_center(3), K_center(3), L_center(3) double precision :: I_center(3), J_center(3), K_center(3), L_center(3)
double precision :: ff, gg, cx, cy, cz double precision :: ff, gg, cx, cy, cz
double precision :: j1b_gauss_2e_j1_schwartz double precision :: env_gauss_2e_j1_schwartz
if( ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then if( ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
j1b_gauss_2e_j1 = j1b_gauss_2e_j1_schwartz(i, j, k, l) env_gauss_2e_j1 = env_gauss_2e_j1_schwartz(i, j, k, l)
return return
endif endif
@ -59,7 +59,7 @@ double precision function j1b_gauss_2e_j1(i, j, k, l)
L_center(p) = nucl_coord(num_l,p) L_center(p) = nucl_coord(num_l,p)
enddo enddo
j1b_gauss_2e_j1 = 0.d0 env_gauss_2e_j1 = 0.d0
do p = 1, ao_prim_num(i) do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p, i) coef1 = ao_coef_normalized_ordered_transp(p, i)
@ -89,18 +89,18 @@ double precision function j1b_gauss_2e_j1(i, j, k, l)
, P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p & , P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p &
, Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q ) , Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q )
j1b_gauss_2e_j1 = j1b_gauss_2e_j1 + coef4 * ( cx + cy + cz ) env_gauss_2e_j1 = env_gauss_2e_j1 + coef4 * ( cx + cy + cz )
enddo ! s enddo ! s
enddo ! r enddo ! r
enddo ! q enddo ! q
enddo ! p enddo ! p
return return
end function j1b_gauss_2e_j1 end
! --- ! ---
double precision function j1b_gauss_2e_j1_schwartz(i, j, k, l) double precision function env_gauss_2e_j1_schwartz(i, j, k, l)
BEGIN_DOC BEGIN_DOC
! !
@ -137,8 +137,6 @@ double precision function j1b_gauss_2e_j1_schwartz(i, j, k, l)
double precision :: schwartz_ij, thr double precision :: schwartz_ij, thr
double precision, allocatable :: schwartz_kl(:,:) double precision, allocatable :: schwartz_kl(:,:)
PROVIDE j1b_pen
dim1 = n_pt_max_integrals dim1 = n_pt_max_integrals
thr = ao_integrals_threshold * ao_integrals_threshold thr = ao_integrals_threshold * ao_integrals_threshold
@ -186,8 +184,7 @@ double precision function j1b_gauss_2e_j1_schwartz(i, j, k, l)
schwartz_kl(0,0) = max( schwartz_kl(0,r) , schwartz_kl(0,0) ) schwartz_kl(0,0) = max( schwartz_kl(0,r) , schwartz_kl(0,0) )
enddo enddo
env_gauss_2e_j1_schwartz = 0.d0
j1b_gauss_2e_j1_schwartz = 0.d0
do p = 1, ao_prim_num(i) do p = 1, ao_prim_num(i)
expo1 = ao_expo_ordered_transp(p, i) expo1 = ao_expo_ordered_transp(p, i)
@ -226,7 +223,7 @@ double precision function j1b_gauss_2e_j1_schwartz(i, j, k, l)
, P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p & , P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p &
, Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q ) , Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q )
j1b_gauss_2e_j1_schwartz = j1b_gauss_2e_j1_schwartz + coef4 * ( cx + cy + cz ) env_gauss_2e_j1_schwartz = env_gauss_2e_j1_schwartz + coef4 * ( cx + cy + cz )
enddo ! s enddo ! s
enddo ! r enddo ! r
enddo ! q enddo ! q
@ -235,7 +232,7 @@ double precision function j1b_gauss_2e_j1_schwartz(i, j, k, l)
deallocate( schwartz_kl ) deallocate( schwartz_kl )
return return
end function j1b_gauss_2e_j1_schwartz end
! --- ! ---
@ -263,14 +260,12 @@ subroutine get_cxcycz_j1( dim1, cx, cy, cz &
double precision :: general_primitive_integral_erf_shifted double precision :: general_primitive_integral_erf_shifted
double precision :: general_primitive_integral_coul_shifted double precision :: general_primitive_integral_coul_shifted
PROVIDE j1b_pen
cx = 0.d0 cx = 0.d0
cy = 0.d0 cy = 0.d0
cz = 0.d0 cz = 0.d0
do ii = 1, nucl_num do ii = 1, nucl_num
expoii = j1b_pen(ii) expoii = env_expo(ii)
Centerii(1:3) = nucl_coord(ii, 1:3) Centerii(1:3) = nucl_coord(ii, 1:3)
call gaussian_product(pp1, P1_center, expoii, Centerii, factii, pp2, P2_center) call gaussian_product(pp1, P1_center, expoii, Centerii, factii, pp2, P2_center)

26
src/ao_tc_eff_map/two_e_1bgauss_j2.irp.f → plugins/local/ao_tc_eff_map/two_e_1bgauss_j2.irp.f
View File

@ -1,6 +1,6 @@
! --- ! ---
double precision function j1b_gauss_2e_j2(i, j, k, l) double precision function env_gauss_2e_j2(i, j, k, l)
BEGIN_DOC BEGIN_DOC
! !
@ -36,12 +36,12 @@ double precision function j1b_gauss_2e_j2(i, j, k, l)
double precision :: I_center(3), J_center(3), K_center(3), L_center(3) double precision :: I_center(3), J_center(3), K_center(3), L_center(3)
double precision :: ff, gg, cx, cy, cz double precision :: ff, gg, cx, cy, cz
double precision :: j1b_gauss_2e_j2_schwartz double precision :: env_gauss_2e_j2_schwartz
dim1 = n_pt_max_integrals dim1 = n_pt_max_integrals
if( ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then if( ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
j1b_gauss_2e_j2 = j1b_gauss_2e_j2_schwartz(i, j, k, l) env_gauss_2e_j2 = env_gauss_2e_j2_schwartz(i, j, k, l)
return return
endif endif
@ -61,7 +61,7 @@ double precision function j1b_gauss_2e_j2(i, j, k, l)
L_center(p) = nucl_coord(num_l,p) L_center(p) = nucl_coord(num_l,p)
enddo enddo
j1b_gauss_2e_j2 = 0.d0 env_gauss_2e_j2 = 0.d0
do p = 1, ao_prim_num(i) do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p, i) coef1 = ao_coef_normalized_ordered_transp(p, i)
@ -91,18 +91,18 @@ double precision function j1b_gauss_2e_j2(i, j, k, l)
, P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p & , P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p &
, Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q ) , Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q )
j1b_gauss_2e_j2 = j1b_gauss_2e_j2 + coef4 * ( cx + cy + cz ) env_gauss_2e_j2 = env_gauss_2e_j2 + coef4 * ( cx + cy + cz )
enddo ! s enddo ! s
enddo ! r enddo ! r
enddo ! q enddo ! q
enddo ! p enddo ! p
return return
end function j1b_gauss_2e_j2 end
! --- ! ---
double precision function j1b_gauss_2e_j2_schwartz(i, j, k, l) double precision function env_gauss_2e_j2_schwartz(i, j, k, l)
BEGIN_DOC BEGIN_DOC
! !
@ -187,7 +187,7 @@ double precision function j1b_gauss_2e_j2_schwartz(i, j, k, l)
enddo enddo
j1b_gauss_2e_j2_schwartz = 0.d0 env_gauss_2e_j2_schwartz = 0.d0
do p = 1, ao_prim_num(i) do p = 1, ao_prim_num(i)
expo1 = ao_expo_ordered_transp(p, i) expo1 = ao_expo_ordered_transp(p, i)
@ -226,7 +226,7 @@ double precision function j1b_gauss_2e_j2_schwartz(i, j, k, l)
, P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p & , P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p &
, Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q ) , Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q )
j1b_gauss_2e_j2_schwartz = j1b_gauss_2e_j2_schwartz + coef4 * ( cx + cy + cz ) env_gauss_2e_j2_schwartz = env_gauss_2e_j2_schwartz + coef4 * ( cx + cy + cz )
enddo ! s enddo ! s
enddo ! r enddo ! r
enddo ! q enddo ! q
@ -235,7 +235,7 @@ double precision function j1b_gauss_2e_j2_schwartz(i, j, k, l)
deallocate( schwartz_kl ) deallocate( schwartz_kl )
return return
end function j1b_gauss_2e_j2_schwartz end
! --- ! ---
@ -263,15 +263,13 @@ subroutine get_cxcycz_j2( dim1, cx, cy, cz &
double precision :: general_primitive_integral_erf_shifted double precision :: general_primitive_integral_erf_shifted
double precision :: general_primitive_integral_coul_shifted double precision :: general_primitive_integral_coul_shifted
PROVIDE j1b_pen j1b_coeff
cx = 0.d0 cx = 0.d0
cy = 0.d0 cy = 0.d0
cz = 0.d0 cz = 0.d0
do ii = 1, nucl_num do ii = 1, nucl_num
expoii = j1b_pen (ii) expoii = env_expo(ii)
coefii = j1b_coeff(ii) coefii = env_coef(ii)
Centerii(1:3) = nucl_coord(ii, 1:3) Centerii(1:3) = nucl_coord(ii, 1:3)
call gaussian_product(pp1, P1_center, expoii, Centerii, factii, pp2, P2_center) call gaussian_product(pp1, P1_center, expoii, Centerii, factii, pp2, P2_center)

0
src/ao_tc_eff_map/two_e_ints_gauss.irp.f → plugins/local/ao_tc_eff_map/two_e_ints_gauss.irp.f
View File

49
src/ao_tc_eff_map/useful_sub.irp.f → plugins/local/ao_tc_eff_map/useful_sub.irp.f
View File

@ -174,7 +174,7 @@ double precision function general_primitive_integral_coul_shifted( dim
general_primitive_integral_coul_shifted = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / dsqrt(p_plus_q) general_primitive_integral_coul_shifted = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / dsqrt(p_plus_q)
return return
end function general_primitive_integral_coul_shifted end
!______________________________________________________________________________________________________________________ !______________________________________________________________________________________________________________________
!______________________________________________________________________________________________________________________ !______________________________________________________________________________________________________________________
@ -354,7 +354,7 @@ double precision function general_primitive_integral_erf_shifted( dim
general_primitive_integral_erf_shifted = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / dsqrt(p_plus_q) general_primitive_integral_erf_shifted = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / dsqrt(p_plus_q)
return return
end function general_primitive_integral_erf_shifted end
!______________________________________________________________________________________________________________________ !______________________________________________________________________________________________________________________
!______________________________________________________________________________________________________________________ !______________________________________________________________________________________________________________________
@ -362,3 +362,48 @@ end function general_primitive_integral_erf_shifted
! ---
subroutine inv_r_times_poly(r, dist_r, dist_vec, poly)
BEGIN_DOC
!
! returns
!
! poly(1) = x / sqrt(x^2+y^2+z^2), poly(2) = y / sqrt(x^2+y^2+z^2), poly(3) = z / sqrt(x^2+y^2+z^2)
!
! with the arguments
!
! r(1) = x, r(2) = y, r(3) = z, dist_r = sqrt(x^2+y^2+z^2)
!
! dist_vec(1) = sqrt(y^2+z^2), dist_vec(2) = sqrt(x^2+z^2), dist_vec(3) = sqrt(x^2+y^2)
!
END_DOC
implicit none
double precision, intent(in) :: r(3), dist_r, dist_vec(3)
double precision, intent(out) :: poly(3)
integer :: i
double precision :: inv_dist
if (dist_r .gt. 1.d-8)then
inv_dist = 1.d0/dist_r
do i = 1, 3
poly(i) = r(i) * inv_dist
enddo
else
do i = 1, 3
if(dabs(r(i)).lt.dist_vec(i)) then
inv_dist = 1.d0/dist_r
poly(i) = r(i) * inv_dist
else
poly(i) = 1.d0
endif
enddo
endif
end
! ---

4
src/basis_correction/51.basis_c.bats → plugins/local/basis_correction/51.basis_c.bats
View File

@ -10,8 +10,8 @@ function run() {
qp set perturbation do_pt2 False qp set perturbation do_pt2 False
qp set determinants n_det_max 8000 qp set determinants n_det_max 8000
qp set determinants n_states 1 qp set determinants n_states 1
qp set davidson threshold_davidson 1.e-10 qp set davidson_keywords threshold_davidson 1.e-10
qp set davidson n_states_diag 8 qp set davidson_keywords n_states_diag 8
qp run fci qp run fci
energy1="$(ezfio get fci energy | tr '[]' ' ' | cut -d ',' -f 1)" energy1="$(ezfio get fci energy | tr '[]' ' ' | cut -d ',' -f 1)"
eq $energy1 $1 $thresh eq $energy1 $1 $thresh

0
src/basis_correction/NEED → plugins/local/basis_correction/NEED
View File

2
src/basis_correction/README.rst → plugins/local/basis_correction/README.rst
View File

@ -12,7 +12,7 @@ This basis set correction relies mainy on :
When HF is a qualitative representation of the electron pairs (i.e. weakly correlated systems), such an approach for \mu(r) is OK. When HF is a qualitative representation of the electron pairs (i.e. weakly correlated systems), such an approach for \mu(r) is OK.
See for instance JPCL, 10, 2931-2937 (2019) for typical flavours of the results. See for instance JPCL, 10, 2931-2937 (2019) for typical flavours of the results.
Thanks to the trivial nature of such a two-body rdm, the equation (22) of J. Chem. Phys. 149, 194301 (2018) can be rewritten in a very efficient way, and therefore the limiting factor of such an approach is the AO->MO four-index transformation of the two-electron integrals. Thanks to the trivial nature of such a two-body rdm, the equation (22) of J. Chem. Phys. 149, 194301 (2018) can be rewritten in a very efficient way, and therefore the limiting factor of such an approach is the AO->MO four-index transformation of the two-electron integrals.
b) "mu_of_r_potential = cas_ful" uses the two-body rdm of CAS-like wave function (i.e. linear combination of Slater determinants developped in an active space with the MOs stored in the EZFIO folder). b) "mu_of_r_potential = cas_full" uses the two-body rdm of CAS-like wave function (i.e. linear combination of Slater determinants developped in an active space with the MOs stored in the EZFIO folder).
If the CAS is properly chosen (i.e. the CAS-like wave function qualitatively represents the wave function of the systems), then such an approach is OK for \mu(r) even in the case of strong correlation. If the CAS is properly chosen (i.e. the CAS-like wave function qualitatively represents the wave function of the systems), then such an approach is OK for \mu(r) even in the case of strong correlation.
+) The use of DFT correlation functionals with multi-determinant reference (Ecmd). These functionals are originally defined in the RS-DFT framework (see for instance Theor. Chem. Acc.114, 305(2005)) and design to capture short-range correlation effects. A important quantity arising in the Ecmd is the exact on-top pair density of the system, and the main differences of approximated Ecmd relies on different approximations for the exact on-top pair density. +) The use of DFT correlation functionals with multi-determinant reference (Ecmd). These functionals are originally defined in the RS-DFT framework (see for instance Theor. Chem. Acc.114, 305(2005)) and design to capture short-range correlation effects. A important quantity arising in the Ecmd is the exact on-top pair density of the system, and the main differences of approximated Ecmd relies on different approximations for the exact on-top pair density.

0
src/basis_correction/TODO → plugins/local/basis_correction/TODO
View File

0
src/basis_correction/basis_correction.irp.f → plugins/local/basis_correction/basis_correction.irp.f
View File

0
src/basis_correction/eff_xi_based_func.irp.f → plugins/local/basis_correction/eff_xi_based_func.irp.f
View File

6
src/basis_correction/pbe_on_top.irp.f → plugins/local/basis_correction/pbe_on_top.irp.f
View File

@ -39,7 +39,7 @@
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate) grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate) grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
if(mu_of_r_potential == "cas_ful")then if(mu_of_r_potential == "cas_full")then
! You take the on-top of the CAS wave function which is computed with mu(r) ! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate) on_top = on_top_cas_mu_r(ipoint,istate)
else else
@ -101,7 +101,7 @@
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate) grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate) grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
if(mu_of_r_potential == "cas_ful")then if(mu_of_r_potential == "cas_full")then
! You take the on-top of the CAS wave function which is computed with mu(r) ! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate) on_top = on_top_cas_mu_r(ipoint,istate)
else else
@ -163,7 +163,7 @@
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate) grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate) grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
if(mu_of_r_potential == "cas_ful")then if(mu_of_r_potential == "cas_full")then
! You take the on-top of the CAS wave function which is computed with mu(r) ! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate) on_top = on_top_cas_mu_r(ipoint,istate)
else else

14
src/basis_correction/print_routine.irp.f → plugins/local/basis_correction/print_routine.irp.f
View File

@ -4,8 +4,8 @@ subroutine print_basis_correction
provide mu_average_prov provide mu_average_prov
if(mu_of_r_potential.EQ."hf")then if(mu_of_r_potential.EQ."hf")then
provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r
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_full".or.mu_of_r_potential.EQ."cas_truncated")then
provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r
provide ecmd_pbe_on_top_mu_of_r ecmd_pbe_on_top_su_mu_of_r provide ecmd_pbe_on_top_mu_of_r ecmd_pbe_on_top_su_mu_of_r
endif endif
@ -25,7 +25,7 @@ subroutine print_basis_correction
if(mu_of_r_potential.EQ."hf")then if(mu_of_r_potential.EQ."hf")then
print*, '' print*, ''
print*,'Using a HF-like two-body density to define mu(r)' print*,'Using a HF-like two-body density to define mu(r)'
print*,'This assumes that HF is a qualitative representation of the wave function ' print*,'This assumes that HF is a qualitative representation of the wave function '
print*,'********************************************' print*,'********************************************'
print*,'Functionals more suited for weak correlation' print*,'Functionals more suited for weak correlation'
print*,'********************************************' print*,'********************************************'
@ -38,10 +38,10 @@ 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".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then else if(mu_of_r_potential.EQ."cas_full".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 '
print*,'********************************************' print*,'********************************************'
print*,'Functionals more suited for weak correlation' print*,'Functionals more suited for weak correlation'
print*,'********************************************' print*,'********************************************'
@ -56,14 +56,14 @@ subroutine print_basis_correction
print*,'' print*,''
print*,'********************************************' print*,'********************************************'
print*,'********************************************' print*,'********************************************'
print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) ' print*,'+) PBE-on-top Ecmd functional : JCP, 152, 174104 (2020) '
print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization' print*,'PBE at mu=0, extrapolated ontop pair density at large mu, usual spin-polarization'
do istate = 1, N_states do istate = 1, N_states
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate) write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_mu_of_r(istate)
enddo enddo
print*,'' print*,''
print*,'********************************************' print*,'********************************************'
print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)' print*,'+) PBE-on-top no spin polarization Ecmd functional : JCP, 152, 174104 (2020)'
print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION' print*,'PBE at mu=0, extrapolated ontop pair density at large mu, and ZERO SPIN POLARIZATION'
do istate = 1, N_states 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) write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD SU-PBE-OT , state ',istate,' = ',ecmd_pbe_on_top_su_mu_of_r(istate)

0
src/basis_correction/print_su_pbe_ot.irp.f → plugins/local/basis_correction/print_su_pbe_ot.irp.f
View File

0
src/basis_correction/weak_corr_func.irp.f → plugins/local/basis_correction/weak_corr_func.irp.f
View File

0
src/bi_ort_ints/NEED → plugins/local/bi_ort_ints/NEED
View File

0
src/bi_ort_ints/README.rst → plugins/local/bi_ort_ints/README.rst
View File

0
src/bi_ort_ints/bi_ort_ints.irp.f → plugins/local/bi_ort_ints/bi_ort_ints.irp.f
View File

37
src/bi_ort_ints/biorthog_mo_for_h.irp.f → plugins/local/bi_ort_ints/biorthog_mo_for_h.irp.f
View File

@ -1,4 +1,39 @@
! ---
BEGIN_PROVIDER [double precision, ao_two_e_coul, (ao_num, ao_num, ao_num, ao_num) ]
BEGIN_DOC
!
! ao_two_e_coul(k,i,l,j) = ( k i | 1/r12 | l j ) = < l k | 1/r12 | j i >
!
END_DOC
integer :: i, j, k, l
double precision, external :: get_ao_two_e_integral
PROVIDE ao_integrals_map
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(ao_num, ao_two_e_coul, ao_integrals_map) &
!$OMP PRIVATE(i, j, k, l)
!$OMP DO
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
! < 1:k, 2:l | 1:i, 2:j >
ao_two_e_coul(k,i,l,j) = get_ao_two_e_integral(i, j, k, l, ao_integrals_map)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
! --- ! ---
double precision function bi_ortho_mo_coul_ints(l, k, j, i) double precision function bi_ortho_mo_coul_ints(l, k, j, i)
@ -25,7 +60,7 @@ double precision function bi_ortho_mo_coul_ints(l, k, j, i)
enddo enddo
enddo enddo
end function bi_ortho_mo_coul_ints end
! --- ! ---

0
src/bi_ort_ints/no_dressing.irp.f → plugins/local/bi_ort_ints/no_dressing.irp.f
View File

0
src/bi_ort_ints/no_dressing_energy.irp.f → plugins/local/bi_ort_ints/no_dressing_energy.irp.f
View File

0
src/bi_ort_ints/no_dressing_naive.irp.f → plugins/local/bi_ort_ints/no_dressing_naive.irp.f
View File

17
src/bi_ort_ints/one_e_bi_ort.irp.f → plugins/local/bi_ort_ints/one_e_bi_ort.irp.f
View File

@ -8,23 +8,6 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)]
ao_one_e_integrals_tc_tot = ao_one_e_integrals ao_one_e_integrals_tc_tot = ao_one_e_integrals
!provide j1b_type
!if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then
!
! print *, ' do things properly !'
! stop
! !do i = 1, ao_num
! ! do j = 1, ao_num
! ! ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) &
! ! + j1b_gauss_hermII (j,i) &
! ! + j1b_gauss_nonherm(j,i) )
! ! enddo
! !enddo
!endif
END_PROVIDER END_PROVIDER
! --- ! ---

0
src/bi_ort_ints/semi_num_ints_mo.irp.f → plugins/local/bi_ort_ints/semi_num_ints_mo.irp.f
View File

0
src/bi_ort_ints/three_body_ijm.irp.f → plugins/local/bi_ort_ints/three_body_ijm.irp.f
View File

0
src/bi_ort_ints/three_body_ijmk.irp.f → plugins/local/bi_ort_ints/three_body_ijmk.irp.f
View File

0
src/bi_ort_ints/three_body_ijmk_n4.irp.f → plugins/local/bi_ort_ints/three_body_ijmk_n4.irp.f
View File

0
src/bi_ort_ints/three_body_ijmk_old.irp.f → plugins/local/bi_ort_ints/three_body_ijmk_old.irp.f
View File

0
src/bi_ort_ints/three_body_ijmkl.irp.f → plugins/local/bi_ort_ints/three_body_ijmkl.irp.f
View File

0
src/bi_ort_ints/three_body_ijmkl_old.irp.f → plugins/local/bi_ort_ints/three_body_ijmkl_old.irp.f
View File

0
src/bi_ort_ints/three_body_ints_bi_ort.irp.f → plugins/local/bi_ort_ints/three_body_ints_bi_ort.irp.f
View File

90
src/bi_ort_ints/total_twoe_pot.irp.f → plugins/local/bi_ort_ints/total_twoe_pot.irp.f
View File

@ -1,91 +1,4 @@
! ---
BEGIN_PROVIDER [double precision, ao_two_e_vartc_tot, (ao_num, ao_num, ao_num, ao_num) ]
integer :: i, j, k, l
provide j1b_type
provide mo_r_coef mo_l_coef
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
ao_two_e_vartc_tot(k,i,l,j) = ao_vartc_int_chemist(k,i,l,j)
enddo
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_num) ]
BEGIN_DOC
!
! ao_two_e_tc_tot(k,i,l,j) = (ki|V^TC(r_12)|lj) = <lk| V^TC(r_12) |ji> where V^TC(r_12) is the total TC operator
!
! including both hermitian and non hermitian parts. THIS IS IN CHEMIST NOTATION.
!
! WARNING :: non hermitian ! acts on "the right functions" (i,j)
!
END_DOC
integer :: i, j, k, l
double precision :: integral_sym, integral_nsym
double precision, external :: get_ao_tc_sym_two_e_pot
provide j1b_type
if(j1b_type .eq. 0) then
PROVIDE ao_tc_sym_two_e_pot_in_map
!!! TODO :: OPENMP
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
integral_sym = get_ao_tc_sym_two_e_pot(i, j, k, l, ao_tc_sym_two_e_pot_map)
! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis
integral_nsym = ao_non_hermit_term_chemist(k,i,l,j)
!print *, ' sym integ = ', integral_sym
!print *, ' non-sym integ = ', integral_nsym
ao_two_e_tc_tot(k,i,l,j) = integral_sym + integral_nsym
!write(111,*) ao_two_e_tc_tot(k,i,l,j)
enddo
enddo
enddo
enddo
else
PROVIDE ao_tc_int_chemist
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
ao_two_e_tc_tot(k,i,l,j) = ao_tc_int_chemist(k,i,l,j)
!write(222,*) ao_two_e_tc_tot(k,i,l,j)
enddo
enddo
enddo
enddo
FREE ao_tc_int_chemist
endif
END_PROVIDER
! --- ! ---
double precision function bi_ortho_mo_ints(l, k, j, i) double precision function bi_ortho_mo_ints(l, k, j, i)
@ -118,8 +31,6 @@ end function bi_ortho_mo_ints
! --- ! ---
! TODO :: transform into DEGEMM
BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num, mo_num, mo_num)] BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC BEGIN_DOC
@ -267,7 +178,6 @@ END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj, (mo_num,mo_num)] BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj, (mo_num,mo_num)]
&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_exchange, (mo_num,mo_num)] &BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_exchange, (mo_num,mo_num)]
&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_anti, (mo_num,mo_num)] &BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_anti, (mo_num,mo_num)]

0
src/bi_ortho_mos/EZFIO.cfg → plugins/local/bi_ortho_mos/EZFIO.cfg
View File

0
src/bi_ortho_mos/NEED → plugins/local/bi_ortho_mos/NEED
View File

0
src/bi_ortho_mos/bi_density.irp.f → plugins/local/bi_ortho_mos/bi_density.irp.f
View File

0
src/bi_ortho_mos/bi_ort_mos_in_r.irp.f → plugins/local/bi_ortho_mos/bi_ort_mos_in_r.irp.f
View File

0
src/bi_ortho_mos/grad_bi_ort_mos_in_r.irp.f → plugins/local/bi_ortho_mos/grad_bi_ort_mos_in_r.irp.f
View File

0
src/bi_ortho_mos/mos_rl.irp.f → plugins/local/bi_ortho_mos/mos_rl.irp.f
View File

0
src/bi_ortho_mos/overlap.irp.f → plugins/local/bi_ortho_mos/overlap.irp.f
View File

0
src/casscf_tc_bi/NEED → plugins/local/casscf_tc_bi/NEED
View File

0
src/casscf_tc_bi/det_manip.irp.f → plugins/local/casscf_tc_bi/det_manip.irp.f
View File

0
src/casscf_tc_bi/grad_dm.irp.f → plugins/local/casscf_tc_bi/grad_dm.irp.f
View File

0
src/casscf_tc_bi/grad_old.irp.f → plugins/local/casscf_tc_bi/grad_old.irp.f
View File

0
src/casscf_tc_bi/gradient.irp.f → plugins/local/casscf_tc_bi/gradient.irp.f
View File

0
src/casscf_tc_bi/test_tc_casscf.irp.f → plugins/local/casscf_tc_bi/test_tc_casscf.irp.f
View File

0
src/cipsi_tc_bi_ortho/EZFIO.cfg → plugins/local/cipsi_tc_bi_ortho/EZFIO.cfg
View File

0
src/cipsi_tc_bi_ortho/NEED → plugins/local/cipsi_tc_bi_ortho/NEED
View File

0
src/cipsi_tc_bi_ortho/cipsi.irp.f → plugins/local/cipsi_tc_bi_ortho/cipsi.irp.f
View File

0
src/cipsi_tc_bi_ortho/energy.irp.f → plugins/local/cipsi_tc_bi_ortho/energy.irp.f
View File

0
src/cipsi_tc_bi_ortho/environment.irp.f → plugins/local/cipsi_tc_bi_ortho/environment.irp.f
View File

0
src/cipsi_tc_bi_ortho/fock_diag.irp.f → plugins/local/cipsi_tc_bi_ortho/fock_diag.irp.f
View File

0
src/cipsi_tc_bi_ortho/get_d.irp.f → plugins/local/cipsi_tc_bi_ortho/get_d.irp.f
View File

0
src/cipsi_tc_bi_ortho/get_d0_good.irp.f → plugins/local/cipsi_tc_bi_ortho/get_d0_good.irp.f
View File

0
src/cipsi_tc_bi_ortho/get_d1_good.irp.f → plugins/local/cipsi_tc_bi_ortho/get_d1_good.irp.f
View File

Some files were not shown because too many files have changed in this diff Show More