Merge pull request #45 from QuantumPackage/dev-lcpq

Dev lcpq

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
 

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

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

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

data/basis/cc-pv5z_ecp_ncsu

@@ -1,5 +1,5 @@
 BORON
- S 9 1.00
+ S 9
 1 11.76050    -0.0036757
 2 5.150520     0.0250517 
 3 2.578276    -0.1249228
@@ -9,15 +9,15 @@ BORON
 7 0.161898     0.4308431
 8 0.081044     0.2486558 
 9 0.040569     0.0317295
-S 1 1.00
+S 1
 1 0.614105     1.000000
-S 1 1.00
+S 1
 1 0.375720     1.000000
-S 1 1.00
+S 1
 1 0.170896     1.000000
-S 1 1.00
+S 1
 1 0.070664     1.000000
-P 9 1.00
+P 9
 1 7.470701     0.0047397
 2 3.735743     0.0376009
 3 1.868068     0.0510600
@@ -27,37 +27,37 @@ P 9 1.00
 7 0.116803     0.2850185
 8 0.058408     0.1448808
 9 0.029207     0.0176962
-P 1 1.00
+P 1 
 1 0.566611     1.000000
-P 1 1.00
+P 1
 1 0.436327     1.000000
-P 1 1.00
+P 1
 1 0.143772     1.000000
-P 1 1.00
+P 1
 1 0.057917     1.000000
-D 1 1.00
+D 1
 1 1.022256     1.000000
-D 1 1.00
+D 1
 1 0.808233     1.000000
-D 1 1.00
+D 1
 1 0.380163     1.000000
-D 1 1.00
+D 1
 1 0.134838     1.000000
-F 1 1.00
+F 1
 1 1.002171     1.000000
-F 1 1.00
+F 1
 1 0.799174     1.000000
-F 1 1.00
+F 1
 1 0.272717     1.000000
-G 1 1.00
+G 1
 1 0.824366     1.000000
-G 1 1.00
+G 1
 1 0.486131     1.000000
-H 1 1.00
+H 1
 1 0.632779     1.000000
  
 CARBON
- S 9 1.00
+ S 9
 1 13.073594    0.0051583
 2 6.541187     0.0603424
 3 4.573411    -0.1978471
@@ -67,15 +67,15 @@ CARBON
 7 0.231300     0.4336405
 8 0.102619     0.2131940
 9 0.051344     0.0049848
-S 1 1.00
+S 1
 1 0.098302     1.000000
-S 1 1.00
+S 1
 1 0.232034     1.000000
-S 1 1.00
+S 1
 1 0.744448     1.000000
-S 1 1.00
+S 1
 1 1.009914     1.000000
-P 9 1.00
+P 9
 1 9.934169     0.0209076
 2 3.886955     0.0572698
 3 1.871016     0.1122682
@@ -85,33 +85,33 @@ P 9 1.00
 7 0.117063     0.2016934
 8 0.058547     0.0453575
 9 0.029281     0.0029775
-P 1 1.00
+P 1
 1 0.084047     1.000000
-P 1 1.00
+P 1
 1 0.216618     1.000000
-P 1 1.00
+P 1
 1 0.576869     1.000000
-P 1 1.00
+P 1
 1 1.006252     1.000000
-D 1 1.00
+D 1
 1 0.206619     1.000000
-D 1 1.00
+D 1
 1 0.606933     1.000000
-D 1 1.00
+D 1
 1 1.001526     1.000000
-D 1 1.00
+D 1
 1 1.504882     1.000000
-F 1 1.00
+F 1
 1 0.400573     1.000000
-F 1 1.00
+F 1
 1 1.099564     1.000000
-F 1 1.00
+F 1
 1 1.501091     1.000000
-G 1 1.00
+G 1
 1 0.797648     1.000000
-G 1 1.00
+G 1
 1 1.401343     1.000000
