Merge pull request #26 from QuantumPackage/dev-stable

Dev stable
This commit is contained in:
AbdAmmar 2023-10-18 23:53:14 +02:00 committed by GitHub
commit 6fba691382
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
93 changed files with 2963 additions and 727 deletions

View File

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

62
config/flang_avx.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 : flang -ffree-line-length-none -I . -mavx -g -fPIC
LAPACK_LIB : -llapack -lblas
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DSET_NESTED
# 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 -mavx
# 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

62
config/gfortran10.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-10 -g -ffree-line-length-none -I . -fPIC
LAPACK_LIB : -lblas -llapack
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --assert -DSET_NESTED
# Global options
################
#
# 1 : Activate
# 0 : Deactivate
#
[OPTION]
MODE : DEBUG ; [ 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
# Profiling flags
#################
#
[PROFILE]
FC : -p -g
FCFLAGS : -Ofast
# Debugging flags
#################
#
# -fcheck=all : Checks uninitialized variables, array subscripts, etc...
# -g : Extra debugging information
#
[DEBUG]
FCFLAGS : -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant -Wuninitialized -fbacktrace -ffpe-trap=zero,overflow,underflow -finit-real=nan
# OpenMP flags
#################
#
[OPENMP]
FC : -fopenmp
IRPF90_FLAGS : --openmp

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 : mpif90 -ffree-line-length-none -I . -g -fPIC -std=legacy
LAPACK_LIB : -lblas -llapack
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 -DMPI -DSET_NESTED
# 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 -msse4.2
# Profiling flags
#################
#
[PROFILE]
FC : -p -g
FCFLAGS : -Ofast -msse4.2
# 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

View File

@ -0,0 +1,63 @@
# Common flags
##############
#
# -mkl=[parallel|sequential] : Use the MKL library
# --ninja : Allow the utilisation of ninja. It is mandatory !
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : ifort -fpic
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --define=WITHOUT_TRAILZ --define=WITHOUT_SHIFTRL -DSET_NESTED
# 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
####################
#
# -xHost : Compile a binary optimized for the current architecture
# -O2 : O3 not better than O2.
# -ip : Inter-procedural optimizations
# -ftz : Flushes denormal results to zero
#
[OPT]
FC : -traceback
FCFLAGS : -xAVX -O2 -ip -ftz -g
# Profiling flags
#################
#
[PROFILE]
FC : -p -g
FCFLAGS : -xSSE4.2 -O2 -ip -ftz
# Debugging flags
#################
#
# -traceback : Activate backtrace on runtime
# -fpe0 : All floating point exaceptions
# -C : Checks uninitialized variables, array subscripts, etc...
# -g : Extra debugging information
# -xSSE2 : Valgrind needs a very simple x86 executable
#
[DEBUG]
FC : -g -traceback
FCFLAGS : -xSSE2 -C -fpe0 -implicitnone
# OpenMP flags
#################
#
[OPENMP]
FC : -qopenmp
IRPF90_FLAGS : --openmp

View File

@ -7,7 +7,7 @@
#
[COMMON]
FC : ifort -fpic
LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --assert -DINTEL -DSET_NESTED

View File

@ -0,0 +1,63 @@
# Common flags
##############
#
# -mkl=[parallel|sequential] : Use the MKL library
# --ninja : Allow the utilisation of ninja. It is mandatory !
# --align=32 : Align all provided arrays on a 32-byte boundary
#
[COMMON]
FC : ifort -fpic
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --define=WITHOUT_TRAILZ --define=WITHOUT_SHIFTRL
# 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
####################
#
# -xHost : Compile a binary optimized for the current architecture
# -O2 : O3 not better than O2.
# -ip : Inter-procedural optimizations
# -ftz : Flushes denormal results to zero
#
[OPT]
FC : -traceback
FCFLAGS : -xAVX -O2 -ip -ftz -g
# Profiling flags
#################
#
[PROFILE]
FC : -p -g
FCFLAGS : -xSSE4.2 -O2 -ip -ftz
# Debugging flags
#################
#
# -traceback : Activate backtrace on runtime
# -fpe0 : All floating point exaceptions
# -C : Checks uninitialized variables, array subscripts, etc...
# -g : Extra debugging information
# -xSSE2 : Valgrind needs a very simple x86 executable
#
[DEBUG]
FC : -g -traceback
FCFLAGS : -xSSE2 -C -fpe0 -implicitnone
# OpenMP flags
#################
#
[OPENMP]
FC : -qopenmp
IRPF90_FLAGS : --openmp

View File

@ -7,7 +7,7 @@
#
[COMMON]
FC : ifort -fpic
LAPACK_LIB : -mkl=parallel -lirc -lsvml -limf -lipps
LAPACK_LIB : -mkl=parallel
IRPF90 : irpf90
IRPF90_FLAGS : --ninja --align=32 --assert -DINTEL

27
configure vendored
View File

@ -211,9 +211,10 @@ EOF
execute << EOF
cd "\${QP_ROOT}"/external
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
cd trexio-${VERSION}
./configure --prefix=\${QP_ROOT} --without-hdf5
./configure --prefix=\${QP_ROOT} --without-hdf5 CFLAGS='-g'
make -j 8 && make -j 8 check && make -j 8 install
tar -zxvf "\${QP_ROOT}"/external/qp2-dependencies/${ARCHITECTURE}/ninja.tar.gz
mv ninja "\${QP_ROOT}"/bin/
@ -224,20 +225,34 @@ EOF
execute << EOF
cd "\${QP_ROOT}"/external
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
cd trexio-${VERSION}
./configure --prefix=\${QP_ROOT}
./configure --prefix=\${QP_ROOT} CFLAGS="-g"
make -j 8 && make -j 8 check && make -j 8 install
EOF
elif [[ ${PACKAGE} = qmckl ]] ; then
VERSION=0.5.2
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
./configure --prefix=\${QP_ROOT} --enable-hpc --disable-doc CFLAGS='-g'
make && make -j 4 check && make install
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
@ -378,13 +393,13 @@ fi
TREXIO=$(find_lib -ltrexio)
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
fi
QMCKL=$(find_lib -lqmckl)
if [[ ${QMCKL} = $(not_found) ]] ; then
error "QMCkl (qmckl) is not installed."
error "QMCkl (qmckl | qmckl-intel) is not installed."
fail
fi

File diff suppressed because it is too large Load Diff

920
data/pseudo/def2 Normal file
View File

@ -0,0 +1,920 @@
$ECP
RB-ECP GEN 28 3
1 ----- f-ul potential -----
-12.3169000 2 3.8431140
3 ----- s-f potential -----
89.5001980 2 5.0365510
0.4937610 2 1.9708490
12.3169000 2 3.8431140
3 ----- p-f potential -----
58.5689740 2 4.2583410
0.4317910 2 1.4707090
12.3169000 2 3.8431140
3 ----- d-f potential -----
26.2248980 2 3.0231270
0.9628390 2 0.6503830
12.3169000 2 3.8431140
SR-ECP GEN 28 3
1 ----- f-ul potential -----
-15.8059920 2 4.6339750
3 ----- s-f potential -----
135.4794300 2 7.4000740
17.5344630 2 3.6063790
15.8059920 2 4.6339750
3 ----- p-f potential -----
88.3597090 2 6.4848680
15.3943720 2 3.2880530
15.8059920 2 4.6339750
3 ----- d-f potential -----
29.8889870 2 4.6228410
6.6594140 2 2.2469040
15.8059920 2 4.6339750
Y-ECP GEN 28 3
2 ----- f-ul potential -----
-19.12219811 2 6.5842120
-2.43637543 2 3.2921060
4 ----- s-f potential -----
135.15384412 2 7.4880494
15.55244130 2 3.7440247
19.12219811 2 6.5842120
2.43637543 2 3.2921060
4 ----- p-f potential -----
87.78499167 2 6.4453772
11.56406599 2 3.2226886
19.12219811 2 6.5842120
2.43637543 2 3.2921060
4 ----- d-f potential -----
29.70100072 2 4.6584472
5.53996847 2 2.3292236
19.12219811 2 6.5842120
2.43637543 2 3.2921060
ZR-ECP GEN 28 3
2 ----- f-ul potential -----
-21.09377605 2 7.5400000
-3.08069427 2 3.7700000
4 ----- s-f potential -----
150.26759106 2 8.2000000
18.97621650 2 4.0897278
21.09377605 2 7.5400000
3.08069427 2 3.7700000
4 ----- p-f potential -----
99.62212372 2 7.1100000
14.16873329 2 3.5967980
21.09377605 2 7.5400000
3.08069427 2 3.7700000
4 ----- d-f potential -----
35.04512355 2 5.3500000
6.11125948 2 2.4918215
21.09377605 2 7.5400000
3.08069427 2 3.7700000
NB-ECP GEN 28 3
2 ----- f-ul potential -----
-22.92954996 2 8.4900000
-3.66630986 2 4.2500000
4 ----- s-f potential -----
165.17914349 2 8.9000000
21.99297437 2 4.4300000
22.92954996 2 8.4900000
3.66630986 2 4.2500000
4 ----- p-f potential -----
111.79441445 2 7.7700000
16.63348326 2 3.9600000
22.92954996 2 8.4900000
3.66630986 2 4.2500000
4 ----- d-f potential -----
38.11224880 2 6.0500000
8.03916727 2 2.8400000
22.92954996 2 8.4900000
3.66630986 2 4.2500000
MO-ECP GEN 28 3
2 ----- f-ul potential -----
-24.80517707 2 9.4500000
-4.15378155 2 4.7200000
4 ----- s-f potential -----
180.10310850 2 9.7145938
24.99722791 2 4.6805004
24.80517707 2 9.4500000
4.15378155 2 4.7200000
4 ----- p-f potential -----
123.77275231 2 8.1421366
19.53022800 2 4.6259863
24.80517707 2 9.4500000
4.15378155 2 4.7200000
4 ----- d-f potential -----
48.37502229 2 6.6184148
8.89205274 2 3.2487516
24.80517707 2 9.4500000
4.15378155 2 4.7200000
TC-ECP GEN 28 3
2 ----- f-ul potential -----
-26.56244747 2 10.4000000
-4.58568054 2 5.2000000
4 ----- s-f potential -----
195.15916591 2 10.4223462
28.09260333 2 5.0365160
26.56244747 2 10.4000000
4.58568054 2 5.2000000
4 ----- p-f potential -----
135.28456622 2 8.9504494
21.80650430 2 4.8544394
26.56244747 2 10.4000000
4.58568054 2 5.2000000
4 ----- d-f potential -----
54.32972942 2 6.9456968
11.15506795 2 3.9705849
26.56244747 2 10.4000000
4.58568054 2 5.2000000
RU-ECP GEN 28 3
2 ----- f-ul potential -----
-28.34061627 2 11.3600000
-4.94462923 2 5.6800000
4 ----- s-f potential -----
209.82297122 2 11.1052693
30.65472642 2 5.4147454
28.34061627 2 11.3600000
4.94462923 2 5.6800000
4 ----- p-f potential -----
146.33618228 2 9.7712707
24.12787723 2 5.0739908
28.34061627 2 11.3600000
4.94462923 2 5.6800000
4 ----- d-f potential -----
67.51589667 2 7.6714231
9.87010415 2 4.1365647
28.34061627 2 11.3600000
4.94462923 2 5.6800000
RH-ECP GEN 28 3
2 ----- f-ul potential -----
-30.09345572 2 12.3100000
-5.21848192 2 6.1600000
4 ----- s-f potential -----
225.34775353 2 11.7200000
32.82318898 2 5.8200000
30.09345572 2 12.3100000
5.21848192 2 6.1600000
4 ----- p-f potential -----
158.70941159 2 10.4200000
26.44410049 2 5.4500000
30.09345572 2 12.3100000
5.21848192 2 6.1600000
4 ----- d-f potential -----
62.75862572 2 8.8200000
10.97871947 2 3.8700000
30.09345572 2 12.3100000
5.21848192 2 6.1600000
PD-ECP GEN 28 3
2 ----- f-ul potential -----
-31.92955431 2 13.2700000
-5.39821694 2 6.6300000
4 ----- s-f potential -----
240.22904033 2 12.4300000
35.17194347 2 6.1707594
31.92955431 2 13.2700000
5.39821694 2 6.6300000
4 ----- p-f potential -----
170.41727605 2 11.0800000
28.47213287 2 5.8295541
31.92955431 2 13.2700000
5.39821694 2 6.6300000
4 ----- d-f potential -----
69.01384488 2 9.5100000
11.75086158 2 4.1397811
31.92955431 2 13.2700000
5.39821694 2 6.6300000
AG-ECP GEN 28 3
2 ----- f-ul potential -----
-33.68992012 2 14.2200000
-5.53112021 2 7.1100000
4 ----- s-f potential -----
255.13936452 2 13.1300000
36.86612154 2 6.5100000
33.68992012 2 14.2200000
5.53112021 2 7.1100000
4 ----- p-f potential -----
182.18186871 2 11.7400000
30.35775148 2 6.2000000
33.68992012 2 14.2200000
5.53112021 2 7.1100000
4 ----- d-f potential -----
73.71926087 2 10.2100000
12.50211712 2 4.3800000
33.68992012 2 14.2200000
5.53112021 2 7.1100000
CD-ECP GEN 28 3
2 ----- f-ul potential -----
-35.47662555 2 15.1847957
-5.61767685 2 7.5923978
4 ----- s-f potential -----
270.00948324 2 13.8358689
38.76730798 2 6.8572704
35.47662555 2 15.1847957
5.61767685 2 7.5923978
4 ----- p-f potential -----
193.82962939 2 12.4049710
31.89652523 2 6.5677995
35.47662555 2 15.1847957
5.61767685 2 7.5923978
4 ----- d-f potential -----
79.19364700 2 10.8969253
13.23082674 2 4.6411649
35.47662555 2 15.1847957
5.61767685 2 7.5923978
IN-ECP GEN 28 3
2 ----- f-ul potential -----
-13.72807800 2 12.53905600
-18.20686600 2 12.55256100
4 ----- s-f potential -----
281.12235000 2 15.39282200
61.90147000 2 8.05586400
13.72807800 2 12.53905600
18.20686600 2 12.55256100
6 ----- p-f potential -----
67.46215400 2 13.92867200
134.94925000 2 13.34723400
14.74614000 2 7.61413200
29.63926200 2 7.31836500
13.72807800 2 12.53905600
18.20686600 2 12.55256100
6 ----- d-f potential -----
35.49325400 2 14.03471500
53.17877300 2 14.51161600
9.17728100 2 5.55055000
12.39241000 2 5.05941500
13.72807800 2 12.53905600
18.20686600 2 12.55256100
SN-ECP GEN 28 3
2 ----- f-ul potential -----
-12.57633300 2 12.28234800
-16.59594400 2 12.27215000
4 ----- s-f potential -----
279.98868200 2 17.42041400
62.37781000 2 7.63115500
12.57633300 2 12.28234800
16.59594400 2 12.27215000
6 ----- p-f potential -----
66.16252300 2 16.13102400
132.17439600 2 15.62807700
16.33941700 2 7.32560800
32.48895900 2 6.94251900
12.57633300 2 12.28234800
16.59594400 2 12.27215000
6 ----- d-f potential -----
36.38744100 2 15.51497600
54.50784100 2 15.18816000
8.69682300 2 5.45602400
12.84020800 2 5.36310500
12.57633300 2 12.28234800
16.59594400 2 12.27215000
SB-ECP GEN 28 3
2 ----- f-ul potential -----
-15.36680100 2 14.44497800
-20.29613800 2 14.44929500
4 ----- s-f potential -----
281.07158100 2 16.33086500
61.71660400 2 8.55654200
15.36680100 2 14.44497800
20.29613800 2 14.44929500
6 ----- p-f potential -----
67.45738000 2 14.47033700
134.93350300 2 13.81619400
14.71634400 2 8.42492400
29.51851200 2 8.09272800
15.36680100 2 14.44497800
20.29613800 2 14.44929500
6 ----- d-f potential -----
35.44781500 2 14.88633100
53.14346600 2 15.14631900
9.17922300 2 5.90826700
13.24025300 2 5.59432200
15.36680100 2 14.44497800
20.29613800 2 14.44929500
TE-ECP GEN 28 3
2 ----- f-ul potential -----
-15.74545000 2 15.20616800
-20.74244800 2 15.20170200
4 ----- s-f potential -----
281.04584300 2 16.81447300
61.62065600 2 8.79352600
15.74545000 2 15.20616800
20.74244800 2 15.20170200
6 ----- p-f potential -----
67.44946400 2 14.87780100
134.90430400 2 14.26973100
14.68954700 2 8.72443500
29.41506300 2 8.29151500
15.74545000 2 15.20616800
20.74244800 2 15.20170200
6 ----- d-f potential -----
35.43205700 2 15.20500800
53.13568700 2 15.22584800
9.06980200 2 6.07176900
13.12230400 2 5.80476000
15.74545000 2 15.20616800
20.74244800 2 15.20170200
I-ECP GEN 28 3
4 ----- f-ul potential -----
-21.84204000 2 19.45860900
-28.46819100 2 19.34926000
-0.24371300 2 4.82376700
-0.32080400 2 4.88431500
7 ----- s-f potential -----
49.99429300 2 40.01583500
281.02531700 2 17.42974700
61.57332600 2 9.00548400
21.84204000 2 19.45860900
28.46819100 2 19.34926000
0.24371300 2 4.82376700
0.32080400 2 4.88431500
8 ----- p-f potential -----
67.44284100 2 15.35546600
134.88113700 2 14.97183300
14.67505100 2 8.96016400
29.37566600 2 8.25909600
21.84204000 2 19.45860900
28.46819100 2 19.34926000
0.24371300 2 4.82376700
0.32080400 2 4.88431500
10 ----- d-f potential -----
35.43952900 2 15.06890800
53.17605700 2 14.55532200
9.06719500 2 6.71864700
13.20693700 2 6.45639300
0.08933500 2 1.19177900
0.05238000 2 1.29115700
21.84204000 2 19.45860900
28.46819100 2 19.34926000
0.24371300 2 4.82376700
0.32080400 2 4.88431500
XE-ECP GEN 28 3
4 ----- f-ul potential -----
-23.08929500 2 20.88155700
-30.07447500 2 20.78344300
-0.28822700 2 5.25338900
-0.38692400 2 5.36118800
7 ----- s-f potential -----
49.99796200 2 40.00518400
281.01330300 2 17.81221400
61.53825500 2 9.30415000
23.08929500 2 20.88155700
30.07447500 2 20.78344300
0.28822700 2 5.25338900
0.38692400 2 5.36118800
8 ----- p-f potential -----
67.43914200 2 15.70177200
134.87471100 2 15.25860800
14.66330000 2 9.29218400
29.35473000 2 8.55900300
23.08929500 2 20.88155700
30.07447500 2 20.78344300
0.28822700 2 5.25338900
0.38692400 2 5.36118800
10 ----- d-f potential -----
35.43690800 2 15.18560000
53.19577200 2 14.28450000
9.04623200 2 7.12188900
13.22368100 2 6.99196300
0.08485300 2 0.62394600
0.04415500 2 0.64728400
23.08929500 2 20.88155700
30.07447500 2 20.78344300
0.28822700 2 5.25338900
0.38692400 2 5.36118800
CS-ECP GEN 46 3
1 ----- f-ul potential -----
-28.8843090 2 3.1232690
3 ----- s-f potential -----
84.5477300 2 4.0797500
16.6541730 2 2.4174060
28.8843090 2 3.1232690
3 ----- p-f potential -----
157.0490590 2 5.5140800
26.4233070 2 2.1603160
28.8843090 2 3.1232690
3 ----- d-f potential -----
13.1727530 2 1.8074100
3.3428330 2 0.8581820
28.8843090 2 3.1232690
BA-ECP GEN 46 3
1 ----- f-ul potential -----
-33.4731740 2 3.5894650
3 ----- s-f potential -----
427.8458160 2 9.5269860
204.4175300 2 4.4875100
33.4731740 2 3.5894650
3 ----- p-f potential -----
293.6058640 2 8.3159300
294.1933160 2 4.2922170
33.4731740 2 3.5894650
3 ----- d-f potential -----
112.5504020 2 5.9161080
181.7826210 2 2.8748420
33.4731740 2 3.5894650
LA-ECP GEN 46 3
1 ----- f-ul potential -----
-36.0100160 2 4.0286000
3 ----- s-f potential -----
91.9321770 2 3.3099000
-3.7887640 2 1.6550000
36.0100160 2 4.0286000
3 ----- p-f potential -----
63.7594860 2 2.8368000
-0.6479580 2 1.4184000
36.0100160 2 4.0286000
3 ----- d-f potential -----
36.1161730 2 2.0213000
0.2191140 2 1.0107000
36.0100160 2 4.0286000
CE-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
580.08345700 2 20.13782900
1 ----- p-h potential -----
310.30283300 2 15.99848200
1 ----- d-h potential -----
167.81394400 2 14.97418700
1 ----- f-h potential -----
-49.39022900 2 23.40245500
1 ----- g-h potential -----
-21.33187900 2 16.57055300
PR-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
577.57312200 2 20.76627800
1 ----- p-h potential -----
295.78584600 2 16.07844800
1 ----- d-h potential -----
150.86705500 2 14.70508900
1 ----- f-h potential -----
-48.73676600 2 23.37896900
1 ----- g-h potential -----
-22.32948800 2 17.44713800
ND-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
574.37098000 2 21.35226700
1 ----- p-h potential -----
280.94644000 2 16.11926500
1 ----- d-h potential -----
138.67062700 2 14.49410300
1 ----- f-h potential -----
-47.52266800 2 23.18386000
1 ----- g-h potential -----
-23.34458700 2 18.34417400
PM-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
575.39574900 2 21.94286500
1 ----- p-h potential -----
281.70451400 2 16.55516100
1 ----- d-h potential -----
123.52473700 2 13.96030800
1 ----- f-h potential -----
-50.74151100 2 24.03354600
1 ----- g-h potential -----
-24.37251000 2 19.26024500
SM-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
572.98533200 2 22.34447100
1 ----- p-h potential -----
272.35914500 2 16.69459000
1 ----- d-h potential -----
115.29390000 2 13.72770500
1 ----- f-h potential -----
-51.10839200 2 24.05909200
1 ----- g-h potential -----
-25.42188500 2 20.19724900
EU-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
607.65933100 2 23.47138400
1 ----- p-h potential -----
264.38547600 2 16.77247900
1 ----- d-h potential -----
115.38137500 2 13.98134300
1 ----- f-h potential -----
-49.40079400 2 23.96288800
1 ----- g-h potential -----
-26.74827300 2 21.23245800
GD-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
637.20086900 2 24.60215100
1 ----- p-h potential -----
261.68960100 2 16.88925000
1 ----- d-h potential -----
106.85653300 2 13.64335800
1 ----- f-h potential -----
-50.68359000 2 24.12691700
1 ----- g-h potential -----
-27.57963000 2 22.13188700
TB-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
668.59715500 2 24.95295600
1 ----- p-h potential -----
266.98047500 2 17.61089900
1 ----- d-h potential -----
97.50659600 2 12.97600900
1 ----- f-h potential -----
-52.17575700 2 24.24886900
1 ----- g-h potential -----
-28.69426800 2 23.13067200
DY-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
705.67122100 2 26.42958600
1 ----- p-h potential -----
254.86698900 2 17.31703400
1 ----- d-h potential -----
95.04518700 2 12.91359900
1 ----- f-h potential -----
-54.57409300 2 24.90787800
1 ----- g-h potential -----
-29.82827700 2 24.14875300
HO-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
755.70313600 2 28.39725700
1 ----- p-h potential -----
253.55199800 2 17.43863300
1 ----- d-h potential -----
89.63567700 2 12.43421200
1 ----- f-h potential -----
-55.48203600 2 25.38701000
1 ----- g-h potential -----
-30.99112500 2 25.18850100
ER-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
800.95287600 2 29.79859200
1 ----- p-h potential -----
262.01986900 2 18.11423700
1 ----- d-h potential -----
80.17055200 2 11.36958700
1 ----- f-h potential -----
-42.33628500 2 21.82123300
1 ----- g-h potential -----
-32.18527800 2 26.25073500
TM-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
845.51074300 2 31.14412200
1 ----- p-h potential -----
258.58523900 2 18.09235300
1 ----- d-h potential -----
80.72905900 2 11.46915900
1 ----- f-h potential -----
-48.70126600 2 23.60554400
1 ----- g-h potential -----
-33.39549600 2 27.32978100
YB-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
891.01377700 2 32.42448400
1 ----- p-h potential -----
264.03695300 2 18.65623200
1 ----- d-h potential -----
73.92391900 2 10.49022200
1 ----- f-h potential -----
-39.59217300 2 20.77418300
1 ----- g-h potential -----
-34.63863800 2 28.43102800
LU-ECP GEN 28 5
1 ----- h-ul potential -----
0.00000000 2 1.00000000
1 ----- s-h potential -----
989.99558400 2 35.16209700
1 ----- p-h potential -----
278.86565200 2 19.46440200
1 ----- d-h potential -----
71.00917800 2 10.00686500
1 ----- f-h potential -----
-47.40589000 2 23.51793200
1 ----- g-h potential -----
-35.55714600 2 29.41223800
HF-ECP GEN 60 3
1 ----- f-ul potential -----
10.04672251 2 1.78576984
3 ----- s-f potential -----
1499.28471073 2 14.76995900
40.28210136 2 7.38497940
-10.04672251 2 1.78576984
3 ----- p-f potential -----
397.73300533 2 9.84948950
19.31640586 2 4.92474450
-10.04672251 2 1.78576984
3 ----- d-f potential -----
101.32980526 2 6.09675640
5.87343821 2 3.04837820
-10.04672251 2 1.78576984
TA-ECP GEN 60 3
1 ----- f-ul potential -----
12.01796094 2 2.01788111
3 ----- s-f potential -----
1345.88064703 2 14.54640770
36.76680620 2 7.27320380
-12.01796094 2 2.01788111
3 ----- p-f potential -----
378.42530145 2 9.93556529
22.29309086 2 4.96778243
-12.01796094 2 2.01788111
3 ----- d-f potential -----
104.88395571 2 6.34737691
8.75584805 2 3.17368846
-12.01796094 2 2.01788111
W-ECP GEN 60 3
1 ----- f-ul potential -----
14.15257947 2 2.25888846
3 ----- s-f potential -----
1192.39588226 2 14.32285640
32.52293315 2 7.16142810
-14.15257947 2 2.25888846
3 ----- p-f potential -----
359.03196711 2 10.02164110
24.03038019 2 5.01082040
-14.15257947 2 2.25888846
3 ----- d-f potential -----
108.30134897 2 6.59799743
10.98252827 2 3.29899871
-14.15257947 2 2.25888846
RE-ECP GEN 60 3
1 ----- f-ul potential -----
16.44985227 2 2.50865059
3 ----- s-f potential -----
1038.95157226 2 14.09930510
29.56173830 2 7.04965250
-16.44985227 2 2.50865059
3 ----- p-f potential -----
339.54350965 2 10.10771690
24.91369646 2 5.05385830
-16.44985227 2 2.50865059
3 ----- d-f potential -----
111.69965275 2 6.84861794
12.62432927 2 3.42430897
-16.44985227 2 2.50865059
OS-ECP GEN 60 3
1 ----- f-ul potential -----
18.90945701 2 2.76707510
3 ----- s-f potential -----
885.40571914 2 13.87575390
25.96704014 2 6.93787690
-18.90945701 2 2.76707510
3 ----- p-f potential -----
320.08390185 2 10.19379260
26.14876493 2 5.09689620
-18.90945701 2 2.76707510
3 ----- d-f potential -----
115.04484313 2 7.09923846
13.62257457 2 3.54961923
-18.90945701 2 2.76707510
IR-ECP GEN 60 3
1 ----- f-ul potential -----
21.53103107 2 3.03407192
3 ----- s-f potential -----
732.26919978 2 13.65220260
26.48472087 2 6.82610130
-21.53103107 2 3.03407192
3 ----- p-f potential -----
299.48947357 2 10.27986840
26.46623354 2 5.13993410
-21.53103107 2 3.03407192
3 ----- d-f potential -----
124.45759451 2 7.34985897
14.03599518 2 3.67492949
-21.53103107 2 3.03407192
PT-ECP GEN 60 3
1 ----- f-ul potential -----
24.31437573 2 3.30956857
3 ----- s-f potential -----
579.22386092 2 13.42865130
29.66949062 2 6.71432560
-24.31437573 2 3.30956857
3 ----- p-f potential -----
280.86077422 2 10.36594420
26.74538204 2 5.18297210
-24.31437573 2 3.30956857
3 ----- d-f potential -----
120.39644429 2 7.60047949
15.81092058 2 3.80023974
-24.31437573 2 3.30956857
AU-ECP GEN 60 3
2 ----- f-ul potential -----
30.49008890 2 4.78982000
5.17107381 2 2.39491000
4 ----- s-f potential -----
426.84667920 2 13.20510000
37.00708285 2 6.60255000
-30.49008890 2 4.78982000
-5.17107381 2 2.39491000
4 ----- p-f potential -----
261.19958038 2 10.45202000
26.96249604 2 5.22601000
-30.49008890 2 4.78982000
-5.17107381 2 2.39491000
4 ----- d-f potential -----
124.79066561 2 7.85110000
16.30072573 2 3.92555000
-30.49008890 2 4.78982000
-5.17107381 2 2.39491000
HG-ECP GEN 60 3
1 ----- f-ul potential -----
30.36499643 2 3.88579112
3 ----- s-f potential -----
275.73721174 2 12.98154870
49.08921249 2 6.49077440
-30.36499643 2 3.88579112
3 ----- p-f potential -----
241.54007398 2 10.53809580
27.39659081 2 5.26904790
-30.36499643 2 3.88579112
3 ----- d-f potential -----
127.86700761 2 8.10172051
16.60831151 2 4.05086026
-30.36499643 2 3.88579112
TL-ECP GEN 60 3
4 ----- f-ul potential -----
15.82548800 2 5.62639900
21.10402100 2 5.54895200
2.91512700 2 2.87494600
3.89690300 2 2.82145100
6 ----- s-f potential -----
281.28466300 2 12.16780500
62.43425100 2 8.29490900
-15.82548800 2 5.62639900
-21.10402100 2 5.54895200
-2.91512700 2 2.87494600
-3.89690300 2 2.82145100
8 ----- p-f potential -----
4.63340800 2 7.15149200
9.34175600 2 5.17286500
72.29925300 2 9.89107200
144.55803700 2 9.00339100
-15.82548800 2 5.62639900
-21.10402100 2 5.54895200
-2.91512700 2 2.87494600
-3.89690300 2 2.82145100
8 ----- d-f potential -----
35.94303900 2 7.13021800
53.90959300 2 6.92690600
10.38193900 2 5.41757000
15.58382200 2 5.13868100
-15.82548800 2 5.62639900
-21.10402100 2 5.54895200
-2.91512700 2 2.87494600
-3.89690300 2 2.82145100
PB-ECP GEN 60 3
2 ----- f-ul potential -----
12.20989200 2 3.88751200
16.19029100 2 3.81196300
4 ----- s-f potential -----
281.28549900 2 12.29630300
62.52021700 2 8.63263400
-12.20989200 2 3.88751200
-16.19029100 2 3.81196300
6 ----- p-f potential -----
72.27689700 2 10.24179000
144.59108300 2 8.92417600
4.75869300 2 6.58134200
9.94062100 2 6.25540300
-12.20989200 2 3.88751200
-16.19029100 2 3.81196300
6 ----- d-f potential -----
35.84850700 2 7.75433600
53.72434200 2 7.72028100
10.11525600 2 4.97026400
14.83373100 2 4.56378900
-12.20989200 2 3.88751200
-16.19029100 2 3.81196300
BI-ECP GEN 60 3
2 ----- f-ul potential -----
13.71338300 2 4.21454600
18.19430800 2 4.13340000
4 ----- s-f potential -----
283.26422700 2 13.04309000
62.47195900 2 8.22168200
-13.71338300 2 4.21454600
-18.19430800 2 4.13340000
6 ----- p-f potential -----
72.00149900 2 10.46777700
144.00227700 2 9.11890100
5.00794500 2 6.75479100
9.99155000 2 6.25259200
-13.71338300 2 4.21454600
-18.19430800 2 4.13340000
6 ----- d-f potential -----
36.39625900 2 8.08147400
54.59766400 2 7.89059500
9.98429400 2 4.95555600
14.98148500 2 4.70455900
-13.71338300 2 4.21454600
-18.19430800 2 4.13340000
PO-ECP GEN 60 3
4 ----- f-ul potential -----
17.42829500 2 5.01327000
23.38035300 2 4.98464000
0.16339200 2 1.32676000
0.32456600 2 1.52875800
6 ----- s-f potential -----
283.24470600 2 13.27722700
62.39646100 2 8.39951800
-17.42829500 2 5.01327000
-23.38035300 2 4.98464000
-0.16339200 2 1.32676000
-0.32456600 2 1.52875800
8 ----- p-f potential -----
71.99171600 2 10.66568200
143.97187100 2 9.28375300
4.94961500 2 6.87274900
9.74049900 2 6.32615000
-17.42829500 2 5.01327000
-23.38035300 2 4.98464000
-0.16339200 2 1.32676000
-0.32456600 2 1.52875800
8 ----- d-f potential -----
36.37838300 2 8.21486600
54.56271500 2 8.00869600
9.88949900 2 5.05522700
14.69387700 2 4.78255300
-17.42829500 2 5.01327000
-23.38035300 2 4.98464000
-0.16339200 2 1.32676000
-0.32456600 2 1.52875800
AT-ECP GEN 60 3
4 ----- f-ul potential -----
19.87019800 2 5.81216300
26.41645200 2 5.75371500
0.99497000 2 2.51347200
1.49070100 2 2.53626100
7 ----- s-f potential -----
49.95715800 2 30.20083200
283.21037100 2 13.61230600
62.28105200 2 8.52934000
-19.87019800 2 5.81216300
-26.41645200 2 5.75371500
-0.99497000 2 2.51347200
-1.49070100 2 2.53626100
8 ----- p-f potential -----
71.98237100 2 10.85406500
143.90353200 2 9.46822900
4.87175900 2 7.03111400
8.98305900 2 6.14385800
-19.87019800 2 5.81216300
-26.41645200 2 5.75371500
-0.99497000 2 2.51347200
-1.49070100 2 2.53626100
8 ----- d-f potential -----
36.36323700 2 8.31351500
54.54897000 2 7.99896500
9.77628500 2 5.17996600
14.26475500 2 4.94222600
-19.87019800 2 5.81216300
-26.41645200 2 5.75371500
-0.99497000 2 2.51347200
-1.49070100 2 2.53626100
RN-ECP GEN 60 3
4 ----- f-ul potential -----
21.79729000 2 6.34857100
28.94680500 2 6.29594900
1.44736500 2 2.88211800
2.17796400 2 2.90804800
7 ----- s-f potential -----
49.96555100 2 30.15124200
283.07000000 2 14.52124100
62.00287000 2 8.05203800
-21.79729000 2 6.34857100
-28.94680500 2 6.29594900
-1.44736500 2 2.88211800
-2.17796400 2 2.90804800
8 ----- p-f potential -----
71.96911900 2 11.00994200
143.86055900 2 9.61762500
4.71476100 2 7.33600800
9.01306500 2 6.40625300
-21.79729000 2 6.34857100
-28.94680500 2 6.29594900
-1.44736500 2 2.88211800
-2.17796400 2 2.90804800
8 ----- d-f potential -----
36.36836500 2 8.36922000
54.55176100 2 8.11697500
9.63448700 2 5.35365600
14.38790200 2 5.09721200
-21.79729000 2 6.34857100
-28.94680500 2 6.29594900
-1.44736500 2 2.88211800
-2.17796400 2 2.90804800
$END

2
external/ezfio vendored

@ -1 +1 @@
Subproject commit d5805497fa0ef30e70e055cde1ecec2963303e93
Subproject commit dba01c4fe0ff7b84c5ecfb1c7c77ec68781311b3

View File

@ -38,9 +38,8 @@ def comp_path(path):
from qp_path import QP_ROOT, QP_SRC, QP_EZFIO
LIB = " -lz -ltrexio"
LIB = " -lz"
EZFIO_LIB = join("$QP_ROOT", "lib", "libezfio_irp.a")
ZMQ_LIB = join("$QP_ROOT", "lib", "libf77zmq.a") + " " + join("$QP_ROOT", "lib", "libzmq.a") + " -lstdc++ -lrt -ldl"
ROOT_BUILD_NINJA = join("$QP_ROOT", "config", "build.ninja")
ROOT_BUILD_NINJA_EXP = join(QP_ROOT, "config", "build.ninja")
ROOT_BUILD_NINJA_EXP_tmp = join(QP_ROOT, "config", "build.ninja.tmp")
@ -118,7 +117,7 @@ def ninja_create_env_variable(pwd_config_file):
lib_lapack = get_compilation_option(pwd_config_file, "LAPACK_LIB")
lib_usr = get_compilation_option(pwd_config_file, "LIB")
str_lib = " ".join([lib_lapack, EZFIO_LIB, ZMQ_LIB, LIB, lib_usr])
str_lib = " ".join([lib_lapack, EZFIO_LIB, LIB, lib_usr])
# Read all LIB files in modules
for directory in [real_join(QP_SRC, m) for m in sorted(os.listdir(QP_SRC))]:

View File

@ -829,4 +829,8 @@ if __name__ == "__main__":
# _|
for (m, dict_ezfio_cfg) in l_dict_ezfio_cfg:
if dict_ezfio_cfg == {}:
print("Error: Empty EZFIO.cfg in ", arguments["--path_module"])
sys.exit(-1)
code_generation(arguments, dict_ezfio_cfg, m)

69
scripts/import_champ_jastrow.py Executable file
View File

@ -0,0 +1,69 @@
#!/usr/bin/env python3
conv = [ 0, 0, 2 , 6 , 13 , 23 , 37 , 55 , 78 , 106 , 140 ]
def import_jastrow(jastrow_filename):
with open(jastrow_filename,'r') as jastrow_file:
lines = [ line.strip() for line in jastrow_file.readlines() ]
lines = [ line for line in lines if line != "" ]
start = 0
end = len(lines)
for i,line in enumerate(lines):
if line.startswith("jastrow_parameter"):
start = i
elif line.startswith("end"):
end = i
lines = lines[start:end]
type_num = (len(lines)-4)//2
nord_a,nord_b,nord_c = [ int(i) for i in lines[1].split()[:3] ]
scale_k = float(lines[2].split()[0])
vec_a = []
for j in range(type_num):
vec_a += [ float(i) for i in lines[3+j].split()[:nord_a+1] ]
vec_b = [ float(i) for i in lines[3+type_num].split()[:nord_b+1] ]
vec_c = []
for j in range(type_num):
vec_c += [ float(i) for i in lines[4+type_num+j].split()[:conv[nord_c]] ]
return {
'type_num' : type_num,
'scale_k' : scale_k,
'nord_a' : nord_a,
'nord_b' : nord_b,
'nord_c' : nord_c,
'vec_a' : vec_a,
'vec_b' : vec_b,
'vec_c' : vec_c,
}
if __name__ == '__main__':
import sys
from ezfio import ezfio
ezfio.set_file(sys.argv[1])
jastrow_file = sys.argv[2]
jastrow = import_jastrow(jastrow_file)
print (jastrow)
ezfio.set_jastrow_jast_type("Qmckl")
ezfio.set_jastrow_jast_qmckl_type_nucl_num(jastrow['type_num'])
charges = ezfio.get_nuclei_nucl_charge()
types = {}
k = 1
for c in charges:
if c not in types:
types[c] = k
k += 1
type_nucl_vector = [types[c] for c in charges]
print(type_nucl_vector)
ezfio.set_jastrow_jast_qmckl_type_nucl_vector(type_nucl_vector)
ezfio.set_jastrow_jast_qmckl_rescale_ee(jastrow['scale_k'])
ezfio.set_jastrow_jast_qmckl_rescale_en([jastrow['scale_k'] for i in type_nucl_vector])
ezfio.set_jastrow_jast_qmckl_aord_num(jastrow['nord_a'])
ezfio.set_jastrow_jast_qmckl_bord_num(jastrow['nord_b'])
ezfio.set_jastrow_jast_qmckl_cord_num(jastrow['nord_c'])
ezfio.set_jastrow_jast_qmckl_c_vector_size(len(jastrow['vec_c']))
ezfio.set_jastrow_jast_qmckl_a_vector(jastrow['vec_a'])
ezfio.set_jastrow_jast_qmckl_b_vector(jastrow['vec_b'])
ezfio.set_jastrow_jast_qmckl_c_vector(jastrow['vec_c'])

View File

@ -115,9 +115,7 @@ def get_l_module_descendant(d_child, l_module):
except KeyError:
print("Error: ", file=sys.stderr)
print("`{0}` is not a submodule".format(module), file=sys.stderr)
print("Check the typo (spelling, case, '/', etc.) ", file=sys.stderr)
# pass
sys.exit(1)
raise
return list(set(l))

View File

@ -38,6 +38,15 @@ else:
QP_ROOT + "/install",
QP_ROOT + "/scripts"] + sys.path
def uint64_to_int64(u):
# Check if the most significant bit is set
if u & (1 << 63):
# Calculate the two's complement
result = -int(np.bitwise_not(np.uint64(u))+1)
else:
# The number is already positive
result = u
return result
def generate_xyz(l):
@ -454,19 +463,40 @@ def write_ezfio(trexio_filename, filename):
else:
print("None")
print("Determinant\t\t...\t", end=' ')
print("Determinant\t...\t", end=' ')
alpha = [ i for i in range(num_alpha) ]
beta = [ i for i in range(num_beta) ]
if trexio.has_mo_spin(trexio_file):
spin = trexio.read_mo_spin(trexio_file)
if max(spin) == 1:
beta = [ i for i in range(mo_num) if spin[i] == 1 ]
alpha = [ i for i in range(len(spin)) if spin[i] == 0 ]
alpha = [ alpha[i] for i in range(num_alpha) ]
beta = [ i for i in range(len(spin)) if spin[i] == 1 ]
beta = [ beta[i] for i in range(num_beta) ]
print("Warning -- UHF orbitals --", end=' ')
alpha_s = ['0']*mo_num
beta_s = ['0']*mo_num
for i in alpha:
alpha_s[i] = '1'
for i in beta:
beta_s[i] = '1'
alpha_s = ''.join(alpha_s)[::-1]
beta_s = ''.join(beta_s)[::-1]
def conv(i):
try:
result = np.int64(i)
except:
result = np.int64(i-2**63-1)
return result
alpha = qp_bitmasks.BitMask(alpha)
beta = qp_bitmasks.BitMask(beta )
print(alpha)
print(beta)
alpha = [ uint64_to_int64(int(i,2)) for i in qp_bitmasks.string_to_bitmask(alpha_s) ][::-1]
beta = [ uint64_to_int64(int(i,2)) for i in qp_bitmasks.string_to_bitmask(beta_s ) ][::-1]
ezfio.set_determinants_bit_kind(8)
ezfio.set_determinants_n_int(1+mo_num//64)
ezfio.set_determinants_n_det(1)
ezfio.set_determinants_n_states(1)
ezfio.set_determinants_psi_det(alpha+beta)
ezfio.set_determinants_psi_coef([[1.0]])
print("OK")

View File

@ -52,7 +52,7 @@
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,c,d_a_2,d_2,deriv_tmp, &
!$OMP alpha, beta, n, l, i,j,c,d_a_2,d_2,deriv_tmp, &
!$OMP overlap_x0,overlap_y0,overlap_z0) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_deriv2_x,ao_deriv2_y,ao_deriv2_z,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &

View File

@ -1,4 +1,4 @@
ao_two_e_erf_ints
ao_two_e_ints
mo_one_e_ints
ao_many_one_e_ints
dft_utils_in_r

View File

@ -1,13 +0,0 @@
[io_ao_two_e_integrals_erf]
type: Disk_access
doc: Read/Write |AO| integrals with the long range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[mu_erf]
type: double precision
doc: cutting of the interaction in the range separated model
interface: ezfio,provider,ocaml
default: 0.5
ezfio_name: mu_erf

View File

@ -1 +0,0 @@
ao_two_e_ints

View File

@ -1,19 +0,0 @@
======================
ao_two_e_erf_ints
======================
Here, all two-electron integrals (:math:`erf(\mu r_{12})/r_{12}`) are computed.
As they have 4 indices and many are zero, they are stored in a map, as defined
in :file:`utils/map_module.f90`.
The main parameter of this module is :option:`ao_two_e_erf_ints mu_erf` which is the range-separation parameter.
To fetch an |AO| integral, use the
`get_ao_two_e_integral_erf(i,j,k,l,ao_integrals_erf_map)` function.
The conventions are:
* For |AO| integrals : (ij|kl) = (11|22) = <ik|jl> = <12|12>

View File

@ -35,3 +35,15 @@ type: logical
doc: Perform Cholesky decomposition of AO integrals
interface: ezfio,provider,ocaml
default: False
[io_ao_two_e_integrals_erf]
type: Disk_access
doc: Read/Write |AO| erf integrals from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None
[use_only_lr]
type: logical
doc: If true, use only the long range part of the two-electron integrals instead of 1/r12
interface: ezfio, provider, ocaml
default: False

