Iterative sbk, broken

plugins/Psiref_Utils/psi_ref_utils.irp.f

@@ -9,7 +9,7 @@ use bitmasks
  ! function.
  END_DOC
  call sort_dets_by_det_search_key(N_det_ref, psi_ref, psi_ref_coef, &
-     psi_ref_sorted_bit, psi_ref_coef_sorted_bit)
+     psi_ref_sorted_bit, psi_ref_coef_sorted_bit, N_states)
 
 END_PROVIDER
 
@@ -152,7 +152,7 @@ END_PROVIDER
  ! function.
  END_DOC
  call sort_dets_by_det_search_key(N_det_ref, psi_non_ref, psi_non_ref_coef, &
-     psi_non_ref_sorted_bit, psi_non_ref_coef_sorted_bit)
+     psi_non_ref_sorted_bit, psi_non_ref_coef_sorted_bit, N_states)
 
 END_PROVIDER
 

plugins/dress_zmq/dress_general.irp.f

@@ -38,10 +38,10 @@ subroutine run_dressing(N_st,energy)
       E_old = sum(psi_energy(:))
       print *,  'Variational energy <Psi|H|Psi>'
       do i=1,N_st
-        print *,  i, psi_energy(i)
+        print *,  i, psi_energy(i)+nuclear_repulsion
       enddo
       !print *, "DELTA IJ", delta_ij(1,1,1)
-      if(.true.) dummy = delta_ij_tmp(1,1,1)
+      PROVIDE delta_ij_tmp
       if(.true.) call delta_ij_done()
       print *,  'Dressed energy <Psi|H+Delta|Psi>'
       do i=1,N_st

plugins/dress_zmq/dress_stoch_routines.irp.f

@@ -229,10 +229,8 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
 
   pullLoop : do while (loop)
     call pull_dress_results(zmq_socket_pull, ind, cur_cp, delta_loc, int_buf, double_buf, det_buf, N_buf, task_id, dress_mwen)
-    !print *, cur_cp, ind
     if(floop) then
       call wall_time(time)
-      print *, "first_pull", time-time0
       time0 = time
       floop = .false.
     end if
@@ -240,8 +238,6 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
       call wall_time(time)
     end if
     
-    print *,  cur_cp, ind
-    
     if(cur_cp == -1) then
       call dress_pulled(ind, int_buf, double_buf, det_buf, N_buf) 
       if (zmq_delete_tasks(zmq_to_qp_run_socket,zmq_socket_pull,task_id,1,more) == -1) then
@@ -295,7 +291,7 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
 
       !print '(I2X, F16.7, 2X, G16.3, 2X, F16.4, A20)', avg+E(istate)+E0, eqt, time-time0, ''
        print '(G10.3, 2X, F16.10, 2X, G16.3, 2X, F16.4, A20)', cps_N(cur_cp), avg+E0+E(istate), eqt, time-time0, ''
-      if ((dabs(eqt) < relative_error .and. cps_N(cur_cp) >= 30)) then
+      if ((dabs(eqt/(avg+E0+E(istate))) < relative_error .and. cps_N(cur_cp) >= 10)) then
         ! Termination
         print *, "TERMINATE"
         if (zmq_abort(zmq_to_qp_run_socket) == -1) then
@@ -511,12 +507,6 @@ END_PROVIDER
   end do
 
 
-print *,  'needed_by_cp'
-do i=1,cur_cp
-  print *,  i, needed_by_cp(i)
-enddo
-  
-
   under_det = 0
   cp_first_tooth = 0
   do i=1,N_dress_jobs

plugins/dress_zmq/dressing.irp.f

@@ -96,12 +96,13 @@ BEGIN_PROVIDER [ double precision, delta_ij_tmp, (N_states,N_det_delta_ij,2) ]
     E_CI_before(:) = psi_energy(:) + nuclear_repulsion
     threshold_selectors = 1.d0
     threshold_generators = 1.d0 
+    SOFT_TOUCH threshold_selectors threshold_generators
   !  if(errr /= 0d0) then
   !    errr = errr / 2d0 
   !  else
   !    errr = 1d-4
   !  end if
-    relative_error = 5.d-5
+    relative_error = 1.d-2
 
     call write_double(6,relative_error,"Convergence of the stochastic algorithm")
     

plugins/shiftedbk/EZFIO.cfg

@@ -3,3 +3,42 @@ type: Perturbation
 doc: Type of zeroth-order Hamiltonian [ EN | Barycentric ]
 interface: ezfio,provider,ocaml
 default: EN
