Merge pull request #45 from QuantumPackage/dev-lcpq

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

--- a/README.md
+++ b/README.md
@ -1,12 +1,13 @@
 # Quantum Package 2.0


-*Quantum package 2.0: an open-source determinant-driven suite of programs*\
+[*Quantum package 2.0: an open-source determinant-driven suite of programs*](https://pubs.acs.org/doi/10.1021/acs.jctc.9b00176)\
 Y. Garniron, K. Gasperich, T. Applencourt, A. Benali, A. Ferté, J. Paquier, B. Pradines, R. Assaraf, P. Reinhardt, J. Toulouse, P. Barbaresco, N. Renon, G. David, J. P. Malrieu, M. Véril, M. Caffarel, P. F. Loos, E. Giner and A. Scemama\
+J. Chem. Theory Comput. (2019)\
 https://arxiv.org/abs/1902.08154


-![QP](https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png)     
+<img src="https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png" width="250">

 # Getting started

--- a/config/ifort_avx.cfg
+++ b/config/ifort_avx.cfg
@ -32,7 +32,7 @@ OPENMP  : 1          ; Append OpenMP flags
 #
 [OPT]
 FC       : -traceback
-FCFLAGS  : -xAVX -O2 -ip -ftz -g 
+FCFLAGS  : -march=corei7-avx -O2 -ip -ftz -g 

 # Profiling flags
 #################
--- a/config/ifort_avx_mpi.cfg
+++ b/config/ifort_avx_mpi.cfg
@ -31,14 +31,14 @@ OPENMP  : 1          ; Append OpenMP flags
 # -ftz                       : Flushes denormal results to zero
 #
 [OPT]
-FCFLAGS  : -xAVX -O2 -ip -ftz -g -traceback  
+FCFLAGS  : -march=corei7-avx -O2 -ip -ftz -g -traceback  

 # Profiling flags
 #################
 #
 [PROFILE]
 FC       : -p -g
-FCFLAGS  : -xSSE4.2 -O2 -ip -ftz 
+FCFLAGS  : -march=corei7 -O2 -ip -ftz 


 # Debugging flags
--- a/44
+++ b/44
@ -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."
--- a/data/basis/cc-pv5z_ecp_ncsu
+++ b/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
--- a/data/basis/cc-pvdz_ecp_ncsu
+++ b/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
 
+
--- a/data/basis/cc-pvqz_ecp_ncsu
+++ b/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
 
+
--- a/data/basis/cc-pvtz_ecp_ncsu
+++ b/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
 
+
--- a/ocaml/Input_determinants_by_hand.ml
+++ b/ocaml/Input_determinants_by_hand.ml
@ -81,9 +81,6 @@ end = struct
  ;;

  let write_n_det n =
-    let n_det_old =
-      Ezfio.get_determinants_n_det ()
-    in
    Det_number.to_int n
    |> Ezfio.set_determinants_n_det
  ;;
--- a/src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f
+++ b/src/ao_two_e_erf_ints/two_e_integrals_erf.irp.f
@ -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
--- a/src/ao_two_e_ints/map_integrals.irp.f
+++ b/src/ao_two_e_ints/map_integrals.irp.f
@ -279,6 +279,100 @@ subroutine get_ao_two_e_integrals_non_zero(j,k,l,sze,out_val,out_val_index,non_z
 end


+subroutine get_ao_two_e_integrals_non_zero_jl(j,l,thresh,sze_max,sze,out_val,out_val_index,non_zero_int)
+  use map_module
+  implicit none
+  BEGIN_DOC
+  ! Gets multiple AO bi-electronic integral from the AO map .
+  ! All non-zero i are retrieved for j,k,l fixed.
+  END_DOC
+  double precision, intent(in)   :: thresh
+  integer, intent(in)            :: j,l, sze,sze_max
+  real(integral_kind), intent(out) :: out_val(sze_max)
+  integer, intent(out)           :: out_val_index(2,sze_max),non_zero_int
+
+  integer                        :: i,k
+  integer(key_kind)              :: hash
+  double precision               :: tmp
+
+  PROVIDE ao_two_e_integrals_in_map
+  non_zero_int = 0
+  if (ao_overlap_abs(j,l) < thresh) then
+    out_val = 0.d0
+    return
+  endif
+
+  non_zero_int = 0
+  do k = 1, sze
+   do i = 1, sze
+     integer, external :: ao_l4
+     double precision, external :: ao_two_e_integral
+     !DIR$ FORCEINLINE
+     if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then
+       cycle
+     endif
+     call two_e_integrals_index(i,j,k,l,hash)
+     call map_get(ao_integrals_map, hash,tmp)
+     if (dabs(tmp) < thresh ) cycle
+     non_zero_int = non_zero_int+1
+     out_val_index(1,non_zero_int) = i
+     out_val_index(2,non_zero_int) = k
+     out_val(non_zero_int) = tmp
+   enddo
+  enddo
+
+end
+
+
+subroutine get_ao_two_e_integrals_non_zero_jl_from_list(j,l,thresh,list,n_list,sze_max,out_val,out_val_index,non_zero_int)
+  use map_module
+  implicit none
+  BEGIN_DOC
+  ! Gets multiple AO two-electron integrals from the AO map .
+  ! All non-zero i are retrieved for j,k,l fixed.
+  END_DOC
+  double precision, intent(in)   :: thresh
+  integer, intent(in)            :: sze_max
+  integer, intent(in)            :: j,l, n_list,list(2,sze_max)
+  real(integral_kind), intent(out) :: out_val(sze_max)
+  integer, intent(out)           :: out_val_index(2,sze_max),non_zero_int
+
+  integer                        :: i,k
+  integer(key_kind)              :: hash
+  double precision               :: tmp
+
+  PROVIDE ao_two_e_integrals_in_map
+  non_zero_int = 0
+  if (ao_overlap_abs(j,l) < thresh) then
+    out_val = 0.d0
+    return
+  endif
+
+  non_zero_int = 0
+ integer :: kk
+  do kk = 1, n_list
+   k = list(1,kk)
+   i = list(2,kk)
+   integer, external :: ao_l4
+   double precision, external :: ao_two_e_integral
+   !DIR$ FORCEINLINE
+   if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < thresh) then
+     cycle
+   endif
+   call two_e_integrals_index(i,j,k,l,hash)
+   call map_get(ao_integrals_map, hash,tmp)
+   if (dabs(tmp) < thresh ) cycle
+   non_zero_int = non_zero_int+1
+   out_val_index(1,non_zero_int) = i
+   out_val_index(2,non_zero_int) = k
+   out_val(non_zero_int) = tmp
+  enddo
+
+end
+
+
+
+
 function get_ao_map_size()
  implicit none
  integer (map_size_kind) :: get_ao_map_size
--- a/src/becke_numerical_grid/EZFIO.cfg
+++ b/src/becke_numerical_grid/EZFIO.cfg
@ -8,3 +8,9 @@ default: 2
 type: integer
 doc: Total number of grid points 
 interface: ezfio
+
+[thresh_grid]
+type: double precision
+doc: threshold on the weight of a given grid point
+interface: ezfio,provider,ocaml
+default: 1.e-20
--- a/src/becke_numerical_grid/atomic_number.irp.f
+++ b/src/becke_numerical_grid/atomic_number.irp.f
@ -0,0 +1,9 @@
+BEGIN_PROVIDER [ integer, grid_atomic_number, (nucl_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Atomic number used to adjust the grid
+ END_DOC
+ grid_atomic_number(:) = max(1,int(nucl_charge(:)))
+
+END_PROVIDER
+
--- a/src/becke_numerical_grid/grid_becke.irp.f
+++ b/src/becke_numerical_grid/grid_becke.irp.f
@ -146,7 +146,7 @@ BEGIN_PROVIDER [double precision, grid_points_per_atom, (3,n_points_integration_
      x = grid_points_radial(j)

      ! value of the radial coordinate for the integration
-      r = knowles_function(alpha_knowles(int(nucl_charge(i))),m_knowles,x)
+      r = knowles_function(alpha_knowles(grid_atomic_number(i)),m_knowles,x)

      ! explicit values of the grid points centered around each atom
      do k = 1, n_points_integration_angular
@ -232,8 +232,8 @@ BEGIN_PROVIDER [double precision, final_weight_at_r, (n_points_integration_angul
    do i = 1, n_points_radial_grid -1 !for each radial grid attached to the "jth" atom
      x = grid_points_radial(i) ! x value for the mapping of the [0, +\infty] to [0,1]
      do k = 1, n_points_integration_angular  ! for each angular point attached to the "jth" atom
-        contrib_integration = derivative_knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)&
-            *knowles_function(alpha_knowles(int(nucl_charge(j))),m_knowles,x)**2
+        contrib_integration = derivative_knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)&
+            *knowles_function(alpha_knowles(grid_atomic_number(j)),m_knowles,x)**2
        final_weight_at_r(k,i,j) = weights_angular_points(k)  * weight_at_r(k,i,j) * contrib_integration * dr_radial_integral
        if(isnan(final_weight_at_r(k,i,j)))then
         print*,'isnan(final_weight_at_r(k,i,j))' 
--- a/src/becke_numerical_grid/grid_becke_per_atom.irp.f
+++ b/src/becke_numerical_grid/grid_becke_per_atom.irp.f
@ -0,0 +1,53 @@
+
+
+ BEGIN_PROVIDER [integer, n_pts_per_atom, (nucl_num)]
+&BEGIN_PROVIDER [integer, n_pts_max_per_atom]
+  BEGIN_DOC
+  ! Number of points which are non zero
+  END_DOC
+  integer                        :: i,j,k,l
+  n_pts_per_atom = 0
+  do j = 1, nucl_num
+    do i = 1, n_points_radial_grid -1
+      do k = 1, n_points_integration_angular
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+          cycle
+        endif
+        n_pts_per_atom(j) += 1
+      enddo
+    enddo
+  enddo
+  n_pts_max_per_atom = maxval(n_pts_per_atom)
+END_PROVIDER
+
+ BEGIN_PROVIDER [double precision, final_grid_points_per_atom, (3,n_pts_max_per_atom,nucl_num)]
+&BEGIN_PROVIDER [double precision, final_weight_at_r_vector_per_atom, (n_pts_max_per_atom,nucl_num) ]
+&BEGIN_PROVIDER [integer, index_final_points_per_atom, (3,n_pts_max_per_atom,nucl_num) ]
+&BEGIN_PROVIDER [integer, index_final_points_per_atom_reverse, (n_points_integration_angular,n_points_radial_grid,nucl_num) ]
+ implicit none
+ integer                        :: i,j,k,l,i_count(nucl_num)
+ double precision               :: r(3)
+ i_count = 0
+  do j = 1, nucl_num
+    do i = 1, n_points_radial_grid -1
+      do k = 1, n_points_integration_angular
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+          cycle
+        endif
+        i_count(j) += 1
+        final_grid_points_per_atom(1,i_count(j),j) = grid_points_per_atom(1,k,i,j)
+        final_grid_points_per_atom(2,i_count(j),j) = grid_points_per_atom(2,k,i,j)
+        final_grid_points_per_atom(3,i_count(j),j) = grid_points_per_atom(3,k,i,j)
+        final_weight_at_r_vector_per_atom(i_count(j),j) = final_weight_at_r(k,i,j)
+        index_final_points_per_atom(1,i_count(j),j) = k
+        index_final_points_per_atom(2,i_count(j),j) = i
+        index_final_points_per_atom(3,i_count(j),j) = j
+        index_final_points_per_atom_reverse(k,i,j) = i_count(j)
+      enddo
+    enddo
+  enddo
+
+
+
+
+END_PROVIDER 
--- a/src/becke_numerical_grid/grid_becke_vector.irp.f
+++ b/src/becke_numerical_grid/grid_becke_vector.irp.f
@ -1,5 +1,6 @@

 BEGIN_PROVIDER [integer, n_points_final_grid]
+  implicit none
  BEGIN_DOC
  ! Number of points which are non zero
  END_DOC
@ -8,9 +9,9 @@ BEGIN_PROVIDER [integer, n_points_final_grid]
  do j = 1, nucl_num
    do i = 1, n_points_radial_grid -1
      do k = 1, n_points_integration_angular
-!       if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then
-!         cycle
-!       endif
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+          cycle
+        endif
        n_points_final_grid += 1
      enddo
    enddo
@ -39,9 +40,9 @@ END_PROVIDER
  do j = 1, nucl_num
    do i = 1, n_points_radial_grid -1
      do k = 1, n_points_integration_angular
-       !if(dabs(final_weight_at_r(k,i,j)) < 1.d-30)then
-       !  cycle
-       !endif
+        if(dabs(final_weight_at_r(k,i,j)) < thresh_grid)then
+          cycle
+        endif
        i_count += 1
        final_grid_points(1,i_count) = grid_points_per_atom(1,k,i,j)
        final_grid_points(2,i_count) = grid_points_per_atom(2,k,i,j)
--- a/src/becke_numerical_grid/step_function_becke.irp.f
+++ b/src/becke_numerical_grid/step_function_becke.irp.f
@ -31,10 +31,6 @@ double precision function cell_function_becke(r,atom_number)
  double precision               :: mu_ij,nu_ij
  double precision               :: distance_i,distance_j,step_function_becke
  integer                        :: j
-  if(int(nucl_charge(atom_number))==0)then
-   cell_function_becke = 0.d0
-   return
-  endif
  distance_i = (r(1) - nucl_coord_transp(1,atom_number) ) * (r(1) - nucl_coord_transp(1,atom_number))
  distance_i += (r(2) - nucl_coord_transp(2,atom_number) ) * (r(2) - nucl_coord_transp(2,atom_number))
  distance_i += (r(3) - nucl_coord_transp(3,atom_number) ) * (r(3) - nucl_coord_transp(3,atom_number))
@ -42,7 +38,6 @@ double precision function cell_function_becke(r,atom_number)
  cell_function_becke = 1.d0
  do j = 1, nucl_num
    if(j==atom_number)cycle
-    if(int(nucl_charge(j))==0)cycle
    distance_j = (r(1) - nucl_coord_transp(1,j) ) * (r(1) - nucl_coord_transp(1,j))
    distance_j+= (r(2) - nucl_coord_transp(2,j) ) * (r(2) - nucl_coord_transp(2,j))
    distance_j+= (r(3) - nucl_coord_transp(3,j) ) * (r(3) - nucl_coord_transp(3,j))
--- a/src/cipsi/cipsi.irp.f
+++ b/src/cipsi/cipsi.irp.f
@ -5,7 +5,7 @@ subroutine run_cipsi
 ! stochastic PT2.
  END_DOC
  integer                        :: i,j,k
-  double precision, allocatable  :: pt2(:), variance(:), norm(:), rpt2(:)
+  double precision, allocatable  :: pt2(:), variance(:), norm(:), rpt2(:), zeros(:)
  integer                        :: n_det_before, to_select

  double precision :: rss
@ -13,7 +13,7 @@ subroutine run_cipsi
  rss = memory_of_double(N_states)*4.d0
  call check_mem(rss,irp_here)

-  allocate (pt2(N_states), rpt2(N_states), norm(N_states), variance(N_states))
+  allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states))

  double precision               :: hf_energy_ref
  logical                        :: has
@ -23,10 +23,11 @@ subroutine run_cipsi

  relative_error=PT2_relative_error

+  zeros = 0.d0
  pt2 = -huge(1.e0)
  rpt2 = -huge(1.e0)
  norm = 0.d0
-  variance = 0.d0
+  variance = huge(1.e0)

  if (s2_eig) then
    call make_s2_eigenfunction
@ -65,7 +66,8 @@ subroutine run_cipsi

  do while (                                                         &
        (N_det < N_det_max) .and.                                    &
-        (maxval(abs(pt2(1:N_states))) > pt2_max) .and.               &
+        (maxval(abs(rpt2(1:N_states))) > pt2_max) .and.               &
+        (maxval(variance(1:N_states)) > variance_max) .and.               &
        (correlation_energy_ratio <= correlation_energy_ratio_max)    &
        )
      write(*,'(A)')  '--------------------------------------------------------------------------------'
@ -83,17 +85,17 @@ subroutine run_cipsi
      SOFT_TOUCH threshold_generators
    endif

+    do k=1,N_states
+      rpt2(k) = pt2(k)/(1.d0 + norm(k))
+    enddo

    correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref)  /     &
-                    (psi_energy_with_nucl_rep(1) + pt2(1) - hf_energy_ref)
+                    (psi_energy_with_nucl_rep(1) + rpt2(1) - hf_energy_ref)
    correlation_energy_ratio = min(1.d0,correlation_energy_ratio)

    call write_double(6,correlation_energy_ratio, 'Correlation ratio')
    call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)

-    do k=1,N_states
-      rpt2(:) = pt2(:)/(1.d0 + norm(k))
-    enddo
    call save_energy(psi_energy_with_nucl_rep, rpt2)

    call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)
