diff --git a/CITATION.cff b/CITATION.cff
index 9cd186d2..a391b8fa 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -1,32 +1,100 @@
# YAML 1.2
# Metadata for citation of this software according to the CFF format (https://citation-file-format.github.io/)
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
-doi: 10.5281/zenodo.825872
+doi: 10.1021/acs.jctc.9b00176
authors:
-- given-names: Anthony
- family-names: Scemama
- affiliation: Laboratoire de Chimie et Physique Quantiques / CNRS
- given-names: Yann
family-names: Garniron
- affiliation: Laboratoire de Chimie et Physique Quantiques / CNRS
-- given-names: Michel
- family-names: Caffarel
- 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: Thomas
family-names: Applencourt
- affiliation: Argonne National Lab
+ affiliation: Computational Science Division, Argonne National Laboratory, Argonne, Illinois 60439, United States
- given-names: Kevin
family-names: Gasperich
- affiliation: Argonne National Lab
+ affiliation: Department of Chemistry, University of Pittsburgh, Pittsburgh, Pennsylvania 15260, United States
- given-names: Anouar
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
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'
-date-released: 2019-02-11
+url: https://quantumpackage.github.io/qp2/
+date-released: 2019-05-13
repository-code: https://github.com/QuantumPackage/qp2
+keywords: [ "computational chemistry", "configuration interaction", "cipsi", "perturbation theory" ]
license: AGPL-3.0-or-later
diff --git a/INSTALL.rst b/INSTALL.rst
index cc5b512a..baf7f051 100644
--- a/INSTALL.rst
+++ b/INSTALL.rst
@@ -45,6 +45,8 @@ Requirements
- |ZeroMQ| : networking library
- `GMP `_ : Gnu Multiple Precision Arithmetic Library
- |OCaml| compiler with |OPAM| package manager
+- `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
@@ -86,6 +88,8 @@ The following packages are supported by the :command:`configure` installer:
* zeromq
* f77zmq
* gmp
+* libcap
+* bwrap
* ocaml ( :math:`\approx` 10 minutes)
* ezfio
* docopt
@@ -243,6 +247,55 @@ With Debian or Ubuntu, you can use
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:
+ ``_
+ 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:
+ ``_
+ 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
-----
diff --git a/README.md b/README.md
index c94adb69..3e85b4cc 100644
--- a/README.md
+++ b/README.md
@@ -1,12 +1,13 @@
# 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\
+J. Chem. Theory Comput. (2019)\
https://arxiv.org/abs/1902.08154
-![QP](https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png)
+
# Getting started
diff --git a/config/ifort_avx.cfg b/config/ifort_avx.cfg
index d3fcd1f0..7a13348c 100644
--- a/config/ifort_avx.cfg
+++ b/config/ifort_avx.cfg
@@ -32,7 +32,7 @@ OPENMP : 1 ; Append OpenMP flags
#
[OPT]
FC : -traceback
-FCFLAGS : -xAVX -O2 -ip -ftz -g
+FCFLAGS : -march=corei7-avx -O2 -ip -ftz -g
# Profiling flags
#################
diff --git a/config/ifort_avx_mpi.cfg b/config/ifort_avx_mpi.cfg
index a6784058..550dbf91 100644
--- a/config/ifort_avx_mpi.cfg
+++ b/config/ifort_avx_mpi.cfg
@@ -31,14 +31,14 @@ OPENMP : 1 ; Append OpenMP flags
# -ftz : Flushes denormal results to zero
#
[OPT]
-FCFLAGS : -xAVX -O2 -ip -ftz -g -traceback
+FCFLAGS : -march=corei7-avx -O2 -ip -ftz -g -traceback
# Profiling flags
#################
#
[PROFILE]
FC : -p -g
-FCFLAGS : -xSSE4.2 -O2 -ip -ftz
+FCFLAGS : -march=corei7 -O2 -ip -ftz
# Debugging flags
diff --git a/configure b/configure
index ffda5016..c343c243 100755
--- a/configure
+++ b/configure
@@ -175,7 +175,7 @@ if [[ "${PACKAGES}.x" != ".x" ]] ; then
fi
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
@@ -206,6 +206,32 @@ EOF
make install
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
@@ -276,7 +302,7 @@ EOF
rm ${QP_ROOT}/external/opam_installer.sh
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)
opam install -y ${OCAML_PACKAGES} || exit 1
@@ -290,7 +316,7 @@ EOF
| sh \${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
- \${QP_ROOT}/bin/opam init --disable-sandboxing --verbose --yes
+ \${QP_ROOT}/bin/opam init --verbose --yes
eval \$(\${QP_ROOT}/bin/opam env)
opam install -y \${OCAML_PACKAGES} || exit 1
EOF
@@ -399,6 +425,18 @@ if [[ ${ZLIB} = $(not_found) ]] ; then
fail
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)
if [[ ${OPAM} = $(not_found) ]] ; then
error "OPAM (ocaml) package manager is not installed."
diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml
index ab6fb2ca..e4c6ff2a 100644
--- a/ocaml/Input_determinants_by_hand.ml
+++ b/ocaml/Input_determinants_by_hand.ml
@@ -81,9 +81,6 @@ end = struct
;;
let write_n_det n =
- let n_det_old =
- Ezfio.get_determinants_n_det ()
- in
Det_number.to_int n
|> Ezfio.set_determinants_n_det
;;
diff --git a/src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f b/src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f
index 4433690f..97debfab 100644
--- a/src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f
+++ b/src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f
@@ -15,6 +15,8 @@ double precision function ao_two_e_integral_erf(i,j,k,l)
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
integer :: iorder_p(3), iorder_q(3)
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
ao_two_e_integral_erf = ao_two_e_integral_schwartz_accel_erf(i,j,k,l)
diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f
index 304fec49..9e729cd4 100644
--- a/src/ao_two_e_ints/map_integrals.irp.f
+++ b/src/ao_two_e_ints/map_integrals.irp.f
@@ -279,6 +279,100 @@ subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_z
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()
implicit none
integer (map_size_kind) :: get_ao_map_size
diff --git a/src/becke_numerical_grid/EZFIO.cfg b/src/becke_numerical_grid/EZFIO.cfg
index ed89428c..ca2100a1 100644
--- a/src/becke_numerical_grid/EZFIO.cfg
+++ b/src/becke_numerical_grid/EZFIO.cfg
@@ -8,3 +8,9 @@ default: 2
type: integer
doc: Total number of grid points
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
diff --git a/src/becke_numerical_grid/atomic_number.irp.f b/src/becke_numerical_grid/atomic_number.irp.f
new file mode 100644
index 00000000..eea1fad7
--- /dev/null
+++ b/src/becke_numerical_grid/atomic_number.irp.f
@@ -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
+
diff --git a/src/becke_numerical_grid/grid_becke.irp.f b/src/becke_numerical_grid/grid_becke.irp.f
index 38d4053f..e72f6460 100644
--- a/src/becke_numerical_grid/grid_becke.irp.f
+++ b/src/becke_numerical_grid/grid_becke.irp.f
@@ -146,7 +146,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_
x = grid_points_radial(j)
! 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
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
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
- contrib_integration = derivative_knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)&
- *knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)**2
+ contrib_integration = derivative_knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)&
+ *knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)**2
final_weight_at_r(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
print*,'isnan(final_weight_at_r(k,i,j))'
diff --git a/src/becke_numerical_grid/grid_becke_per_atom.irp.f b/src/becke_numerical_grid/grid_becke_per_atom.irp.f
new file mode 100644
index 00000000..6026350b
--- /dev/null
+++ b/src/becke_numerical_grid/grid_becke_per_atom.irp.f
@@ -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
diff --git a/src/becke_numerical_grid/grid_becke_vector.irp.f b/src/becke_numerical_grid/grid_becke_vector.irp.f
index a595cd0b..3c2a6b91 100644
--- a/src/becke_numerical_grid/grid_becke_vector.irp.f
+++ b/src/becke_numerical_grid/grid_becke_vector.irp.f
@@ -1,5 +1,6 @@
BEGIN_PROVIDER [integer, n_points_final_grid]
+ implicit none
BEGIN_DOC
! Number of points which are non zero
END_DOC
@@ -8,9 +9,9 @@ BEGIN_PROVIDER [integer, n_points_final_grid]
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)) < 1.d-30)then
-! cycle
-! endif
+ if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+ cycle
+ endif
n_points_final_grid += 1
enddo
enddo
@@ -39,9 +40,9 @@ END_PROVIDER
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)) < 1.d-30)then
- ! cycle
- !endif
+ if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+ cycle
+ endif
i_count += 1
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)
diff --git a/src/becke_numerical_grid/step_function_becke.irp.f b/src/becke_numerical_grid/step_function_becke.irp.f
index b6335c3d..2905c6c0 100644
--- a/src/becke_numerical_grid/step_function_becke.irp.f
+++ b/src/becke_numerical_grid/step_function_becke.irp.f
@@ -31,10 +31,6 @@ double precision function cell_function_becke(r,atom_number)
double precision :: mu_ij,nu_ij
double precision :: distance_i,distance_j,step_function_becke
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(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))
@@ -42,7 +38,6 @@ double precision function cell_function_becke(r,atom_number)
cell_function_becke = 1.d0
do j = 1, nucl_num
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(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))
diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f
index f0cab384..7e292d6e 100644
--- a/src/cipsi/cipsi.irp.f
+++ b/src/cipsi/cipsi.irp.f
@@ -5,7 +5,7 @@ subroutine run_cipsi
! stochastic PT2.
END_DOC
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
double precision :: rss
@@ -13,7 +13,7 @@ subroutine run_cipsi
rss = memory_of_double(N_states)*4.d0
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
logical :: has
@@ -23,10 +23,11 @@ subroutine run_cipsi
relative_error=PT2_relative_error
+ zeros = 0.d0
pt2 = -huge(1.e0)
rpt2 = -huge(1.e0)
norm = 0.d0
- variance = 0.d0
+ variance = huge(1.e0)
if (s2_eig) then
call make_s2_eigenfunction
@@ -65,7 +66,8 @@ subroutine run_cipsi
do while ( &
(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) &
)
write(*,'(A)') '--------------------------------------------------------------------------------'
@@ -83,17 +85,17 @@ subroutine run_cipsi
SOFT_TOUCH threshold_generators
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) / &
- (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)
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)
- do k=1,N_states
- rpt2(:) = pt2(:)/(1.d0 + norm(k))
- enddo
call save_energy(psi_energy_with_nucl_rep, rpt2)
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
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 = min(to_select, N_det_max-n_det_before)
call ZMQ_selection(to_select, pt2, variance, norm)
PROVIDE psi_coef
@@ -114,32 +115,30 @@ subroutine run_cipsi
call diagonalize_CI
call save_wavefunction
- rpt2(:) = 0.d0
- call save_energy(psi_energy_with_nucl_rep, rpt2)
+ call save_energy(psi_energy_with_nucl_rep, zeros)
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
if (.not.qp_stop()) then
if (N_det < N_det_max) then
call diagonalize_CI
call save_wavefunction
- rpt2(:) = 0.d0
- call save_energy(psi_energy_with_nucl_rep, rpt2)
+ call save_energy(psi_energy_with_nucl_rep, zeros)
endif
if (do_pt2) then
- pt2 = 0.d0
- variance = 0.d0
- norm = 0.d0
+ pt2(:) = 0.d0
+ variance(:) = 0.d0
+ norm(:) = 0.d0
threshold_generators = 1d0
SOFT_TOUCH threshold_generators
call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, &
norm,0) ! Stochastic PT2
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
print *, 'N_det = ', N_det
print *, 'N_sop = ', N_occ_pattern
@@ -148,12 +147,11 @@ subroutine run_cipsi
do k=1,N_states
- rpt2(:) = pt2(:)/(1.d0 + norm(k))
+ rpt2(k) = pt2(k)/(1.d0 + norm(k))
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 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 print_extrapolated_energy()
endif
diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f
index 0e143a9c..9f891320 100644
--- a/src/cipsi/pt2_stoch_routines.irp.f
+++ b/src/cipsi/pt2_stoch_routines.irp.f
@@ -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_transp_rows_loc psi_bilinear_matrix_transp_columns
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
PROVIDE psi_occ_pattern_hii det_to_occ_pattern
endif
- if (N_det < max(1000,N_states)) then
+ if (N_det < max(4,N_states)) then
pt2=0.d0
variance=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
stop 'Unable to put state_average_weight on ZMQ server'
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
stop 'Unable to put pt2_stoch_istate on ZMQ server'
endif
@@ -333,13 +336,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
pt2(k) = 0.d0
enddo
- ! Adjust PT2 weights for next selection
- 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
+ call update_pt2_and_variance_weights(pt2, variance, norm, N_states)
end subroutine
diff --git a/src/cipsi/run_selection_slave.irp.f b/src/cipsi/run_selection_slave.irp.f
index 480ef12b..c1542445 100644
--- a/src/cipsi/run_selection_slave.irp.f
+++ b/src/cipsi/run_selection_slave.irp.f
@@ -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_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_order N_int pt2_F
- PROVIDE psi_selectors_coef_transp psi_det_sorted
+ PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F pseudo_sym
+ PROVIDE psi_selectors_coef_transp psi_det_sorted weight_selection
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
else
pt2(:) = 0.d0
+ variance(:) = 0.d0
+ norm(:) = 0.d0
endif
rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)
diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f
index eedcd67f..df31bc39 100644
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@@ -6,15 +6,108 @@ BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ]
! Weights adjusted along the selection to make the PT2 contributions
! of each state coincide.
END_DOC
- pt2_match_weight = 1.d0
+ pt2_match_weight(:) = 1.d0
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) ]
implicit none
BEGIN_DOC
! Weights used in the selection criterion
END_DOC
- selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states)
+ 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)
+
+ 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
@@ -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
norm(istate) = norm(istate) + coef * coef
-! if (h0_type == "Variance") then
-! sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
-! else
+ if (weight_selection /= 5) then
+ ! Energy selection
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
if(pseudo_sym)then
if(dabs(mat(1, p1, p2)).lt.thresh_sym)then
diff --git a/src/cipsi/slave_cipsi.irp.f b/src/cipsi/slave_cipsi.irp.f
index 58b53f94..c9710c18 100644
--- a/src/cipsi/slave_cipsi.irp.f
+++ b/src/cipsi/slave_cipsi.irp.f
@@ -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 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 N_det_selectors pt2_stoch_istate N_det
+ PROVIDE N_det_selectors pt2_stoch_istate N_det selection_weight pseudo_sym
end
subroutine run_slave_main
@@ -220,8 +220,12 @@ subroutine run_slave_main
call mpi_print('zmq_get_dvector state_average_weight')
IRP_ENDIF
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)
- 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 write_double(6,(t1-t0),'Broadcast time')
diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f
index 157479d9..ae2b7519 100644
--- a/src/cipsi/stochastic_cipsi.irp.f
+++ b/src/cipsi/stochastic_cipsi.irp.f
@@ -4,7 +4,7 @@ subroutine run_stochastic_cipsi
! Selected Full Configuration Interaction with Stochastic selection and PT2.
END_DOC
integer :: i,j,k
- double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:)
+ double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:), zeros(:)
integer :: to_select
logical, external :: qp_stop
@@ -18,7 +18,7 @@ subroutine run_stochastic_cipsi
rss = memory_of_double(N_states)*4.d0
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
logical :: has
@@ -26,6 +26,7 @@ subroutine run_stochastic_cipsi
relative_error=PT2_relative_error
+ zeros = 0.d0
pt2 = -huge(1.e0)
rpt2 = -huge(1.e0)
norm = 0.d0
@@ -63,14 +64,14 @@ subroutine run_stochastic_cipsi
do while ( &
(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. &
(correlation_energy_ratio <= correlation_energy_ratio_max) &
)
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)
pt2 = 0.d0
@@ -79,17 +80,17 @@ subroutine run_stochastic_cipsi
call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, &
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) / &
- (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)
- call save_energy(psi_energy_with_nucl_rep, rpt2)
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)
- do k=1,N_states
- rpt2(:) = pt2(:)/(1.d0 + norm(k))
- enddo
call save_energy(psi_energy_with_nucl_rep, rpt2)
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 save_wavefunction
- rpt2(:) = 0.d0
- call save_energy(psi_energy_with_nucl_rep, rpt2)
+ call save_energy(psi_energy_with_nucl_rep, zeros)
if (qp_stop()) exit
enddo
@@ -117,20 +117,18 @@ subroutine run_stochastic_cipsi
if (N_det < N_det_max) then
call diagonalize_CI
call save_wavefunction
- rpt2(:) = 0.d0
- call save_energy(psi_energy_with_nucl_rep, rpt2)
+ call save_energy(psi_energy_with_nucl_rep, zeros)
endif
- pt2 = 0.d0
- variance = 0.d0
- norm = 0.d0
+ pt2(:) = 0.d0
+ variance(:) = 0.d0
+ norm(:) = 0.d0
call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, &
norm,0) ! Stochastic PT2
do k=1,N_states
- rpt2(:) = pt2(:)/(1.d0 + norm(k))
+ rpt2(k) = pt2(k)/(1.d0 + norm(k))
enddo
- 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)
diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f
index c7c11eec..081d998f 100644
--- a/src/cipsi/zmq_selection.irp.f
+++ b/src/cipsi/zmq_selection.irp.f
@@ -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_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_order
+ PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym
+
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
stop 'Unable to put state_average_weight on ZMQ server'
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
stop 'Unable to put threshold_generators on ZMQ server'
endif
@@ -85,7 +89,11 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
endif
integer :: nproc_target
- nproc_target = nproc
+ if (N_det < 3*nproc) then
+ nproc_target = N_det/4
+ else
+ nproc_target = nproc
+ endif
double precision :: mem
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)')
@@ -131,13 +139,7 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
norm(k) = norm(k) * f(k)
enddo
- ! Adjust PT2 weights for next selection
- 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
+ call update_pt2_and_variance_weights(pt2, variance, norm, N_states)
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
type(selection_buffer), intent(inout) :: b
integer, intent(in) :: N
- double precision, intent(inout) :: pt2(N_states)
- double precision, intent(inout) :: variance(N_states)
- double precision, intent(inout) :: norm(N_states)
+ double precision, intent(out) :: pt2(N_states)
+ double precision, intent(out) :: variance(N_states)
+ double precision, intent(out) :: norm(N_states)
double precision :: pt2_mwen(N_states)
double precision :: variance_mwen(N_states)
double precision :: norm_mwen(N_states)
diff --git a/src/davidson/print_e_components.irp.f b/src/davidson/print_e_components.irp.f
new file mode 100644
index 00000000..ddf83474
--- /dev/null
+++ b/src/davidson/print_e_components.irp.f
@@ -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
diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg
index b8fc7406..95128969 100644
--- a/src/determinants/EZFIO.cfg
+++ b/src/determinants/EZFIO.cfg
@@ -28,12 +28,18 @@ doc: Force the wave function to be an eigenfunction of |S^2|
interface: ezfio,provider,ocaml
default: True
-[used_weight]
+[weight_one_e_dm]
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))
interface: ezfio,provider,ocaml
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]
type: Threshold
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
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]
type: Threshold
diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f
index a9630977..e4f76bca 100644
--- a/src/determinants/density_matrix.irp.f
+++ b/src/determinants/density_matrix.irp.f
@@ -305,9 +305,9 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
logical :: exists
state_average_weight(:) = 1.d0
- if (used_weight == 0) then
+ if (weight_one_e_dm == 0) then
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
else
call ezfio_has_determinants_state_average_weight(exists)
diff --git a/src/determinants/utils.irp.f b/src/determinants/utils.irp.f
index 20d9e1e5..3aec16f9 100644
--- a/src/determinants/utils.irp.f
+++ b/src/determinants/utils.irp.f
@@ -43,4 +43,3 @@ BEGIN_PROVIDER [ double precision, S2_matrix_all_dets,(N_det,N_det) ]
!$OMP END PARALLEL DO
END_PROVIDER
-
diff --git a/src/dft_utils_in_r/ao_in_r.irp.f b/src/dft_utils_in_r/ao_in_r.irp.f
index 17892832..767f329c 100644
--- a/src/dft_utils_in_r/ao_in_r.irp.f
+++ b/src/dft_utils_in_r/ao_in_r.irp.f
@@ -121,3 +121,26 @@
enddo
END_PROVIDER
+
+ BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom, (ao_num,n_pts_max_per_atom,nucl_num)]
+&BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom_transp, (n_pts_max_per_atom,ao_num,nucl_num)]
+ implicit none
+ 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
+
diff --git a/src/fci/pt2.irp.f b/src/fci/pt2.irp.f
index 45860fb5..6135864f 100644
--- a/src/fci/pt2.irp.f
+++ b/src/fci/pt2.irp.f
@@ -46,7 +46,7 @@ subroutine run
call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, &
norm,0) ! Stochastic PT2
do k=1,N_states
- rpt2(:) = pt2(:)/(1.d0 + norm(k))
+ rpt2(k) = pt2(k)/(1.d0 + norm(k))
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)
diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f
index a8037982..ad87bc8e 100644
--- a/src/iterations/print_summary.irp.f
+++ b/src/iterations/print_summary.irp.f
@@ -31,18 +31,19 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
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,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
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)
if (N_states_p > 1) then
write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1)
write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0
endif
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(*,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)
@@ -97,5 +98,7 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_
enddo
endif
+ call print_energy_components()
+
end subroutine
diff --git a/src/nuclei/atomic_radii.irp.f b/src/nuclei/atomic_radii.irp.f
index 82487b9d..439b5cec 100644
--- a/src/nuclei/atomic_radii.irp.f
+++ b/src/nuclei/atomic_radii.irp.f
@@ -66,7 +66,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom, (nucl_num)]
implicit none
integer :: i
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
END_PROVIDER
@@ -74,7 +74,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom_ua, (nucl_num)]
implicit none
integer :: i
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
END_PROVIDER