Merge branch 'dev-lcpq' into dev-lct

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
+  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
   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

INSTALL.rst

@@ -45,6 +45,8 @@ Requirements
 - |ZeroMQ| : networking library
 - `GMP <https://gmplib.org/>`_ : Gnu Multiple Precision Arithmetic Library
 - |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
 
 
@@ -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:
+  `<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
 -----
 

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)     
+<img src="https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png" width="250">
 
 # Getting started
 

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."

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
   ;;

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, zeros)
-    call save_energy(psi_energy_with_nucl_rep, rpt2)
     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, zeros)
-        call save_energy(psi_energy_with_nucl_rep, rpt2)
     endif
 
     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
       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

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
+  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
 

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_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()
@@ -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)

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
+   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
+        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)
-!        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 

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')

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, zeros)
-    call save_energy(psi_energy_with_nucl_rep, rpt2)
     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, zeros)
-        call save_energy(psi_energy_with_nucl_rep, rpt2)
     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, &
       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)

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
+  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
+  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
 
@@ -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(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                   :: variance_mwen(N_states)
   double precision                   :: norm_mwen(N_states)

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

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

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)

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
 
-

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)

src/iterations/print_summary.irp.f

@@ -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)
-  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)
@@ -43,6 +43,7 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_
   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) '# 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