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

Merge pull request #45 from QuantumPackage/dev-lcpq

Browse Source 
Dev lcpq
This commit is contained in:
Anthony Scemama 2019-06-08 16:24:33 +02:00 committed by GitHub
commit 231b071d78
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
30 changed files with 625 additions and 122 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

96
CITATION.cff
View File

@ -1,32 +1,100 @@
# YAML 1.2 # YAML 1.2
# Metadata for citation of this software according to the CFF format (https://citation-file-format.github.io/) # Metadata for citation of this software according to the CFF format (https://citation-file-format.github.io/)
cff-version: 1.0.3 cff-version: 1.0.3
message: If you use this software, please cite it using these metadata. message: "If you use this software, please cite it using these metadata."
title: Quantum Package title: Quantum Package
doi: 10.5281/zenodo.825872 doi: 10.1021/acs.jctc.9b00176
authors: authors:
- given-names: Anthony
family-names: Scemama
affiliation: Laboratoire de Chimie et Physique Quantiques / CNRS
- given-names: Yann - given-names: Yann
family-names: Garniron family-names: Garniron
affiliation: Laboratoire de Chimie et Physique Quantiques / CNRS affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France
- given-names: Michel
family-names: Caffarel
affiliation: Laboratoire de Chimie et Physique Quantiques / CNRS
- given-names: Thomas - given-names: Thomas
family-names: Applencourt family-names: Applencourt
affiliation: Argonne National Lab affiliation: Computational Science Division, Argonne National Laboratory, Argonne, Illinois 60439, United States
- given-names: Kevin - given-names: Kevin
family-names: Gasperich family-names: Gasperich
affiliation: Argonne National Lab affiliation: Department of Chemistry, University of Pittsburgh, Pittsburgh, Pennsylvania 15260, United States
- given-names: Anouar - given-names: Anouar
family-names: Benali family-names: Benali
affiliation: Argonne National Lab affiliation: Computational Science Division, Argonne National Laboratory, Argonne, Illinois 60439, United States
- given-names: Anthony
family-names: Ferté
affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France
- given-names: Julien
family-names: Paquier
affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France
- given-names: Barthélémy
family-names: Pradines
affiliation: Institut des Sciences du Calcul et des Données, Sorbonne Université, F-75005 Paris, France
- given-names: Roland
family-names: Assaraf
affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France
- given-names: Peter
family-names: Reinhardt
affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France
- given-names: Julien
family-names: Toulouse
affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France
- given-names: Pierrette
family-names: Barbaresco
affiliation: CALMIP, Université de Toulouse, CNRS, INPT, INSA, UPS, UMS 3667, Toulouse, France
- given-names: Nicolas
family-names: Renon
affiliation: CALMIP, Université de Toulouse, CNRS, INPT, INSA, UPS, UMS 3667, Toulouse, France
- given-names: Grégoire
family-names: David
affiliation: Aix-Marseille Univ, CNRS, ICR, Marseille, France
- given-names: Jean-Paul
family-names: Malrieu
affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France
- given-names: Mickaël
family-names: Véril
affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France
- given-names: Michel
family-names: Caffarel
affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France
- given-names: Pierre-François
family-names: Loos
affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France
- given-names: Emmanuel - given-names: Emmanuel
family-names: Giner family-names: Giner
affiliation: Laboratoire de Chimie Theorique / CNRS affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France
- given-names: Anthony
family-names: Scemama
affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France
abstract: "Quantum chemistry is a discipline which relies heavily on very
expensive numerical computations. The scaling of correlated wave function
methods lies, in their standard implementation, between O(N^5) and O(exp(N)),
where N is proportional to the system size. Therefore, performing accurate
calculations on chemically meaningful systems requires (i) approximations that
can lower the computational scaling and (ii) efficient implementations that
take advantage of modern massively parallel architectures. Quantum Package is
an open-source programming environment for quantum chemistry specially designed
for wave function methods. Its main goal is the development of
determinant-driven selected configuration interaction (sCI) methods and
multireference second-order perturbation theory (PT2). The determinant-driven
framework allows the programmer to include any arbitrary set of determinants in
the reference space, hence providing greater methodological freedom. The sCI
method implemented in Quantum Package is based on the CIPSI (Configuration
Interaction using a Perturbative Selection made Iteratively) algorithm which
complements the variational sCI energy with a PT2 correction. Additional
external plugins have been recently added to perform calculations with
multireference coupled cluster theory and range-separated density-functional
theory. All the programs are developed with the IRPF90 code generator, which
simplifies collaborative work and the development of new features. Quantum
Package strives to allow easy implementation and experimentation of new
methods, while making parallel computation as simple and efficient as possible
on modern supercomputer architectures. Currently, the code enables, routinely,
to realize runs on roughly 2 000 CPU cores, with tens of millions of
determinants in the reference space. Moreover, we have been able to push up to
12 288 cores in order to test its parallel efficiency. In the present
manuscript, we also introduce some key new developments: (i) a renormalized
second-order perturbative correction for efficient extrapolation to the full CI
limit and (ii) a stochastic version of the CIPSI selection performed
simultaneously to the PT2 calculation at no extra cost."
version: '2.0' version: '2.0'
date-released: 2019-02-11 url: https://quantumpackage.github.io/qp2/
date-released: 2019-05-13
repository-code: https://github.com/QuantumPackage/qp2 repository-code: https://github.com/QuantumPackage/qp2
keywords: [ "computational chemistry", "configuration interaction", "cipsi", "perturbation theory" ]
license: AGPL-3.0-or-later license: AGPL-3.0-or-later

53
INSTALL.rst
View File

@ -45,6 +45,8 @@ Requirements
- |ZeroMQ| : networking library - |ZeroMQ| : networking library
- `GMP <https://gmplib.org/>`_ : Gnu Multiple Precision Arithmetic Library - `GMP <https://gmplib.org/>`_ : Gnu Multiple Precision Arithmetic Library
- |OCaml| compiler with |OPAM| package manager - |OCaml| compiler with |OPAM| package manager
- `Bubblewrap <https://github.com/projectatomic/bubblewrap>`_ : Sandboxing tool required by Opam
- `libcap https://git.kernel.org/pub/scm/linux/kernel/git/morgan/libcap.git`_ : POSIX capabilities required by Bubblewrap
- |Ninja| : a parallel build system - |Ninja| : a parallel build system
@ -86,6 +88,8 @@ The following packages are supported by the :command:`configure` installer:
* zeromq * zeromq
* f77zmq * f77zmq
* gmp * gmp
* libcap
* bwrap
* ocaml ( :math:`\approx` 10 minutes) * ocaml ( :math:`\approx` 10 minutes)
* ezfio * ezfio
* docopt * docopt
@ -243,6 +247,55 @@ With Debian or Ubuntu, you can use
sudo apt install libgmp-dev sudo apt install libgmp-dev
libcap
------
Libcap is a library for getting and setting POSIX.1e draft 15 capabilities.
* Download the latest version of libcap here:
`<https://git.kernel.org/pub/scm/linux/kernel/git/morgan/libcap.git/snapshot/libcap-2.25.tar.gz>`_
and move it in the :file:`${QP_ROOT}/external` directory
* Extract the archive, go into the :file:`libcap-*/libcap` directory and run
the following command
.. code:: bash
prefix=$QP_ROOT make install
With Debian or Ubuntu, you can use
.. code:: bash
sudo apt install libcap-dev
Bubblewrap
----------
Bubblewrap is an unprivileged sandboxing tool.
* Download Bubblewrap here:
`<https://github.com/projectatomic/bubblewrap/releases/download/v0.3.3/bubblewrap-0.3.3.tar.xz>`_
and move it in the :file:`${QP_ROOT}/external` directory
* Extract the archive, go into the :file:`bubblewrap-*` directory and run
the following commands
.. code:: bash
./configure --prefix=$QP_ROOT && make -j 8
make install-exec-am
With Debian or Ubuntu, you can use
.. code:: bash
sudo apt install bubblewrap
OCaml OCaml
----- -----