+[energy]
+type: double precision
+doc: Calculated Selected FCI energy
+interface: ezfio
+
+[energy_pt2]
+type: double precision
+doc: Calculated FCI energy + PT2
+interface: ezfio
+
+[iterative_save]
+type: integer
+doc: Save data at each iteration :  1(Append) | 2(Overwrite) | 3(NoSave)
+interface: ezfio,ocaml
+default: 2
+
+[n_iter]
+interface: ezfio
+doc: number of iterations
+type:integer
+
+[n_det_iter]
+interface: ezfio
+doc: number of determinants at iteration
+type: integer
+size: (full_ci_zmq.n_iter)
+
+[energy_iter]
+interface: ezfio
+doc: The energy without a pt2 correction for n_det
+type: double precision
+size: (determinants.n_states,full_ci_zmq.n_iter)
+
+[pt2_iter]
+interface: ezfio
+doc: The pt2 correction for n_det
+type: double precision
+size: (determinants.n_states,full_ci_zmq.n_iter)
+

plugins/shiftedbk/selection_buffer.irp.f

@@ -156,3 +156,89 @@ end subroutine
 
 
 
+subroutine unique_selection_buffer(b)
+  use selection_types
+  implicit none
+  BEGIN_DOC
+! Removes duplicate determinants in the selection buffer
+  END_DOC
+  type(selection_buffer), intent(inout) :: b
+  integer, allocatable :: iorder(:)
+  integer(bit_kind), pointer :: detmp(:,:,:)
+  double precision, pointer      :: val(:)
+  integer                        :: i,j,k
+  integer(bit_kind), allocatable :: bit_tmp(:)
+  logical,allocatable            :: duplicate(:)
+
+  logical :: found_duplicates
+  integer*8, external :: det_search_key
+
+  if (b%N == 0 .or. b%cur == 0) return
+  allocate (duplicate(b%cur), val(size(b%val)), detmp(N_int, 2, size(b%val)), bit_tmp(b%cur))
+  call sort_dets_by_det_search_key(b%cur, b%det, b%val, detmp, val, 1) 
+
+  deallocate(b%det, b%val)
+  do i=b%cur+1,b%N
+    val(i) = 0.d0
+    detmp(1:N_int,1:2,i) = 0_bit_kind
+  enddo
+  b%det => detmp
+  b%val => val
+
+  do i=1,b%cur
+    bit_tmp(i) = det_search_key(b%det(1,1,i),N_int)
+    duplicate(i) = .False.
+  enddo
+
+  do i=1,b%cur-1
+    if (duplicate(i)) then
+      cycle
+    endif
+    j = i+1
+    do while (bit_tmp(j)==bit_tmp(i))
+      if (duplicate(j)) then
+        j += 1
+        if (j > b%cur) then
+          exit
+        else
+          cycle
+        endif
+      endif
+      duplicate(j) = .True.
+      do k=1,N_int
+        if ( (b%det(k,1,i) /= b%det(k,1,j) ) &
+        .or. (b%det(k,2,i) /= b%det(k,2,j) ) ) then
+          duplicate(j) = .False.
+          exit
+        endif
+      enddo
+      j += 1
+      if (j > b%cur) then
+        exit
+      endif
+    enddo
+  enddo
+
+  found_duplicates = .False.
+  do i=1,b%cur
+    if (duplicate(i)) then
+      found_duplicates = .True.
+      exit
+    endif
+  enddo
+
+  if (found_duplicates) then
+    k=0
+    do i=1,N_det
+      if (.not.duplicate(i)) then
+        k += 1
+        b%det(:,:,k) = b%det(:,:,i)
+        b%val(k)  = b%val(i)
+      endif
+    enddo
+    b%cur = k
+  endif
+  deallocate (duplicate,bit_tmp)
+end
+
+

plugins/shiftedbk/shifted_bk_iter.irp.f