-H 1 1.00
+H 1
 1 1.001703     1.000000
  
 NITROGEN
@@ -320,7 +320,7 @@ G   1
 H   1
   1       0.249975    1.00000000
  
-MG
+MAGNESIUM
  S   12
   1      63.931893   -0.00079400
   2      31.602596    0.00747900
@@ -500,7 +500,7 @@ G   1
 H   1
   1       0.441969    1.00000000
  
-SI
+SILICON
  S   12
   1      96.651837   -0.00044000
   2      48.652547    0.01867100
@@ -625,7 +625,7 @@ S   1
   1       0.216983    1.00000000
 S   1
   1       0.094992    1.00000000
-PHOSPHORUS   12
+P   12
   1      48.154282    0.01288400
   2      25.406431    0.09709500
   3      13.404555    0.17821500
@@ -638,7 +638,7 @@ PHOSPHORUS   12
  10       0.152550   -0.00006600
  11       0.080486    0.00012900
  12       0.042465   -0.00002900
-PHOSPHORUS   12
+P   12
   1      48.154282   -0.00315200
   2      25.406431   -0.02300600
   3      13.404555   -0.04239800
@@ -651,13 +651,13 @@ PHOSPHORUS   12
  10       0.152550    0.34480600
  11       0.080486    0.18173100
  12       0.042465    0.03664900
-PHOSPHORUS   1
+P   1
   1       2.103768    1.00000000
-PHOSPHORUS   1
+P   1
   1       0.988246    1.00000000
-PHOSPHORUS   1
+P   1
   1       0.464229    1.00000000
-PHOSPHORUS   1
+P   1
   1       0.218072    1.00000000
 D   1
   1       2.461614    1.00000000
@@ -681,7 +681,7 @@ H   1
   1       0.737477    1.00000000
  
 SULFUR
- SULFUR   12
+ S   12
   1     306.317903    0.00006400
   2     146.602801   -0.00078500
   3      70.163647    0.02247100
@@ -694,7 +694,7 @@ SULFUR
  10       0.403554    0.00443600
  11       0.193140   -0.00101500
  12       0.092436    0.00050700
-SULFUR   12
+S   12
   1     306.317903    0.00002100
   2     146.602801   -0.00000400
   3      70.163647   -0.00611900
@@ -707,13 +707,13 @@ SULFUR   12
  10       0.403554    0.52928200
  11       0.193140    0.46625400
  12       0.092436    0.12580000
-SULFUR   1
+S   1
   1       1.476765    1.00000000
-SULFUR   1
+S   1
   1       0.646261    1.00000000
-SULFUR   1
+S   1
   1       0.282816    1.00000000
-SULFUR   1
+S   1
   1       0.123766    1.00000000
 P   12
   1      55.148271    0.01344700
@@ -1206,7 +1206,7 @@ H  1
 I  1
   1       2.384691  1.00000000
  
-TI
+TITANIUM
  S  13
   1      68.910511  0.00061600
   2      33.720700  -0.00750100
@@ -2742,7 +2742,7 @@ H  1
 I  1
   1       3.729940  1.00000000
  
-NI
+NICKEL
  S  13
   1      97.161835  0.00070900
   2      51.187866  -0.01239900

data/basis/cc-pvdz_ecp_ncsu

@@ -1,5 +1,5 @@
 BORON
- S 9 1.00
+ S 9
 1 11.76050    -0.0036757
 2 5.150520     0.0250517 
 3 2.578276    -0.1249228
@@ -9,9 +9,9 @@ BORON
 7 0.161898     0.4308431
 8 0.081044     0.2486558 
 9 0.040569     0.0317295
-S 1 1.00
+S 1
 1 0.082513     1.000000
-P 9 1.00
+P 9
 1 7.470701     0.0047397
 2 3.735743     0.0376009
 3 1.868068     0.0510600