5
README.md
View File

@ -1,12 +1,13 @@
# Quantum Package 2.0 # Quantum Package 2.0
*Quantum package 2.0: an open-source determinant-driven suite of programs*\ [*Quantum package 2.0: an open-source determinant-driven suite of programs*](https://pubs.acs.org/doi/10.1021/acs.jctc.9b00176)\
Y. Garniron, K. Gasperich, T. Applencourt, A. Benali, A. Ferté, J. Paquier, B. Pradines, R. Assaraf, P. Reinhardt, J. Toulouse, P. Barbaresco, N. Renon, G. David, J. P. Malrieu, M. Véril, M. Caffarel, P. F. Loos, E. Giner and A. Scemama\ Y. Garniron, K. Gasperich, T. Applencourt, A. Benali, A. Ferté, J. Paquier, B. Pradines, R. Assaraf, P. Reinhardt, J. Toulouse, P. Barbaresco, N. Renon, G. David, J. P. Malrieu, M. Véril, M. Caffarel, P. F. Loos, E. Giner and A. Scemama\
J. Chem. Theory Comput. (2019)\
https://arxiv.org/abs/1902.08154 https://arxiv.org/abs/1902.08154
![QP](https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png) <img src="https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png" width="250">
# Getting started # Getting started

2
config/ifort_avx.cfg
View File

@ -32,7 +32,7 @@ OPENMP : 1 ; Append OpenMP flags
# #
[OPT] [OPT]
FC : -traceback FC : -traceback
FCFLAGS : -xAVX -O2 -ip -ftz -g FCFLAGS : -march=corei7-avx -O2 -ip -ftz -g
# Profiling flags # Profiling flags
################# #################

4
config/ifort_avx_mpi.cfg
View File

@ -31,14 +31,14 @@ OPENMP : 1 ; Append OpenMP flags
# -ftz : Flushes denormal results to zero # -ftz : Flushes denormal results to zero
# #
[OPT] [OPT]
FCFLAGS : -xAVX -O2 -ip -ftz -g -traceback FCFLAGS : -march=corei7-avx -O2 -ip -ftz -g -traceback
# Profiling flags # Profiling flags
################# #################
# #
[PROFILE] [PROFILE]
FC : -p -g FC : -p -g
FCFLAGS : -xSSE4.2 -O2 -ip -ftz FCFLAGS : -march=corei7 -O2 -ip -ftz
# Debugging flags # Debugging flags

44
configure vendored
View File

@ -175,7 +175,7 @@ if [[ "${PACKAGES}.x" != ".x" ]] ; then
fi fi
if [[ ${PACKAGES} = all ]] ; then if [[ ${PACKAGES} = all ]] ; then
PACKAGES="zlib ninja irpf90 zeromq f77zmq gmp ocaml ezfio docopt resultsFile bats" PACKAGES="zlib ninja irpf90 zeromq f77zmq gmp libcap bwrap ocaml ezfio docopt resultsFile bats"
fi fi
@ -206,6 +206,32 @@ EOF
make install make install
EOF EOF
elif [[ ${PACKAGE} = libcap ]] ; then
download \
"https://git.kernel.org/pub/scm/linux/kernel/git/morgan/libcap.git/snapshot/libcap-2.25.tar.gz" \
"${QP_ROOT}"/external/libcap.tar.gz
execute << EOF
cd "\${QP_ROOT}"/external
tar --gunzip --extract --file libcap.tar.gz
rm libcap.tar.gz
cd libcap-*/libcap
prefix=$QP_ROOT make install
EOF
elif [[ ${PACKAGE} = bwrap ]] ; then
download \
"https://github.com/projectatomic/bubblewrap/releases/download/v0.3.3/bubblewrap-0.3.3.tar.xz" \
"${QP_ROOT}"/external/bwrap.tar.xz
execute << EOF
cd "\${QP_ROOT}"/external
tar --xz --extract --file bwrap.tar.xz
rm bwrap.tar.xz
cd bubblewrap*
./configure --prefix=$QP_ROOT && make -j 8
make install-exec-am
EOF
elif [[ ${PACKAGE} = irpf90 ]] ; then elif [[ ${PACKAGE} = irpf90 ]] ; then
@ -276,7 +302,7 @@ EOF
rm ${QP_ROOT}/external/opam_installer.sh rm ${QP_ROOT}/external/opam_installer.sh
source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true
${QP_ROOT}/bin/opam init --disable-sandboxing --verbose --yes ${QP_ROOT}/bin/opam init --verbose --yes
eval $(${QP_ROOT}/bin/opam env) eval $(${QP_ROOT}/bin/opam env)
opam install -y ${OCAML_PACKAGES} || exit 1 opam install -y ${OCAML_PACKAGES} || exit 1
@ -290,7 +316,7 @@ EOF
| sh \${QP_ROOT}/external/opam_installer.sh | sh \${QP_ROOT}/external/opam_installer.sh
rm \${QP_ROOT}/external/opam_installer.sh rm \${QP_ROOT}/external/opam_installer.sh
source \${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true source \${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true
\${QP_ROOT}/bin/opam init --disable-sandboxing --verbose --yes \${QP_ROOT}/bin/opam init --verbose --yes
eval \$(\${QP_ROOT}/bin/opam env) eval \$(\${QP_ROOT}/bin/opam env)
opam install -y \${OCAML_PACKAGES} || exit 1 opam install -y \${OCAML_PACKAGES} || exit 1
EOF EOF
@ -399,6 +425,18 @@ if [[ ${ZLIB} = $(not_found) ]] ; then
fail fail
fi fi
BWRAP=$(find_exe bwrap)
if [[ ${BWRAP} = $(not_found) ]] ; then
error "Bubblewrap (bwrap) is not installed."
fail
fi
LIBCAP=$(find_lib -lcap)
if [[ ${LIBCAP} = $(not_found) ]] ; then
error "Libcap (libcap) is not installed."
fail
fi
OPAM=$(find_exe opam) OPAM=$(find_exe opam)
if [[ ${OPAM} = $(not_found) ]] ; then if [[ ${OPAM} = $(not_found) ]] ; then
error "OPAM (ocaml) package manager is not installed." error "OPAM (ocaml) package manager is not installed."

3
ocaml/Input_determinants_by_hand.ml
View File

@ -81,9 +81,6 @@ end = struct
;; ;;
let write_n_det n = let write_n_det n =
let n_det_old =
Ezfio.get_determinants_n_det ()
in
Det_number.to_int n Det_number.to_int n
|> Ezfio.set_determinants_n_det |> Ezfio.set_determinants_n_det
;; ;;

2
src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f
View File

@ -16,6 +16,8 @@ double precision function ao_two_e_integral_erf(i,j,k,l)
integer :: iorder_p(3), iorder_q(3) integer :: iorder_p(3), iorder_q(3)
double precision :: ao_two_e_integral_schwartz_accel_erf double precision :: ao_two_e_integral_schwartz_accel_erf
provide mu_erf
if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then if (ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
ao_two_e_integral_erf = ao_two_e_integral_schwartz_accel_erf(i,j,k,l) ao_two_e_integral_erf = ao_two_e_integral_schwartz_accel_erf(i,j,k,l)
return return

94
src/ao_two_e_ints/map_integrals.irp.f
View File

@ -279,6 +279,100 @@ subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_z
end end
subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out_val_index,non_zero_int)
use map_module
implicit none
BEGIN_DOC
! Gets multiple AO bi-electronic integral from the AO map .
! All non-zero i are retrieved for j,k,l fixed.
END_DOC
double precision, intent(in) :: thresh
integer, intent(in) :: j,l, sze,sze_max
real(integral_kind), intent(out) :: out_val(sze_max)
integer, intent(out) :: out_val_index(2,sze_max),non_zero_int
integer :: i,k
integer(key_kind) :: hash
double precision :: tmp
PROVIDE ao_two_e_integrals_in_map
non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
non_zero_int = 0
do k = 1, sze
do i = 1, sze
integer, external :: ao_l4
double precision, external :: ao_two_e_integral
!DIR$ FORCEINLINE
if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then
cycle
endif
call two_e_integrals_index(i,j,k,l,hash)
call map_get(ao_integrals_map, hash,tmp)
if (dabs(tmp) < thresh ) cycle
non_zero_int = non_zero_int+1
out_val_index(1,non_zero_int) = i
out_val_index(2,non_zero_int) = k
out_val(non_zero_int) = tmp
enddo
enddo
end
subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,sze_max,out_val,out_val_index,non_zero_int)
use map_module
implicit none
BEGIN_DOC
! Gets multiple AO two-electron integrals from the AO map .
! All non-zero i are retrieved for j,k,l fixed.
END_DOC
double precision, intent(in) :: thresh
integer, intent(in) :: sze_max
integer, intent(in) :: j,l, n_list,list(2,sze_max)
real(integral_kind), intent(out) :: out_val(sze_max)
integer, intent(out) :: out_val_index(2,sze_max),non_zero_int
integer :: i,k
integer(key_kind) :: hash
double precision :: tmp
PROVIDE ao_two_e_integrals_in_map
non_zero_int = 0
if (ao_overlap_abs(j,l) < thresh) then
out_val = 0.d0
return
endif
non_zero_int = 0
integer :: kk
do kk = 1, n_list
k = list(1,kk)
i = list(2,kk)
integer, external :: ao_l4
double precision, external :: ao_two_e_integral
!DIR$ FORCEINLINE
if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then
cycle
endif
call two_e_integrals_index(i,j,k,l,hash)
call map_get(ao_integrals_map, hash,tmp)
if (dabs(tmp) < thresh ) cycle
non_zero_int = non_zero_int+1
out_val_index(1,non_zero_int) = i
out_val_index(2,non_zero_int) = k
out_val(non_zero_int) = tmp
enddo
end
function get_ao_map_size() function get_ao_map_size()
implicit none implicit none
integer (map_size_kind) :: get_ao_map_size integer (map_size_kind) :: get_ao_map_size

6
src/becke_numerical_grid/EZFIO.cfg
View File

@ -8,3 +8,9 @@ default: 2
type: integer type: integer
doc: Total number of grid points doc: Total number of grid points
interface: ezfio interface: ezfio
[thresh_grid]
type: double precision
doc: threshold on the weight of a given grid point
interface: ezfio,provider,ocaml
default: 1.e-20

9
src/becke_numerical_grid/atomic_number.irp.f Normal file
View File

@ -0,0 +1,9 @@
BEGIN_PROVIDER [ integer, grid_atomic_number, (nucl_num) ]
implicit none
BEGIN_DOC
! Atomic number used to adjust the grid
END_DOC
grid_atomic_number(:) = max(1,int(nucl_charge(:)))
END_PROVIDER

6
src/becke_numerical_grid/grid_becke.irp.f
View File

@ -146,7 +146,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_
x = grid_points_radial(j) x = grid_points_radial(j)
! value of the radial coordinate for the integration ! value of the radial coordinate for the integration
r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x) r = knowles_function(alpha_knowles(grid_atomic_number(i)),m_knowles,x)
! explicit values of the grid points centered around each atom ! explicit values of the grid points centered around each atom
do k = 1, n_points_integration_angular do k = 1, n_points_integration_angular
@ -232,8 +232,8 @@ BEGIN_PROVIDER [double precision, final_weight_at_r, (n_points_integration_angul
do i = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom do i = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom
x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1] x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1]
do k = 1, n_points_integration_angular ! for each angular point attached to the "jth" atom do k = 1, n_points_integration_angular ! for each angular point attached to the "jth" atom
contrib_integration = derivative_knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)& contrib_integration = derivative_knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)&
*knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)**2 *knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)**2
final_weight_at_r(k,i,j) = weights_angular_points(k) * weight_at_r(k,i,j) * contrib_integration * dr_radial_integral final_weight_at_r(k,i,j) = weights_angular_points(k) * weight_at_r(k,i,j) * contrib_integration * dr_radial_integral
if(isnan(final_weight_at_r(k,i,j)))then if(isnan(final_weight_at_r(k,i,j)))then
print*,'isnan(final_weight_at_r(k,i,j))' print*,'isnan(final_weight_at_r(k,i,j))'

