From c7aa11c14cbc972747a85e4a256ee66ecf20887a Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 5 Nov 2018 16:26:37 +0100 Subject: [PATCH] CITATION.cff --- CITATION.cff | 32 ++++ src/Determinants/EZFIO.cfg | 19 +- src/Determinants/determinants.irp.f | 2 +- src/Determinants/slater_rules.irp.f | 285 +++++++++++++++------------- src/Iterations/io.irp.f | 11 +- src/Iterations/print_summary.irp.f | 16 +- 6 files changed, 211 insertions(+), 154 deletions(-) create mode 100644 CITATION.cff diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 00000000..3eb00f5b --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,32 @@ +# 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. +title: Quantum Package +doi: 10.5281/zenodo.825872 +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 +- given-names: Thomas + family-names: Applencourt + affiliation: Argonne National Lab +- given-names: Kevin + family-names: Gasperich + affiliation: Argonne National Lab +- given-names: Anouar + family-names: Benali + affiliation: Argonne National Lab +- given-names: Emmanuel + family-names: Giner + affiliation: Laboratoire de Chimie Theorique / CNRS +version: '2.0' +date-released: "Not yet released" +repository-code: https://github.com/LCPQ/quantum_package +license: GPL-3.0-or-later diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg index 0dd61828..2baa13b2 100644 --- a/src/Determinants/EZFIO.cfg +++ b/src/Determinants/EZFIO.cfg @@ -4,17 +4,11 @@ doc: Max number of determinants in the wave function interface: ezfio,provider,ocaml default: 1000000 -[N_det_max_property] +[N_det_max_full] type: Det_number_max -doc: Max number of determinants in the wave function when you select for a given property +doc: Maximum number of determinants where H is fully diagonalized interface: ezfio,provider,ocaml -default: 10000 - -[N_det_max_jacobi] -type: Det_number_max -doc: Maximum number of determinants diagonalized by Jacobi -interface: ezfio,provider,ocaml -default: 2000 +default: 1000 [N_states] type: States_number @@ -46,7 +40,6 @@ doc: 0: 1/(c_0^2), 1: 1/N_states, 2: input state-average weight, 3: 1/(Norm_L3(P interface: ezfio,provider,ocaml default: 0 - [threshold_generators] type: Threshold doc: Thresholds on generators (fraction of the square of the norm) @@ -132,3 +125,9 @@ doc: Weight of the states in state-average calculations. interface: ezfio size: (determinants.n_states) +[zero_doubles] +type: logical +doc: If true, set to zero all the matrix elements for two determinants differing by a double excitation. +interface: ezfio,provider,ocaml +default: False + diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 3d8f6b2c..668a4cef 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -5,7 +5,7 @@ BEGIN_PROVIDER [ character*(64), diag_algorithm ] BEGIN_DOC ! Diagonalization algorithm (Davidson or Lapack) END_DOC - if (N_det > N_det_max_jacobi) then + if (N_det > N_det_max_full) then diag_algorithm = "Davidson" else diag_algorithm = "Lapack" diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 9986ee6e..a337fc29 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -532,47 +532,51 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2) integer :: spin select case (degree) case (2) - call get_double_excitation(key_i,key_j,exc,phase,Nint) - ! Mono alpha, mono beta - if (exc(0,1,1) == 1) then - if ( (exc(1,1,1) == exc(1,2,2)).and.(exc(1,1,2) == exc(1,2,1)) ) then - s2 = -phase + if (zero_doubles) then + hij = 0.d0 + else + call get_double_excitation(key_i,key_j,exc,phase,Nint) + ! Mono alpha, mono beta + if (exc(0,1,1) == 1) then + if ( (exc(1,1,1) == exc(1,2,2)).and.(exc(1,1,2) == exc(1,2,1)) ) then + s2 = -phase + endif + if(exc(1,1,1) == exc(1,2,2) )then + hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) ==exc(1,1,2))then + hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif + ! Double alpha + else if (exc(0,1,1) == 2) then + hij = phase*(get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + ! Double beta + else if (exc(0,1,2) == 2) then + hij = phase*(get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) endif - if(exc(1,1,1) == exc(1,2,2) )then - hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) - else if (exc(1,2,1) ==exc(1,1,2))then - hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) - else - hij = phase*get_mo_bielec_integral( & - exc(1,1,1), & - exc(1,1,2), & - exc(1,2,1), & - exc(1,2,2) ,mo_integrals_map) - endif - ! Double alpha - else if (exc(0,1,1) == 2) then - hij = phase*(get_mo_bielec_integral( & - exc(1,1,1), & - exc(2,1,1), & - exc(1,2,1), & - exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral( & - exc(1,1,1), & - exc(2,1,1), & - exc(2,2,1), & - exc(1,2,1) ,mo_integrals_map) ) - ! Double beta - else if (exc(0,1,2) == 2) then - hij = phase*(get_mo_bielec_integral( & - exc(1,1,2), & - exc(2,1,2), & - exc(1,2,2), & - exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral( & - exc(1,1,2), & - exc(2,1,2), & - exc(2,2,2), & - exc(1,2,2) ,mo_integrals_map) ) endif case (1) call get_mono_excitation(key_i,key_j,exc,phase,Nint) @@ -634,44 +638,48 @@ subroutine i_H_j(key_i,key_j,Nint,hij) integer :: spin select case (degree) case (2) - call get_double_excitation(key_i,key_j,exc,phase,Nint) - if (exc(0,1,1) == 1) then - ! Mono alpha, mono beta - if(exc(1,1,1) == exc(1,2,2) )then - hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) - else if (exc(1,2,1) ==exc(1,1,2))then - hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) - else - hij = phase*get_mo_bielec_integral( & - exc(1,1,1), & - exc(1,1,2), & - exc(1,2,1), & - exc(1,2,2) ,mo_integrals_map) + if (zero_doubles) then + hij = 0.d0 + else + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + if(exc(1,1,1) == exc(1,2,2) )then + hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) ==exc(1,1,2))then + hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) endif - else if (exc(0,1,1) == 2) then - ! Double alpha - hij = phase*(get_mo_bielec_integral( & - exc(1,1,1), & - exc(2,1,1), & - exc(1,2,1), & - exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral( & - exc(1,1,1), & - exc(2,1,1), & - exc(2,2,1), & - exc(1,2,1) ,mo_integrals_map) ) - else if (exc(0,1,2) == 2) then - ! Double beta - hij = phase*(get_mo_bielec_integral( & - exc(1,1,2), & - exc(2,1,2), & - exc(1,2,2), & - exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral( & - exc(1,1,2), & - exc(2,1,2), & - exc(2,2,2), & - exc(1,2,2) ,mo_integrals_map) ) endif case (1) call get_mono_excitation(key_i,key_j,exc,phase,Nint) @@ -731,38 +739,42 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree) call get_excitation_degree(key_i,key_j,degree,Nint) select case (degree) case (2) - call get_double_excitation(key_i,key_j,exc,phase,Nint) - if (exc(0,1,1) == 1) then - ! Mono alpha, mono beta - hij = phase*get_mo_bielec_integral( & - exc(1,1,1), & - exc(1,1,2), & - exc(1,2,1), & - exc(1,2,2) ,mo_integrals_map) - else if (exc(0,1,1) == 2) then - ! Double alpha - hij = phase*(get_mo_bielec_integral( & - exc(1,1,1), & - exc(2,1,1), & - exc(1,2,1), & - exc(2,2,1) ,mo_integrals_map) - & - get_mo_bielec_integral( & - exc(1,1,1), & - exc(2,1,1), & - exc(2,2,1), & - exc(1,2,1) ,mo_integrals_map) ) - else if (exc(0,1,2) == 2) then - ! Double beta - hij = phase*(get_mo_bielec_integral( & - exc(1,1,2), & - exc(2,1,2), & - exc(1,2,2), & - exc(2,2,2) ,mo_integrals_map) - & - get_mo_bielec_integral( & - exc(1,1,2), & - exc(2,1,2), & - exc(2,2,2), & - exc(1,2,2) ,mo_integrals_map) ) + if (zero_doubles) then + hij = 0.d0 + else + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) + endif endif case (1) call get_mono_excitation(key_i,key_j,exc,phase,Nint) @@ -3146,18 +3158,21 @@ subroutine i_H_j_double_spin(key_i,key_j,Nint,hij) double precision, external :: get_mo_bielec_integral PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map - - call get_double_excitation_spin(key_i,key_j,exc,phase,Nint) - hij = phase*(get_mo_bielec_integral( & - exc(1,1), & - exc(2,1), & - exc(1,2), & - exc(2,2), mo_integrals_map) - & - get_mo_bielec_integral( & - exc(1,1), & - exc(2,1), & - exc(2,2), & - exc(1,2), mo_integrals_map) ) + if (zero_doubles) then + hij = 0.d0 + else + call get_double_excitation_spin(key_i,key_j,exc,phase,Nint) + hij = phase*(get_mo_bielec_integral( & + exc(1,1), & + exc(2,1), & + exc(1,2), & + exc(2,2), mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1), & + exc(2,1), & + exc(2,2), & + exc(1,2), mo_integrals_map) ) + endif end subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij) @@ -3176,19 +3191,23 @@ subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij) PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map - call get_mono_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint) - call get_mono_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase2,Nint) - phase = phase*phase2 - if (exc(1,1,1) == exc(1,2,2)) then - hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) - else if (exc(1,2,1) == exc(1,1,2)) then - hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + if (zero_doubles) then + hij = 0.d0 else - hij = phase*get_mo_bielec_integral( & - exc(1,1,1), & - exc(1,1,2), & - exc(1,2,1), & - exc(1,2,2) ,mo_integrals_map) + call get_mono_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint) + call get_mono_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase2,Nint) + phase = phase*phase2 + if (exc(1,1,1) == exc(1,2,2)) then + hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) == exc(1,1,2)) then + hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif endif end diff --git a/src/Iterations/io.irp.f b/src/Iterations/io.irp.f index 3a145996..821f5e84 100644 --- a/src/Iterations/io.irp.f +++ b/src/Iterations/io.irp.f @@ -8,20 +8,15 @@ BEGIN_PROVIDER [ integer, n_iter ] PROVIDE ezfio_filename if (mpi_master) then - call ezfio_has_iterations_n_iter(has) - if (has) then - write(6, *) '.. >>>>> [ IO Read ] n_iter <<<<< ..' - call ezfio_get_iterations_n_iter(n_iter) - else - n_iter = 1 double precision :: zeros(N_states,100) integer :: izeros(100) zeros = 0.d0 izeros = 0 + call ezfio_set_iterations_n_iter(0) call ezfio_set_iterations_energy_iterations(zeros) call ezfio_set_iterations_pt2_iterations(zeros) - call ezfio_set_iterations_n_det_iterations(zeros) - endif + call ezfio_set_iterations_n_det_iterations(izeros) + n_iter = 1 endif IRP_IF MPI_DEBUG print *, irp_here, mpi_rank diff --git a/src/Iterations/print_summary.irp.f b/src/Iterations/print_summary.irp.f index 95f01ce9..923a5ce2 100644 --- a/src/Iterations/print_summary.irp.f +++ b/src/Iterations/print_summary.irp.f @@ -9,7 +9,7 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_) integer :: N_states_p character*(9) :: pt2_string character*(512) :: fmt - + double precision :: f(N_states) if (do_pt2) then pt2_string = ' ' @@ -19,6 +19,10 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_) N_states_p = min(N_det,N_states) + do i=1,N_states_p + f(i) = 1.d0/(1+norm_(i)) + enddo + print *, '' print '(A,I12)', 'Summary at N_det = ', N_det print '(A)', '-----------------------------------' @@ -40,6 +44,7 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_) write(*,fmt) '# PT2'//pt2_string, (pt2_(k), error_(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) if (N_states_p > 1) then write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_(k)-e_(1)-pt2_(1)), & dsqrt(error_(k)*error_(k)+error_(1)*error_(1)), k=1,N_states_p) @@ -59,7 +64,8 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_) print *, 'PT norm = ', dsqrt(norm_(k)) print *, 'PT2 = ', pt2_(k) print *, 'E = ', e_(k) - print *, 'E+PT2'//pt2_string//' = ', e_(k)+pt2_(k), ' +/- ', error_(k) + print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_(k), ' +/- ', error_(k) + print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_(k)*f(k), ' +/- ', error_(k)*f(k) enddo print *, '-----' @@ -75,6 +81,12 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_) print*,'Delta E = ', (e_(i)+ pt2_(i) - (e_(1) + pt2_(1))), & (e_(i)+ pt2_(i) - (e_(1) + pt2_(1))) * 27.211396641308d0 enddo + print *, '-----' + print*, 'Variational + perturbative Energy difference (au | eV)' + do i=2, N_states_p + print*,'Delta E = ', (e_(i)+ pt2_(i)*f(i) - (e_(1) + pt2_(1)*f(1))), & + (e_(i)+ pt2_(i)*f(i) - (e_(1) + pt2_(1)*f(1))) * 27.211396641308d0 + enddo endif end subroutine