@@ -21,13 +21,13 @@ P 9 1.00
 7 0.116803     0.2850185
 8 0.058408     0.1448808
 9 0.029207     0.0176962
-P 1 1.00
+P 1
 1 0.086803     1.000000
-D 1 1.00
+D 1
 1 0.349879     1.000000
  
 CARBON
- S 9 1.00
+ S 9
 1 13.073594    0.0051583
 2 6.541187     0.0603424
 3 4.573411    -0.1978471
@@ -37,9 +37,9 @@ CARBON
 7 0.231300     0.4336405
 8 0.102619     0.2131940
 9 0.051344     0.0049848
-S 1 1.00
+S 1
 1 0.127852     1.000000
-P 9 1.00
+P 9
 1 9.934169     0.0209076
 2 3.886955     0.0572698
 3 1.871016     0.1122682
@@ -49,9 +49,9 @@ P 9 1.00
 7 0.117063     0.2016934
 8 0.058547     0.0453575
 9 0.029281     0.0029775
-P 1 1.00
+P 1
 1 0.149161     1.000000
-D 1 1.00
+D 1
 1 0.561160     1.000000
  
 NITROGEN
@@ -170,7 +170,7 @@ P   1
 D   1
   1       0.093145    1.00000000
  
-MG
+MAGNESIUM
  S   12
   1      63.931893   -0.00079400
   2      31.602596    0.00747900
@@ -290,7 +290,7 @@ P   1
 D   1
   1       0.189387    1.00000000
  
-SI
+SILICON
  S   12
   1      96.651837   -0.00044000
   2      48.652547    0.01867100
@@ -379,7 +379,7 @@ S   12
  12       0.072398    0.12409900
 S   1
   1       0.111116    1.00000000
-PHOSPHORUS   12
+P   12
   1      48.154282    0.01288400
   2      25.406431    0.09709500
   3      13.404555    0.17821500
@@ -392,7 +392,7 @@ PHOSPHORUS   12
  10       0.152550   -0.00006600
  11       0.080486    0.00012900
  12       0.042465   -0.00002900
-PHOSPHORUS   12
+P   12
   1      48.154282   -0.00315200
   2      25.406431   -0.02300600
   3      13.404555   -0.04239800
@@ -405,13 +405,13 @@ PHOSPHORUS   12
  10       0.152550    0.34480600
  11       0.080486    0.18173100
  12       0.042465    0.03664900
-PHOSPHORUS   1
+P   1
   1       0.110006    1.00000000
 D   1
   1       0.373518    1.00000000
  
 SULFUR
- SULFUR   12
+ S   12
   1     306.317903    0.00006400
   2     146.602801   -0.00078500
   3      70.163647    0.02247100
@@ -424,7 +424,7 @@ SULFUR
  10       0.403554    0.00443600
  11       0.193140   -0.00101500
  12       0.092436    0.00050700
-SULFUR   12
+S   12
   1     306.317903    0.00002100
   2     146.602801   -0.00000400
   3      70.163647   -0.00611900
@@ -437,7 +437,7 @@ SULFUR   12
  10       0.403554    0.52928200
  11       0.193140    0.46625400
  12       0.092436    0.12580000
-SULFUR   1
+S   1
   1       0.138193    1.00000000
 P   12
   1      55.148271    0.01344700
@@ -708,7 +708,7 @@ D  1
 F  1
   1       0.280673  1.00000000
  
-TI
+TITANIUM
  S  13
   1      68.910511  0.00061600
   2      33.720700  -0.00750100
@@ -1770,3 +1770,4 @@ D  1
 F  1
   1       3.171936  1.00000000
  
+

data/basis/cc-pvqz_ecp_ncsu

@@ -1,5 +1,5 @@
 BORON
- S 9 1.00
+ S 9
 1 11.76050    -0.0036757
 2 5.150520     0.0250517 
 3 2.578276    -0.1249228