View File

@ -1,3 +1,4 @@
hamiltonian
ao_one_e_ints
pseudo
bitmask

View File

@ -90,7 +90,7 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_erf_in_map ]
if (write_ao_two_e_integrals_erf) then
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
call ezfio_set_ao_two_e_erf_ints_io_ao_two_e_integrals_erf("Read")
call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals_erf('Read')
endif
END_PROVIDER

View File

@ -4,7 +4,7 @@ subroutine save_erf_two_e_integrals_ao
PROVIDE ao_two_e_integrals_erf_in_map
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
call ezfio_set_ao_two_e_erf_ints_io_ao_two_e_integrals_erf('Read')
call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals_erf('Read')
end
subroutine save_erf_two_e_ints_ao_into_ints_ao

View File

@ -21,9 +21,9 @@ double precision function ao_two_e_integral(i, j, k, l)
double precision :: P_new(0:max_dim,3),P_center(3),fact_p,pp
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
double precision :: ao_two_e_integral_schwartz_accel
double precision :: ao_two_e_integral_cosgtos
double precision, external :: ao_two_e_integral_erf
double precision, external :: ao_two_e_integral_cosgtos
double precision, external :: ao_two_e_integral_schwartz_accel
if(use_cosgtos) then
@ -31,13 +31,15 @@ double precision function ao_two_e_integral(i, j, k, l)
ao_two_e_integral = ao_two_e_integral_cosgtos(i, j, k, l)
else
else if (use_only_lr) then
if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
ao_two_e_integral = ao_two_e_integral_erf(i, j, k, l)
else if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l)
else
else
dim1 = n_pt_max_integrals
@ -117,8 +119,6 @@ double precision function ao_two_e_integral(i, j, k, l)
enddo ! q
enddo ! p
endif
endif
endif

View File

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

View File