@ -103,9 +105,8 @@ subroutine run_cipsi
    if (qp_stop()) exit 

    n_det_before = N_det
-    to_select = N_det
+    to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
    to_select = max(N_states_diag, to_select)
-!    to_select = min(to_select, N_det_max-n_det_before)
    call ZMQ_selection(to_select, pt2, variance, norm)

    PROVIDE  psi_coef
@ -114,32 +115,30 @@ subroutine run_cipsi

    call diagonalize_CI
    call save_wavefunction
-    rpt2(:) = 0.d0
-    call save_energy(psi_energy_with_nucl_rep, rpt2)
+    call save_energy(psi_energy_with_nucl_rep, zeros)
    if (qp_stop()) exit 
+print *,  (N_det < N_det_max)
+print *,  (maxval(abs(rpt2(1:N_states))) > pt2_max)
+print *,  (maxval(variance(1:N_states)) > variance_max)
+print *,  (correlation_energy_ratio <= correlation_energy_ratio_max)
  enddo

  if (.not.qp_stop()) then
    if (N_det < N_det_max) then
        call diagonalize_CI
        call save_wavefunction
-        rpt2(:) = 0.d0
-        call save_energy(psi_energy_with_nucl_rep, rpt2)
+        call save_energy(psi_energy_with_nucl_rep, zeros)
    endif

    if (do_pt2) then