53
src/becke_numerical_grid/grid_becke_per_atom.irp.f Normal file
View File

@ -0,0 +1,53 @@
BEGIN_PROVIDER [integer, n_pts_per_atom, (nucl_num)]
&BEGIN_PROVIDER [integer, n_pts_max_per_atom]
BEGIN_DOC
! Number of points which are non zero
END_DOC
integer :: i,j,k,l
n_pts_per_atom = 0
do j = 1, nucl_num
do i = 1, n_points_radial_grid -1
do k = 1, n_points_integration_angular
if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
cycle
endif
n_pts_per_atom(j) += 1
enddo
enddo
enddo
n_pts_max_per_atom = maxval(n_pts_per_atom)
END_PROVIDER
BEGIN_PROVIDER [double precision, final_grid_points_per_atom, (3,n_pts_max_per_atom,nucl_num)]
&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_per_atom, (n_pts_max_per_atom,nucl_num) ]
&BEGIN_PROVIDER [integer, index_final_points_per_atom, (3,n_pts_max_per_atom,nucl_num) ]
&BEGIN_PROVIDER [integer, index_final_points_per_atom_reverse, (n_points_integration_angular,n_points_radial_grid,nucl_num) ]
implicit none
integer :: i,j,k,l,i_count(nucl_num)
double precision :: r(3)
i_count = 0
do j = 1, nucl_num
do i = 1, n_points_radial_grid -1
do k = 1, n_points_integration_angular
if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
cycle
endif
i_count(j) += 1
final_grid_points_per_atom(1,i_count(j),j) = grid_points_per_atom(1,k,i,j)
final_grid_points_per_atom(2,i_count(j),j) = grid_points_per_atom(2,k,i,j)
final_grid_points_per_atom(3,i_count(j),j) = grid_points_per_atom(3,k,i,j)
final_weight_at_r_vector_per_atom(i_count(j),j) = final_weight_at_r(k,i,j)
index_final_points_per_atom(1,i_count(j),j) = k
index_final_points_per_atom(2,i_count(j),j) = i
index_final_points_per_atom(3,i_count(j),j) = j
index_final_points_per_atom_reverse(k,i,j) = i_count(j)
enddo
enddo
enddo
END_PROVIDER