@@ -0,0 +1,159 @@
+program shifted_bk
+  implicit none
+  integer                        :: i,j,k
+  double precision, allocatable  :: pt2(:)
+  integer                        :: degree
+  integer                        :: n_det_before
+  double precision               :: threshold_davidson_in
+  
+  allocate (pt2(N_states))
+
+  double precision               :: hf_energy_ref
+  logical                        :: has
+  double precision               :: relative_error, absolute_error
+  integer                        :: N_states_p
+  character*(512)                :: fmt
+
+  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
+
+
+  pt2 = -huge(1.e0)
+  threshold_davidson_in = threshold_davidson
+  threshold_davidson = threshold_davidson_in * 100.d0
+  SOFT_TOUCH threshold_davidson
+
+  call diagonalize_CI_dressed
+  call save_wavefunction
+  
+  call ezfio_has_hartree_fock_energy(has)
+  if (has) then
+    call ezfio_get_hartree_fock_energy(hf_energy_ref)
+  else
+    hf_energy_ref = ref_bitmask_energy
+  endif
+
+  if (N_det > N_det_max) then
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    N_det = N_det_max
+    soft_touch N_det psi_det psi_coef
+    call diagonalize_CI_dressed
+    call save_wavefunction
+    N_states_p = min(N_det,N_states)
+  endif
+  
+  n_det_before = 0
+
+  character*(8) :: pt2_string
+  double precision :: threshold_selectors_save, threshold_generators_save
+  threshold_selectors_save  = threshold_selectors
+  threshold_generators_save = threshold_generators
+  double precision :: error(N_states), energy(N_states)
+  error = 0.d0
+
+  threshold_selectors = 1.d0
+  threshold_generators = 1d0 
+
+  if (.True.) then 
+    pt2_string = '(sh-Bk) '
+    do while ( (N_det < N_det_max) )
+      write(*,'(A)')  '--------------------------------------------------------------------------------'
+
+      N_det_delta_ij = N_det
+
+      do i=1,N_states
+        energy(i) = psi_energy(i)+nuclear_repulsion
+      enddo
+
+      PROVIDE delta_ij_tmp
+      call delta_ij_done()
+
+      call diagonalize_ci_dressed
+      do i=1,N_states
+        pt2(i) = ci_energy_dressed(i) - energy(i)
+      enddo
+
+      N_states_p = min(N_det,N_states)
+
+      print *, ''
+      print '(A,I12)',  'Summary at N_det = ', N_det
+      print '(A)',      '-----------------------------------'
+      print *, ''
+      print *, ''
+
+      write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
+      write(*,fmt)
+      write(fmt,*) '(12X,', 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) '# E          ', energy(1:N_states_p)
+      if (N_states_p > 1) then
+        write(*,fmt) '# Excit. (au)', energy(1:N_states_p)-energy(1)
+        write(*,fmt) '# Excit. (eV)', (energy(1:N_states_p)-energy(1))*27.211396641308d0
+      endif
+      write(fmt,*) '(A12,', 2*N_states_p, '(1X,F14.8))'
+      write(*,fmt) '# PT2'//pt2_string, (pt2(k), error(k), k=1,N_states_p)
+      write(*,'(A)') '#'
+      write(*,fmt) '# E+PT2      ', (energy(k)+pt2(k),error(k), k=1,N_states_p)
+      if (N_states_p > 1) then
+        write(*,fmt) '# Excit. (au)', ( (energy(k)+pt2(k)-energy(1)-pt2(1)), &
+          dsqrt(error(k)*error(k)+error(1)*error(1)), k=1,N_states_p)
+        write(*,fmt) '# Excit. (eV)', ( (energy(k)+pt2(k)-energy(1)-pt2(1))*27.211396641308d0, &
+          dsqrt(error(k)*error(k)+error(1)*error(1))*27.211396641308d0, k=1,N_states_p)
+      endif
+      write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
+      write(*,fmt)
+      print *,  ''
+
+      print *,  'N_det             = ', N_det
+      print *,  'N_states          = ', N_states
+
+      do k=1, N_states_p
+        print*,'State ',k
+        print *,  'PT2             = ', pt2(k)
+        print *,  'E               = ', energy(k)
+        print *,  'E+PT2'//pt2_string//'   = ', energy(k)+pt2(k)
+      enddo
+
+      print *,  '-----'
+      if(N_states.gt.1)then
+        print *, 'Variational Energy difference (au | eV)'
+        do i=2, N_states_p
+          print*,'Delta E = ', (energy(i) - energy(1)), &
+            (energy(i) - energy(1)) * 27.211396641308d0
+        enddo
+        print *,  '-----'
+        print*, 'Variational + perturbative Energy difference (au | eV)'
+        do i=2, N_states_p
+          print*,'Delta E = ', (energy(i)+ pt2(i) - (energy(1) + pt2(1))), &
+            (energy(i)+ pt2(i) - (energy(1) + pt2(1))) * 27.211396641308d0
+        enddo
+      endif
+      call ezfio_set_shiftedbk_energy_pt2(energy(1)+pt2(1))
+!      call dump_fci_iterations_value(N_det,energy,pt2) 
+
+      n_det_before = N_det
+      
+      PROVIDE  psi_coef
+      PROVIDE  psi_det
+      PROVIDE  psi_det_sorted
+
+      if (N_det >= N_det_max) then
+        threshold_davidson = threshold_davidson_in
+      end if
+      call save_wavefunction
+      call ezfio_set_shiftedbk_energy(energy(1))
+      call ezfio_set_shiftedbk_energy_pt2(ci_energy_dressed(1))
+    enddo
+  endif
+
+
+
+
+end
+

plugins/shiftedbk/shifted_bk_routines.irp.f