@@ -9,13 +9,13 @@ BORON
 7 0.161898     0.4308431
 8 0.081044     0.2486558 
 9 0.040569     0.0317295
-S 1 1.00
+S 1
 1 0.422217     1.000000 
-S 1 1.00
+S 1
 1 0.305133     1.000000
-S 1 1.00
+S 1
 1 0.082968     1.000000
-P 9 1.00
+P 9
 1 7.470701     0.0047397
 2 3.735743     0.0376009
 3 1.868068     0.0510600
@@ -25,27 +25,27 @@ P 9 1.00
 7 0.116803     0.2850185
 8 0.058408     0.1448808
 9 0.029207     0.0176962
-P 1 1.00
+P 1 
 1 0.447031     1.000000
-P 1 1.00
+P 1 
 1 0.196614     1.000000
-P 1 1.00
+P 1 
 1 0.066445     1.000000
-D 1 1.00
+D 1 
 1 1.142614     1.000000
-D 1 1.00
+D 1 
 1 0.410733     1.000000
-D 1 1.00
+D 1 
 1 0.149100     1.000000
-F 1 1.00
+F 1 
 1 0.870011     1.000000
-F 1 1.00
+F 1 
 1 0.315902     1.000000
-G 1 1.00
+G 1 
 1 0.710746     1.000000
  
 CARBON
- S 9 1.00
+ S 9
 1 13.073594    0.0051583
 2 6.541187     0.0603424
 3 4.573411    -0.1978471
@@ -55,13 +55,13 @@ CARBON
 7 0.231300     0.4336405
 8 0.102619     0.2131940
 9 0.051344     0.0049848
-S 1 1.00
+S 1 
 1 0.846879     1.000000
-S 1 1.00
+S 1 
 1 0.269659     1.000000
-S 1 1.00
+S 1 
 1 0.109576     1.000000
-P 9 1.00
+P 9 
 1 9.934169     0.0209076
 2 3.886955     0.0572698
 3 1.871016     0.1122682
@@ -71,23 +71,23 @@ P 9 1.00
 7 0.117063     0.2016934
 8 0.058547     0.0453575
 9 0.029281     0.0029775
-P 1 1.00
+P 1
 1 0.804681     1.000000
-P 1 1.00
+P 1
 1 0.313254     1.000000
-P 1 1.00
+P 1
 1 0.105389     1.000000
-D 1 1.00
+D 1 
 1 2.013760     1.000000
-D 1 1.00
+D 1
 1 0.684884     1.000000
-D 1 1.00
+D 1
 1 0.240171     1.000000
-F 1 1.00
+F 1
 1 0.457302     1.000000
-F 1 1.00
+F 1
 1 1.324930     1.000000
-G 1 1.00
+G 1
 1 1.034180     1.000000
  
 NITROGEN
@@ -260,7 +260,7 @@ F   1
 G   1
   1       0.169674    1.00000000
  
-MG
+MAGNESIUM
  S   12
   1      63.931893   -0.00079400
   2      31.602596    0.00747900
@@ -416,7 +416,7 @@ F   1
 G   1
   1       0.352315    1.00000000
  
-SI
+SILICON
  S   12
   1      96.651837   -0.00044000
   2      48.652547    0.01867100
@@ -527,7 +527,7 @@ S   1
   1       0.276109    1.00000000
 S   1
   1       0.110986    1.00000000
-PHOSPHORUS   12
+P   12
   1      48.154282    0.01288400
   2      25.406431    0.09709500
   3      13.404555    0.17821500
@@ -540,7 +540,7 @@ PHOSPHORUS   12
  10       0.152550   -0.00006600
  11       0.080486    0.00012900
  12       0.042465   -0.00002900
-PHOSPHORUS   12
+P   12
   1      48.154282   -0.00315200
   2      25.406431   -0.02300600
   3      13.404555   -0.04239800
