CITATION.cff

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

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
+default: 1000
-
-[N_det_max_jacobi]
-type: Det_number_max
-doc: Maximum number of determinants diagonalized by Jacobi
-interface: ezfio,provider,ocaml
-default: 2000
 
 [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
+

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"

src/Determinants/slater_rules.irp.f

@@ -532,6 +532,9 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
   integer :: spin
   select case (degree)
     case (2)
+      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
@@ -574,6 +577,7 @@ subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
               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)
       !DIR$ FORCEINLINE
@@ -634,6 +638,9 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
   integer :: spin
   select case (degree)
     case (2)
+      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
@@ -673,6 +680,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
               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)
       !DIR$ FORCEINLINE
@@ -731,6 +739,9 @@ 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)
+      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
@@ -764,6 +775,7 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
               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)
       !DIR$ FORCEINLINE
@@ -3146,7 +3158,9 @@ 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
-
+  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),                                                    &
@@ -3158,6 +3172,7 @@ subroutine i_H_j_double_spin(key_i,key_j,Nint,hij)
         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,6 +3191,9 @@ subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij)
 
   PROVIDE big_array_exchange_integrals mo_bielec_integrals_in_map
 
+  if (zero_doubles) then
+    hij = 0.d0
+  else
     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
@@ -3190,6 +3208,7 @@ subroutine i_H_j_double_alpha_beta(key_i,key_j,Nint,hij)
           exc(1,2,1),                                                  &
           exc(1,2,2) ,mo_integrals_map)
     endif
+  endif
 end
 
 

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)
+      call ezfio_set_iterations_n_det_iterations(izeros)
-    endif
+      n_iter = 1
   endif
   IRP_IF MPI_DEBUG
     print *,  irp_here, mpi_rank

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