13
src/becke_numerical_grid/grid_becke_vector.irp.f
View File

@ -1,5 +1,6 @@
BEGIN_PROVIDER [integer, n_points_final_grid] BEGIN_PROVIDER [integer, n_points_final_grid]
implicit none
BEGIN_DOC BEGIN_DOC
! Number of points which are non zero ! Number of points which are non zero
END_DOC END_DOC
@ -8,9 +9,9 @@ BEGIN_PROVIDER [integer, n_points_final_grid]
do j = 1, nucl_num do j = 1, nucl_num
do i = 1, n_points_radial_grid -1 do i = 1, n_points_radial_grid -1
do k = 1, n_points_integration_angular do k = 1, n_points_integration_angular
! if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
! cycle cycle
! endif endif
n_points_final_grid += 1 n_points_final_grid += 1
enddo enddo
enddo enddo
@ -39,9 +40,9 @@ END_PROVIDER
do j = 1, nucl_num do j = 1, nucl_num
do i = 1, n_points_radial_grid -1 do i = 1, n_points_radial_grid -1
do k = 1, n_points_integration_angular do k = 1, n_points_integration_angular
!if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
! cycle cycle
!endif endif
i_count += 1 i_count += 1
final_grid_points(1,i_count) = grid_points_per_atom(1,k,i,j) final_grid_points(1,i_count) = grid_points_per_atom(1,k,i,j)
final_grid_points(2,i_count) = grid_points_per_atom(2,k,i,j) final_grid_points(2,i_count) = grid_points_per_atom(2,k,i,j)

5
src/becke_numerical_grid/step_function_becke.irp.f
View File

@ -31,10 +31,6 @@ double precision function cell_function_becke(r,atom_number)
double precision :: mu_ij,nu_ij double precision :: mu_ij,nu_ij
double precision :: distance_i,distance_j,step_function_becke double precision :: distance_i,distance_j,step_function_becke
integer :: j integer :: j
if(int(nucl_charge(atom_number))==0)then
cell_function_becke = 0.d0
return
endif
distance_i = (r(1) - nucl_coord_transp(1,atom_number) ) * (r(1) - nucl_coord_transp(1,atom_number)) distance_i = (r(1) - nucl_coord_transp(1,atom_number) ) * (r(1) - nucl_coord_transp(1,atom_number))
distance_i += (r(2) - nucl_coord_transp(2,atom_number) ) * (r(2) - nucl_coord_transp(2,atom_number)) distance_i += (r(2) - nucl_coord_transp(2,atom_number) ) * (r(2) - nucl_coord_transp(2,atom_number))
distance_i += (r(3) - nucl_coord_transp(3,atom_number) ) * (r(3) - nucl_coord_transp(3,atom_number)) distance_i += (r(3) - nucl_coord_transp(3,atom_number) ) * (r(3) - nucl_coord_transp(3,atom_number))
@ -42,7 +38,6 @@ double precision function cell_function_becke(r,atom_number)
cell_function_becke = 1.d0 cell_function_becke = 1.d0
do j = 1, nucl_num do j = 1, nucl_num
if(j==atom_number)cycle if(j==atom_number)cycle
if(int(nucl_charge(j))==0)cycle
distance_j = (r(1) - nucl_coord_transp(1,j) ) * (r(1) - nucl_coord_transp(1,j)) distance_j = (r(1) - nucl_coord_transp(1,j) ) * (r(1) - nucl_coord_transp(1,j))
distance_j+= (r(2) - nucl_coord_transp(2,j) ) * (r(2) - nucl_coord_transp(2,j)) distance_j+= (r(2) - nucl_coord_transp(2,j) ) * (r(2) - nucl_coord_transp(2,j))
distance_j+= (r(3) - nucl_coord_transp(3,j) ) * (r(3) - nucl_coord_transp(3,j)) distance_j+= (r(3) - nucl_coord_transp(3,j) ) * (r(3) - nucl_coord_transp(3,j))

46
src/cipsi/cipsi.irp.f
View File