@@ -20,7 +20,7 @@ END_PROVIDER
    fock_diag_tmp_(:,:,:) = 0.d0
    integer                        :: i
    
-   N_det_increase_factor = 1d0
+   N_det_increase_factor = dble(N_states)
    
    
    n_det_add = max(1, int(float(N_det) * N_det_increase_factor))
@@ -120,17 +120,16 @@ subroutine delta_ij_done()
   old_det_gen = N_det_generators
  
 
-  if (dress_stoch_istate == N_states) then
+  ! Add buffer only when the last state is computed
-    ! Add buffer only when the last state is computed
+  call unique_selection_buffer(global_sb)
-    call sort_selection_buffer(global_sb)
+  call sort_selection_buffer(global_sb)
-    call fill_H_apply_buffer_no_selection(global_sb%cur,global_sb%det,N_int,0) 
+  call fill_H_apply_buffer_no_selection(global_sb%cur,global_sb%det,N_int,0) 
-    call copy_H_apply_buffer_to_wf()
+  call copy_H_apply_buffer_to_wf()
-    if (s2_eig.or.(N_states > 1) ) then
+  if (s2_eig.or.(N_states > 1) ) then
-      call make_s2_eigenfunction
+    call make_s2_eigenfunction
-    endif
-    call undress_with_alpha(old_generators, old_det_gen, psi_det(1,1,N_det_delta_ij+1), N_det-N_det_delta_ij)
-    call save_wavefunction
   endif
+  call undress_with_alpha(old_generators, old_det_gen, psi_det(1,1,N_det_delta_ij+1), N_det-N_det_delta_ij)
+  call save_wavefunction
 
 end subroutine
 

src/Determinants/H_apply.irp.f

@@ -192,8 +192,8 @@ subroutine copy_H_apply_buffer_to_wf
   call normalize(psi_coef,N_det)
   SOFT_TOUCH N_det psi_det psi_coef
   
-  logical :: found_duplicates
+!  logical :: found_duplicates
-  call remove_duplicates_in_psi_det(found_duplicates)
+!  call remove_duplicates_in_psi_det(found_duplicates)
 
 end
 

src/Determinants/determinants.irp.f

@@ -321,21 +321,24 @@ end subroutine
    END_DOC
    
    call sort_dets_by_det_search_key(N_det, psi_det, psi_coef,        &
-       psi_det_sorted_bit, psi_coef_sorted_bit)
+       psi_det_sorted_bit, psi_coef_sorted_bit, N_states)
 END_PROVIDER
  
-subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out)
+subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out, N_st)
    use bitmasks
    implicit none
-   integer, intent(in)            :: Ndet
+   integer, intent(in)            :: Ndet, N_st
    integer(bit_kind), intent(in)  :: det_in  (N_int,2,psi_det_size)
-   double precision , intent(in)  :: coef_in(psi_det_size,N_states)
+   double precision , intent(in)  :: coef_in(psi_det_size,N_st)
    integer(bit_kind), intent(out) :: det_out (N_int,2,psi_det_size)
-   double precision , intent(out) :: coef_out(psi_det_size,N_states)
+   double precision , intent(out) :: coef_out(psi_det_size,N_st)
    BEGIN_DOC
    ! Determinants are sorted are sorted according to their det_search_key.
    ! Useful to accelerate the search of a random determinant in the wave
    ! function.
+   !
+   ! /!\ The first dimension of coef_out and coef_in need to be psi_det_size
+   !
    END_DOC
    integer                        :: i,j,k
    integer, allocatable           :: iorder(:)
@@ -356,7 +359,7 @@ subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out)
        det_out(j,1,i) = det_in(j,1,iorder(i))
        det_out(j,2,i) = det_in(j,2,iorder(i))
      enddo
-     do k=1,N_states
+     do k=1,N_st
        coef_out(i,k) = coef_in(iorder(i),k)
      enddo
    enddo

src/Determinants/psi_cas.irp.f

@@ -54,7 +54,7 @@ END_PROVIDER
  ! function.
  END_DOC
  call sort_dets_by_det_search_key(N_det_cas, psi_cas, psi_cas_coef, &
-     psi_cas_sorted_bit, psi_cas_coef_sorted_bit)
+     psi_cas_sorted_bit, psi_cas_coef_sorted_bit, N_states)
 
 END_PROVIDER
 
@@ -107,7 +107,7 @@ END_PROVIDER
  ! function.
  END_DOC
  call sort_dets_by_det_search_key(N_det_cas, psi_non_cas, psi_non_cas_coef, &
-     psi_non_cas_sorted_bit, psi_non_cas_coef_sorted_bit)
+     psi_non_cas_sorted_bit, psi_non_cas_coef_sorted_bit, N_states)
 
 END_PROVIDER