-      pt2 = 0.d0
-      variance = 0.d0
-      norm = 0.d0
+      pt2(:) = 0.d0
+      variance(:) = 0.d0
+      norm(:) = 0.d0
      threshold_generators = 1d0
      SOFT_TOUCH threshold_generators
      call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, &
        norm,0) ! Stochastic PT2
      SOFT_TOUCH threshold_generators
-      do k=1,N_states
-        rpt2(:) = pt2(:)/(1.d0 + norm(k))
-      enddo
-      call save_energy(psi_energy_with_nucl_rep, pt2)
    endif
    print *,  'N_det             = ', N_det
    print *,  'N_sop             = ', N_occ_pattern
@ -148,12 +147,11 @@ subroutine run_cipsi


    do k=1,N_states
-      rpt2(:) = pt2(:)/(1.d0 + norm(k))
+      rpt2(k) = pt2(k)/(1.d0 + norm(k))
    enddo
-    call save_energy(psi_energy_with_nucl_rep, rpt2)

    call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
-    call save_energy(psi_energy_with_nucl_rep, pt2)
+    call save_energy(psi_energy_with_nucl_rep, rpt2)
    call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)
    call print_extrapolated_energy()
  endif
--- a/src/cipsi/pt2_stoch_routines.irp.f
+++ b/src/cipsi/pt2_stoch_routines.irp.f
@ -129,13 +129,13 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
  PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
  PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
  PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted
-  PROVIDE psi_det_hii N_generators_bitmask
+  PROVIDE psi_det_hii N_generators_bitmask selection_weight pseudo_sym

  if (h0_type == 'SOP') then
    PROVIDE psi_occ_pattern_hii det_to_occ_pattern
  endif

-  if (N_det < max(1000,N_states)) then
+  if (N_det < max(4,N_states)) then
    pt2=0.d0
    variance=0.d0
    norm=0.d0
@ -182,6 +182,9 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
      if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then
        stop 'Unable to put state_average_weight on ZMQ server'
      endif
+      if (zmq_put_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) then
+        stop 'Unable to put selection_weight on ZMQ server'
+      endif
      if (zmq_put_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) then
        stop 'Unable to put pt2_stoch_istate on ZMQ server'
      endif
@ -333,13 +336,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in)
    pt2(k) = 0.d0
  enddo

-  ! Adjust PT2 weights for next selection
-  double precision :: pt2_avg
-  pt2_avg = sum(pt2) / dble(N_states)
-  do k=1,N_states
-    pt2_match_weight(k) *= (pt2(k)/pt2_avg)**2
-  enddo
-  SOFT_TOUCH pt2_match_weight
+  call update_pt2_and_variance_weights(pt2, variance, norm, N_states)

 end subroutine

--- a/src/cipsi/run_selection_slave.irp.f
+++ b/src/cipsi/run_selection_slave.irp.f
@ -25,8 +25,8 @@ subroutine run_selection_slave(thread,iproc,energy)
  PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
  PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
  PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
-  PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F
-  PROVIDE psi_selectors_coef_transp psi_det_sorted
+  PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F pseudo_sym
+  PROVIDE psi_selectors_coef_transp psi_det_sorted weight_selection


  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
@ -230,6 +230,8 @@ subroutine pull_selection_results(zmq_socket_pull, pt2, variance, norm, val, det
      endif
  else
      pt2(:) = 0.d0
+      variance(:) = 0.d0
+      norm(:) = 0.d0
  endif

  rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0)
--- a/src/cipsi/selection.irp.f
+++ b/src/cipsi/selection.irp.f
@ -6,15 +6,108 @@ BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ]
 ! Weights adjusted along the selection to make the PT2 contributions
 ! of each state coincide.
 END_DOC
- pt2_match_weight = 1.d0
+ pt2_match_weight(:) = 1.d0
 END_PROVIDER

+BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Weights adjusted along the selection to make the variances 
+ ! of each state coincide.
+ END_DOC
+ variance_match_weight(:) = 1.d0
+END_PROVIDER
+
+subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st)
+  implicit none
+  BEGIN_DOC
+! Updates the rPT2- and Variance- matching weights.
+  END_DOC
+  integer, intent(in)          :: N_st
+  double precision, intent(in) :: pt2(N_st)
+  double precision, intent(in) :: variance(N_st)
+  double precision, intent(in) :: norm(N_st)
+
+  double precision :: avg, rpt2(N_st), element, dt, x
+  integer          :: k
+  integer, save    :: i_iter=0
+  integer, parameter :: i_itermax = 3
+  double precision, allocatable, save :: memo_variance(:,:), memo_pt2(:,:)
+
+  if (i_iter == 0) then
+    allocate(memo_variance(N_st,i_itermax), memo_pt2(N_st,i_itermax))
+    memo_pt2(:,:) = 1.d0
+    memo_variance(:,:) = 1.d0
+  endif
+
+  i_iter = i_iter+1
+  if (i_iter > i_itermax) then
+    i_iter = 1
+  endif
+
+  dt = 4.d0 
+
+  do k=1,N_st
+    rpt2(k) = pt2(k)/(1.d0 + norm(k))                                                     
+  enddo
+
+  avg = sum(rpt2(1:N_st)) / dble(N_st)
+  do k=1,N_st
+    element = exp(dt*(rpt2(k)/avg -1.d0))
+    element = min(1.5d0 , element)
+    element = max(0.5d0 , element)
+    memo_pt2(k,i_iter) = element
+    pt2_match_weight(k) = product(memo_pt2(k,:))
+  enddo
+
+  avg = sum(variance(1:N_st)) / dble(N_st)
+  do k=1,N_st
+    element = exp(dt*(variance(k)/avg -1.d0))
+    element = min(1.5d0 , element)
+    element = max(0.5d0 , element)
+    memo_variance(k,i_iter) = element
+    variance_match_weight(k) = product(memo_variance(k,:))
+  enddo
+
+  print *, '# PT2 weight ', real(pt2_match_weight(:),4)
+  print *, '# var weight ', real(variance_match_weight(:),4)
+  SOFT_TOUCH pt2_match_weight variance_match_weight
+end
+
+
 BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ]
   implicit none
   BEGIN_DOC
   ! Weights used in the selection criterion
   END_DOC