@ -5,7 +5,7 @@ subroutine run_cipsi
! stochastic PT2. ! stochastic PT2.
END_DOC END_DOC
integer :: i,j,k integer :: i,j,k
double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:) double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:), zeros(:)
integer :: n_det_before, to_select integer :: n_det_before, to_select
double precision :: rss double precision :: rss
@ -13,7 +13,7 @@ subroutine run_cipsi
rss = memory_of_double(N_states)*4.d0 rss = memory_of_double(N_states)*4.d0
call check_mem(rss,irp_here) call check_mem(rss,irp_here)
allocate (pt2(N_states), rpt2(N_states), norm(N_states), variance(N_states)) allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states))
double precision :: hf_energy_ref double precision :: hf_energy_ref
logical :: has logical :: has
@ -23,10 +23,11 @@ subroutine run_cipsi
relative_error=PT2_relative_error relative_error=PT2_relative_error
zeros = 0.d0
pt2 = -huge(1.e0) pt2 = -huge(1.e0)
rpt2 = -huge(1.e0) rpt2 = -huge(1.e0)
norm = 0.d0 norm = 0.d0
variance = 0.d0 variance = huge(1.e0)
if (s2_eig) then if (s2_eig) then
call make_s2_eigenfunction call make_s2_eigenfunction
@ -65,7 +66,8 @@ subroutine run_cipsi
do while ( & do while ( &
(N_det < N_det_max) .and. & (N_det < N_det_max) .and. &
(maxval(abs(pt2(1:N_states))) > pt2_max) .and. & (maxval(abs(rpt2(1:N_states))) > pt2_max) .and. &
(maxval(variance(1:N_states)) > variance_max) .and. &
(correlation_energy_ratio <= correlation_energy_ratio_max) & (correlation_energy_ratio <= correlation_energy_ratio_max) &
) )
write(*,'(A)') '--------------------------------------------------------------------------------' write(*,'(A)') '--------------------------------------------------------------------------------'
@ -83,17 +85,17 @@ subroutine run_cipsi
SOFT_TOUCH threshold_generators SOFT_TOUCH threshold_generators
endif endif
do k=1,N_states
rpt2(k) = pt2(k)/(1.d0 + norm(k))
enddo
correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / & correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / &
(psi_energy_with_nucl_rep(1) + pt2(1) - hf_energy_ref) (psi_energy_with_nucl_rep(1) + rpt2(1) - hf_energy_ref)
correlation_energy_ratio = min(1.d0,correlation_energy_ratio) correlation_energy_ratio = min(1.d0,correlation_energy_ratio)
call write_double(6,correlation_energy_ratio, 'Correlation ratio') call write_double(6,correlation_energy_ratio, 'Correlation ratio')
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
do k=1,N_states
rpt2(:) = pt2(:)/(1.d0 + norm(k))
enddo
call save_energy(psi_energy_with_nucl_rep, rpt2) call save_energy(psi_energy_with_nucl_rep, rpt2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)
@ -103,9 +105,8 @@ subroutine run_cipsi
if (qp_stop()) exit if (qp_stop()) exit
n_det_before = N_det n_det_before = N_det
to_select = N_det to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
to_select = max(N_states_diag, to_select) to_select = max(N_states_diag, to_select)
! to_select = min(to_select, N_det_max-n_det_before)
call ZMQ_selection(to_select, pt2, variance, norm) call ZMQ_selection(to_select, pt2, variance, norm)
PROVIDE psi_coef PROVIDE psi_coef
@ -114,32 +115,30 @@ subroutine run_cipsi
call diagonalize_CI call diagonalize_CI
call save_wavefunction call save_wavefunction
rpt2(:) = 0.d0 call save_energy(psi_energy_with_nucl_rep, zeros)
call save_energy(psi_energy_with_nucl_rep, rpt2)
if (qp_stop()) exit if (qp_stop()) exit
print *, (N_det < N_det_max)
print *, (maxval(abs(rpt2(1:N_states))) > pt2_max)
print *, (maxval(variance(1:N_states)) > variance_max)
print *, (correlation_energy_ratio <= correlation_energy_ratio_max)
enddo enddo
if (.not.qp_stop()) then if (.not.qp_stop()) then
if (N_det < N_det_max) then if (N_det < N_det_max) then
call diagonalize_CI call diagonalize_CI
call save_wavefunction call save_wavefunction
rpt2(:) = 0.d0 call save_energy(psi_energy_with_nucl_rep, zeros)
call save_energy(psi_energy_with_nucl_rep, rpt2)
endif endif
if (do_pt2) then if (do_pt2) then
pt2 = 0.d0 pt2(:) = 0.d0
variance = 0.d0 variance(:) = 0.d0
norm = 0.d0 norm(:) = 0.d0
threshold_generators = 1d0 threshold_generators = 1d0
SOFT_TOUCH threshold_generators SOFT_TOUCH threshold_generators
call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, & call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, &
norm,0) ! Stochastic PT2 norm,0) ! Stochastic PT2
SOFT_TOUCH threshold_generators SOFT_TOUCH threshold_generators
do k=1,N_states
rpt2(:) = pt2(:)/(1.d0 + norm(k))
enddo
call save_energy(psi_energy_with_nucl_rep, pt2)
endif endif
print *, 'N_det = ', N_det print *, 'N_det = ', N_det
print *, 'N_sop = ', N_occ_pattern print *, 'N_sop = ', N_occ_pattern
@ -148,12 +147,11 @@ subroutine run_cipsi
do k=1,N_states do k=1,N_states
rpt2(:) = pt2(:)/(1.d0 + norm(k)) rpt2(k) = pt2(k)/(1.d0 + norm(k))
enddo enddo
call save_energy(psi_energy_with_nucl_rep, rpt2)
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
call save_energy(psi_energy_with_nucl_rep, pt2) call save_energy(psi_energy_with_nucl_rep, rpt2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)
call print_extrapolated_energy() call print_extrapolated_energy()
endif endif

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

@ -129,13 +129,13 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted
PROVIDE psi_det_hii N_generators_bitmask PROVIDE psi_det_hii N_generators_bitmask selection_weight pseudo_sym
if (h0_type == 'SOP') then if (h0_type == 'SOP') then
PROVIDE psi_occ_pattern_hii det_to_occ_pattern PROVIDE psi_occ_pattern_hii det_to_occ_pattern
endif endif
if (N_det < max(1000,N_states)) then if (N_det < max(4,N_states)) then
pt2=0.d0 pt2=0.d0
variance=0.d0 variance=0.d0
norm=0.d0 norm=0.d0
@ -182,6 +182,9 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then
stop 'Unable to put state_average_weight on ZMQ server' stop 'Unable to put state_average_weight on ZMQ server'
endif endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) then
stop 'Unable to put selection_weight on ZMQ server'
endif
if (zmq_put_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) then if (zmq_put_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) then
stop 'Unable to put pt2_stoch_istate on ZMQ server' stop 'Unable to put pt2_stoch_istate on ZMQ server'
endif endif
@ -333,13 +336,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
pt2(k) = 0.d0 pt2(k) = 0.d0
enddo enddo
! Adjust PT2 weights for next selection call update_pt2_and_variance_weights(pt2, variance, norm, N_states)
double precision :: pt2_avg
pt2_avg = sum(pt2) / dble(N_states)
do k=1,N_states
pt2_match_weight(k) *= (pt2(k)/pt2_avg)**2
enddo
SOFT_TOUCH pt2_match_weight
end subroutine end subroutine

6
src/cipsi/run_selection_slave.irp.f
View File

@ -25,8 +25,8 @@ subroutine run_selection_slave(thread,iproc,energy)
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F pseudo_sym
PROVIDE psi_selectors_coef_transp psi_det_sorted PROVIDE psi_selectors_coef_transp psi_det_sorted weight_selection
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
@ -230,6 +230,8 @@ subroutine pull_selection_results(zmq_socket_pull, pt2, variance, norm, val, det
endif endif
else else
pt2(:) = 0.d0 pt2(:) = 0.d0
variance(:) = 0.d0
norm(:) = 0.d0
endif endif
rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0) rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)