@ -18,12 +18,13 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_direct_bi_ort, (mo_num, mo_num,
double precision :: integral, wall1, wall0
PROVIDE mo_l_coef mo_r_coef
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
three_e_3_idx_direct_bi_ort = 0.d0
print *, ' Providing the three_e_3_idx_direct_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
@ -79,6 +80,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_1_bi_ort, (mo_num, mo_num
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
@ -135,6 +137,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_cycle_2_bi_ort, (mo_num, mo_num
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
@ -191,6 +194,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch23_bi_ort, (mo_num, mo_num,
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
@ -247,6 +251,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch13_bi_ort, (mo_num, mo_num,
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
@ -303,6 +308,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort, (mo_num, mo_num,
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &
@ -349,6 +355,7 @@ BEGIN_PROVIDER [ double precision, three_e_3_idx_exch12_bi_ort_new, (mo_num, mo_
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
call give_integrals_3_body_bi_ort(1, 1, 1, 1, 1, 1, integral)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,m,integral) &

View File

@ -17,10 +17,10 @@
!
! notice the -1 sign: in this way three_e_4_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
! three_e_4_idx_direct_bi_ort (m,j,k,i) : Lk Ri Imm Ijj + Lj Rj Imm Iki + Lm Rm Ijj Iki
! three_e_4_idx_exch13_bi_ort (m,j,k,i) : Lk Rm Imi Ijj + Lj Rj Imi Ikm + Lm Ri Ijj Ikm
! three_e_4_idx_direct_bi_ort (m,j,k,i) : Lk Ri Imm Ijj + Lj Rj Imm Iki + Lm Rm Ijj Iki
! three_e_4_idx_exch13_bi_ort (m,j,k,i) : Lk Rm Imi Ijj + Lj Rj Imi Ikm + Lm Ri Ijj Ikm
! three_e_4_idx_exch23_bi_ort (m,j,k,i) : Lk Ri Imj Ijm + Lj Rm Imj Iki + Lm Rj Ijm Iki
! three_e_4_idx_cycle_1_bi_ort(m,j,k,i) : Lk Rm Imj Iji + Lj Ri Imj Ikm + Lm Rj Iji Ikm
! three_e_4_idx_cycle_1_bi_ort(m,j,k,i) : Lk Rm Imj Iji + Lj Ri Imj Ikm + Lm Rj Iji Ikm
!
END_DOC
@ -74,8 +74,8 @@
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (n, ipoint, tmp_loc_1, tmp_loc_2, tmp_2d, tmp1, tmp2) &
!$OMP SHARED (mo_num, n_points_final_grid, i, k, &
!$OMP PRIVATE (k, i, j, m, n, ipoint, tmp_loc_1, tmp_loc_2, tmp_2d, tmp1, tmp2) &
!$OMP SHARED (mo_num, n_points_final_grid, &
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
!$OMP tmp_aux_1, tmp_aux_2, &
@ -125,17 +125,17 @@
do n = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp_loc_1 = mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,n)
tmp_loc_2 = mos_l_in_r_array_transp(ipoint,n) * mos_r_in_r_array_transp(ipoint,i)
tmp1(ipoint,1,n) = int2_grad1_u12_bimo_t(ipoint,1,n,i) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,1,k,n) * tmp_loc_2
tmp1(ipoint,2,n) = int2_grad1_u12_bimo_t(ipoint,2,n,i) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,2,k,n) * tmp_loc_2
tmp1(ipoint,3,n) = int2_grad1_u12_bimo_t(ipoint,3,n,i) * tmp_loc_1 + int2_grad1_u12_bimo_t(ipoint,3,k,n) * tmp_loc_2
tmp1(ipoint,4,n) = int2_grad1_u12_bimo_t(ipoint,1,n,i) * int2_grad1_u12_bimo_t(ipoint,1,k,n) &
+ int2_grad1_u12_bimo_t(ipoint,2,n,i) * int2_grad1_u12_bimo_t(ipoint,2,k,n) &
+ int2_grad1_u12_bimo_t(ipoint,3,n,i) * int2_grad1_u12_bimo_t(ipoint,3,k,n)
tmp2(ipoint,1,n) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,n)
tmp2(ipoint,2,n) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,n)
tmp2(ipoint,3,n) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,n)
@ -225,7 +225,7 @@
print *, ' wall time for three_e_4_idx_bi_ort', wall1 - wall0
call print_memory_usage()
END_PROVIDER
END_PROVIDER
! ---

View File

@ -29,6 +29,9 @@ BEGIN_PROVIDER [ double precision, three_body_ints_bi_ort, (mo_num, mo_num, mo_n
!provide x_W_ki_bi_ortho_erf_rk
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
provide int2_grad1_u12_ao_transp final_grid_points int2_grad1_u12_bimo_t
provide mo_l_coef mo_r_coef mos_l_in_r_array_transp mos_r_in_r_array_transp n_points_final_grid
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &

View File

@ -9,8 +9,8 @@ function run_stoch() {
test_exe casscf || skip
qp set perturbation do_pt2 True
qp set determinants n_det_max $3
qp set davidson threshold_davidson 1.e-10
qp set davidson n_states_diag 4
qp set davidson_keywords threshold_davidson 1.e-10
qp set davidson_keywords n_states_diag 4
qp run casscf | tee casscf.out
energy1="$(ezfio get casscf energy_pt2 | tr '[]' ' ' | cut -d ',' -f 1)"
eq $energy1 $1 $thresh

View File

@ -73,3 +73,9 @@ type: logical
doc: If |true|, the pt2_max value in the CIPSI iterations will automatically adapt, otherwise it is fixed at the value given in the EZFIO folder
interface: ezfio,provider,ocaml
default: True
[small_active_space]
type: logical
doc: If |true|, the pt2_max value in the CIPSI is set to 10-10 and will not change
interface: ezfio,provider,ocaml
default: False

View File

@ -3,3 +3,45 @@ casscf
======
|CASSCF| program with the CIPSI algorithm.
Example of inputs
-----------------
a) Small active space : standard CASSCF
---------------------------------------
Let's do O2 (triplet) in aug-cc-pvdz with the following geometry (xyz format, Bohr units)
3
O 0.0000000000 0.0000000000 -1.1408000000
O 0.0000000000 0.0000000000 1.1408000000
# Create the ezfio folder
qp create_ezfio -b aug-cc-pvdz O2.xyz -m 3 -a -o O2_avdz
# Start with an ROHF guess
qp run scf | tee ${EZFIO_FILE}.rohf.out
# Get the ROHF energy for check
qp get hartree_fock energy # should be -149.4684509
# Define the full valence active space: the two 1s are doubly occupied, the other 8 valence orbitals are active
# CASSCF(12e,10orb)
qp set_mo_class -c "[1-2]" -a "[3-10]" -v "[11-46]"
# Specify that you want an near exact CASSCF, i.e. the CIPSI selection will stop at pt2_max = 10^-10
qp set casscf_cipsi small_active_space True
# RUN THE CASSCF
qp run casscf | tee ${EZFIO_FILE}.casscf.out
# you should find around -149.7243542
b) Large active space : Exploit the selected CI in the active space
-------------------------------------------------------------------
#Let us start from the small active space calculation orbitals and add another 10 virtuals: CASSCF(12e,20orb)
qp set_mo_class -c "[1-2]" -a "[3-20]" -v "[21-46]"
# As this active space is larger, you unset the small_active_space feature
qp set casscf_cipsi small_active_space False
# As it is a large active space, the energy convergence thereshold is set to be 0.0001
qp run casscf | tee ${EZFIO_FILE}.casscf_large.out
# you should find around -149.9046

View File

@ -8,17 +8,23 @@ program casscf
! touch no_vvvv_integrals
n_det_max_full = 500
touch n_det_max_full
pt2_relative_error = 0.04
if(small_active_space)then
pt2_relative_error = 0.00001
else
thresh_scf = 1.d-4
pt2_relative_error = 0.04
endif
touch pt2_relative_error
! call run_stochastic_cipsi
call run
end
subroutine run
implicit none
double precision :: energy_old, energy, pt2_max_before, ept2_before,delta_E
double precision :: energy_old, energy, pt2_max_before,delta_E
logical :: converged,state_following_casscf_cipsi_save
integer :: iteration
integer :: iteration,istate
double precision, allocatable :: E_PT2(:), PT2(:), Ev(:), ept2_before(:)
allocate(E_PT2(N_states), PT2(N_states), Ev(N_states), ept2_before(N_states))
converged = .False.
energy = 0.d0
@ -28,13 +34,20 @@ subroutine run
state_following_casscf = .True.
touch state_following_casscf
ept2_before = 0.d0
if(adaptive_pt2_max)then
pt2_max = 0.005
if(small_active_space)then
pt2_max = 1.d-10
SOFT_TOUCH pt2_max
else
if(adaptive_pt2_max)then
pt2_max = 0.005
SOFT_TOUCH pt2_max
endif
endif
do while (.not.converged)
print*,'pt2_max = ',pt2_max
call run_stochastic_cipsi
call run_stochastic_cipsi(Ev,PT2)
print*,'Ev,PT2',Ev(1),PT2(1)
E_PT2(1:N_states) = Ev(1:N_states) + PT2(1:N_states)
energy_old = energy
energy = eone+etwo+ecore
pt2_max_before = pt2_max
@ -42,15 +55,13 @@ subroutine run
call write_time(6)
call write_int(6,iteration,'CAS-SCF iteration = ')
call write_double(6,energy,'CAS-SCF energy = ')
if(n_states == 1)then
double precision :: E_PT2, PT2
call ezfio_get_casscf_cipsi_energy_pt2(E_PT2)
call ezfio_get_casscf_cipsi_energy(PT2)
PT2 -= E_PT2
call write_double(6,E_PT2,'E + PT2 energy = ')
call write_double(6,PT2,' PT2 = ')
! if(n_states == 1)then
! call ezfio_get_casscf_cipsi_energy_pt2(E_PT2)
! call ezfio_get_casscf_cipsi_energy(PT2)
call write_double(6,E_PT2(1:N_states),'E + PT2 energy = ')
call write_double(6,PT2(1:N_states),' PT2 = ')
call write_double(6,pt2_max,' PT2_MAX = ')
endif
! endif
print*,''
call write_double(6,norm_grad_vec2,'Norm of gradients = ')
@ -65,15 +76,20 @@ subroutine run
else if (criterion_casscf == "gradients")then
converged = norm_grad_vec2 < thresh_scf
else if (criterion_casscf == "e_pt2")then
delta_E = dabs(E_PT2 - ept2_before)
delta_E = 0.d0
do istate = 1, N_states
delta_E += dabs(E_PT2(istate) - ept2_before(istate))
enddo
converged = dabs(delta_E) < thresh_casscf
endif
ept2_before = E_PT2
if(adaptive_pt2_max)then
pt2_max = dabs(energy_improvement / (pt2_relative_error))
pt2_max = min(pt2_max, pt2_max_before)
if(n_act_orb.ge.n_big_act_orb)then
pt2_max = max(pt2_max,pt2_min_casscf)
if(.not.small_active_space)then
if(adaptive_pt2_max)then
pt2_max = dabs(energy_improvement / (pt2_relative_error))
pt2_max = min(pt2_max, pt2_max_before)
if(n_act_orb.ge.n_big_act_orb)then
pt2_max = max(pt2_max,pt2_min_casscf)
endif
endif
endif
print*,''
@ -94,8 +110,10 @@ subroutine run
read_wf = .True.
call clear_mo_map
SOFT_TOUCH mo_coef N_det psi_det psi_coef
if(adaptive_pt2_max)then
SOFT_TOUCH pt2_max
if(.not.small_active_space)then
if(adaptive_pt2_max)then
SOFT_TOUCH pt2_max
endif
endif
if(iteration .gt. 3)then
state_following_casscf = state_following_casscf_cipsi_save
@ -104,6 +122,25 @@ subroutine run
endif
enddo
integer :: i
print*,'Converged CASSCF '
print*,'--------------------------'
write(6,*) ' occupation numbers of orbitals '
do i=1,mo_num
write(6,*) i,occnum(i)
end do
print*,'--------------'
!
! write(6,*)
! write(6,*) ' the diagonal of the inactive effective Fock matrix '
! write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num)
! write(6,*)
print*,'Fock MCSCF'
do i = 1, mo_num
write(*,*)i,mcscf_fock_diag_mo(i)
! write(*,*)mcscf_fock_alpha_mo(i,i)
enddo
end

View File

@ -17,6 +17,35 @@ BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ]
END_PROVIDER
BEGIN_PROVIDER [double precision, D0tu_alpha_ao, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, D0tu_beta_ao, (ao_num, ao_num)]
implicit none
integer :: i,ii,j,u,t,uu,tt
double precision, allocatable :: D0_tmp_alpha(:,:),D0_tmp_beta(:,:)
allocate(D0_tmp_alpha(mo_num, mo_num),D0_tmp_beta(mo_num, mo_num))
D0_tmp_beta = 0.d0
D0_tmp_alpha = 0.d0
do i = 1, n_core_inact_orb
ii = list_core_inact(i)
D0_tmp_alpha(ii,ii) = 1.d0
D0_tmp_beta(ii,ii) = 1.d0
enddo
print*,'Diagonal elements of the 1RDM in the active space'
do u=1,n_act_orb
uu = list_act(u)
print*,uu,one_e_dm_mo_alpha_average(uu,uu),one_e_dm_mo_beta_average(uu,uu)
do t=1,n_act_orb
tt = list_act(t)
D0_tmp_alpha(tt,uu) = one_e_dm_mo_alpha_average(tt,uu)
D0_tmp_beta(tt,uu) = one_e_dm_mo_beta_average(tt,uu)
enddo
enddo
call mo_to_ao_no_overlap(D0_tmp_alpha,mo_num,D0tu_alpha_ao,ao_num)
call mo_to_ao_no_overlap(D0_tmp_beta,mo_num,D0tu_beta_ao,ao_num)
END_PROVIDER
BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ]
BEGIN_DOC
! The second-order density matrix in the basis of the starting MOs ONLY IN THE RANGE OF ACTIVE MOS

View File

@ -77,4 +77,119 @@ BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ]
END_PROVIDER
BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_ao, (ao_num, ao_num)]
&BEGIN_PROVIDER [ double precision, mcscf_fock_beta_ao, (ao_num, ao_num)]
implicit none
BEGIN_DOC
! mcscf_fock_alpha_ao are set to usual Fock like operator but computed with the MCSCF densities on the AO basis
END_DOC
SCF_density_matrix_ao_alpha = D0tu_alpha_ao
SCF_density_matrix_ao_beta = D0tu_beta_ao
soft_touch SCF_density_matrix_ao_alpha SCF_density_matrix_ao_beta
mcscf_fock_beta_ao = fock_matrix_ao_beta
mcscf_fock_alpha_ao = fock_matrix_ao_alpha
END_PROVIDER
BEGIN_PROVIDER [ double precision, mcscf_fock_alpha_mo, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, mcscf_fock_beta_mo, (mo_num, mo_num)]
implicit none
BEGIN_DOC
! Mo_mcscf_fock_alpha are set to usual Fock like operator but computed with the MCSCF densities on the MO basis
END_DOC
call ao_to_mo(mcscf_fock_alpha_ao,ao_num,mcscf_fock_alpha_mo,mo_num)
call ao_to_mo(mcscf_fock_beta_ao,ao_num,mcscf_fock_beta_mo,mo_num)
END_PROVIDER
BEGIN_PROVIDER [ double precision, mcscf_fock_mo, (mo_num,mo_num) ]
&BEGIN_PROVIDER [ double precision, mcscf_fock_diag_mo, (mo_num)]
implicit none
BEGIN_DOC
! MCSF Fock matrix on the MO basis.
! For open shells, the ROHF Fock Matrix is ::
!
! | Rcc | F^b | Fcv |
! |-----------------------|
! | F^b | Roo | F^a |
! |-----------------------|
! | Fcv | F^a | Rvv |
!
! C: Core, O: Open, V: Virtual
!
! Rcc = Acc Fcc^a + Bcc Fcc^b
! Roo = Aoo Foo^a + Boo Foo^b
! Rvv = Avv Fvv^a + Bvv Fvv^b
! Fcv = (F^a + F^b)/2
!
! F^a: Fock matrix alpha (MO), F^b: Fock matrix beta (MO)
! A,B: Coupling parameters
!
! J. Chem. Phys. 133, 141102 (2010), https://doi.org/10.1063/1.3503173
! Coupling parameters from J. Chem. Phys. 125, 204110 (2006); https://doi.org/10.1063/1.2393223.
! cc oo vv
! A -0.5 0.5 1.5
! B 1.5 0.5 -0.5
!
END_DOC
integer :: i,j,n
if (elec_alpha_num == elec_beta_num) then
mcscf_fock_mo = mcscf_fock_alpha_mo
else
! Core
do j = 1, elec_beta_num
! Core
do i = 1, elec_beta_num
mcscf_fock_mo(i,j) = - 0.5d0 * mcscf_fock_alpha_mo(i,j) &
+ 1.5d0 * mcscf_fock_beta_mo(i,j)
enddo
! Open
do i = elec_beta_num+1, elec_alpha_num
mcscf_fock_mo(i,j) = mcscf_fock_beta_mo(i,j)
enddo
! Virtual
do i = elec_alpha_num+1, mo_num
mcscf_fock_mo(i,j) = 0.5d0 * mcscf_fock_alpha_mo(i,j) &
+ 0.5d0 * mcscf_fock_beta_mo(i,j)
enddo
enddo
! Open
do j = elec_beta_num+1, elec_alpha_num
! Core
do i = 1, elec_beta_num
mcscf_fock_mo(i,j) = mcscf_fock_beta_mo(i,j)
enddo
! Open
do i = elec_beta_num+1, elec_alpha_num
mcscf_fock_mo(i,j) = 0.5d0 * mcscf_fock_alpha_mo(i,j) &
+ 0.5d0 * mcscf_fock_beta_mo(i,j)
enddo
! Virtual
do i = elec_alpha_num+1, mo_num
mcscf_fock_mo(i,j) = mcscf_fock_alpha_mo(i,j)
enddo
enddo
! Virtual
do j = elec_alpha_num+1, mo_num
! Core
do i = 1, elec_beta_num
mcscf_fock_mo(i,j) = 0.5d0 * mcscf_fock_alpha_mo(i,j) &
+ 0.5d0 * mcscf_fock_beta_mo(i,j)
enddo
! Open
do i = elec_beta_num+1, elec_alpha_num
mcscf_fock_mo(i,j) = mcscf_fock_alpha_mo(i,j)
enddo
! Virtual
do i = elec_alpha_num+1, mo_num
mcscf_fock_mo(i,j) = 1.5d0 * mcscf_fock_alpha_mo(i,j) &
- 0.5d0 * mcscf_fock_beta_mo(i,j)
enddo
enddo
endif
do i = 1, mo_num
mcscf_fock_diag_mo(i) = mcscf_fock_mo(i,i)
enddo
END_PROVIDER

View File

@ -1,10 +1,11 @@
subroutine run_stochastic_cipsi
subroutine run_stochastic_cipsi(Ev,PT2)
use selection_types
implicit none
BEGIN_DOC
! Selected Full Configuration Interaction with Stochastic selection and PT2.
END_DOC
integer :: i,j,k
double precision, intent(out) :: Ev(N_states), PT2(N_states)
double precision, allocatable :: zeros(:)
integer :: to_select
type(pt2_type) :: pt2_data, pt2_data_err
@ -79,12 +80,14 @@ subroutine run_stochastic_cipsi
to_select = max(N_states_diag, to_select)
Ev(1:N_states) = psi_energy_with_nucl_rep(1:N_states)
call pt2_dealloc(pt2_data)
call pt2_dealloc(pt2_data_err)
call pt2_alloc(pt2_data, N_states)
call pt2_alloc(pt2_data_err, N_states)
call ZMQ_pt2(psi_energy_with_nucl_rep,pt2_data,pt2_data_err,relative_error,to_select) ! Stochastic PT2 and selection
PT2(1:N_states) = pt2_data % pt2(1:N_states)
correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / &
(psi_energy_with_nucl_rep(1) + pt2_data % rpt2(1) - hf_energy_ref)
correlation_energy_ratio = min(1.d0,correlation_energy_ratio)

View File