-   selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states)
+   select case (weight_selection)
+
+     case (0)
+      print *,  'Using input weights in selection'
+      selection_weight(1:N_states) = state_average_weight(1:N_states)
+
+     case (1)
+      print *,  'Using 1/c_max^2 weight in selection'
+      selection_weight(1:N_states) = c0_weight(1:N_states) 
+
+     case (2)
+      print *,  'Using pt2-matching weight in selection'
+      selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states)
+
+     case (3)
+      print *,  'Using variance-matching weight in selection'
+      selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states)
+
+     case (4)
+      print *,  'Using variance- and pt2-matching weights in selection'
+      selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) * pt2_match_weight(1:N_states)
+
+     case (5)
+      print *,  'Using variance-matching weight in selection'
+      selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states)
+
+    end select
+
 END_PROVIDER


@ -621,11 +714,13 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
        variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi
        norm(istate) = norm(istate) + coef * coef

-!        if (h0_type == "Variance") then
-!          sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
-!        else
+        if (weight_selection /= 5) then
+          ! Energy selection
          sum_e_pert = sum_e_pert + e_pert * selection_weight(istate)
-!        endif
+        else
+          ! Variance selection
+          sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
+        endif
      end do
      if(pseudo_sym)then
       if(dabs(mat(1, p1, p2)).lt.thresh_sym)then 
--- a/src/cipsi/slave_cipsi.irp.f
+++ b/src/cipsi/slave_cipsi.irp.f
@ -17,7 +17,7 @@ subroutine provide_everything
  PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag
  PROVIDE pt2_e0_denominator mo_num N_int ci_energy mpi_master zmq_state zmq_context
  PROVIDE psi_det psi_coef threshold_generators state_average_weight
-  PROVIDE N_det_selectors pt2_stoch_istate N_det
+  PROVIDE N_det_selectors pt2_stoch_istate N_det selection_weight pseudo_sym
 end

 subroutine run_slave_main
@ -220,8 +220,12 @@ subroutine run_slave_main
        call mpi_print('zmq_get_dvector state_average_weight')
      IRP_ENDIF
      if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle
+      IRP_IF MPI_DEBUG
+        call mpi_print('zmq_get_dvector selection_weight')
+      IRP_ENDIF
+      if (zmq_get_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) cycle
      pt2_e0_denominator(1:N_states) = energy(1:N_states)
-      SOFT_TOUCH pt2_e0_denominator state_average_weight pt2_stoch_istate threshold_generators
+      SOFT_TOUCH pt2_e0_denominator state_average_weight pt2_stoch_istate threshold_generators selection_weight

      call wall_time(t1)
      call write_double(6,(t1-t0),'Broadcast time')
--- a/src/cipsi/stochastic_cipsi.irp.f
+++ b/src/cipsi/stochastic_cipsi.irp.f
@ -4,7 +4,7 @@ subroutine run_stochastic_cipsi
 ! Selected Full Configuration Interaction with Stochastic selection and PT2.
  END_DOC
  integer                        :: i,j,k
-  double precision, allocatable  :: pt2(:), variance(:), norm(:), rpt2(:)
+  double precision, allocatable  :: pt2(:), variance(:), norm(:), rpt2(:), zeros(:)
  integer                        :: to_select
  logical, external :: qp_stop

@ -18,7 +18,7 @@ subroutine run_stochastic_cipsi
  rss = memory_of_double(N_states)*4.d0
  call check_mem(rss,irp_here)

-  allocate (pt2(N_states), rpt2(N_states), norm(N_states), variance(N_states))
+  allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states))

  double precision               :: hf_energy_ref
  logical                        :: has
@ -26,6 +26,7 @@ subroutine run_stochastic_cipsi

  relative_error=PT2_relative_error

+  zeros = 0.d0
  pt2 = -huge(1.e0)
  rpt2 = -huge(1.e0)
  norm = 0.d0
@ -63,14 +64,14 @@ subroutine run_stochastic_cipsi

  do while (                                                         &
        (N_det < N_det_max) .and.                                    &
-        (maxval(abs(pt2(1:N_states))) > pt2_max) .and.               &
+        (maxval(abs(rpt2(1:N_states))) > pt2_max) .and.               &
        (maxval(abs(variance(1:N_states))) > variance_max) .and.     &
        (correlation_energy_ratio <= correlation_energy_ratio_max)   &
        )
      write(*,'(A)')  '--------------------------------------------------------------------------------'


-    to_select = N_det*int(sqrt(dble(N_states)))
+    to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
    to_select = max(N_states_diag, to_select)

    pt2 = 0.d0
@ -79,17 +80,17 @@ subroutine run_stochastic_cipsi
    call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, &
      norm, to_select) ! Stochastic PT2 and selection