105
src/cipsi/selection.irp.f
View File

@ -6,15 +6,108 @@ BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ]
! Weights adjusted along the selection to make the PT2 contributions ! Weights adjusted along the selection to make the PT2 contributions
! of each state coincide. ! of each state coincide.
END_DOC END_DOC
pt2_match_weight = 1.d0 pt2_match_weight(:) = 1.d0
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ]
implicit none
BEGIN_DOC
! Weights adjusted along the selection to make the variances
! of each state coincide.
END_DOC
variance_match_weight(:) = 1.d0
END_PROVIDER
subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st)
implicit none
BEGIN_DOC
! Updates the rPT2- and Variance- matching weights.
END_DOC
integer, intent(in) :: N_st
double precision, intent(in) :: pt2(N_st)
double precision, intent(in) :: variance(N_st)
double precision, intent(in) :: norm(N_st)
double precision :: avg, rpt2(N_st), element, dt, x
integer :: k
integer, save :: i_iter=0
integer, parameter :: i_itermax = 3
double precision, allocatable, save :: memo_variance(:,:), memo_pt2(:,:)
if (i_iter == 0) then
allocate(memo_variance(N_st,i_itermax), memo_pt2(N_st,i_itermax))
memo_pt2(:,:) = 1.d0
memo_variance(:,:) = 1.d0
endif
i_iter = i_iter+1
if (i_iter > i_itermax) then
i_iter = 1
endif
dt = 4.d0
do k=1,N_st
rpt2(k) = pt2(k)/(1.d0 + norm(k))
enddo
avg = sum(rpt2(1:N_st)) / dble(N_st)
do k=1,N_st
element = exp(dt*(rpt2(k)/avg -1.d0))
element = min(1.5d0 , element)
element = max(0.5d0 , element)
memo_pt2(k,i_iter) = element
pt2_match_weight(k) = product(memo_pt2(k,:))
enddo
avg = sum(variance(1:N_st)) / dble(N_st)
do k=1,N_st
element = exp(dt*(variance(k)/avg -1.d0))
element = min(1.5d0 , element)
element = max(0.5d0 , element)
memo_variance(k,i_iter) = element
variance_match_weight(k) = product(memo_variance(k,:))
enddo
print *, '# PT2 weight ', real(pt2_match_weight(:),4)
print *, '# var weight ', real(variance_match_weight(:),4)
SOFT_TOUCH pt2_match_weight variance_match_weight
end
BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ]
implicit none implicit none
BEGIN_DOC BEGIN_DOC
! Weights used in the selection criterion ! Weights used in the selection criterion
END_DOC END_DOC
select case (weight_selection)
case (0)
print *, 'Using input weights in selection'
selection_weight(1:N_states) = state_average_weight(1:N_states)
case (1)
print *, 'Using 1/c_max^2 weight in selection'
selection_weight(1:N_states) = c0_weight(1:N_states)
case (2)
print *, 'Using pt2-matching weight in selection'
selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states)
case (3)
print *, 'Using variance-matching weight in selection'
selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states)
case (4)
print *, 'Using variance- and pt2-matching weights in selection'
selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)
case (5)
print *, 'Using variance-matching weight in selection'
selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states)
end select
END_PROVIDER END_PROVIDER
@ -621,11 +714,13 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi
norm(istate) = norm(istate) + coef * coef norm(istate) = norm(istate) + coef * coef
! if (h0_type == "Variance") then if (weight_selection /= 5) then
! sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate) ! Energy selection
! else
sum_e_pert = sum_e_pert + e_pert * selection_weight(istate) sum_e_pert = sum_e_pert + e_pert * selection_weight(istate)
! endif else
! Variance selection
sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
endif
end do end do
if(pseudo_sym)then if(pseudo_sym)then
if(dabs(mat(1, p1, p2)).lt.thresh_sym)then if(dabs(mat(1, p1, p2)).lt.thresh_sym)then

8
src/cipsi/slave_cipsi.irp.f
View File

@ -17,7 +17,7 @@ subroutine provide_everything
PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag
PROVIDE pt2_e0_denominator mo_num N_int ci_energy mpi_master zmq_state zmq_context PROVIDE pt2_e0_denominator mo_num N_int ci_energy mpi_master zmq_state zmq_context
PROVIDE psi_det psi_coef threshold_generators state_average_weight PROVIDE psi_det psi_coef threshold_generators state_average_weight
PROVIDE N_det_selectors pt2_stoch_istate N_det PROVIDE N_det_selectors pt2_stoch_istate N_det selection_weight pseudo_sym
end end
subroutine run_slave_main subroutine run_slave_main
@ -220,8 +220,12 @@ subroutine run_slave_main
call mpi_print('zmq_get_dvector state_average_weight') call mpi_print('zmq_get_dvector state_average_weight')
IRP_ENDIF IRP_ENDIF
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle
IRP_IF MPI_DEBUG
call mpi_print('zmq_get_dvector selection_weight')
IRP_ENDIF
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) cycle
pt2_e0_denominator(1:N_states) = energy(1:N_states) pt2_e0_denominator(1:N_states) = energy(1:N_states)
SOFT_TOUCH pt2_e0_denominator state_average_weight pt2_stoch_istate threshold_generators SOFT_TOUCH pt2_e0_denominator state_average_weight pt2_stoch_istate threshold_generators selection_weight
call wall_time(t1) call wall_time(t1)
call write_double(6,(t1-t0),'Broadcast time') call write_double(6,(t1-t0),'Broadcast time')

34
src/cipsi/stochastic_cipsi.irp.f
View File

