diff --git a/CITATION.cff b/CITATION.cff index 9cd186d2..a391b8fa 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -1,32 +1,100 @@ # YAML 1.2 # Metadata for citation of this software according to the CFF format (https://citation-file-format.github.io/) cff-version: 1.0.3 -message: If you use this software, please cite it using these metadata. +message: "If you use this software, please cite it using these metadata." title: Quantum Package -doi: 10.5281/zenodo.825872 +doi: 10.1021/acs.jctc.9b00176 authors: -- given-names: Anthony - family-names: Scemama - affiliation: Laboratoire de Chimie et Physique Quantiques / CNRS - given-names: Yann family-names: Garniron - affiliation: Laboratoire de Chimie et Physique Quantiques / CNRS -- given-names: Michel - family-names: Caffarel - affiliation: Laboratoire de Chimie et Physique Quantiques / CNRS + affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France - given-names: Thomas family-names: Applencourt - affiliation: Argonne National Lab + affiliation: Computational Science Division, Argonne National Laboratory, Argonne, Illinois 60439, United States - given-names: Kevin family-names: Gasperich - affiliation: Argonne National Lab + affiliation: Department of Chemistry, University of Pittsburgh, Pittsburgh, Pennsylvania 15260, United States - given-names: Anouar family-names: Benali - affiliation: Argonne National Lab + affiliation: Computational Science Division, Argonne National Laboratory, Argonne, Illinois 60439, United States +- given-names: Anthony + family-names: Ferté + affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France +- given-names: Julien + family-names: Paquier + affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France +- given-names: Barthélémy + family-names: Pradines + affiliation: Institut des Sciences du Calcul et des Données, Sorbonne Université, F-75005 Paris, France +- given-names: Roland + family-names: Assaraf + affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France +- given-names: Peter + family-names: Reinhardt + affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France +- given-names: Julien + family-names: Toulouse + affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France +- given-names: Pierrette + family-names: Barbaresco + affiliation: CALMIP, Université de Toulouse, CNRS, INPT, INSA, UPS, UMS 3667, Toulouse, France +- given-names: Nicolas + family-names: Renon + affiliation: CALMIP, Université de Toulouse, CNRS, INPT, INSA, UPS, UMS 3667, Toulouse, France +- given-names: Grégoire + family-names: David + affiliation: Aix-Marseille Univ, CNRS, ICR, Marseille, France +- given-names: Jean-Paul + family-names: Malrieu + affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France +- given-names: Mickaël + family-names: Véril + affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France +- given-names: Michel + family-names: Caffarel + affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France +- given-names: Pierre-François + family-names: Loos + affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France - given-names: Emmanuel family-names: Giner - affiliation: Laboratoire de Chimie Theorique / CNRS + affiliation: Laboratoire de Chimie Théorique, Sorbonne Université, CNRS, Paris, France +- given-names: Anthony + family-names: Scemama + affiliation: Laboratoire de Chimie et Physique Quantiques (UMR 5626), Université de Toulouse, CNRS, UPS, Toulouse, France +abstract: "Quantum chemistry is a discipline which relies heavily on very +expensive numerical computations. The scaling of correlated wave function +methods lies, in their standard implementation, between O(N^5) and O(exp(N)), +where N is proportional to the system size. Therefore, performing accurate +calculations on chemically meaningful systems requires (i) approximations that +can lower the computational scaling and (ii) efficient implementations that +take advantage of modern massively parallel architectures. Quantum Package is +an open-source programming environment for quantum chemistry specially designed +for wave function methods. Its main goal is the development of +determinant-driven selected configuration interaction (sCI) methods and +multireference second-order perturbation theory (PT2). The determinant-driven +framework allows the programmer to include any arbitrary set of determinants in +the reference space, hence providing greater methodological freedom. The sCI +method implemented in Quantum Package is based on the CIPSI (Configuration +Interaction using a Perturbative Selection made Iteratively) algorithm which +complements the variational sCI energy with a PT2 correction. Additional +external plugins have been recently added to perform calculations with +multireference coupled cluster theory and range-separated density-functional +theory. All the programs are developed with the IRPF90 code generator, which +simplifies collaborative work and the development of new features. Quantum +Package strives to allow easy implementation and experimentation of new +methods, while making parallel computation as simple and efficient as possible +on modern supercomputer architectures. Currently, the code enables, routinely, +to realize runs on roughly 2 000 CPU cores, with tens of millions of +determinants in the reference space. Moreover, we have been able to push up to +12 288 cores in order to test its parallel efficiency. In the present +manuscript, we also introduce some key new developments: (i) a renormalized +second-order perturbative correction for efficient extrapolation to the full CI +limit and (ii) a stochastic version of the CIPSI selection performed +simultaneously to the PT2 calculation at no extra cost." version: '2.0' -date-released: 2019-02-11 +url: https://quantumpackage.github.io/qp2/ +date-released: 2019-05-13 repository-code: https://github.com/QuantumPackage/qp2 +keywords: [ "computational chemistry", "configuration interaction", "cipsi", "perturbation theory" ] license: AGPL-3.0-or-later diff --git a/INSTALL.rst b/INSTALL.rst index cc5b512a..baf7f051 100644 --- a/INSTALL.rst +++ b/INSTALL.rst @@ -45,6 +45,8 @@ Requirements - |ZeroMQ| : networking library - `GMP `_ : Gnu Multiple Precision Arithmetic Library - |OCaml| compiler with |OPAM| package manager +- `Bubblewrap `_ : Sandboxing tool required by Opam +- `libcap https://git.kernel.org/pub/scm/linux/kernel/git/morgan/libcap.git`_ : POSIX capabilities required by Bubblewrap - |Ninja| : a parallel build system @@ -86,6 +88,8 @@ The following packages are supported by the :command:`configure` installer: * zeromq * f77zmq * gmp +* libcap +* bwrap * ocaml ( :math:`\approx` 10 minutes) * ezfio * docopt @@ -243,6 +247,55 @@ With Debian or Ubuntu, you can use sudo apt install libgmp-dev +libcap +------ + +Libcap is a library for getting and setting POSIX.1e draft 15 capabilities. + +* Download the latest version of libcap here: + ``_ + and move it in the :file:`${QP_ROOT}/external` directory + +* Extract the archive, go into the :file:`libcap-*/libcap` directory and run + the following command + +.. code:: bash + + prefix=$QP_ROOT make install + +With Debian or Ubuntu, you can use + +.. code:: bash + + sudo apt install libcap-dev + + +Bubblewrap +---------- + +Bubblewrap is an unprivileged sandboxing tool. + +* Download Bubblewrap here: + ``_ + and move it in the :file:`${QP_ROOT}/external` directory + +* Extract the archive, go into the :file:`bubblewrap-*` directory and run + the following commands + +.. code:: bash + + ./configure --prefix=$QP_ROOT && make -j 8 + make install-exec-am + + +With Debian or Ubuntu, you can use + +.. code:: bash + + sudo apt install bubblewrap + + + OCaml ----- diff --git a/README.md b/README.md index c94adb69..3e85b4cc 100644 --- a/README.md +++ b/README.md @@ -1,12 +1,13 @@ # Quantum Package 2.0 -*Quantum package 2.0: an open-source determinant-driven suite of programs*\ +[*Quantum package 2.0: an open-source determinant-driven suite of programs*](https://pubs.acs.org/doi/10.1021/acs.jctc.9b00176)\ Y. Garniron, K. Gasperich, T. Applencourt, A. Benali, A. Ferté, J. Paquier, B. Pradines, R. Assaraf, P. Reinhardt, J. Toulouse, P. Barbaresco, N. Renon, G. David, J. P. Malrieu, M. Véril, M. Caffarel, P. F. Loos, E. Giner and A. Scemama\ +J. Chem. Theory Comput. (2019)\ https://arxiv.org/abs/1902.08154 -![QP](https://raw.githubusercontent.com/QuantumPackage/qp2/master/data/qp2.png) + # Getting started diff --git a/configure b/configure index ffda5016..c343c243 100755 --- a/configure +++ b/configure @@ -175,7 +175,7 @@ if [[ "${PACKAGES}.x" != ".x" ]] ; then fi if [[ ${PACKAGES} = all ]] ; then - PACKAGES="zlib ninja irpf90 zeromq f77zmq gmp ocaml ezfio docopt resultsFile bats" + PACKAGES="zlib ninja irpf90 zeromq f77zmq gmp libcap bwrap ocaml ezfio docopt resultsFile bats" fi @@ -206,6 +206,32 @@ EOF make install EOF + elif [[ ${PACKAGE} = libcap ]] ; then + + download \ + "https://git.kernel.org/pub/scm/linux/kernel/git/morgan/libcap.git/snapshot/libcap-2.25.tar.gz" \ + "${QP_ROOT}"/external/libcap.tar.gz + execute << EOF + cd "\${QP_ROOT}"/external + tar --gunzip --extract --file libcap.tar.gz + rm libcap.tar.gz + cd libcap-*/libcap + prefix=$QP_ROOT make install +EOF + + elif [[ ${PACKAGE} = bwrap ]] ; then + + download \ + "https://github.com/projectatomic/bubblewrap/releases/download/v0.3.3/bubblewrap-0.3.3.tar.xz" \ + "${QP_ROOT}"/external/bwrap.tar.xz + execute << EOF + cd "\${QP_ROOT}"/external + tar --xz --extract --file bwrap.tar.xz + rm bwrap.tar.xz + cd bubblewrap* + ./configure --prefix=$QP_ROOT && make -j 8 + make install-exec-am +EOF elif [[ ${PACKAGE} = irpf90 ]] ; then @@ -276,7 +302,7 @@ EOF rm ${QP_ROOT}/external/opam_installer.sh source ${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true - ${QP_ROOT}/bin/opam init --disable-sandboxing --verbose --yes + ${QP_ROOT}/bin/opam init --verbose --yes eval $(${QP_ROOT}/bin/opam env) opam install -y ${OCAML_PACKAGES} || exit 1 @@ -290,7 +316,7 @@ EOF | sh \${QP_ROOT}/external/opam_installer.sh rm \${QP_ROOT}/external/opam_installer.sh source \${OPAMROOT}/opam-init/init.sh > /dev/null 2> /dev/null || true - \${QP_ROOT}/bin/opam init --disable-sandboxing --verbose --yes + \${QP_ROOT}/bin/opam init --verbose --yes eval \$(\${QP_ROOT}/bin/opam env) opam install -y \${OCAML_PACKAGES} || exit 1 EOF @@ -399,6 +425,18 @@ if [[ ${ZLIB} = $(not_found) ]] ; then fail fi +BWRAP=$(find_exe bwrap) +if [[ ${BWRAP} = $(not_found) ]] ; then + error "Bubblewrap (bwrap) is not installed." + fail +fi + +LIBCAP=$(find_lib -lcap) +if [[ ${LIBCAP} = $(not_found) ]] ; then + error "Libcap (libcap) is not installed." + fail +fi + OPAM=$(find_exe opam) if [[ ${OPAM} = $(not_found) ]] ; then error "OPAM (ocaml) package manager is not installed." diff --git a/ocaml/Input_determinants_by_hand.ml b/ocaml/Input_determinants_by_hand.ml index ab6fb2ca..e4c6ff2a 100644 --- a/ocaml/Input_determinants_by_hand.ml +++ b/ocaml/Input_determinants_by_hand.ml @@ -81,9 +81,6 @@ end = struct ;; let write_n_det n = - let n_det_old = - Ezfio.get_determinants_n_det () - in Det_number.to_int n |> Ezfio.set_determinants_n_det ;; diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index f0cab384..7e292d6e 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -5,7 +5,7 @@ subroutine run_cipsi ! stochastic PT2. END_DOC integer :: i,j,k - double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:) + double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:), zeros(:) integer :: n_det_before, to_select double precision :: rss @@ -13,7 +13,7 @@ subroutine run_cipsi rss = memory_of_double(N_states)*4.d0 call check_mem(rss,irp_here) - allocate (pt2(N_states), rpt2(N_states), norm(N_states), variance(N_states)) + allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states)) double precision :: hf_energy_ref logical :: has @@ -23,10 +23,11 @@ subroutine run_cipsi relative_error=PT2_relative_error + zeros = 0.d0 pt2 = -huge(1.e0) rpt2 = -huge(1.e0) norm = 0.d0 - variance = 0.d0 + variance = huge(1.e0) if (s2_eig) then call make_s2_eigenfunction @@ -65,7 +66,8 @@ subroutine run_cipsi do while ( & (N_det < N_det_max) .and. & - (maxval(abs(pt2(1:N_states))) > pt2_max) .and. & + (maxval(abs(rpt2(1:N_states))) > pt2_max) .and. & + (maxval(variance(1:N_states)) > variance_max) .and. & (correlation_energy_ratio <= correlation_energy_ratio_max) & ) write(*,'(A)') '--------------------------------------------------------------------------------' @@ -83,17 +85,17 @@ subroutine run_cipsi SOFT_TOUCH threshold_generators endif + do k=1,N_states + rpt2(k) = pt2(k)/(1.d0 + norm(k)) + enddo correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / & - (psi_energy_with_nucl_rep(1) + pt2(1) - hf_energy_ref) + (psi_energy_with_nucl_rep(1) + rpt2(1) - hf_energy_ref) correlation_energy_ratio = min(1.d0,correlation_energy_ratio) call write_double(6,correlation_energy_ratio, 'Correlation ratio') call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) - do k=1,N_states - rpt2(:) = pt2(:)/(1.d0 + norm(k)) - enddo call save_energy(psi_energy_with_nucl_rep, rpt2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) @@ -103,9 +105,8 @@ subroutine run_cipsi if (qp_stop()) exit n_det_before = N_det - to_select = N_det + to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor) to_select = max(N_states_diag, to_select) -! to_select = min(to_select, N_det_max-n_det_before) call ZMQ_selection(to_select, pt2, variance, norm) PROVIDE psi_coef @@ -114,32 +115,30 @@ subroutine run_cipsi call diagonalize_CI call save_wavefunction - rpt2(:) = 0.d0 - call save_energy(psi_energy_with_nucl_rep, rpt2) + call save_energy(psi_energy_with_nucl_rep, zeros) if (qp_stop()) exit +print *, (N_det < N_det_max) +print *, (maxval(abs(rpt2(1:N_states))) > pt2_max) +print *, (maxval(variance(1:N_states)) > variance_max) +print *, (correlation_energy_ratio <= correlation_energy_ratio_max) enddo if (.not.qp_stop()) then if (N_det < N_det_max) then call diagonalize_CI call save_wavefunction - rpt2(:) = 0.d0 - call save_energy(psi_energy_with_nucl_rep, rpt2) + call save_energy(psi_energy_with_nucl_rep, zeros) endif if (do_pt2) then - pt2 = 0.d0 - variance = 0.d0 - norm = 0.d0 + pt2(:) = 0.d0 + variance(:) = 0.d0 + norm(:) = 0.d0 threshold_generators = 1d0 SOFT_TOUCH threshold_generators call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, & norm,0) ! Stochastic PT2 SOFT_TOUCH threshold_generators - do k=1,N_states - rpt2(:) = pt2(:)/(1.d0 + norm(k)) - enddo - call save_energy(psi_energy_with_nucl_rep, pt2) endif print *, 'N_det = ', N_det print *, 'N_sop = ', N_occ_pattern @@ -148,12 +147,11 @@ subroutine run_cipsi do k=1,N_states - rpt2(:) = pt2(:)/(1.d0 + norm(k)) + rpt2(k) = pt2(k)/(1.d0 + norm(k)) enddo - call save_energy(psi_energy_with_nucl_rep, rpt2) call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) - call save_energy(psi_energy_with_nucl_rep, pt2) + call save_energy(psi_energy_with_nucl_rep, rpt2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) call print_extrapolated_energy() endif diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 0e143a9c..9f891320 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -129,13 +129,13 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in) PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted - PROVIDE psi_det_hii N_generators_bitmask + PROVIDE psi_det_hii N_generators_bitmask selection_weight pseudo_sym if (h0_type == 'SOP') then PROVIDE psi_occ_pattern_hii det_to_occ_pattern endif - if (N_det < max(1000,N_states)) then + if (N_det < max(4,N_states)) then pt2=0.d0 variance=0.d0 norm=0.d0 @@ -182,6 +182,9 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in) if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then stop 'Unable to put state_average_weight on ZMQ server' endif + if (zmq_put_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) then + stop 'Unable to put selection_weight on ZMQ server' + endif if (zmq_put_ivector(zmq_to_qp_run_socket,1,'pt2_stoch_istate',pt2_stoch_istate,1) == -1) then stop 'Unable to put pt2_stoch_istate on ZMQ server' endif @@ -333,13 +336,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in) pt2(k) = 0.d0 enddo - ! Adjust PT2 weights for next selection - double precision :: pt2_avg - pt2_avg = sum(pt2) / dble(N_states) - do k=1,N_states - pt2_match_weight(k) *= (pt2(k)/pt2_avg)**2 - enddo - SOFT_TOUCH pt2_match_weight + call update_pt2_and_variance_weights(pt2, variance, norm, N_states) end subroutine diff --git a/src/cipsi/run_selection_slave.irp.f b/src/cipsi/run_selection_slave.irp.f index 480ef12b..c1542445 100644 --- a/src/cipsi/run_selection_slave.irp.f +++ b/src/cipsi/run_selection_slave.irp.f @@ -25,8 +25,8 @@ subroutine run_selection_slave(thread,iproc,energy) PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F - PROVIDE psi_selectors_coef_transp psi_det_sorted + PROVIDE psi_bilinear_matrix_transp_order N_int pt2_F pseudo_sym + PROVIDE psi_selectors_coef_transp psi_det_sorted weight_selection zmq_to_qp_run_socket = new_zmq_to_qp_run_socket() @@ -230,6 +230,8 @@ subroutine pull_selection_results(zmq_socket_pull, pt2, variance, norm, val, det endif else pt2(:) = 0.d0 + variance(:) = 0.d0 + norm(:) = 0.d0 endif rc = f77_zmq_recv( zmq_socket_pull, ntask, 4, 0) diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index eedcd67f..df31bc39 100644 --- a/src/cipsi/selection.irp.f +++ b/src/cipsi/selection.irp.f @@ -6,15 +6,108 @@ BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ] ! Weights adjusted along the selection to make the PT2 contributions ! of each state coincide. END_DOC - pt2_match_weight = 1.d0 + pt2_match_weight(:) = 1.d0 END_PROVIDER +BEGIN_PROVIDER [ double precision, variance_match_weight, (N_states) ] + implicit none + BEGIN_DOC + ! Weights adjusted along the selection to make the variances + ! of each state coincide. + END_DOC + variance_match_weight(:) = 1.d0 +END_PROVIDER + +subroutine update_pt2_and_variance_weights(pt2, variance, norm, N_st) + implicit none + BEGIN_DOC +! Updates the rPT2- and Variance- matching weights. + END_DOC + integer, intent(in) :: N_st + double precision, intent(in) :: pt2(N_st) + double precision, intent(in) :: variance(N_st) + double precision, intent(in) :: norm(N_st) + + double precision :: avg, rpt2(N_st), element, dt, x + integer :: k + integer, save :: i_iter=0 + integer, parameter :: i_itermax = 3 + double precision, allocatable, save :: memo_variance(:,:), memo_pt2(:,:) + + if (i_iter == 0) then + allocate(memo_variance(N_st,i_itermax), memo_pt2(N_st,i_itermax)) + memo_pt2(:,:) = 1.d0 + memo_variance(:,:) = 1.d0 + endif + + i_iter = i_iter+1 + if (i_iter > i_itermax) then + i_iter = 1 + endif + + dt = 4.d0 + + do k=1,N_st + rpt2(k) = pt2(k)/(1.d0 + norm(k)) + enddo + + avg = sum(rpt2(1:N_st)) / dble(N_st) + do k=1,N_st + element = exp(dt*(rpt2(k)/avg -1.d0)) + element = min(1.5d0 , element) + element = max(0.5d0 , element) + memo_pt2(k,i_iter) = element + pt2_match_weight(k) = product(memo_pt2(k,:)) + enddo + + avg = sum(variance(1:N_st)) / dble(N_st) + do k=1,N_st + element = exp(dt*(variance(k)/avg -1.d0)) + element = min(1.5d0 , element) + element = max(0.5d0 , element) + memo_variance(k,i_iter) = element + variance_match_weight(k) = product(memo_variance(k,:)) + enddo + + print *, '# PT2 weight ', real(pt2_match_weight(:),4) + print *, '# var weight ', real(variance_match_weight(:),4) + SOFT_TOUCH pt2_match_weight variance_match_weight +end + + BEGIN_PROVIDER [ double precision, selection_weight, (N_states) ] implicit none BEGIN_DOC ! Weights used in the selection criterion END_DOC - selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) + select case (weight_selection) + + case (0) + print *, 'Using input weights in selection' + selection_weight(1:N_states) = state_average_weight(1:N_states) + + case (1) + print *, 'Using 1/c_max^2 weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) + + case (2) + print *, 'Using pt2-matching weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * pt2_match_weight(1:N_states) + + case (3) + print *, 'Using variance-matching weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) + + case (4) + print *, 'Using variance- and pt2-matching weights in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) * pt2_match_weight(1:N_states) + + case (5) + print *, 'Using variance-matching weight in selection' + selection_weight(1:N_states) = c0_weight(1:N_states) * variance_match_weight(1:N_states) + + end select + END_PROVIDER @@ -621,11 +714,13 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi norm(istate) = norm(istate) + coef * coef -! if (h0_type == "Variance") then -! sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate) -! else + if (weight_selection /= 5) then + ! Energy selection sum_e_pert = sum_e_pert + e_pert * selection_weight(istate) -! endif + else + ! Variance selection + sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate) + endif end do if(pseudo_sym)then if(dabs(mat(1, p1, p2)).lt.thresh_sym)then diff --git a/src/cipsi/slave_cipsi.irp.f b/src/cipsi/slave_cipsi.irp.f index 58b53f94..c9710c18 100644 --- a/src/cipsi/slave_cipsi.irp.f +++ b/src/cipsi/slave_cipsi.irp.f @@ -17,7 +17,7 @@ subroutine provide_everything PROVIDE H_apply_buffer_allocated mo_two_e_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag PROVIDE pt2_e0_denominator mo_num N_int ci_energy mpi_master zmq_state zmq_context PROVIDE psi_det psi_coef threshold_generators state_average_weight - PROVIDE N_det_selectors pt2_stoch_istate N_det + PROVIDE N_det_selectors pt2_stoch_istate N_det selection_weight pseudo_sym end subroutine run_slave_main @@ -220,8 +220,12 @@ subroutine run_slave_main call mpi_print('zmq_get_dvector state_average_weight') IRP_ENDIF if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle + IRP_IF MPI_DEBUG + call mpi_print('zmq_get_dvector selection_weight') + IRP_ENDIF + if (zmq_get_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) cycle pt2_e0_denominator(1:N_states) = energy(1:N_states) - SOFT_TOUCH pt2_e0_denominator state_average_weight pt2_stoch_istate threshold_generators + SOFT_TOUCH pt2_e0_denominator state_average_weight pt2_stoch_istate threshold_generators selection_weight call wall_time(t1) call write_double(6,(t1-t0),'Broadcast time') diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 157479d9..ae2b7519 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -4,7 +4,7 @@ subroutine run_stochastic_cipsi ! Selected Full Configuration Interaction with Stochastic selection and PT2. END_DOC integer :: i,j,k - double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:) + double precision, allocatable :: pt2(:), variance(:), norm(:), rpt2(:), zeros(:) integer :: to_select logical, external :: qp_stop @@ -18,7 +18,7 @@ subroutine run_stochastic_cipsi rss = memory_of_double(N_states)*4.d0 call check_mem(rss,irp_here) - allocate (pt2(N_states), rpt2(N_states), norm(N_states), variance(N_states)) + allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm(N_states), variance(N_states)) double precision :: hf_energy_ref logical :: has @@ -26,6 +26,7 @@ subroutine run_stochastic_cipsi relative_error=PT2_relative_error + zeros = 0.d0 pt2 = -huge(1.e0) rpt2 = -huge(1.e0) norm = 0.d0 @@ -63,14 +64,14 @@ subroutine run_stochastic_cipsi do while ( & (N_det < N_det_max) .and. & - (maxval(abs(pt2(1:N_states))) > pt2_max) .and. & + (maxval(abs(rpt2(1:N_states))) > pt2_max) .and. & (maxval(abs(variance(1:N_states))) > variance_max) .and. & (correlation_energy_ratio <= correlation_energy_ratio_max) & ) write(*,'(A)') '--------------------------------------------------------------------------------' - to_select = N_det*int(sqrt(dble(N_states))) + to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor) to_select = max(N_states_diag, to_select) pt2 = 0.d0 @@ -79,17 +80,17 @@ subroutine run_stochastic_cipsi call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, & norm, to_select) ! Stochastic PT2 and selection + do k=1,N_states + rpt2(k) = pt2(k)/(1.d0 + norm(k)) + enddo + correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / & - (psi_energy_with_nucl_rep(1) + pt2(1) - hf_energy_ref) + (psi_energy_with_nucl_rep(1) + rpt2(1) - hf_energy_ref) correlation_energy_ratio = min(1.d0,correlation_energy_ratio) - call save_energy(psi_energy_with_nucl_rep, rpt2) call write_double(6,correlation_energy_ratio, 'Correlation ratio') call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) - do k=1,N_states - rpt2(:) = pt2(:)/(1.d0 + norm(k)) - enddo call save_energy(psi_energy_with_nucl_rep, rpt2) call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det) @@ -108,8 +109,7 @@ subroutine run_stochastic_cipsi call diagonalize_CI call save_wavefunction - rpt2(:) = 0.d0 - call save_energy(psi_energy_with_nucl_rep, rpt2) + call save_energy(psi_energy_with_nucl_rep, zeros) if (qp_stop()) exit enddo @@ -117,20 +117,18 @@ subroutine run_stochastic_cipsi if (N_det < N_det_max) then call diagonalize_CI call save_wavefunction - rpt2(:) = 0.d0 - call save_energy(psi_energy_with_nucl_rep, rpt2) + call save_energy(psi_energy_with_nucl_rep, zeros) endif - pt2 = 0.d0 - variance = 0.d0 - norm = 0.d0 + pt2(:) = 0.d0 + variance(:) = 0.d0 + norm(:) = 0.d0 call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, & norm,0) ! Stochastic PT2 do k=1,N_states - rpt2(:) = pt2(:)/(1.d0 + norm(k)) + rpt2(k) = pt2(k)/(1.d0 + norm(k)) enddo - call save_energy(psi_energy_with_nucl_rep, rpt2) call save_energy(psi_energy_with_nucl_rep, rpt2) call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) diff --git a/src/cipsi/zmq_selection.irp.f b/src/cipsi/zmq_selection.irp.f index c7c11eec..081d998f 100644 --- a/src/cipsi/zmq_selection.irp.f +++ b/src/cipsi/zmq_selection.irp.f @@ -21,7 +21,8 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm) PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns - PROVIDE psi_bilinear_matrix_transp_order + PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym + call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection') @@ -45,6 +46,9 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm) if (zmq_put_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) then stop 'Unable to put state_average_weight on ZMQ server' endif + if (zmq_put_dvector(zmq_to_qp_run_socket,1,'selection_weight',selection_weight,N_states) == -1) then + stop 'Unable to put selection_weight on ZMQ server' + endif if (zmq_put_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) then stop 'Unable to put threshold_generators on ZMQ server' endif @@ -85,7 +89,11 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm) endif integer :: nproc_target - nproc_target = nproc + if (N_det < 3*nproc) then + nproc_target = N_det/4 + else + nproc_target = nproc + endif double precision :: mem mem = 8.d0 * N_det * (N_int * 2.d0 * 3.d0 + 3.d0 + 5.d0) / (1024.d0**3) call write_double(6,mem,'Estimated memory/thread (Gb)') @@ -131,13 +139,7 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm) norm(k) = norm(k) * f(k) enddo - ! Adjust PT2 weights for next selection - double precision :: pt2_avg - pt2_avg = sum(pt2) / dble(N_states) - do k=1,N_states - pt2_match_weight(k) *= (pt2(k)/pt2_avg)**2 - enddo - SOFT_TOUCH pt2_match_weight + call update_pt2_and_variance_weights(pt2, variance, norm, N_states) end subroutine @@ -159,9 +161,9 @@ subroutine selection_collector(zmq_socket_pull, b, N, pt2, variance, norm) integer(ZMQ_PTR), intent(in) :: zmq_socket_pull type(selection_buffer), intent(inout) :: b integer, intent(in) :: N - double precision, intent(inout) :: pt2(N_states) - double precision, intent(inout) :: variance(N_states) - double precision, intent(inout) :: norm(N_states) + double precision, intent(out) :: pt2(N_states) + double precision, intent(out) :: variance(N_states) + double precision, intent(out) :: norm(N_states) double precision :: pt2_mwen(N_states) double precision :: variance_mwen(N_states) double precision :: norm_mwen(N_states) diff --git a/src/davidson/print_e_components.irp.f b/src/davidson/print_e_components.irp.f new file mode 100644 index 00000000..ddf83474 --- /dev/null +++ b/src/davidson/print_e_components.irp.f @@ -0,0 +1,54 @@ +subroutine print_energy_components() + implicit none + BEGIN_DOC +! Prints the different components of the energy. + END_DOC + integer, save :: ifirst = 0 + double precision :: Vee, Ven, Vnn, Vecp, T, f + integer :: i,j,k + + Vnn = nuclear_repulsion + + print *, 'Energy components' + print *, '=================' + print *, '' + do k=1,N_states + + Ven = 0.d0 + Vecp = 0.d0 + T = 0.d0 + + do j=1,mo_num + do i=1,mo_num + f = one_e_dm_mo_alpha(i,j,k) + one_e_dm_mo_beta(i,j,k) + Ven = Ven + f * mo_integrals_n_e(i,j) + Vecp = Vecp + f * mo_pseudo_integrals(i,j) + T = T + f * mo_kinetic_integrals(i,j) + enddo + enddo + Vee = psi_energy(k) - Ven - Vecp - T + + if (ifirst == 0) then + ifirst = 1 + print *, 'Vnn : Nucleus-Nucleus potential energy' + print *, 'Ven : Electron-Nucleus potential energy' + print *, 'Vee : Electron-Electron potential energy' + print *, 'Vecp : Potential energy of the pseudo-potentials' + print *, 'T : Electronic kinetic energy' + print *, '' + endif + + print *, 'State ', k + print *, '---------' + print *, '' + print *, 'Vnn = ', Vnn + print *, 'Ven = ', Ven + print *, 'Vee = ', Vee + print *, 'Vecp = ', Vecp + print *, 'T = ', T + print *, '' + enddo + + print *, '' + +end diff --git a/src/determinants/EZFIO.cfg b/src/determinants/EZFIO.cfg index b8fc7406..95128969 100644 --- a/src/determinants/EZFIO.cfg +++ b/src/determinants/EZFIO.cfg @@ -28,12 +28,18 @@ doc: Force the wave function to be an eigenfunction of |S^2| interface: ezfio,provider,ocaml default: True -[used_weight] +[weight_one_e_dm] type: integer doc: Weight used in the calculation of the one-electron density matrix. 0: 1./(c_0^2), 1: 1/N_states, 2: input state-average weight, 3: 1/(Norm_L3(Psi)) interface: ezfio,provider,ocaml default: 1 +[weight_selection] +type: integer +doc: Weight used in the selection. 0: input state-average weight, 1: 1./(c_0^2), 2: rPT2 matching, 3: variance matching, 4: variance and rPT2 matching, 5: variance minimization and matching +interface: ezfio,provider,ocaml +default: 2 + [threshold_generators] type: Threshold doc: Thresholds on generators (fraction of the square of the norm) @@ -89,6 +95,11 @@ doc: Weight of the states in state-average calculations. interface: ezfio size: (determinants.n_states) +[selection_factor] +type: double precision +doc: f such that the number of determinants to add is f * N_det * sqrt(N_states) +interface: ezfio,provider,ocaml +default: 1. [thresh_sym] type: Threshold diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f index a9630977..e4f76bca 100644 --- a/src/determinants/density_matrix.irp.f +++ b/src/determinants/density_matrix.irp.f @@ -305,9 +305,9 @@ BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ] logical :: exists state_average_weight(:) = 1.d0 - if (used_weight == 0) then + if (weight_one_e_dm == 0) then state_average_weight(:) = c0_weight(:) - else if (used_weight == 1) then + else if (weight_one_e_dm == 1) then state_average_weight(:) = 1./N_states else call ezfio_has_determinants_state_average_weight(exists) diff --git a/src/determinants/utils.irp.f b/src/determinants/utils.irp.f index 20d9e1e5..3aec16f9 100644 --- a/src/determinants/utils.irp.f +++ b/src/determinants/utils.irp.f @@ -43,4 +43,3 @@ BEGIN_PROVIDER [ double precision, S2_matrix_all_dets,(N_det,N_det) ] !$OMP END PARALLEL DO END_PROVIDER - diff --git a/src/fci/pt2.irp.f b/src/fci/pt2.irp.f index 45860fb5..6135864f 100644 --- a/src/fci/pt2.irp.f +++ b/src/fci/pt2.irp.f @@ -46,7 +46,7 @@ subroutine run call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, & norm,0) ! Stochastic PT2 do k=1,N_states - rpt2(:) = pt2(:)/(1.d0 + norm(k)) + rpt2(k) = pt2(k)/(1.d0 + norm(k)) enddo call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index a8037982..ad87bc8e 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -31,18 +31,19 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_ write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' write(*,fmt) - write(fmt,*) '(12X,', N_states_p, '(6X,A7,1X,I6,10X))' + write(fmt,*) '(13X,', N_states_p, '(6X,A7,1X,I6,10X))' write(*,fmt) ('State',k, k=1,N_states_p) write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))' write(*,fmt) - write(fmt,*) '(A12,', N_states_p, '(1X,F14.8,15X))' + write(fmt,*) '(A13,', N_states_p, '(1X,F14.8,15X))' write(*,fmt) '# E ', e_(1:N_states_p) if (N_states_p > 1) then write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1) write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0 endif write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))' - write(*,fmt) '# PT2'//pt2_string, (pt2_(k), error_(k), k=1,N_states_p) + write(*,fmt) '# PT2 '//pt2_string, (pt2_(k), error_(k), k=1,N_states_p) + write(*,fmt) '# rPT2'//pt2_string, (pt2_(k)*f(k), error_(k)*f(k), k=1,N_states_p) write(*,'(A)') '#' write(*,fmt) '# E+PT2 ', (e_(k)+pt2_(k),error_(k), k=1,N_states_p) write(*,fmt) '# E+rPT2 ', (e_(k)+pt2_(k)*f(k),error_(k)*f(k), k=1,N_states_p) @@ -97,5 +98,7 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_ enddo endif + call print_energy_components() + end subroutine