+    do k=1,N_states
+      rpt2(k) = pt2(k)/(1.d0 + norm(k))
+    enddo
+
    correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref)  /     &
-                    (psi_energy_with_nucl_rep(1) + pt2(1) - hf_energy_ref)
+                    (psi_energy_with_nucl_rep(1) + rpt2(1) - hf_energy_ref)
    correlation_energy_ratio = min(1.d0,correlation_energy_ratio)

-    call save_energy(psi_energy_with_nucl_rep, rpt2)
    call write_double(6,correlation_energy_ratio, 'Correlation ratio')
    call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)

-    do k=1,N_states
-      rpt2(:) = pt2(:)/(1.d0 + norm(k))
-    enddo
    call save_energy(psi_energy_with_nucl_rep, rpt2)

    call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)
@ -108,8 +109,7 @@ subroutine run_stochastic_cipsi

    call diagonalize_CI
    call save_wavefunction
-    rpt2(:) = 0.d0
-    call save_energy(psi_energy_with_nucl_rep, rpt2)
+    call save_energy(psi_energy_with_nucl_rep, zeros)
    if (qp_stop()) exit 
  enddo

@ -117,20 +117,18 @@ subroutine run_stochastic_cipsi
    if (N_det < N_det_max) then
        call diagonalize_CI
        call save_wavefunction
-        rpt2(:) = 0.d0
-        call save_energy(psi_energy_with_nucl_rep, rpt2)
+        call save_energy(psi_energy_with_nucl_rep, zeros)
    endif

-    pt2 = 0.d0
-    variance = 0.d0
-    norm = 0.d0
+    pt2(:) = 0.d0
+    variance(:) = 0.d0
+    norm(:) = 0.d0
    call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, &
      norm,0) ! Stochastic PT2

    do k=1,N_states
-      rpt2(:) = pt2(:)/(1.d0 + norm(k))
+      rpt2(k) = pt2(k)/(1.d0 + norm(k))
    enddo
-    call save_energy(psi_energy_with_nucl_rep, rpt2)

    call save_energy(psi_energy_with_nucl_rep, rpt2)
    call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
--- a/src/cipsi/zmq_selection.irp.f
+++ b/src/cipsi/zmq_selection.irp.f
@ -21,7 +21,8 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
    PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
    PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
    PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
-    PROVIDE psi_bilinear_matrix_transp_order
+    PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym
+

    call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection')

@ -45,6 +46,9 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
    if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then
      stop 'Unable to put state_average_weight on ZMQ server'
    endif
+    if (zmq_put_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) then
+      stop 'Unable to put selection_weight on ZMQ server'
+    endif
    if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) then
      stop 'Unable to put threshold_generators on ZMQ server'
    endif
@ -85,7 +89,11 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
  endif

  integer :: nproc_target
-  nproc_target = nproc
+  if (N_det < 3*nproc) then
+    nproc_target = N_det/4
+  else
+    nproc_target = nproc
+  endif
  double precision :: mem
  mem = 8.d0 * N_det * (N_int * 2.d0 * 3.d0 +  3.d0 + 5.d0) / (1024.d0**3)
  call write_double(6,mem,'Estimated memory/thread (Gb)')
@ -131,13 +139,7 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
    norm(k) = norm(k) * f(k)
  enddo

-  ! Adjust PT2 weights for next selection
-  double precision :: pt2_avg
-  pt2_avg = sum(pt2) / dble(N_states)
-  do k=1,N_states
-    pt2_match_weight(k) *= (pt2(k)/pt2_avg)**2
-  enddo
-  SOFT_TOUCH pt2_match_weight
+  call update_pt2_and_variance_weights(pt2, variance, norm, N_states)

 end subroutine

@ -159,9 +161,9 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2, variance, norm)
  integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull
  type(selection_buffer), intent(inout) :: b
  integer, intent(in)            :: N
-  double precision, intent(inout)    :: pt2(N_states)
-  double precision, intent(inout)    :: variance(N_states)
-  double precision, intent(inout)    :: norm(N_states)
+  double precision, intent(out)      :: pt2(N_states)
+  double precision, intent(out)      :: variance(N_states)
+  double precision, intent(out)      :: norm(N_states)
  double precision                   :: pt2_mwen(N_states)
  double precision                   :: variance_mwen(N_states)
  double precision                   :: norm_mwen(N_states)
--- a/src/davidson/print_e_components.irp.f
+++ b/src/davidson/print_e_components.irp.f
@ -0,0 +1,54 @@
+subroutine print_energy_components()
+  implicit none
+  BEGIN_DOC
+! Prints the different components of the energy.
+  END_DOC
+  integer, save :: ifirst = 0
+  double precision :: Vee, Ven, Vnn, Vecp, T, f
+  integer  :: i,j,k
+
+  Vnn = nuclear_repulsion
+
+  print *, 'Energy components'
+  print *, '================='
+  print *, ''
+  do k=1,N_states
+
+    Ven  = 0.d0
+    Vecp = 0.d0
+    T    = 0.d0
+
+    do j=1,mo_num
+      do i=1,mo_num
+        f = one_e_dm_mo_alpha(i,j,k) + one_e_dm_mo_beta(i,j,k)
+        Ven  = Ven  + f * mo_integrals_n_e(i,j)
+        Vecp = Vecp + f * mo_pseudo_integrals(i,j)
+        T    = T    + f * mo_kinetic_integrals(i,j)
+      enddo
+    enddo
+    Vee = psi_energy(k) - Ven - Vecp - T
+    
+    if (ifirst == 0) then
+      ifirst = 1
+      print *, 'Vnn  : Nucleus-Nucleus   potential energy'
+      print *, 'Ven  : Electron-Nucleus  potential energy'
+      print *, 'Vee  : Electron-Electron potential energy'
+      print *, 'Vecp : Potential energy of the pseudo-potentials'
+      print *, 'T    : Electronic kinetic energy'
+      print *, ''
+    endif
+
+    print *, 'State ', k
+    print *, '---------'
+    print *, ''
+    print *, 'Vnn  = ', Vnn
+    print *, 'Ven  = ', Ven
+    print *, 'Vee  = ', Vee
+    print *, 'Vecp = ', Vecp
+    print *, 'T    = ', T
+    print *, ''
+  enddo
+
+  print *, ''
+
+end
--- a/src/determinants/EZFIO.cfg
+++ b/src/determinants/EZFIO.cfg
@ -28,12 +28,18 @@ doc: Force the wave function to be an eigenfunction of |S^2|
 interface: ezfio,provider,ocaml
 default: True