@ -4,7 +4,7 @@ subroutine run_stochastic_cipsi
! Selected Full Configuration Interaction with Stochastic selection and PT2. ! Selected Full Configuration Interaction with Stochastic selection and PT2.
END_DOC END_DOC
integer :: i,j,k integer :: i,j,k
double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:) double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:), zeros(:)
integer :: to_select integer :: to_select
logical, external :: qp_stop logical, external :: qp_stop
@ -18,7 +18,7 @@ subroutine run_stochastic_cipsi
rss = memory_of_double(N_states)*4.d0 rss = memory_of_double(N_states)*4.d0
call check_mem(rss,irp_here) call check_mem(rss,irp_here)
allocate (pt2(N_states), rpt2(N_states), norm(N_states), variance(N_states)) allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states))
double precision :: hf_energy_ref double precision :: hf_energy_ref
logical :: has logical :: has
@ -26,6 +26,7 @@ subroutine run_stochastic_cipsi
relative_error=PT2_relative_error relative_error=PT2_relative_error
zeros = 0.d0
pt2 = -huge(1.e0) pt2 = -huge(1.e0)
rpt2 = -huge(1.e0) rpt2 = -huge(1.e0)
norm = 0.d0 norm = 0.d0
@ -63,14 +64,14 @@ subroutine run_stochastic_cipsi
do while ( & do while ( &
(N_det < N_det_max) .and. & (N_det < N_det_max) .and. &
(maxval(abs(pt2(1:N_states))) > pt2_max) .and. & (maxval(abs(rpt2(1:N_states))) > pt2_max) .and. &
(maxval(abs(variance(1:N_states))) > variance_max) .and. & (maxval(abs(variance(1:N_states))) > variance_max) .and. &
(correlation_energy_ratio <= correlation_energy_ratio_max) & (correlation_energy_ratio <= correlation_energy_ratio_max) &
) )
write(*,'(A)') '--------------------------------------------------------------------------------' write(*,'(A)') '--------------------------------------------------------------------------------'
to_select = N_det*int(sqrt(dble(N_states))) to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
to_select = max(N_states_diag, to_select) to_select = max(N_states_diag, to_select)
pt2 = 0.d0 pt2 = 0.d0
@ -79,17 +80,17 @@ subroutine run_stochastic_cipsi
call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, & call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, &
norm, to_select) ! Stochastic PT2 and selection norm, to_select) ! Stochastic PT2 and selection
do k=1,N_states
rpt2(k) = pt2(k)/(1.d0 + norm(k))
enddo
correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / & correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / &
(psi_energy_with_nucl_rep(1) + pt2(1) - hf_energy_ref) (psi_energy_with_nucl_rep(1) + rpt2(1) - hf_energy_ref)
correlation_energy_ratio = min(1.d0,correlation_energy_ratio) correlation_energy_ratio = min(1.d0,correlation_energy_ratio)
call save_energy(psi_energy_with_nucl_rep, rpt2)
call write_double(6,correlation_energy_ratio, 'Correlation ratio') call write_double(6,correlation_energy_ratio, 'Correlation ratio')
call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
do k=1,N_states
rpt2(:) = pt2(:)/(1.d0 + norm(k))
enddo
call save_energy(psi_energy_with_nucl_rep, rpt2) call save_energy(psi_energy_with_nucl_rep, rpt2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)
@ -108,8 +109,7 @@ subroutine run_stochastic_cipsi
call diagonalize_CI call diagonalize_CI
call save_wavefunction call save_wavefunction
rpt2(:) = 0.d0 call save_energy(psi_energy_with_nucl_rep, zeros)
call save_energy(psi_energy_with_nucl_rep, rpt2)
if (qp_stop()) exit if (qp_stop()) exit
enddo enddo
@ -117,20 +117,18 @@ subroutine run_stochastic_cipsi
if (N_det < N_det_max) then if (N_det < N_det_max) then
call diagonalize_CI call diagonalize_CI
call save_wavefunction call save_wavefunction
rpt2(:) = 0.d0 call save_energy(psi_energy_with_nucl_rep, zeros)
call save_energy(psi_energy_with_nucl_rep, rpt2)
endif endif
pt2 = 0.d0 pt2(:) = 0.d0
variance = 0.d0 variance(:) = 0.d0
norm = 0.d0 norm(:) = 0.d0
call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, & call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, &
norm,0) ! Stochastic PT2 norm,0) ! Stochastic PT2
do k=1,N_states do k=1,N_states
rpt2(:) = pt2(:)/(1.d0 + norm(k)) rpt2(k) = pt2(k)/(1.d0 + norm(k))
enddo enddo
call save_energy(psi_energy_with_nucl_rep, rpt2)
call save_energy(psi_energy_with_nucl_rep, rpt2) call save_energy(psi_energy_with_nucl_rep, rpt2)
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)

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

@ -21,7 +21,8 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym
call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection') call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection')
@ -45,6 +46,9 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then
stop 'Unable to put state_average_weight on ZMQ server' stop 'Unable to put state_average_weight on ZMQ server'
endif endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) then
stop 'Unable to put selection_weight on ZMQ server'
endif
if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) then if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) then
stop 'Unable to put threshold_generators on ZMQ server' stop 'Unable to put threshold_generators on ZMQ server'
endif endif
@ -85,7 +89,11 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
endif endif
integer :: nproc_target integer :: nproc_target
if (N_det < 3*nproc) then
nproc_target = N_det/4
else
nproc_target = nproc nproc_target = nproc
endif
double precision :: mem double precision :: mem
mem = 8.d0 * N_det * (N_int * 2.d0 * 3.d0 + 3.d0 + 5.d0) / (1024.d0**3) mem = 8.d0 * N_det * (N_int * 2.d0 * 3.d0 + 3.d0 + 5.d0) / (1024.d0**3)
call write_double(6,mem,'Estimated memory/thread (Gb)') call write_double(6,mem,'Estimated memory/thread (Gb)')
@ -131,13 +139,7 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
norm(k) = norm(k) * f(k) norm(k) = norm(k) * f(k)
enddo enddo
! Adjust PT2 weights for next selection call update_pt2_and_variance_weights(pt2, variance, norm, N_states)
double precision :: pt2_avg
pt2_avg = sum(pt2) / dble(N_states)
do k=1,N_states
pt2_match_weight(k) *= (pt2(k)/pt2_avg)**2
enddo
SOFT_TOUCH pt2_match_weight
end subroutine end subroutine
@ -159,9 +161,9 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2, variance, norm)
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
type(selection_buffer), intent(inout) :: b type(selection_buffer), intent(inout) :: b
integer, intent(in) :: N integer, intent(in) :: N
double precision, intent(inout) :: pt2(N_states) double precision, intent(out) :: pt2(N_states)
double precision, intent(inout) :: variance(N_states) double precision, intent(out) :: variance(N_states)
double precision, intent(inout) :: norm(N_states) double precision, intent(out) :: norm(N_states)
double precision :: pt2_mwen(N_states) double precision :: pt2_mwen(N_states)
double precision :: variance_mwen(N_states) double precision :: variance_mwen(N_states)
double precision :: norm_mwen(N_states) double precision :: norm_mwen(N_states)

54
src/davidson/print_e_components.irp.f Normal file
View File

@ -0,0 +1,54 @@
subroutine print_energy_components()
implicit none
BEGIN_DOC
! Prints the different components of the energy.
END_DOC
integer, save :: ifirst = 0
double precision :: Vee, Ven, Vnn, Vecp, T, f
integer :: i,j,k
Vnn = nuclear_repulsion
print *, 'Energy components'
print *, '================='
print *, ''
do k=1,N_states
Ven = 0.d0
Vecp = 0.d0
T = 0.d0
do j=1,mo_num
do i=1,mo_num
f = one_e_dm_mo_alpha(i,j,k) + one_e_dm_mo_beta(i,j,k)
Ven = Ven + f * mo_integrals_n_e(i,j)
Vecp = Vecp + f * mo_pseudo_integrals(i,j)
T = T + f * mo_kinetic_integrals(i,j)
enddo
enddo
Vee = psi_energy(k) - Ven - Vecp - T
if (ifirst == 0) then
ifirst = 1
print *, 'Vnn : Nucleus-Nucleus potential energy'
print *, 'Ven : Electron-Nucleus potential energy'
print *, 'Vee : Electron-Electron potential energy'
print *, 'Vecp : Potential energy of the pseudo-potentials'
print *, 'T : Electronic kinetic energy'
print *, ''
endif
print *, 'State ', k
print *, '---------'
print *, ''
print *, 'Vnn = ', Vnn
print *, 'Ven = ', Ven
print *, 'Vee = ', Vee
print *, 'Vecp = ', Vecp
print *, 'T = ', T
print *, ''
enddo
print *, ''
end