@@ -553,11 +553,11 @@ PHOSPHORUS   12
  10       0.152550    0.34480600
  11       0.080486    0.18173100
  12       0.042465    0.03664900
-PHOSPHORUS   1
+P   1
   1       1.581855    1.00000000
-PHOSPHORUS   1
+P   1
   1       0.658132    1.00000000
-PHOSPHORUS   1
+P   1
   1       0.273816    1.00000000
 D   1
   1       0.984476    1.00000000
@@ -573,7 +573,7 @@ G   1
   1       0.589583    1.00000000
  
 SULFUR
- SULFUR   12
+ S   12
   1     306.317903    0.00006400
   2     146.602801   -0.00078500
   3      70.163647    0.02247100
@@ -586,7 +586,7 @@ SULFUR
  10       0.403554    0.00443600
  11       0.193140   -0.00101500
  12       0.092436    0.00050700
-SULFUR   12
+S   12
   1     306.317903    0.00002100
   2     146.602801   -0.00000400
   3      70.163647   -0.00611900
@@ -599,11 +599,11 @@ SULFUR   12
  10       0.403554    0.52928200
  11       0.193140    0.46625400
  12       0.092436    0.12580000
-SULFUR   1
+S   1
   1       0.893039    1.00000000
-SULFUR   1
+S   1
   1       0.358094    1.00000000
-SULFUR   1
+S   1
   1       0.143590    1.00000000
 P   12
   1      55.148271    0.01344700
@@ -1014,7 +1014,7 @@ G  1
 H  1
   1       3.057835  1.00000000
  
-TI
+TITANIUM
  S  13
   1      68.910511  0.00061600
   2      33.720700  -0.00750100
@@ -2262,7 +2262,7 @@ G  1
 H  1
   1       2.901601  1.00000000
  
-NI
+NICKEL
  S  13
   1      97.161835  0.00070900
   2      51.187866  -0.01239900
@@ -2886,3 +2886,4 @@ G  1
 H  1
   1       4.405055  1.00000000
  
+

data/basis/cc-pvtz_ecp_ncsu

@@ -1,5 +1,5 @@
 BORON
- S 9 1.00
+ S 9
 1 11.76050    -0.0036757
 2 5.150520     0.0250517 
 3 2.578276    -0.1249228
@@ -9,11 +9,11 @@ BORON
 7 0.161898     0.4308431
 8 0.081044     0.2486558 
 9 0.040569     0.0317295
-S 1 1.00
+S 1
 1 0.626026     1.000000
-S 1 1.00
+S 1
 1 0.092094     1.000000
-P 9 1.00
+P 9
 1 7.470701     0.0047397
 2 3.735743     0.0376009
 3 1.868068     0.0510600
@@ -23,19 +23,19 @@ P 9 1.00
 7 0.116803     0.2850185
 8 0.058408     0.1448808
 9 0.029207     0.0176962
-P 1 1.00
+P 1 
 1 0.235016     1.000000
-P 1 1.00
+P 1
 1 0.082056     1.000000
-D 1 1.00
+D 1
 1 0.699153     1.000000
-D 1 1.00
+D 1
 1 0.207316     1.000000
-F 1 1.00
+F 1
 1 0.478872     1.000000
  
 CARBON
- S 9 1.00
+ S 9 
 1 13.073594    0.0051583
 2 6.541187     0.0603424
 3 4.573411    -0.1978471
@@ -45,11 +45,11 @@ CARBON
 7 0.231300     0.4336405
 8 0.102619     0.2131940
 9 0.051344     0.0049848
-S 1 1.00
+S 1 
 1 0.921552     1.000000
-S 1 1.00
+S 1 
 1 0.132800     1.000000
-P 9 1.00
+P 9 
 1 9.934169     0.0209076
 2 3.886955     0.0572698
 3 1.871016     0.1122682
@@ -59,15 +59,15 @@ P 9 1.00
 7 0.117063     0.2016934
 8 0.058547     0.0453575
 9 0.029281     0.0029775