-[used_weight]
+[weight_one_e_dm]
 type: integer
 doc: Weight used in the calculation of the one-electron density matrix. 0: 1./(c_0^2), 1: 1/N_states, 2: input state-average weight, 3: 1/(Norm_L3(Psi))
 interface: ezfio,provider,ocaml
 default: 1

+[weight_selection]
+type: integer
+doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: rPT2 matching, 3: variance matching, 4: variance and rPT2 matching, 5: variance minimization and matching                                  
+interface: ezfio,provider,ocaml
+default: 2
+
 [threshold_generators]
 type: Threshold
 doc: Thresholds on generators (fraction of the square of the norm) 
@ -89,6 +95,11 @@ doc: Weight of the states in state-average calculations.
 interface: ezfio
 size: (determinants.n_states)

+[selection_factor]
+type: double precision
+doc: f such that the number of determinants to add is f * N_det * sqrt(N_states)
+interface: ezfio,provider,ocaml
+default: 1.

 [thresh_sym]
 type: Threshold
--- a/src/determinants/density_matrix.irp.f
+++ b/src/determinants/density_matrix.irp.f
@ -305,9 +305,9 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
   logical                        :: exists

   state_average_weight(:) = 1.d0
-   if (used_weight == 0) then
+   if (weight_one_e_dm == 0) then
     state_average_weight(:) = c0_weight(:)
-   else if (used_weight == 1) then
+   else if (weight_one_e_dm == 1) then
     state_average_weight(:) = 1./N_states
   else
     call ezfio_has_determinants_state_average_weight(exists)
--- a/src/determinants/utils.irp.f
+++ b/src/determinants/utils.irp.f
@ -43,4 +43,3 @@ BEGIN_PROVIDER [ double precision, S2_matrix_all_dets,(N_det,N_det) ]
 !$OMP END PARALLEL DO
 END_PROVIDER

-
--- a/src/dft_utils_in_r/ao_in_r.irp.f
+++ b/src/dft_utils_in_r/ao_in_r.irp.f
@ -121,3 +121,26 @@
 enddo
 END_PROVIDER

+
+ BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom, (ao_num,n_pts_max_per_atom,nucl_num)]
+&BEGIN_PROVIDER[double precision, aos_in_r_array_per_atom_transp, (n_pts_max_per_atom,ao_num,nucl_num)]
+ implicit none
+ BEGIN_DOC
+ ! aos_in_r_array_per_atom(i,j,k)        = value of the ith ao on the jth grid point attached on the kth atom 
+ END_DOC
+ integer :: i,j,k
+ double precision :: aos_array(ao_num), r(3)
+ do k = 1, nucl_num
+  do i = 1, n_pts_per_atom(k)
+   r(1) = final_grid_points_per_atom(1,i,k)
+   r(2) = final_grid_points_per_atom(2,i,k)
+   r(3) = final_grid_points_per_atom(3,i,k)
+   call give_all_aos_at_r(r,aos_array)
+   do j = 1, ao_num
+    aos_in_r_array_per_atom(j,i,k) = aos_array(j)
+    aos_in_r_array_per_atom_transp(i,j,k) = aos_array(j)
+   enddo
+  enddo
+ enddo
+ END_PROVIDER
+
--- a/src/fci/pt2.irp.f
+++ b/src/fci/pt2.irp.f
@ -46,7 +46,7 @@ subroutine run
  call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, &
     norm,0) ! Stochastic PT2
  do k=1,N_states
-    rpt2(:) = pt2(:)/(1.d0 + norm(k))
+    rpt2(k) = pt2(k)/(1.d0 + norm(k))
  enddo

  call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
--- a/src/iterations/print_summary.irp.f
+++ b/src/iterations/print_summary.irp.f
@ -31,18 +31,19 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_

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

+  call print_energy_components()
+
 end subroutine

--- a/src/nuclei/atomic_radii.irp.f
+++ b/src/nuclei/atomic_radii.irp.f
@ -66,7 +66,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom, (nucl_num)]
 implicit none
 integer :: i
 do i = 1, nucl_num
-  slater_bragg_radii_per_atom(i) = slater_bragg_radii(int(nucl_charge(i)))
+  slater_bragg_radii_per_atom(i) = slater_bragg_radii(max(1,int(nucl_charge(i))))
 enddo
 END_PROVIDER

@ -74,7 +74,7 @@ BEGIN_PROVIDER [double precision, slater_bragg_radii_per_atom_ua, (nucl_num)]
 implicit none
 integer :: i
 do i = 1, nucl_num
-  slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(int(nucl_charge(i)))
+  slater_bragg_radii_per_atom_ua(i) = slater_bragg_radii_ua(max(1,int(nucl_charge(i))))
 enddo
 END_PROVIDER