13
src/determinants/EZFIO.cfg
View File

@ -28,12 +28,18 @@ doc: Force the wave function to be an eigenfunction of |S^2|
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: True default: True
[used_weight] [weight_one_e_dm]
type: integer type: integer
doc: Weight used in the calculation of the one-electron density matrix. 0: 1./(c_0^2), 1: 1/N_states, 2: input state-average weight, 3: 1/(Norm_L3(Psi)) doc: Weight used in the calculation of the one-electron density matrix. 0: 1./(c_0^2), 1: 1/N_states, 2: input state-average weight, 3: 1/(Norm_L3(Psi))
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: 1 default: 1
[weight_selection]
type: integer
doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: rPT2 matching, 3: variance matching, 4: variance and rPT2 matching, 5: variance minimization and matching
interface: ezfio,provider,ocaml
default: 2
[threshold_generators] [threshold_generators]
type: Threshold type: Threshold
doc: Thresholds on generators (fraction of the square of the norm) doc: Thresholds on generators (fraction of the square of the norm)
@ -89,6 +95,11 @@ doc: Weight of the states in state-average calculations.
interface: ezfio interface: ezfio
size: (determinants.n_states) size: (determinants.n_states)
[selection_factor]
type: double precision
doc: f such that the number of determinants to add is f * N_det * sqrt(N_states)
interface: ezfio,provider,ocaml
default: 1.
[thresh_sym] [thresh_sym]
type: Threshold type: Threshold

4
src/determinants/density_matrix.irp.f
View File

@ -305,9 +305,9 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
logical :: exists logical :: exists
state_average_weight(:) = 1.d0 state_average_weight(:) = 1.d0
if (used_weight == 0) then if (weight_one_e_dm == 0) then
state_average_weight(:) = c0_weight(:) state_average_weight(:) = c0_weight(:)
else if (used_weight == 1) then else if (weight_one_e_dm == 1) then
state_average_weight(:) = 1./N_states state_average_weight(:) = 1./N_states
else else
call ezfio_has_determinants_state_average_weight(exists) call ezfio_has_determinants_state_average_weight(exists)

1
src/determinants/utils.irp.f
View File

@ -43,4 +43,3 @@ BEGIN_PROVIDER [ double precision, S2_matrix_all_dets,(N_det,N_det) ]
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
END_PROVIDER END_PROVIDER

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

@ -121,3 +121,26 @@
enddo enddo
END_PROVIDER END_PROVIDER
BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom, (ao_num,n_pts_max_per_atom,nucl_num)]
&BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom_transp, (n_pts_max_per_atom,ao_num,nucl_num)]
implicit none
BEGIN_DOC
! aos_in_r_array_per_atom(i,j,k) = value of the ith ao on the jth grid point attached on the kth atom
END_DOC
integer :: i,j,k
double precision :: aos_array(ao_num), r(3)
do k = 1, nucl_num
do i = 1, n_pts_per_atom(k)
r(1) = final_grid_points_per_atom(1,i,k)
r(2) = final_grid_points_per_atom(2,i,k)
r(3) = final_grid_points_per_atom(3,i,k)
call give_all_aos_at_r(r,aos_array)
do j = 1, ao_num
aos_in_r_array_per_atom(j,i,k) = aos_array(j)
aos_in_r_array_per_atom_transp(i,j,k) = aos_array(j)
enddo
enddo
enddo
END_PROVIDER

2
src/fci/pt2.irp.f
View File

@ -46,7 +46,7 @@ subroutine run
call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, & call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, &
norm,0) ! Stochastic PT2 norm,0) ! Stochastic PT2
do k=1,N_states do k=1,N_states
rpt2(:) = pt2(:)/(1.d0 + norm(k)) rpt2(k) = pt2(k)/(1.d0 + norm(k))
enddo enddo
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)

7
src/iterations/print_summary.irp.f
View File

@ -31,11 +31,11 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt) write(*,fmt)
write(fmt,*) '(12X,', N_states_p, '(6X,A7,1X,I6,10X))' write(fmt,*) '(13X,', N_states_p, '(6X,A7,1X,I6,10X))'
write(*,fmt) ('State',k, k=1,N_states_p) write(*,fmt) ('State',k, k=1,N_states_p)
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt) write(*,fmt)
write(fmt,*) '(A12,', N_states_p, '(1X,F14.8,15X))' write(fmt,*) '(A13,', N_states_p, '(1X,F14.8,15X))'
write(*,fmt) '# E ', e_(1:N_states_p) write(*,fmt) '# E ', e_(1:N_states_p)
if (N_states_p > 1) then if (N_states_p > 1) then
write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1) write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1)
@ -43,6 +43,7 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_
endif endif
write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))' write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))'
write(*,fmt) '# PT2 '//pt2_string, (pt2_(k), error_(k), k=1,N_states_p) write(*,fmt) '# PT2 '//pt2_string, (pt2_(k), error_(k), k=1,N_states_p)
write(*,fmt) '# rPT2'//pt2_string, (pt2_(k)*f(k), error_(k)*f(k), k=1,N_states_p)
write(*,'(A)') '#' write(*,'(A)') '#'
write(*,fmt) '# E+PT2 ', (e_(k)+pt2_(k),error_(k), k=1,N_states_p) write(*,fmt) '# E+PT2 ', (e_(k)+pt2_(k),error_(k), k=1,N_states_p)
write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_(k)*f(k),error_(k)*f(k), k=1,N_states_p) write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_(k)*f(k),error_(k)*f(k), k=1,N_states_p)
@ -97,5 +98,7 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_
enddo enddo
endif endif
call print_energy_components()
end subroutine end subroutine

4
src/nuclei/atomic_radii.irp.f
View File

@ -66,7 +66,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom, (nucl_num)]
implicit none implicit none
integer :: i integer :: i
do i = 1, nucl_num do i = 1, nucl_num
slater_bragg_radii_per_atom(i) = slater_bragg_radii(int(nucl_charge(i))) slater_bragg_radii_per_atom(i) = slater_bragg_radii(max(1,int(nucl_charge(i))))
enddo enddo
END_PROVIDER END_PROVIDER
@ -74,7 +74,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom_ua, (nucl_num)]
implicit none implicit none
integer :: i integer :: i
do i = 1, nucl_num do i = 1, nucl_num
slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(int(nucl_charge(i))) slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(max(1,int(nucl_charge(i))))
enddo enddo
END_PROVIDER END_PROVIDER