@ -9,7 +9,7 @@ function run() {
qp set_file $1
qp edit --check
qp set determinants n_states 3
qp set davidson threshold_davidson 1.e-12
qp set davidson_keywords threshold_davidson 1.e-12
qp set mo_two_e_ints io_mo_two_e_integrals Write
qp set_frozen_core
qp run cis
@ -59,7 +59,7 @@ function run() {
@test "ClO" { # 1.65582s 2.06465s
[[ -n $TRAVIS ]] && skip
run clo.ezfio -534.263560525680 -534.256601571199 -534.062020844428
run clo.ezfio -534.2635737789097 -534.2566081298855 -534.0620070783308
}
@test "SO" { # 1.9667s 2.91234s
@ -69,7 +69,7 @@ function run() {
@test "OH" { # 2.201s 2.65573s
[[ -n $TRAVIS ]] && skip
run oh.ezfio -75.4314648243896 -75.4254639668256 -75.2707675632313
run oh.ezfio -75.4314822573358 -75.4254733392003 -75.2707586997333
}
@test "H2O2" { # 2.27079s 3.07875s
@ -109,7 +109,7 @@ function run() {
@test "DHNO" { # 6.42976s 12.9899s
[[ -n $TRAVIS ]] && skip
run dhno.ezfio -130.4472288472718 -130.3571808164850 -130.2196257046987
run dhno.ezfio -130.447238897118 -130.357186843611 -130.219626716369
}
@test "CH4" { # 6.4969s 10.9157s
@ -129,7 +129,7 @@ function run() {
@test "[Cu(NH3)4]2+" { # 29.7711s 3.45478m
[[ -n ${TRAVIS} ]] && skip
run cu_nh3_4_2plus.ezfio -1862.97958885180 -1862.92457657404 -1862.91134959451
run cu_nh3_4_2plus.ezfio -1862.97958844302 -1862.92454785007 -1862.91130869967
}

View File

@ -8,10 +8,9 @@ function run() {
test_exe cisd || skip
qp edit --check
qp set determinants n_states 2
qp set davidson threshold_davidson 1.e-12
qp set davidson n_states_diag 24
qp run cis
qp run cisd
qp set davidson_keywords threshold_davidson 1.e-12
qp set davidson_keywords n_states_diag 24
qp run cisd
energy1="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 1)"
energy2="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 2)"
eq $energy1 $1 $thresh
@ -19,7 +18,7 @@ function run() {
}
@test "B-B" { #
@test "B-B" { #
qp set_file b2_stretched.ezfio
qp set_frozen_core
run -49.120607088648597 -49.055152453388231
@ -34,7 +33,7 @@ function run() {
@test "HBO" { # 4.42968s 19.6099s
qp set_file hbo.ezfio
qp set_frozen_core
run -100.2019254455993 -99.79484127741013
run -100.2019254455993 -99.79484127741013
}
@test "HCO" { # 6.6077s 28.6801s
@ -46,7 +45,7 @@ function run() {
@test "H2O" { # 7.0651s 30.6642s
qp set_file h2o.ezfio
qp set_frozen_core
run -76.22975602077072 -75.80609108747208
run -76.22975602077072 -75.80609108747208
}
@ -78,7 +77,7 @@ function run() {
[[ -n $TRAVIS ]] && skip
qp set_file oh.ezfio
qp set_frozen_core
run -75.6087472926588 -75.5370393736601
run -75.6088105201621 -75.5370802925698
}
@test "CH4" { # 19.821s 1.38648m
@ -105,8 +104,9 @@ function run() {
@test "DHNO" { # 24.7077s 1.46487m
[[ -n $TRAVIS ]] && skip
qp set_file dhno.ezfio
qp set_mo_class --core="[1-7]" --act="[8-64]"
run -130.458814562403 -130.356308303681
qp set_mo_class --core="[1-7]" --act="[8-64]"
run -130.4659881027444 -130.2692384198501
# run -130.458814562403 -130.356308303681
}
@test "H3COH" { # 24.7248s 1.85043m
@ -120,7 +120,7 @@ function run() {
[[ -n $TRAVIS ]] && skip
qp set_file cu_nh3_4_2plus.ezfio
qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]"
run -1862.98689579931 -1862.6883044626563
run -1862.98310702274 -1862.88506319755
}
@ -135,14 +135,14 @@ function run() {
[[ -n $TRAVIS ]] && skip
qp set_file c2h2.ezfio
qp set_mo_class --act="[1-30]" --del="[31-36]"
run -12.3566731164213 -11.9495394759914
run -12.3566731164213 -11.9495394759914
}
@test "ClO" { # 37.6949s
[[ -n $TRAVIS ]] && skip
qp set_file clo.ezfio
qp set_frozen_core
run -534.5404021326773 -534.3818725793897
run -534.540464615019 -534.381904487587
}
@test "F2" { # 45.2078s
@ -155,7 +155,7 @@ function run() {
@test "SO2" { # 47.6922s
[[ -n $TRAVIS ]] && skip
qp set_file so2.ezfio
qp set_mo_class --core="[1-8]" --act="[9-87]"
qp set_mo_class --core="[1-8]" --act="[9-87]"
run -41.5746738710350 -41.3800467740750
}
@ -177,7 +177,7 @@ function run() {
[[ -n $TRAVIS ]] && skip
qp set_file n2.ezfio
qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-60]"
run -109.275693633982 -108.757794570948
run -109.275693633982 -108.757794570948
}
@test "HCN" { # 133.8696s

View File

@ -286,7 +286,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
! Small
h(N_st_diag*itermax,N_st_diag*itermax), &
h_p(N_st_diag*itermax,N_st_diag*itermax), &
! h_p(N_st_diag*itermax,N_st_diag*itermax), &
y(N_st_diag*itermax,N_st_diag*itermax), &
s_(N_st_diag*itermax,N_st_diag*itermax), &
s_tmp(N_st_diag*itermax,N_st_diag*itermax), &
@ -340,7 +340,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
exit
endif
do iter=1,itermax-1
iter = 0
do while (iter < itermax-1)
iter += 1
! do iter=1,itermax-1
shift = N_st_diag*(iter-1)
shift2 = N_st_diag*iter
@ -430,30 +433,30 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
call dgemm('T','N', shift2, shift2, sze, &
1.d0, U, size(U,1), W, size(W,1), &
0.d0, h, size(h_p,1))
0.d0, h, size(h,1))
call dgemm('T','N', shift2, shift2, sze, &
1.d0, U, size(U,1), U, size(U,1), &
0.d0, s_tmp, size(s_tmp,1))
! Penalty method
! --------------
if (s2_eig) then
h_p = s_
do k=1,shift2
h_p(k,k) = h_p(k,k) - expected_s2
enddo
if (only_expected_s2) then
alpha = 0.1d0
h_p = h + alpha*h_p
else
alpha = 0.0001d0
h_p = h + alpha*h_p
endif
else
h_p = h
alpha = 0.d0
endif
! ! Penalty method
! ! --------------
!
! if (s2_eig) then
! h_p = s_
! do k=1,shift2
! h_p(k,k) = h_p(k,k) - expected_s2
! enddo
! if (only_expected_s2) then
! alpha = 0.1d0
! h_p = h + alpha*h_p
! else
! alpha = 0.0001d0
! h_p = h + alpha*h_p
! endif
! else
! h_p = h
! alpha = 0.d0
! endif
! Diagonalize h_p
! ---------------
@ -473,8 +476,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
call dsygv(1,'V','U',shift2,y,size(y,1), &
s_tmp,size(s_tmp,1), lambda, work,lwork,info)
deallocate(work)
if (info /= 0) then
stop 'DSYGV Diagonalization failed'
if (info > 0) then
! Numerical errors propagate. We need to reduce the number of iterations
itermax = iter-1
exit
endif
! Compute Energy for each eigenvector

View File

@ -493,3 +493,101 @@ subroutine get_occupation_from_dets(istate,occupation)
enddo
end
BEGIN_PROVIDER [double precision, difference_dm, (mo_num, mo_num, N_states)]
implicit none
BEGIN_DOC
! difference_dm(i,j,istate) = dm(i,j,1) - dm(i,j,istate)
END_DOC
integer :: istate
do istate = 1, N_states
difference_dm(:,:,istate) = one_e_dm_mo_alpha(:,:,1) + one_e_dm_mo_beta(:,:,1) &
- (one_e_dm_mo_alpha(:,:,istate) + one_e_dm_mo_beta(:,:,istate))
enddo
END_PROVIDER
BEGIN_PROVIDER [double precision, difference_dm_eigvect, (mo_num, mo_num, N_states) ]
&BEGIN_PROVIDER [double precision, difference_dm_eigval, (mo_num, N_states) ]
implicit none
BEGIN_DOC
! eigenvalues and eigevenctors of the difference_dm
END_DOC
integer :: istate,i
do istate = 2, N_states
call lapack_diag(difference_dm_eigval(1,istate),difference_dm_eigvect(1,1,istate)&
,difference_dm(1,1,istate),mo_num,mo_num)
print*,'Eigenvalues of difference_dm for state ',istate
do i = 1, mo_num
print*,i,difference_dm_eigval(i,istate)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer , n_attachment, (N_states)]
&BEGIN_PROVIDER [ integer , n_dettachment, (N_states)]
&BEGIN_PROVIDER [ integer , list_attachment, (mo_num,N_states)]
&BEGIN_PROVIDER [ integer , list_dettachment, (mo_num,N_states)]
implicit none
integer :: i,istate
integer :: list_attachment_tmp(mo_num)
n_attachment = 0
n_dettachment = 0
do istate = 2, N_states
do i = 1, mo_num
if(difference_dm_eigval(i,istate).lt.0.d0)then ! dettachment_orbitals
n_dettachment(istate) += 1
list_dettachment(n_dettachment(istate),istate) = i ! they are already sorted
else
n_attachment(istate) += 1
list_attachment_tmp(n_attachment(istate)) = i ! they are not sorted
endif
enddo
! sorting the attachment
do i = 0, n_attachment(istate) - 1
list_attachment(i+1,istate) = list_attachment_tmp(n_attachment(istate) - i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, attachment_numbers_sorted, (mo_num, N_states)]
&BEGIN_PROVIDER [ double precision, dettachment_numbers_sorted, (mo_num, N_states)]
implicit none
integer :: i,istate
do istate = 2, N_states
print*,'dettachment'
do i = 1, n_dettachment(istate)
dettachment_numbers_sorted(i,istate) = difference_dm_eigval(list_dettachment(i,istate),istate)
print*,i,list_dettachment(i,istate),dettachment_numbers_sorted(i,istate)
enddo
print*,'attachment'
do i = 1, n_attachment(istate)
attachment_numbers_sorted(i,istate) = difference_dm_eigval(list_attachment(i,istate),istate)
print*,i,list_attachment(i,istate),attachment_numbers_sorted(i,istate)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, attachment_orbitals, (ao_num, mo_num, N_states)]
&BEGIN_PROVIDER [ double precision, dettachment_orbitals, (ao_num, mo_num, N_states)]
implicit none
integer :: i,j,k,istate
attachment_orbitals = 0.d0
dettachment_orbitals = 0.d0
do istate = 2, N_states
do i = 1, n_dettachment(istate)
do j = 1, mo_num
do k = 1, ao_num
dettachment_orbitals(k,list_dettachment(i,istate),istate) += mo_coef(k,j) * difference_dm_eigvect(j,list_dettachment(i,istate),istate)
enddo
enddo
enddo
do i = 1, n_attachment(istate)
do j = 1, mo_num
do k = 1, ao_num
attachment_orbitals(k,i,istate) += mo_coef(k,j) * difference_dm_eigvect(j,list_attachment(i,istate),istate)
enddo
enddo
enddo
enddo
END_PROVIDER

View File

@ -26,10 +26,10 @@
enddo
enddo
! print*,'electron part for z_dipole = ',z_dipole_moment
! print*,'electron part for y_dipole = ',y_dipole_moment
! print*,'electron part for x_dipole = ',x_dipole_moment
!
print*,'electron part for z_dipole = ',z_dipole_moment
print*,'electron part for y_dipole = ',y_dipole_moment
print*,'electron part for x_dipole = ',x_dipole_moment
nuclei_part_z = 0.d0
nuclei_part_y = 0.d0
nuclei_part_x = 0.d0
@ -38,10 +38,10 @@
nuclei_part_y += nucl_charge(i) * nucl_coord(i,2)
nuclei_part_x += nucl_charge(i) * nucl_coord(i,1)
enddo
! print*,'nuclei part for z_dipole = ',nuclei_part_z
! print*,'nuclei part for y_dipole = ',nuclei_part_y
! print*,'nuclei part for x_dipole = ',nuclei_part_x
!
print*,'nuclei part for z_dipole = ',nuclei_part_z
print*,'nuclei part for y_dipole = ',nuclei_part_y
print*,'nuclei part for x_dipole = ',nuclei_part_x
do istate = 1, N_states
z_dipole_moment(istate) += nuclei_part_z
y_dipole_moment(istate) += nuclei_part_y

View File

@ -4,6 +4,4 @@ mo_one_e_ints
mo_two_e_ints
ao_one_e_ints
ao_two_e_ints
mo_two_e_erf_ints
ao_two_e_erf_ints
mu_of_r

View File

@ -3,7 +3,7 @@ BEGIN_PROVIDER [double precision, mu_erf_dft]
BEGIN_DOC
! range separation parameter used in RS-DFT.
!
! It is set to mu_erf in order to be consistent with the module "ao_two_e_erf_ints"
! It is set to mu_erf in order to be consistent with the module "hamiltonian"
END_DOC
mu_erf_dft = mu_erf

View File

@ -146,3 +146,51 @@ end
end subroutine g0_dg0
subroutine g0_dg0_d2g0(rho, rho_a, rho_b, g0, dg0drho, d2g0drho2)
implicit none
BEGIN_DOC
! Give the on-top pair distribution function g0 second derivative according to rho d2g0drho2
END_DOC
double precision, intent (in) :: rho, rho_a, rho_b
double precision, intent (out) :: g0, dg0drho, d2g0drho2
double precision :: pi
double precision :: g0_UEG_mu_inf, dg0drs, d2g0drs2, d2rsdrho2
double precision :: C1, F1, D1, E1, B1, rs
pi = dacos(-1.d0)
C1 = 0.0819306d0
F1 = 0.752411d0
D1 = -0.0127713d0
E1 = 0.00185898d0
B1 = 0.7317d0 - F1
if(dabs(rho).gt.1.d-20)then
rs = (3.d0 / (4.d0*pi*rho))**(1.d0/3.d0)
else
rs = (3.d0 / (4.d0*pi*1.d-20))**(1.d0/3.d0)
endif
g0 = g0_UEG_mu_inf(rho_a, rho_b)
if(dabs(F1*rs).lt.50.d0)then
dg0drs = 0.5d0*((-B1 + 2.d0*C1*rs + 3.d0*D1*rs**2 + 4.d0*E1*rs**3)-F1*(1.d0 - B1*rs + C1*rs**2 + D1*rs**3 + E1*rs**4))*dexp(-F1*rs)
d2g0drs2 = 0.5d0*((2.d0*C1 + 6.d0*D1*rs + 12*E1*rs**2) - 2.d0*F1*(-B1 + 2.d0*C1*rs + 3.d0*D1*rs**2 + 4.d0*E1*rs**3)&
&+ (F1**2)*(1.d0 - B1*rs + C1*rs**2 + D1*rs**3 + E1*rs**4))*dexp(-F1*rs)
else
dg0drs = 0.d0
d2g0drs2 = 0.d0
endif
if(dabs(rho).gt.1.d-20)then
dg0drho = -((6.d0*dsqrt(pi)*rho**2)**(-2.d0/3.d0))*dg0drs
d2rsdrho2 = -8.d0*dsqrt(pi)*rho*(6.d0*dsqrt(pi)*rho**2)**(-5.d0/3.d0)
d2g0drho2 = dg0drho*d2rsdrho2 -((6.d0*dsqrt(pi)*rho**2)**(-4.d0/3.d0))*d2g0drs2
else
dg0drho = -((6.d0*dsqrt(pi)*1.d-40)**(-2.d0/3.d0))*dg0drs
d2rsdrho2 = -8.d0*dsqrt(pi)*(1.d-20)*(6.d0*dsqrt(pi)*1.d-40)**(-5.d0/3.d0)
d2g0drho2 = dg0drho*d2rsdrho2 - ((6.d0*dsqrt(pi)*1.d-40)**(-4.d0/3.d0))*d2g0drs2
endif
end subroutine g0_dg0

View File

@ -1,6 +1,5 @@
ao_basis
ao_one_e_ints
ao_two_e_erf_ints
ao_two_e_ints
aux_quantities
becke_numerical_grid
@ -24,13 +23,13 @@ functionals
generators_cas
generators_full
hartree_fock
hamiltonian
iterations
kohn_sham
kohn_sham_rs
mo_basis
mo_guess
mo_one_e_ints
mo_two_e_erf_ints
mo_two_e_ints
mpi
nuclei

View File

@ -41,8 +41,10 @@ program fci
write(json_unit,json_array_open_fmt) 'fci'
double precision, allocatable :: Ev(:),PT2(:)
allocate(Ev(N_states), PT2(N_states))
if (do_pt2) then
call run_stochastic_cipsi
call run_stochastic_cipsi(Ev,PT2)
else
call run_cipsi
endif

View File

@ -0,0 +1,8 @@
[mu_erf]
type: double precision
doc: cutting of the interaction in the range separated model
interface: ezfio,provider,ocaml
default: 0.5
ezfio_name: mu_erf

0
src/hamiltonian/NEED Normal file
View File

View File

@ -0,0 +1,5 @@
===========
hamiltonian
===========
Parameters of the Hamiltonian.

View File

@ -37,7 +37,7 @@ subroutine print_extrapolated_energy
write(*,*) 'minimum PT2 ', 'Extrapolated energy ', ' Excitation (a.u) ', ' Excitation (eV) '
write(*,*) '=========== ', '=================== ', '=================== ', '==================='
do k=2,N_iter_p
write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,k), extrapolated_energy(k,i), &
write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter_p+1-k), extrapolated_energy(k,i), &
extrapolated_energy(k,i) - extrapolated_energy(k,1), &
(extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
enddo

69
src/jastrow/EZFIO.cfg Normal file
View File

@ -0,0 +1,69 @@
[jast_type]
doc: Type of Jastrow [None| Mu | Qmckl]
type: character*(32)
interface: ezfio, provider, ocaml
default: None
[jast_qmckl_type_nucl_num]
doc: Number of different nuclei types in QMCkl jastrow
type: integer
interface: ezfio, provider
[jast_qmckl_type_nucl_vector]
doc: Nucleus type in QMCkl jastrow
type: integer
size: (nuclei.nucl_num)
interface: ezfio, provider
[jast_qmckl_rescale_ee]
doc: Rescaling factor for electron-electron in QMCkl Jastrow
type: double precision
interface: ezfio, provider
[jast_qmckl_rescale_en]
doc: Rescaling factor for electron-nucleus in QMCkl Jastrow
type: double precision
size: (jastrow.jast_qmckl_type_nucl_num)
interface: ezfio, provider
[jast_qmckl_aord_num]
doc: Order of polynomials in e-n parameters of QMCkl jastrow
type: integer
interface: ezfio, provider
[jast_qmckl_bord_num]
doc: Order of polynomials in e-e parameters of QMCkl jastrow
type: integer
interface: ezfio, provider
[jast_qmckl_cord_num]
doc: Order of polynomials in e-e-n parameters of QMCkl jastrow
type: integer
interface: ezfio, provider
[jast_qmckl_c_vector_size]
doc: Number of parameters for c_vector
type: integer
interface: ezfio, provider
[jast_qmckl_a_vector]
doc: electron-nucleus parameters in QMCkl Jastrow
type: double precision
size: (jastrow.jast_qmckl_type_nucl_num*jastrow.jast_qmckl_aord_num+jastrow.jast_qmckl_type_nucl_num)
interface: ezfio, provider
[jast_qmckl_b_vector]
doc: electron-electron parameters in QMCkl Jastrow
type: double precision
size: (jastrow.jast_qmckl_bord_num+1)
interface: ezfio, provider
[jast_qmckl_c_vector]
doc: electron-electron-nucleus parameters in QMCkl Jastrow
type: double precision
size: (jastrow.jast_qmckl_c_vector_size)
interface: ezfio, provider

2
src/jastrow/NEED Normal file
View File

@ -0,0 +1,2 @@
nuclei
electrons

3
src/jastrow/README.md Normal file
View File

@ -0,0 +1,3 @@
# Jastrow
Information relative to the Jastrow factor in trans-correlated calculations.

View File

@ -13,7 +13,7 @@ function run() {
qp set scf_utils thresh_scf 1.e-10
qp set dft_keywords exchange_functional $functional
qp set dft_keywords correlation_functional $functional
qp set ao_two_e_erf_ints mu_erf 0.5
qp set hamiltonian mu_erf 0.5
qp set becke_numerical_grid grid_type_sgn 1
qp_reset --mos $1
qp run rs_ks_scf

View File

@ -11,11 +11,13 @@ subroutine run_optimization
implicit none
double precision :: e_cipsi, e_opt, delta_e
double precision, allocatable :: Ev(:),PT2(:)
integer :: nb_iter,i
logical :: not_converged
character (len=100) :: filename
PROVIDE psi_det psi_coef mo_two_e_integrals_in_map ao_pseudo_integrals
allocate(Ev(N_states),PT2(N_states))
not_converged = .True.
nb_iter = 0
@ -38,7 +40,7 @@ subroutine run_optimization
print*,''
print*,'********** cipsi step **********'
! cispi calculation
call run_stochastic_cipsi
call run_stochastic_cipsi(Ev,PT2)
! State average energy after the cipsi step
call state_average_energy(e_cipsi)

View File

@ -1,6 +0,0 @@
[io_mo_two_e_integrals_erf]
type: Disk_access
doc: Read/Write MO integrals with the long range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None

View File

@ -1,3 +0,0 @@
ao_two_e_erf_ints
mo_two_e_ints
mo_basis

View File

@ -1,20 +0,0 @@
======================
mo_two_e_erf_ints
======================
Here, all two-electron integrals (:math:`erf({\mu}_{erf} * r_{12})/r_{12}`) are computed.
As they have 4 indices and many are zero, they are stored in a map, as defined
in :file:`Utils/map_module.f90`.
The range separation parameter :math:`{\mu}_{erf}` is the variable :option:`ao_two_e_erf_ints mu_erf`.
To fetch an |MO| integral, use
`get_mo_two_e_integral_erf(i,j,k,l,mo_integrals_map_erf)`
The conventions are:
* For |MO| integrals : <ij|kl> = <12|12>
Be aware that it might not be the same conventions for |MO| and |AO| integrals.

View File

@ -17,3 +17,10 @@ doc: If `True`, computes all integrals except for the integrals having 3 or 4 vi
interface: ezfio,provider,ocaml
default: false
[io_mo_two_e_integrals_erf]
type: Disk_access
doc: Read/Write MO integrals with the long range interaction from/to disk [ Write | Read | None ]
interface: ezfio,provider,ocaml
default: None

View File

@ -55,7 +55,7 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_erf_in_map ]
if (write_mo_two_e_integrals_erf) then
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
call ezfio_set_mo_two_e_erf_ints_io_mo_two_e_integrals_erf("Read")
call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals_erf("Read")
endif
END_PROVIDER

View File

@ -4,7 +4,7 @@ subroutine save_erf_two_e_integrals_mo
PROVIDE mo_two_e_integrals_erf_in_map
call ezfio_set_work_empty(.False.)
call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
call ezfio_set_mo_two_e_erf_ints_io_mo_two_e_integrals_erf('Read')
call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals_erf('Read')
end

View File

@ -1,3 +1,4 @@
qmckl
jastrow
ao_tc_eff_map
bi_ortho_mos

View File

@ -51,7 +51,7 @@
r1(2) = final_grid_points(2,ipoint)
r1(3) = final_grid_points(3,ipoint)
do jpoint = 1, n_points_extra_final_grid ! r2
do jpoint = 1, n_points_extra_final_grid ! r2
r2(1) = final_grid_points_extra(1,jpoint)
r2(2) = final_grid_points_extra(2,jpoint)
@ -63,9 +63,9 @@
dy = grad1_u2b(2)
dz = grad1_u2b(3)
grad1_u12_num(jpoint,ipoint,1) = dx
grad1_u12_num(jpoint,ipoint,2) = dy
grad1_u12_num(jpoint,ipoint,3) = dz
grad1_u12_num(jpoint,ipoint,1) = dx
grad1_u12_num(jpoint,ipoint,2) = dy
grad1_u12_num(jpoint,ipoint,3) = dz
grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
enddo
@ -92,7 +92,7 @@
v1b_r1 = j1b_nucl(r1)
call grad1_j1b_nucl(r1, grad1_v1b)
do jpoint = 1, n_points_extra_final_grid ! r2
do jpoint = 1, n_points_extra_final_grid ! r2
r2(1) = final_grid_points_extra(1,jpoint)
r2(2) = final_grid_points_extra(2,jpoint)
@ -106,9 +106,9 @@
dy = (grad1_u2b(2) * v1b_r1 + u2b_r12 * grad1_v1b(2)) * v1b_r2
dz = (grad1_u2b(3) * v1b_r1 + u2b_r12 * grad1_v1b(3)) * v1b_r2
grad1_u12_num(jpoint,ipoint,1) = dx
grad1_u12_num(jpoint,ipoint,2) = dy
grad1_u12_num(jpoint,ipoint,3) = dz
grad1_u12_num(jpoint,ipoint,1) = dx
grad1_u12_num(jpoint,ipoint,2) = dy
grad1_u12_num(jpoint,ipoint,3) = dz
grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
enddo
@ -121,86 +121,109 @@
double precision :: f
f = 1.d0 / dble(elec_num - 1)
integer*8 :: n_points, n_points_max, k
integer :: ipoint_block, ipoint_end
n_points_max = n_points_extra_final_grid * n_points_final_grid
n_points = 100_8*n_points_extra_final_grid
double precision, allocatable :: rij(:,:,:)
allocate( rij(3, 2, n_points_extra_final_grid) )
allocate( rij(3, 2, n_points) )
use qmckl
integer(qmckl_exit_code) :: rc
integer*8 :: npoints
npoints = n_points_extra_final_grid
double precision, allocatable :: gl(:,:,:)
allocate( gl(2,4,n_points_extra_final_grid) )
do ipoint = 1, n_points_final_grid ! r1
allocate( gl(2,4,n_points) )
do jpoint = 1, n_points_extra_final_grid ! r2
rij(1:3, 1, jpoint) = final_grid_points (1:3, ipoint)
rij(1:3, 2, jpoint) = final_grid_points_extra(1:3, jpoint)
do ipoint_block = 1, n_points_final_grid, 100 ! r1
ipoint_end = min(n_points_final_grid, ipoint_block+100)
k=0
do ipoint = ipoint_block, ipoint_end
do jpoint = 1, n_points_extra_final_grid ! r2
k=k+1
rij(1:3, 1, k) = final_grid_points (1:3, ipoint)
rij(1:3, 2, k) = final_grid_points_extra(1:3, jpoint)
end do
enddo
rc = qmckl_set_electron_coord(qmckl_ctx_jastrow, 'N', npoints, rij, npoints*6_8)
rc = qmckl_set_electron_coord(qmckl_ctx_jastrow, 'N', n_points, rij, n_points*6_8)
if (rc /= QMCKL_SUCCESS) then
print *, irp_here, 'qmckl error in set_electron_coord'
rc = qmckl_check(qmckl_ctx_jastrow, rc)
stop -1
endif
! ---
! e-e term
rc = qmckl_get_jastrow_champ_factor_ee_gl(qmckl_ctx_jastrow, gl, 8_8*npoints)
rc = qmckl_get_jastrow_champ_factor_ee_gl(qmckl_ctx_jastrow, gl, 8_8*n_points)
if (rc /= QMCKL_SUCCESS) then
print *, irp_here, 'qmckl error in fact_ee_gl'
print *, irp_here, ' qmckl error in fact_ee_gl'
rc = qmckl_check(qmckl_ctx_jastrow, rc)
stop -1
endif
do jpoint = 1, n_points_extra_final_grid ! r2
grad1_u12_num(jpoint,ipoint,1) = gl(1,1,jpoint)
grad1_u12_num(jpoint,ipoint,2) = gl(1,2,jpoint)
grad1_u12_num(jpoint,ipoint,3) = gl(1,3,jpoint)
k=0
do ipoint = ipoint_block, ipoint_end
do jpoint = 1, n_points_extra_final_grid ! r2
k=k+1
grad1_u12_num(jpoint,ipoint,1) = gl(1,1,k)
grad1_u12_num(jpoint,ipoint,2) = gl(1,2,k)
grad1_u12_num(jpoint,ipoint,3) = gl(1,3,k)
enddo
enddo
! ---
! e-e-n term
! ---
! e-e-n term
! rc = qmckl_get_jastrow_champ_factor_een_gl(qmckl_ctx_jastrow, gl, 8_8*npoints)
! if (rc /= QMCKL_SUCCESS) then
! print *, irp_here, 'qmckl error in fact_een_gl'
! stop -1
! endif
! rc = qmckl_get_jastrow_champ_factor_een_gl(qmckl_ctx_jastrow, gl, 8_8*n_points)
! if (rc /= QMCKL_SUCCESS) then
! print *, irp_here, 'qmckl error in fact_een_gl'
! rc = qmckl_check(qmckl_ctx_jastrow, rc)
! stop -1
! endif
!
! do jpoint = 1, n_points_extra_final_grid ! r2
! grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + gl(1,1,jpoint)
! grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + gl(1,2,jpoint)
! grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + gl(1,3,jpoint)
! enddo
! k=0
! do ipoint = 1, n_points_final_grid ! r1
! do jpoint = 1, n_points_extra_final_grid ! r2
! k=k+1
! grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + gl(1,1,k)
! grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + gl(1,2,k)
! grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + gl(1,3,k)
! enddo
! enddo
! ---
! e-n term
rc = qmckl_get_jastrow_champ_factor_en_gl(qmckl_ctx_jastrow, gl, 8_8*npoints)
rc = qmckl_get_jastrow_champ_factor_en_gl(qmckl_ctx_jastrow, gl, 8_8*n_points)
if (rc /= QMCKL_SUCCESS) then
print *, irp_here, 'qmckl error in fact_en_gl'
rc = qmckl_check(qmckl_ctx_jastrow, rc)
stop -1
endif
do jpoint = 1, n_points_extra_final_grid ! r2
grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + f * gl(1,1,jpoint)
grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + f * gl(1,2,jpoint)
grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + f * gl(1,3,jpoint)
k=0
do ipoint = ipoint_block, ipoint_end ! r1
do jpoint = 1, n_points_extra_final_grid ! r2
k = k+1
grad1_u12_num(jpoint,ipoint,1) = grad1_u12_num(jpoint,ipoint,1) + f * gl(1,1,k)
grad1_u12_num(jpoint,ipoint,2) = grad1_u12_num(jpoint,ipoint,2) + f * gl(1,2,k)
grad1_u12_num(jpoint,ipoint,3) = grad1_u12_num(jpoint,ipoint,3) + f * gl(1,3,k)
dx = grad1_u12_num(jpoint,ipoint,1)
dy = grad1_u12_num(jpoint,ipoint,2)
dz = grad1_u12_num(jpoint,ipoint,3)
grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
enddo
enddo
do jpoint = 1, n_points_extra_final_grid ! r2
dx = grad1_u12_num(jpoint,ipoint,1)
dy = grad1_u12_num(jpoint,ipoint,2)
dz = grad1_u12_num(jpoint,ipoint,3)
grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
enddo
enddo !ipoint_block
enddo
deallocate(gl, rij)
@ -212,7 +235,7 @@
endif
call wall_time(time1)
print*, ' Wall time for grad1_u12_num & grad1_u12_squared_num (min) =', (time1-time0)/60.d0
print*, ' Wall time for grad1_u12_num & grad1_u12_squared_num (min) =', (time1-time0)/60.d0
END_PROVIDER

View File

@ -1,5 +1,6 @@
BEGIN_PROVIDER [ integer*8, qmckl_ctx_jastrow ]
use qmckl
use iso_c_binding
implicit none
BEGIN_DOC
! Context for the QMCKL library
@ -8,6 +9,10 @@ BEGIN_PROVIDER [ integer*8, qmckl_ctx_jastrow ]
qmckl_ctx_jastrow = qmckl_context_create()
rc = qmckl_set_jastrow_champ_spin_independent(qmckl_ctx_jastrow, 1)
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
rc = qmckl_set_nucleus_num(qmckl_ctx_jastrow, nucl_num*1_8)
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
@ -26,77 +31,47 @@ BEGIN_PROVIDER [ integer*8, qmckl_ctx_jastrow ]
! Jastrow parameters
rc = qmckl_set_jastrow_champ_type_nucl_num (qmckl_ctx_jastrow, 2_8)
rc = qmckl_set_jastrow_champ_type_nucl_num(qmckl_ctx_jastrow, 1_8*jast_qmckl_type_nucl_num)
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
rc = qmckl_set_jastrow_champ_type_nucl_vector (qmckl_ctx_jastrow, (/0_8,1_8,1_8/), 1_8*nucl_num)
rc = qmckl_set_jastrow_champ_type_nucl_vector(qmckl_ctx_jastrow, 1_8*jast_qmckl_type_nucl_vector-1_8, 1_8*nucl_num)
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
rc = qmckl_set_jastrow_champ_rescale_factor_ee (qmckl_ctx_jastrow, 0.6d0)
rc = qmckl_set_jastrow_champ_rescale_factor_ee(qmckl_ctx_jastrow, jast_qmckl_rescale_ee)
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
rc = qmckl_set_jastrow_champ_rescale_factor_en (qmckl_ctx_jastrow, (/0.6d0, 0.6d0 /), 2_8 )
rc = qmckl_set_jastrow_champ_rescale_factor_en(qmckl_ctx_jastrow, jast_qmckl_rescale_en, 1_8*jast_qmckl_type_nucl_num)
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
rc = qmckl_set_jastrow_champ_aord_num (qmckl_ctx_jastrow, 5_8)
rc = qmckl_set_jastrow_champ_aord_num(qmckl_ctx_jastrow, jast_qmckl_aord_num*1_8)
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
rc = qmckl_set_jastrow_champ_bord_num (qmckl_ctx_jastrow, 5_8)
rc = qmckl_set_jastrow_champ_a_vector(qmckl_ctx_jastrow, jast_qmckl_a_vector, 1_8*size(jast_qmckl_a_vector))
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
rc = qmckl_set_jastrow_champ_cord_num (qmckl_ctx_jastrow, 0_8)
rc = qmckl_set_jastrow_champ_bord_num(qmckl_ctx_jastrow, jast_qmckl_bord_num*1_8)
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
! double precision :: a_vector(12) = dble(&
! (/ 0.00000000, 0.00000000, -0.71168405, -0.44415699, -0.13865109, 0.07002267 , &
! 0.00000000, 0.00000000, -0.11379992, 0.04542846, 0.01696997, -0.01809299 /) )
! double precision :: b_vector(6) = dble(&
! (/ 0.00000000, 0.65603311, 0.14581988, 0.03138163, 0.00153156, -0.00447302 /) )
! double precision :: c_vector(46) = &
! (/ 1.06384279d0, -1.44303973d0, -0.92409833d0, 0.11845356d0, -0.02980776d0, &
! 1.07048863d0, 0.06009623d0, -0.01854872d0, -0.00915398d0, 0.01324198d0, &
! -0.00504959d0, -0.01202497d0, -0.00531644d0, 0.15101629d0, -0.00723831d0, &
! -0.00384182d0, -0.00295036d0, -0.00114583d0, 0.00158107d0, -0.00078107d0, &
! -0.00080000d0, -0.14140576d0, -0.00237271d0, -0.03006706d0, 0.01537009d0, &
! -0.02327226d0, 0.16502789d0, -0.01458259d0, -0.09946065d0, 0.00850029d0, &
! -0.02969361d0, -0.01159547d0, 0.00516313d0, 0.00405247d0, -0.02200886d0, &
! 0.03376709d0, 0.01277767d0, -0.01523013d0, -0.00739224d0, -0.00463953d0, &
! 0.00003174d0, -0.01421128d0, 0.00808140d0, 0.00612988d0, -0.00610632d0, &
! 0.01926215d0 /)
! a_vector = 0.d0
! b_vector = 0.d0
! c_vector = 0.d0
double precision :: a_vector(12) = dble(&
(/ 0.00000000 , 0.00000000, -0.45105821, -0.23519218, -0.03825391, 0.10072866, &
0.00000000 , 0.00000000, -0.06930592, -0.02909224, -0.00134650, 0.01477242 /) )
double precision :: b_vector(6) = dble(&
(/ 0.00000000, 0.00000000, 0.29217862, -0.00450671, -0.02925982, -0.01381532 /) )
double precision :: c_vector(46)
c_vector = 0.d0
rc = qmckl_set_jastrow_champ_a_vector(qmckl_ctx_jastrow, a_vector, 12_8)
rc = qmckl_set_jastrow_champ_b_vector(qmckl_ctx_jastrow, jast_qmckl_b_vector, 1_8*size(jast_qmckl_b_vector))
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
rc = qmckl_set_jastrow_champ_b_vector(qmckl_ctx_jastrow, b_vector, 6_8)
rc = qmckl_set_jastrow_champ_cord_num(qmckl_ctx_jastrow, jast_qmckl_cord_num*1_8)
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
! rc = qmckl_set_jastrow_champ_c_vector(qmckl_ctx_jastrow, c_vector, 46_8)
! rc = qmckl_check(qmckl_ctx_jastrow, rc)
! if (rc /= QMCKL_SUCCESS) stop -1
if (jast_qmckl_cord_num > 0) then
rc = qmckl_set_jastrow_champ_c_vector(qmckl_ctx_jastrow, jast_qmckl_c_vector, 1_8*jast_qmckl_c_vector_size)
rc = qmckl_check(qmckl_ctx_jastrow, rc)
if (rc /= QMCKL_SUCCESS) stop -1
endif
END_PROVIDER

View File

@ -270,7 +270,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
integer, intent(out) :: n_real_eigv
double precision, intent(out) :: reigvec(n,n), leigvec(n,n), eigval(n)
integer :: i, j
integer :: i, j,k
integer :: n_good
double precision :: thr, thr_cut, thr_diag, thr_norm
double precision :: accu_d, accu_nd
@ -278,6 +278,8 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
integer, allocatable :: list_good(:), iorder(:)
double precision, allocatable :: WR(:), WI(:), VL(:,:), VR(:,:)
double precision, allocatable :: S(:,:)
double precision, allocatable :: phi_1_tilde(:),phi_2_tilde(:),chi_1_tilde(:),chi_2_tilde(:)
allocate(phi_1_tilde(n),phi_2_tilde(n),chi_1_tilde(n),chi_2_tilde(n))
! -------------------------------------------------------------------------------------
@ -301,11 +303,78 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
call lapack_diag_non_sym(n, A, WR, WI, VL, VR)
!call lapack_diag_non_sym_new(n, A, WR, WI, VL, VR)
!print *, ' '
!print *, ' eigenvalues'
!do i = 1, n
! write(*, '(1000(F16.10,X))') WR(i), WI(i)
!enddo
print *, ' '
print *, ' eigenvalues'
i = 1
do while(i .le. n)
write(*, '(I3,X,1000(F16.10,X))')i, WR(i), WI(i)
if(.false.)then
if(WI(i).ne.0.d0)then
print*,'*****************'
print*,'WARNING ! IMAGINARY EIGENVALUES !!!'
write(*, '(1000(F16.10,X))') WR(i), WI(i+1)
! phi = VR(:,i), psi = VR(:,i+1), |Phi_i> = phi + j psi , |Phi_i+1> = phi - j psi
! chi = VL(:,i), xhi = VL(:,i+1), |Chi_i> = chi + j xhi , |Chi_i+1> = chi - j xhi
!
accu_chi_phi = 0.d0
accu_xhi_psi = 0.d0
accu_chi_psi = 0.d0
accu_xhi_phi = 0.d0
double precision :: accu_chi_phi, accu_xhi_psi, accu_chi_psi, accu_xhi_phi
double precision :: mat_ovlp(2,2),eigval_tmp(2),eigvec(2,2),mat_ovlp_orig(2,2)
do j = 1, n
accu_chi_phi += VL(j,i) * VR(j,i)
accu_xhi_psi += VL(j,i+1) * VR(j,i+1)
accu_chi_psi += VL(j,i) * VR(j,i+1)
accu_xhi_phi += VL(j,i+1) * VR(j,i)
enddo
mat_ovlp_orig(1,1) = accu_chi_phi
mat_ovlp_orig(2,1) = accu_xhi_phi
mat_ovlp_orig(1,2) = accu_chi_psi
mat_ovlp_orig(2,2) = accu_xhi_psi
print*,'old overlap matrix '
write(*,'(100(F16.10,X))')mat_ovlp_orig(1:2,1)
write(*,'(100(F16.10,X))')mat_ovlp_orig(1:2,2)
mat_ovlp(1,1) = accu_xhi_phi
mat_ovlp(2,1) = accu_chi_phi
mat_ovlp(1,2) = accu_xhi_psi
mat_ovlp(2,2) = accu_chi_psi
!print*,'accu_chi_phi = ',accu_chi_phi
!print*,'accu_xhi_psi = ',accu_xhi_psi
!print*,'accu_chi_psi = ',accu_chi_psi
!print*,'accu_xhi_phi = ',accu_xhi_phi
print*,'new overlap matrix '
write(*,'(100(F16.10,X))')mat_ovlp(1:2,1)
write(*,'(100(F16.10,X))')mat_ovlp(1:2,2)
call lapack_diag(eigval_tmp,eigvec,mat_ovlp,2,2)
print*,'eigval_tmp(1) = ',eigval_tmp(1)
print*,'eigvec(1) = ',eigvec(1:2,1)
print*,'eigval_tmp(2) = ',eigval_tmp(2)
print*,'eigvec(2) = ',eigvec(1:2,2)
print*,'*****************'
phi_1_tilde = 0.d0
phi_2_tilde = 0.d0
chi_1_tilde = 0.d0
chi_2_tilde = 0.d0
do j = 1, n
phi_1_tilde(j) += VR(j,i) * eigvec(1,1) + VR(j,i+1) * eigvec(2,1)
phi_2_tilde(j) += VR(j,i) * eigvec(1,2) + VR(j,i+1) * eigvec(2,2)
chi_1_tilde(j) += VL(j,i+1) * eigvec(1,1) + VL(j,i) * eigvec(2,1)
chi_2_tilde(j) += VL(j,i+1) * eigvec(1,2) + VL(j,i) * eigvec(2,2)
enddo
VR(1:n,i) = phi_1_tilde(1:n)
VR(1:n,i+1) = phi_2_tilde(1:n)
! Vl(1:n,i) = -chi_1_tilde(1:n)
! Vl(1:n,i+1) = chi_2_tilde(1:n)
i+=1
endif
endif
i+=1
enddo
!print *, ' right eigenvect bef'
!do i = 1, n
! write(*, '(1000(F16.10,X))') VR(:,i)
@ -331,7 +400,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
!thr = 100d0
thr = Im_thresh_tcscf
do i = 1, n
!print*, 'Re(i) + Im(i)', WR(i), WI(i)
print*, 'Re(i) + Im(i)', WR(i), WI(i)
if(dabs(WI(i)) .lt. thr) then
n_good += 1
else
@ -405,7 +474,7 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(n_real_eigv))/dble(n_real_eigv) .lt. thr_d) ) then
!print *, ' lapack vectors are normalized and bi-orthogonalized'
print *, ' lapack vectors are normalized and bi-orthogonalized'
deallocate(S)
return
@ -422,13 +491,14 @@ subroutine non_hrmt_bieig(n, A, thr_d, thr_nd, leigvec, reigvec, n_real_eigv, ei
else
!print *, ' lapack vectors are not normalized neither bi-orthogonalized'
print *, ' lapack vectors are not normalized neither bi-orthogonalized'
! ---
! call impose_orthog_degen_eigvec(n, eigval, reigvec)
! call impose_orthog_degen_eigvec(n, eigval, leigvec)
call reorder_degen_eigvec(n, eigval, leigvec, reigvec)
call impose_biorthog_degen_eigvec(n, eigval, leigvec, reigvec)

View File

@ -1857,7 +1857,7 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_
integer :: i, j
double precision, allocatable :: SS(:,:)
!print *, ' check bi-orthogonality'
print *, ' check bi-orthogonality'
! ---
@ -1865,10 +1865,10 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_
, Vl, size(Vl, 1), Vr, size(Vr, 1) &
, 0.d0, S, size(S, 1) )
!print *, ' overlap matrix:'
!do i = 1, m
! write(*,'(1000(F16.10,X))') S(i,:)
!enddo
print *, ' overlap matrix:'
do i = 1, m
write(*,'(1000(F16.10,X))') S(i,:)
enddo
accu_d = 0.d0
accu_nd = 0.d0
@ -1883,8 +1883,8 @@ subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_
enddo
accu_nd = dsqrt(accu_nd) / dble(m)
!print *, ' accu_nd = ', accu_nd
!print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m)
print *, ' accu_nd = ', accu_nd
print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m)
! ---
@ -1944,6 +1944,96 @@ subroutine check_orthog(n, m, V, accu_d, accu_nd, S)
end subroutine check_orthog
! ---
subroutine reorder_degen_eigvec(n, e0, L0, R0)
implicit none
integer, intent(in) :: n
double precision, intent(in) :: e0(n)
double precision, intent(inout) :: L0(n,n), R0(n,n)
logical :: complex_root
integer :: i, j, k, m
double precision :: ei, ej, de, de_thr
double precision :: accu_d, accu_nd
integer, allocatable :: deg_num(:)
double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:)
! ---
allocate( deg_num(n) )
do i = 1, n
deg_num(i) = 1
enddo
de_thr = thr_degen_tc
do i = 1, n-1
ei = e0(i)
! already considered in degen vectors
if(deg_num(i).eq.0) cycle
do j = i+1, n
ej = e0(j)
de = dabs(ei - ej)
if(de .lt. de_thr) then
deg_num(i) = deg_num(i) + 1
deg_num(j) = 0
endif
enddo
enddo
do i = 1, n
if(deg_num(i) .gt. 1) then
print *, ' degen on', i, deg_num(i), e0(i)
endif
enddo
! ---
do i = 1, n
m = deg_num(i)
if(m .gt. 1) then
allocate(L(n,m))
allocate(R(n,m),S(m,m))
do j = 1, m
L(1:n,j) = L0(1:n,i+j-1)
R(1:n,j) = R0(1:n,i+j-1)
enddo
call dgemm( 'T', 'N', m, m, n, 1.d0 &
, L, size(L, 1), R, size(R, 1) &
, 0.d0, S, size(S, 1) )
print*,'Overlap matrix '
accu_nd = 0.D0
do j = 1, m
write(*,'(100(F16.10,X))')S(1:m,j)
do k = 1, m
if(j==k)cycle
accu_nd += dabs(S(j,k))
enddo
enddo
print*,'accu_nd = ',accu_nd
! if(accu_nd .gt.1.d-10)then
! stop
! endif
do j = 1, m
L0(1:n,i+j-1) = L(1:n,j)
R0(1:n,i+j-1) = R(1:n,j)
enddo
deallocate(L, R,S)
endif
enddo
end subroutine reorder_degen_eigvec
subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
@ -1987,11 +2077,11 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
enddo
enddo
!do i = 1, n
! if(deg_num(i) .gt. 1) then
! print *, ' degen on', i, deg_num(i), e0(i)
! endif
!enddo
do i = 1, n
if(deg_num(i) .gt. 1) then
print *, ' degen on', i, deg_num(i), e0(i)
endif
enddo
! ---
@ -2010,7 +2100,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
! ---
call impose_orthog_svd(n, m, L)
! call impose_orthog_svd(n, m, L)
call impose_orthog_svd(n, m, R)
!call impose_orthog_GramSchmidt(n, m, L)
!call impose_orthog_GramSchmidt(n, m, R)
@ -2031,6 +2121,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
!deallocate(S, S_inv_half)
call impose_biorthog_svd(n, m, L, R)
! call impose_biorthog_inverse(n, m, L, R)
!call impose_biorthog_qr(n, m, thr_d, thr_nd, L, R)
@ -2045,6 +2136,7 @@ subroutine impose_biorthog_degen_eigvec(n, e0, L0, R0)
endif
enddo
! call impose_biorthog_inverse(n, n, L0, R0)
end subroutine impose_biorthog_degen_eigvec
@ -2420,10 +2512,10 @@ subroutine impose_biorthog_svd(n, m, L, R)
, L, size(L, 1), R, size(R, 1) &
, 0.d0, S, size(S, 1) )
!print *, ' overlap bef SVD: '
!do i = 1, m
! write(*, '(1000(F16.10,X))') S(i,:)
!enddo
print *, ' overlap bef SVD: '
do i = 1, m
write(*, '(1000(F16.10,X))') S(i,:)
enddo
! ---
@ -2495,10 +2587,10 @@ subroutine impose_biorthog_svd(n, m, L, R)
, L, size(L, 1), R, size(R, 1) &
, 0.d0, S, size(S, 1) )
!print *, ' overlap aft SVD: '
!do i = 1, m
! write(*, '(1000(F16.10,X))') S(i,:)
!enddo
print *, ' overlap aft SVD: '
do i = 1, m
write(*, '(1000(F16.10,X))') S(i,:)
enddo
deallocate(S)
@ -2506,6 +2598,50 @@ subroutine impose_biorthog_svd(n, m, L, R)
end subroutine impose_biorthog_svd
subroutine impose_biorthog_inverse(n, m, L, R)
implicit none
integer, intent(in) :: n, m
double precision, intent(inout) :: L(n,m)
double precision, intent(in) :: R(n,m)
double precision, allocatable :: Lt(:,:),S(:,:)
integer :: i,j
allocate(Lt(m,n))
allocate(S(m,m))
call dgemm( 'T', 'N', m, m, n, 1.d0 &
, L, size(L, 1), R, size(R, 1) &
, 0.d0, S, size(S, 1) )
print *, ' overlap bef SVD: '
do i = 1, m
write(*, '(1000(F16.10,X))') S(i,:)
enddo
call get_pseudo_inverse(R,n,n,m,Lt,m,1.d-6)
do i = 1, m
do j = 1, n
L(j,i) = Lt(i,j)
enddo
enddo
! ---
call dgemm( 'T', 'N', m, m, n, 1.d0 &
, L, size(L, 1), R, size(R, 1) &
, 0.d0, S, size(S, 1) )
print *, ' overlap aft SVD: '
do i = 1, m
write(*, '(1000(F16.10,X))') S(i,:)
enddo
deallocate(S,Lt)
end subroutine impose_biorthog_svd
! ---
subroutine impose_weighted_biorthog_qr(m, n, thr_d, thr_nd, Vl, W, Vr)

View File

@ -5,7 +5,7 @@ interface: ezfio, provider
[nucl_label]
doc: Nuclear labels
type: character*(32)
type: character*(32)
size: (nuclei.nucl_num)
interface: ezfio, provider
@ -17,7 +17,7 @@ interface: ezfio, provider
[nucl_coord]
doc: Nuclear coordinates in the format (:, {x,y,z})
type: double precision
type: double precision
size: (nuclei.nucl_num,3)
interface: ezfio
@ -37,11 +37,12 @@ type: logical
doc: If true, the calculation uses periodic boundary conditions
interface: ezfio, provider, ocaml
default: false
[n_pts_charge]
type: integer
doc: Number of point charges to be added to the potential
interface: ezfio
default: 0
default: 0
[pts_charge_z]
type: double precision

View File

@ -15,13 +15,27 @@ program tc_natorb_bi_ortho
PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r
my_n_pt_a_grid = tc_grid1_a
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
if(j1b_type .ge. 100) then
my_extra_grid_becke = .True.
PROVIDE tc_grid2_a tc_grid2_r
my_n_pt_r_extra_grid = tc_grid2_r
my_n_pt_a_extra_grid = tc_grid2_a
touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
call write_int(6, my_n_pt_r_extra_grid, 'radial internal grid over')
call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over')
endif
read_wf = .True.
touch read_wf
call print_energy_and_mos()
call save_tc_natorb()
call print_angles_tc()
!call minimize_tc_orb_angles()
end
@ -35,9 +49,12 @@ subroutine save_tc_natorb()
print*,'Saving the natorbs '
provide natorb_tc_leigvec_ao natorb_tc_reigvec_ao
mo_l_coef = natorb_tc_leigvec_ao
mo_r_coef = natorb_tc_reigvec_ao
touch mo_l_coef mo_r_coef
call ezfio_set_bi_ortho_mos_mo_l_coef(natorb_tc_leigvec_ao)
call ezfio_set_bi_ortho_mos_mo_r_coef(natorb_tc_reigvec_ao)
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
call save_ref_determinant_nstates_1()
call ezfio_set_determinants_read_wf(.False.)

View File

@ -4,10 +4,11 @@
subroutine provide_all_three_ints_bi_ortho()
BEGIN_DOC
! routine that provides all necessary three-electron integrals
! routine that provides all necessary three-electron integrals
END_DOC
implicit none
PROVIDE ao_two_e_integrals_in_map
if(three_body_h_tc) then
@ -17,14 +18,14 @@ subroutine provide_all_three_ints_bi_ortho()
endif
if(three_e_4_idx_term) then
PROVIDE three_e_4_idx_direct_bi_ort three_e_4_idx_cycle_1_bi_ort three_e_4_idx_exch23_bi_ort three_e_4_idx_exch13_bi_ort
PROVIDE three_e_4_idx_direct_bi_ort three_e_4_idx_cycle_1_bi_ort three_e_4_idx_exch23_bi_ort three_e_4_idx_exch13_bi_ort
endif
if(pure_three_body_h_tc)then
provide three_body_ints_bi_ort
endif
if(.not. double_normal_ord .and. three_e_5_idx_term) then
PROVIDE three_e_5_idx_direct_bi_ort
PROVIDE three_e_5_idx_direct_bi_ort
elseif(double_normal_ord .and. (.not. three_e_5_idx_term)) then
PROVIDE normal_two_body_bi_orth
endif
@ -44,9 +45,9 @@ subroutine htilde_mu_mat_opt_bi_ortho_tot(key_j, key_i, Nint, htot)
!
! <key_j |H_tilde | key_i> where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis
!!
! Returns the total matrix element
! Returns the total matrix element
!! WARNING !!
!
!
! Non hermitian !!
!
END_DOC
@ -69,9 +70,9 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
!
! <key_j |H_tilde | key_i> where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis
!!
! Returns the detail of the matrix element in terms of single, two and three electron contribution.
! Returns the detail of the matrix element in terms of single, two and three electron contribution.
!! WARNING !!
!
!
! Non hermitian !!
!
END_DOC
@ -82,7 +83,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: hmono, htwoe, hthree, htot
integer :: degree
integer :: degree
hmono = 0.d0
htwoe = 0.d0
@ -94,7 +95,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
if(.not.pure_three_body_h_tc) then
if(degree .gt. 2) return
if(degree == 0) then
call diag_htilde_mu_mat_fock_bi_ortho (Nint, key_i, hmono, htwoe, hthree, htot)
else if (degree == 1) then
@ -103,7 +104,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
endif
else
else
if(degree .gt. 3) return
@ -122,7 +123,7 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
if(degree==0) then
htot += nuclear_repulsion
endif
end
! ---
@ -133,9 +134,9 @@ subroutine htilde_mu_mat_opt_bi_ortho_no_3e(key_j, key_i, Nint, htot)
!
! <key_j |H_tilde | key_i> where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis
!!
! Returns the detail of the matrix element WITHOUT ANY CONTRIBUTION FROM THE THREE ELECTRON TERMS
! Returns the detail of the matrix element WITHOUT ANY CONTRIBUTION FROM THE THREE ELECTRON TERMS
!! WARNING !!
!
!
! Non hermitian !!
!
END_DOC
@ -146,7 +147,7 @@ subroutine htilde_mu_mat_opt_bi_ortho_no_3e(key_j, key_i, Nint, htot)
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
double precision, intent(out) :: htot
integer :: degree
integer :: degree
htot = 0.d0
@ -164,7 +165,7 @@ subroutine htilde_mu_mat_opt_bi_ortho_no_3e(key_j, key_i, Nint, htot)
if(degree==0) then
htot += nuclear_repulsion
endif
end
! ---

View File

@ -13,6 +13,17 @@ program tc_bi_ortho
my_n_pt_a_grid = tc_grid1_a
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
if(j1b_type .ge. 100) then
my_extra_grid_becke = .True.
PROVIDE tc_grid2_a tc_grid2_r
my_n_pt_r_extra_grid = tc_grid2_r
my_n_pt_a_extra_grid = tc_grid2_a
touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
call write_int(6, my_n_pt_r_extra_grid, 'radial internal grid over')
call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over')
endif
read_wf = .True.
touch read_wf

View File

@ -34,4 +34,19 @@ subroutine test
do i= 1, 3
print*,tc_bi_ortho_dipole(i,1)
enddo
integer, allocatable :: occ(:,:)
integer :: n_occ_ab(2)
allocate(occ(N_int*bit_kind_size,2))
call bitstring_to_list_ab(ref_bitmask, occ, n_occ_ab, N_int)
integer :: ispin,j,jorb
double precision :: accu
accu = 0.d0
do ispin=1, 2
do i = 1, n_occ_ab(ispin)
jorb = occ(i,ispin)
accu += mo_bi_orth_bipole_z(jorb,jorb)
enddo
enddo
print*,'accu = ',accu
end

View File

@ -15,7 +15,7 @@ BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho, (N_det,N_det)]
integer :: i, j
double precision :: htot
PROVIDE N_int
call provide_all_three_ints_bi_ortho
i = 1
j = 1

View File

@ -29,9 +29,22 @@
write(*, '(100(F16.10,X))') -dm_tmp(:,i)
enddo
print *, ' Transition density matrix AO'
do i = 1, ao_num
write(*, '(100(F16.10,X))') tc_transition_matrix_ao(:,i,1,1)
enddo
stop
thr_d = 1.d-6
thr_nd = 1.d-6
thr_deg = 1.d-3
do i = 1, mo_num
do j = 1, mo_num
if(dabs(dm_tmp(j,i)).lt.thr_d)then
dm_tmp(j,i) = 0.d0
endif
enddo
enddo
! if(n_core_orb.ne.0)then
! call diag_mat_per_fock_degen_core( fock_diag, dm_tmp, list_core, n_core_orb, mo_num, thr_d, thr_nd, thr_deg &
! , natorb_tc_leigvec_mo, natorb_tc_reigvec_mo, natorb_tc_eigval)

View File

@ -90,6 +90,7 @@
enddo
enddo
enddo
print*,'tc_bi_ortho_dipole(3) elec = ',tc_bi_ortho_dipole(3,1)
nuclei_part = 0.d0
do m = 1, 3

View File

@ -91,7 +91,7 @@
print *, ' parameters for nuclei jastrow'
print *, ' i, Z, j1b_pen, j1b_pen_coef'
do i = 1, nucl_num
write(*,"(I4, 2x, 3(E15.7, 2X))") i, nucl_charge(i), j1b_pen(i), j1b_pen_coef(i)
write(*,'(I4, 2x, 3(E15.7, 2X))') i, nucl_charge(i), j1b_pen(i), j1b_pen_coef(i)
enddo
END_PROVIDER

View File

@ -8,15 +8,15 @@ function run_Ne() {
rm -rf Ne_tc_scf
echo Ne > Ne.xyz
qp create_ezfio -b cc-pcvdz Ne.xyz -o Ne_tc_scf
qp run scf
qp run scf
qp set ao_two_e_erf_ints mu_erf 0.87
qp set hamiltonian mu_erf 0.87
qp set tc_keywords j1b_type 3
qp set tc_keywords j1b_pen [1.5]
qp set tc_keywords bi_ortho True
qp set tc_keywords j1b_pen [1.5]
qp set tc_keywords bi_ortho True
qp set tc_keywords test_cycle_tc True
qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
eref=-128.552134
energy="$(qp get tc_scf bitc_energy)"
eq $energy $eref 1e-6
@ -24,22 +24,22 @@ function run_Ne() {
@test "Ne" {
run_Ne
run_Ne
}
function run_C() {
rm -rf C_tc_scf
echo C > C.xyz
qp create_ezfio -b cc-pcvdz C.xyz -o C_tc_scf -m 3
qp run scf
qp run scf
qp set ao_two_e_erf_ints mu_erf 0.87
qp set hamiltonian mu_erf 0.87
qp set tc_keywords j1b_type 3
qp set tc_keywords j1b_pen [1.5]
qp set tc_keywords bi_ortho True
qp set tc_keywords j1b_pen [1.5]
qp set tc_keywords bi_ortho True
qp set tc_keywords test_cycle_tc True
qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
eref=-37.691254356408791
energy="$(qp get tc_scf bitc_energy)"
eq $energy $eref 1e-6
@ -47,7 +47,7 @@ function run_C() {
@test "C" {
run_C
run_C
}
@ -55,15 +55,15 @@ function run_O() {
rm -rf O_tc_scf
echo O > O.xyz
qp create_ezfio -b cc-pcvdz O.xyz -o O_tc_scf -m 3
qp run scf
qp run scf
qp set ao_two_e_erf_ints mu_erf 0.87
qp set hamiltonian mu_erf 0.87
qp set tc_keywords j1b_type 3
qp set tc_keywords j1b_pen [1.5]
qp set tc_keywords bi_ortho True
qp set tc_keywords j1b_pen [1.5]
qp set tc_keywords bi_ortho True
qp set tc_keywords test_cycle_tc True
qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
eref=-74.814687229354590
energy="$(qp get tc_scf bitc_energy)"
eq $energy $eref 1e-6
@ -71,7 +71,7 @@ function run_O() {
@test "O" {
run_O
run_O
}
@ -79,16 +79,16 @@ function run_O() {
function run_ch2() {
rm -rf ch2_tc_scf
cp ${QP_ROOT}/tests/input/ch2.xyz .
qp create_ezfio -b "C:cc-pcvdz|H:cc-pvdz" ch2.xyz -o ch2_tc_scf
qp run scf
qp create_ezfio -b "C:cc-pcvdz|H:cc-pvdz" ch2.xyz -o ch2_tc_scf
qp run scf
qp set ao_two_e_erf_ints mu_erf 0.87
qp set hamiltonian mu_erf 0.87
qp set tc_keywords j1b_type 3
qp set tc_keywords j1b_pen '[1.5,10000,10000]'
qp set tc_keywords bi_ortho True
qp set tc_keywords test_cycle_tc True
qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
eref=-38.903247818077737
energy="$(qp get tc_scf bitc_energy)"
eq $energy $eref 1e-6
@ -96,6 +96,6 @@ function run_ch2() {
@test "ch2" {
run_ch2
run_ch2
}

View File

@ -239,7 +239,7 @@ subroutine rh_tcscf_diis()
write(json_unit, json_real_fmt) ' delta Energy ', e_delta
write(json_unit, json_real_fmt) ' DIIS error ', er_DIIS
write(json_unit, json_real_fmt) ' level_shift ', level_shift_tcscf
write(json_unit, json_int_fmt) ' DIIS ', dim_DIIS
write(json_unit, json_int_fmtx) ' DIIS ', dim_DIIS
write(json_unit, json_real_fmt) ' Wall time (min)', (t1-t0)/60.d0
call unlock_io

View File

@ -402,6 +402,7 @@ subroutine print_energy_and_mos(good_angles)
print *, ' TC energy = ', TC_HF_energy
print *, ' TC SCF energy gradient = ', grad_non_hermit
print *, ' Max angle Left/right = ', max_angle_left_right
call print_angles_tc()
if(max_angle_left_right .lt. thresh_lr_angle) then
print *, ' Maximum angle BELOW 45 degrees, everthing is OK !'

View File

@ -30,9 +30,22 @@ subroutine main()
allocate(Sl(mo_num,mo_num), Sr(mo_num,mo_num), Pf(mo_num,mo_num))
call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 &
, mo_l_coef, size(mo_l_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
, 0.d0, Sl, size(Sl, 1) )
call LTxSxR(ao_num, mo_num, mo_l_coef, ao_overlap, mo_r_coef, Sl)
!call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 &
! , mo_l_coef, size(mo_l_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
! , 0.d0, Sl, size(Sl, 1) )
print *, ''
print *, ' left-right orthog matrix:'
do i = 1, mo_num
write(*,'(100(F8.4,X))') Sl(:,i)
enddo
call LTxSxR(ao_num, mo_num, mo_l_coef, ao_overlap, mo_l_coef, Sl)
!call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 &
! , mo_l_coef, size(mo_l_coef, 1), mo_l_coef, size(mo_l_coef, 1) &
! , 0.d0, Sl, size(Sl, 1) )
print *, ''
print *, ' left-orthog matrix:'
@ -40,9 +53,10 @@ subroutine main()
write(*,'(100(F8.4,X))') Sl(:,i)
enddo
call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 &
, mo_r_coef, size(mo_r_coef, 1), mo_r_coef, size(mo_r_coef, 1) &
, 0.d0, Sr, size(Sr, 1) )
call LTxSxR(ao_num, mo_num, mo_r_coef, ao_overlap, mo_r_coef, Sr)
! call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 &
! , mo_r_coef, size(mo_r_coef, 1), mo_r_coef, size(mo_r_coef, 1) &
! , 0.d0, Sr, size(Sr, 1) )
print *, ''
print *, ' right-orthog matrix:'

View File

@ -1,5 +1,4 @@
fci
mo_two_e_erf_ints
aux_quantities
hartree_fock
two_body_rdm

View File

@ -0,0 +1,168 @@
program molden_detachment_attachment
implicit none
read_wf=.True.
touch read_wf
call molden_attachment
end
subroutine molden_attachment
implicit none
BEGIN_DOC
! Produces a Molden file
END_DOC
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
integer :: i,j,k,l
double precision, parameter :: a0 = 0.529177249d0
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'.attachement.mol'
print*,'output = ',trim(output)
i_unit_output = getUnitAndOpen(output,'w')
write(i_unit_output,'(A)') '[Molden Format]'
write(i_unit_output,'(A)') '[Atoms] Angs'
do i = 1, nucl_num
write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))') &
trim(element_name(int(nucl_charge(i)))), &
i, &
int(nucl_charge(i)), &
nucl_coord(i,1)*a0, nucl_coord(i,2)*a0, nucl_coord(i,3)*a0
enddo
write(i_unit_output,'(A)') '[GTO]'
character*(1) :: character_shell
integer :: i_shell,i_prim,i_ao
integer :: iorder(ao_num)
integer :: nsort(ao_num)
i_shell = 0
i_prim = 0
do i=1,nucl_num
write(i_unit_output,*) i, 0
do j=1,nucl_num_shell_aos(i)
i_shell +=1
i_ao = nucl_list_shell_aos(i,j)
character_shell = trim(ao_l_char(i_ao))
write(i_unit_output,*) character_shell, ao_prim_num(i_ao), '1.00'
do k = 1, ao_prim_num(i_ao)
i_prim +=1
write(i_unit_output,'(ES20.10,2X,ES20.10)') ao_expo(i_ao,k), ao_coef(i_ao,k)
enddo
l = i_ao
do while ( ao_l(l) == ao_l(i_ao) )
nsort(l) = i*10000 + j*100
l += 1
if (l > ao_num) exit
enddo
enddo
write(i_unit_output,*)''
enddo
do i=1,ao_num
iorder(i) = i
! p
if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 3
! d
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 6
! f
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 6
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 7
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 8
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 9
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 10
! g
else if ((ao_power(i,1) == 4 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 4 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 4 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 6
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 7
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 8
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 9
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 10
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 11
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 12
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 13
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 14
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 15
endif
enddo
call isort(nsort,iorder,ao_num)
write(i_unit_output,'(A)') '[MO]'
integer :: istate
istate = 2
do i=1,n_dettachment(istate)
write (i_unit_output,*) 'Sym= 1'
write (i_unit_output,*) 'Ene=', dettachment_numbers_sorted(i,istate)
write (i_unit_output,*) 'Spin= Alpha'
write (i_unit_output,*) 'Occup=', dettachment_numbers_sorted(i,istate)
do j=1,ao_num
write(i_unit_output, '(I6,2X,ES20.10)') j, dettachment_orbitals(iorder(j),i,istate)
enddo
enddo
do i=1,n_attachment(istate)
write (i_unit_output,*) 'Sym= 1'
write (i_unit_output,*) 'Ene=', attachment_numbers_sorted(i,istate)
write (i_unit_output,*) 'Spin= Alpha'
write (i_unit_output,*) 'Occup=', attachment_numbers_sorted(i,istate)
do j=1,ao_num
write(i_unit_output, '(I6,2X,ES20.10)') j, attachment_orbitals(iorder(j),i,istate)
enddo
enddo
close(i_unit_output)
end

1
src/trexio/LIB Normal file
View File

@ -0,0 +1 @@
-ltrexio

View File

@ -123,7 +123,7 @@
state_av_act_2_rdm_spin_trace_mo = state_av_act_2_rdm_ab_mo &
+ state_av_act_2_rdm_aa_mo &
+ state_av_act_2_rdm_bb_mo
!
! call orb_range_2_rdm_state_av_openmp(state_av_act_2_rdm_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
call wall_time(wall_2)

View File

@ -2,7 +2,7 @@ program test_2_rdm
implicit none
read_wf = .True.
touch read_wf
! call routine_active_only
call routine_active_only
call routine_full_mos
end

1
src/zmq/LIB Normal file
View File

@ -0,0 +1 @@
-lf77zmq -lzmq