Author	SHA1	Message	Date
Anthony Scemama	231b071d78	Merge pull request #45 from QuantumPackage/dev-lcpq Dev lcpq	2019-06-08 16:24:33 +02:00
Anthony Scemama	22ab996d00	Merge branch 'anbenali-Fix_Bug_ECP'	2019-06-08 15:52:55 +02:00
Anthony Scemama	00375418e5	Merge branch 'Fix_Bug_ECP' of https://github.com/anbenali/qp2 into anbenali-Fix_Bug_ECP	2019-06-08 15:52:33 +02:00
Anthony Scemama	cf4f9990bc	Updated citation.ccf	2019-06-08 15:44:18 +02:00
Anthony Scemama	c7dedf49ed	More efficient multi-state selection	2019-06-08 15:14:57 +02:00
Anthony Scemama	4c5759fe6f	Modified README	2019-06-08 15:14:38 +02:00
Anthony Scemama	c7ef5fc8c7	Fixed tests	2019-06-05 18:57:17 +02:00
Anthony Scemama	3c6b417ae7	Type conversion	2019-06-05 18:00:04 +02:00
Anthony Scemama	9c028127a0	Updated configure	2019-06-05 17:48:19 +02:00
Anthony Scemama	ab10efa313	Merge branch 'dev-lcpq' of github.com:QuantumPackage/qp2 into dev-lcpq	2019-06-05 17:39:13 +02:00
Anthony Scemama	f36629173d	Updated configure	2019-06-05 17:39:01 +02:00
Anthony Scemama	8a127a9bab	added selection_weight to fci slave	2019-06-05 17:34:36 +02:00
Anthony Scemama	8cdcd08579	Merge branch 'dev-lcpq' of github.com:QuantumPackage/qp2 into dev-lcpq	2019-06-05 17:28:46 +02:00
Anthony Scemama	720a63aadf	Bubblewrap and libcap install scripts	2019-06-05 17:28:15 +02:00
Anthony Scemama	b71579ab43	Moved print_e_components	2019-06-05 16:17:53 +02:00
Anthony Scemama	7c285bddf3	Printing of components of energy:	2019-06-05 15:38:05 +02:00
Anthony Scemama	453cfa0b65	Improved weights	2019-06-05 15:07:36 +02:00
Anthony Scemama	04ca07a540	Modified weighted selection for exp	2019-06-04 18:13:52 +02:00
Anthony Scemama	f514ba0acd	Fixed bug in deterministic variance	2019-06-04 18:11:50 +02:00
Anthony Scemama	bd9816d6e3	Fixed ZMQ stalling because of pseudo_sym	2019-06-04 18:10:56 +02:00
Anthony Scemama	6455b281ff	Compiles	2019-06-04 11:56:06 +02:00
Anthony Scemama	a1c9ec2503	Merge branch 'dev-lcpq' of github.com:QuantumPackage/qp2 into dev-lcpq	2019-06-04 11:39:11 +02:00
Anthony Scemama	4a72ca6b12	Added switch for multiple selection weights, including variance	2019-06-04 11:20:00 +02:00
Anthony Scemama	ce0a5f4e70	Added selection factor	2019-06-04 11:16:20 +02:00
Anthony Scemama	e53e7585f9	Fixed rPT2 small bug	2019-06-04 11:15:43 +02:00
Anthony Scemama	fedc20dc31	rPT2 matching instead of pt2 matching	2019-06-04 11:14:42 +02:00
anbenali	70536970b0	Fix small format issue in ECP	2019-05-29 12:03:06 -05:00
Anthony Scemama	6d15bf2c6d	Merge branch 'dev-lcpq' of github.com:QuantumPackage/qp2 into dev-lcpq	2019-05-28 18:49:34 +02:00
Anthony Scemama	96c17686b4	fixed bugs with dummy atom and becke grid	2019-05-28 18:49:21 +02:00
Anthony Scemama	982855eeb5	AVX work with AMD	2019-05-28 15:39:11 +02:00
Anthony Scemama	d50cca6f88	Updated atomic data	2019-05-28 15:39:11 +02:00
Anthony Scemama	714d53363e	Update documentation and qp_set_frozne_core	2019-05-28 15:39:11 +02:00
Anthony Scemama	7bc6c7e709	Fixed tests	2019-05-28 15:39:11 +02:00
Anthony Scemama	926378c1bc	Fixed tests	2019-05-28 15:39:11 +02:00
Anthony Scemama	f74e57ddef	Fixing travis	2019-05-28 15:39:11 +02:00
Anthony Scemama	b86a5ba963	research.bib	2019-05-28 15:39:09 +02:00
Anthony Scemama	3d70b0f9c2	Deterministic PT2 for < 1000 dets	2019-05-28 15:39:09 +02:00
Anouar Benali	f133cb1e9f	NCSU ECP and Basis set (#41 ) * Adding Basis Sets and ECP from NCSU. When available He core ECP are chosen instead of large core. This applies to Na Mg Al Si P S Cl and Ar * ECP from ncsu * Update Readme file with description of Basis Set	2019-05-28 15:39:09 +02:00
Anthony Scemama	a4e328d40d	Added EPYC config file	2019-05-28 15:39:09 +02:00
Anthony Scemama	e98f745dba	Removed install and uninstall from create_executables_list.sh	2019-05-28 15:39:09 +02:00
Barthelemy Pradines LCT	da51f0960b	Merge branch 'dev-lct' of https://github.com/QuantumPackage/qp2 into dev-lct	2019-05-27 15:17:26 +02:00
Barthelemy Pradines LCT	6584a0c707	added missing dependency inao_two_e_erf_ints/two_e_integrals_erf.irp.f	2019-05-27 15:17:10 +02:00
anbenali	5d951efd51	Update Readme file with description of Basis Set	2019-05-24 13:55:43 -05:00
anbenali	0b3cd3fdef	ECP from ncsu	2019-05-24 12:19:23 -05:00
anbenali	83a9ce67db	Adding Basis Sets and ECP from NCSU. When available He core ECP are chosen instead of large core. This applies to Na Mg Al Si P S Cl and Ar	2019-05-24 12:18:48 -05:00
Pierre Loos	5e732b0f61	modified stuffs	2019-05-23 16:06:37 +02:00
Emmanuel Giner LCT	35cdb13bd4	added some stuffs for getting the bielec integrals	2019-05-20 08:54:39 +02:00
Emmanuel Giner LCT	a7c1a04840	added get integral ao with two index only	2019-05-16 12:18:43 +02:00