-P 1 1.00
+P 1 
 1 0.376742     1.000000
-P 1 1.00
+P 1 
 1 0.126772     1.000000
-D 1 1.00
+D 1 
 1 1.141611     1.000000
-D 1 1.00
+D 1 
 1 0.329486     1.000000
-F 1 1.00
+F 1 
 1 0.773485     1.000000
  
 NITROGEN
@@ -210,7 +210,7 @@ D   1
 F   1
   1       0.134692    1.00000000
  
-MG
+MAGNESIUM
  S   12
   1      63.931893   -0.00079400
   2      31.602596    0.00747900
@@ -346,7 +346,7 @@ D   1
 F   1
   1       0.254342    1.00000000
  
-SI
+SILICON
  S   12
   1      96.651837   -0.00044000
   2      48.652547    0.01867100
@@ -445,7 +445,7 @@ S   1
   1       0.115599    1.00000000
 S   1
   1       0.725938    1.00000000
-PHOSPHORUS   12
+P   12
   1      48.154282    0.01288400
   2      25.406431    0.09709500
   3      13.404555    0.17821500
@@ -458,7 +458,7 @@ PHOSPHORUS   12
  10       0.152550   -0.00006600
  11       0.080486    0.00012900
  12       0.042465   -0.00002900
-PHOSPHORUS   12
+P   12
   1      48.154282   -0.00315200
   2      25.406431   -0.02300600
   3      13.404555   -0.04239800
@@ -471,9 +471,9 @@ PHOSPHORUS   12
  10       0.152550    0.34480600
  11       0.080486    0.18173100
  12       0.042465    0.03664900
-PHOSPHORUS   1
+P   1
   1       0.093233    1.00000000
-PHOSPHORUS   1
+P   1
   1       0.625438    1.00000000
 D   1
   1       0.214392    1.00000000
@@ -483,7 +483,7 @@ F   1
   1       0.457026    1.00000000
  
 SULFUR
- SULFUR   12
+ S   12
   1     306.317903    0.00006400
   2     146.602801   -0.00078500
   3      70.163647    0.02247100
@@ -496,7 +496,7 @@ SULFUR
  10       0.403554    0.00443600
  11       0.193140   -0.00101500
  12       0.092436    0.00050700
-SULFUR   12
+S   12
   1     306.317903    0.00002100
   2     146.602801   -0.00000400
   3      70.163647   -0.00611900
@@ -509,9 +509,9 @@ SULFUR   12
  10       0.403554    0.52928200
  11       0.193140    0.46625400
  12       0.092436    0.12580000
-SULFUR   1
+S   1
   1       0.898845    1.00000000
-SULFUR   1
+S   1
   1       0.146364    1.00000000
 P   12
   1      55.148271    0.01344700
@@ -848,7 +848,7 @@ F  1
 G  1
   1       0.413265  1.00000000
  
-TI
+TITANIUM
  S  13
   1      68.910511  0.00061600
   2      33.720700  -0.00750100
@@ -1820,7 +1820,7 @@ F  1
 G  1
   1       2.525364  1.00000000
  
-NI
+NICKEL
  S  13
   1      97.161835  0.00070900
   2      51.187866  -0.01239900
@@ -2306,3 +2306,4 @@ F  1
 G  1
   1       4.099533  1.00000000
  
+

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/ao_two_e_erf_ints/two_e_integrals_erf.irp.f

@@ -16,6 +16,8 @@ double precision function ao_two_e_integral_erf(i,j,k,l)
   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)
      return

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

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

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
+

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)&
+        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
         if(isnan(final_weight_at_r(k,i,j)))then
          print*,'isnan(final_weight_at_r(k,i,j))' 

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 

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
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
-!         cycle
+          cycle
-!       endif
+        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
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
-       !  cycle
+          cycle
-       !endif
+        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)

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

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

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
 

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