corrected bitmask in Full_CI_ZMQ - seemingly working shiftedbk (#68)

* corrected bitmask in Full_CI_ZMQ - seemingly working shiftedbk

* removed set_generators_bitmasks_as_holes_and_particles from dress_zmq

* removed dress_zmq dependency to MRCC_Utils

plugins/Full_CI_ZMQ/pt2_stoch.irp.f

@@ -19,6 +19,7 @@ subroutine run
   double precision               :: E_CI_before, relative_error, absolute_error, eqt
 
   allocate (pt2(N_states))
+  call diagonalize_CI()
   pt2 = 0.d0
   
   E_CI_before = pt2_E0_denominator(1) + nuclear_repulsion

plugins/Full_CI_ZMQ/run_pt2_slave.irp.f

@@ -41,7 +41,7 @@ subroutine run_pt2_slave(thread,iproc,energy)
 
   buf%N = 0
   n_tasks = 0
-  call create_selection_buffer(1, 2, buf)
+  call create_selection_buffer(0, 0, buf)
 
   done = .False.
   do while (.not.done)

plugins/Full_CI_ZMQ/selection.irp.f

@@ -253,6 +253,7 @@ subroutine get_m0(gen, phasemask, bannedOrb, vect, mask, h, p, sp, coefs)
   deallocate(lbanned)
 end 
 
+
 subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,buf,subset)
   use bitmasks
   use selection_types
@@ -598,7 +599,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
   
   logical, external :: detEq
   
-  
   if(sp == 3) then
     s1 = 1
     s2 = 2
@@ -616,15 +616,13 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
     do p2=ib,mo_tot_num
       if(bannedOrb(p2, s2)) cycle
       if(banned(p1,p2)) cycle
+
       if( sum(abs(mat(1:N_states, p1, p2))) == 0d0) cycle
       call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
       
       Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
       min_e_pert = 0d0
       
-!     double precision :: hij       
-!     call i_h_j(psi_det_generators(1,1,i_generator), det, N_int, hij)
-
       do istate=1,N_states
         delta_E = E0(istate) - Hii
         val = mat(istate, p1, p2) + mat(istate, p1, p2) 
@@ -635,7 +633,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
         e_pert = 0.5d0 * (tmp - delta_E)
         pt2(istate) = pt2(istate) + e_pert
         min_e_pert = min(e_pert,min_e_pert)
-!        ci(istate) = e_pert / hij
       end do
       
       if(min_e_pert <= buf%mini) then

plugins/Full_CI_ZMQ/selection_buffer.irp.f

@@ -39,7 +39,7 @@ subroutine add_to_selection_buffer(b, det, val)
   double precision, intent(in) :: val
   integer :: i
 
-  if(val <= b%mini) then
+  if(b%N > 0 .and. val <= b%mini) then
     b%cur += 1
     b%det(1:N_int,1:2,b%cur) = det(1:N_int,1:2)
     b%val(b%cur) = val
@@ -119,7 +119,7 @@ subroutine sort_selection_buffer(b)
   integer(bit_kind), pointer :: detmp(:,:,:)
   integer :: i, nmwen
   logical, external :: detEq
-  if (b%cur == 0) return
+  if (b%N == 0 .or. b%cur == 0) return
   nmwen = min(b%N, b%cur)
 
   allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)))

plugins/dress_zmq/EZFIO.cfg

@@ -0,0 +1,17 @@
+[energy]
+type: double precision
+doc: Calculated energy
+interface: ezfio
+
+[thresh_dressed_ci]
+type: Threshold
+doc: Threshold on the convergence of the dressed CI energy
+interface: ezfio,provider,ocaml
+default: 1.e-5
+
+[n_it_max_dressed_ci]
+type: Strictly_positive_int
+doc: Maximum number of dressed CI iterations
+interface: ezfio,provider,ocaml
+default: 10
+

plugins/dress_zmq/alpha_factory.irp.f

@@ -75,26 +75,19 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
   monoAdo = .true.
   monoBdo = .true.
   
+
   do k=1,N_int
-    !hole    (k,1) = iand(psi_det_generators(k,1,i_generator), generators_bitmask(k,1,s_hole,bitmask_index))
+    hole    (k,1) = iand(psi_det_generators(k,1,i_generator), generators_bitmask(k,1,s_hole,bitmask_index))
-    !hole    (k,2) = iand(psi_det_generators(k,2,i_generator), generators_bitmask(k,2,s_hole,bitmask_index))
+    hole    (k,2) = iand(psi_det_generators(k,2,i_generator), generators_bitmask(k,2,s_hole,bitmask_index))
-    !particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), generators_bitmask(k,1,s_part,bitmask_index))
+    particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), generators_bitmask(k,1,s_part,bitmask_index))
-    !particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), generators_bitmask(k,2,s_part,bitmask_index))
+    particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), generators_bitmask(k,2,s_part,bitmask_index))
-    hole    (k,1) = iand(psi_det_generators(k,1,i_generator), full_ijkl_bitmask(k))
+    !hole    (k,1) = iand(psi_det_generators(k,1,i_generator), full_ijkl_bitmask(k))
-    hole    (k,2) = iand(psi_det_generators(k,2,i_generator), full_ijkl_bitmask(k))
+    !hole    (k,2) = iand(psi_det_generators(k,2,i_generator), full_ijkl_bitmask(k))
-    particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), full_ijkl_bitmask(k))
+    !particle(k,1) = iand(not(psi_det_generators(k,1,i_generator)), full_ijkl_bitmask(k))
-    particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), full_ijkl_bitmask(k))
+    !particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), full_ijkl_bitmask(k))
  
   enddo
 
-  !if(i_generator == 34) then
-  !  call debug_det(psi_det_generators(1,1,34), N_int)
-  !  call debug_det(generators_bitmask(1,1,s_part,bitmask_index), N_int)
-  !  call debug_det(particle, N_int)
-  !  print *, "dddddddddddd"
-  !  stop
-  !end if
-
   integer                        :: N_holes(2), N_particles(2)
   integer                        :: hole_list(N_int*bit_kind_size,2)
   integer                        :: particle_list(N_int*bit_kind_size,2)

plugins/dress_zmq/dress_stoch_routines.irp.f

@@ -286,7 +286,7 @@ subroutine dress_collector(zmq_socket_pull, E, relative_error, delta, delta_s2,
 
 
       call wall_time(time)
-      if ((dabs(eqt) < relative_error .and. cps_N(cur_cp) >= 5)  .or. total_computed == N_det_generators) then
+      if ((dabs(eqt) < relative_error .and. cps_N(cur_cp) >= 30)  .or. total_computed == N_det_generators) then
         ! Termination
         print '(2X, F16.7, 2X, G16.3, 2X, F16.4, A20)', avg+E(istate)+E0, eqt, time-time0, ''
         if (zmq_abort(zmq_to_qp_run_socket) == -1) then

plugins/dress_zmq/dress_zmq_routines.irp.f

@@ -5,7 +5,7 @@ subroutine dress_zmq()
 
   read_wf = .True.
   SOFT_TOUCH read_wf
-  call set_generators_bitmasks_as_holes_and_particles
+  
   if (.True.) then
     integer :: i,j
     do j=1,N_states

plugins/dress_zmq/dressing.irp.f

@@ -86,7 +86,7 @@ BEGIN_PROVIDER [ double precision, delta_ij, (N_states,N_det,2) ]
 !  else
 !    errr = 1d-4
 !  end if
-  relative_error = 1.d-3
+  relative_error = 1.d-4
   call write_double(6,relative_error,"Convergence of the stochastic algorithm")
 
   call ZMQ_dress(E_CI_before, dress, del, del_s2, abs(relative_error))

plugins/shiftedbk/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-dress_zmq DavidsonDressed MRCC_Utils
+dress_zmq DavidsonDressed 

plugins/shiftedbk/shifted_bk.irp.f

@@ -1,17 +1,23 @@
 
-program mrcc_sto
+program shifted_bk
   implicit none
   BEGIN_DOC
 ! TODO
   END_DOC
+
+  call diagonalize_CI()
   call dress_zmq()
 end
 
 
  BEGIN_PROVIDER [ double precision, fock_diag_tmp_, (2,mo_tot_num+1,Nproc) ]
 &BEGIN_PROVIDER [ integer, current_generator_, (Nproc) ]
+&BEGIN_PROVIDER [ double precision, a_h_i, (N_det, Nproc) ]
+&BEGIN_PROVIDER [ double precision, a_s2_i, (N_det, Nproc) ]
   implicit none
   current_generator_(:) = 0
+  a_h_i = 0d0
+  a_s2_i = 0d0
  END_PROVIDER
 
 
@@ -31,40 +37,48 @@ subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minili
   integer(bit_kind), intent(in)   :: alpha(Nint,2), det_minilist(Nint, 2, n_minilist)
   integer,intent(in)              :: minilist(n_minilist)
   double precision, intent(inout) :: delta_ij_loc(Nstates,N_det,2)
-  double precision :: hii, hij, sij, delta_e
+  double precision :: haa, hij, sij
   double precision, external :: diag_H_mat_elem_fock
   integer                        :: i,j,k,l,m, l_sd
-  double precision, save :: tot = 0d0
+  double precision :: hdress, sdress
-  double precision :: de(N_states), val, tmp
+  double precision :: de, a_h_psi(Nstates), c_alpha
   
 
+  a_h_psi = 0d0
+  
   if(current_generator_(iproc) /= i_gen) then
     current_generator_(iproc) = i_gen
     call build_fock_tmp(fock_diag_tmp_(1,1,iproc),psi_det_generators(1,1,i_gen),N_int)
   end if
 
-  hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_gen),alpha,fock_diag_tmp_(1,1,iproc),N_int)
+  haa = diag_H_mat_elem_fock(psi_det_generators(1,1,i_gen),alpha,fock_diag_tmp_(1,1,iproc),N_int)
-  do i=1,N_states
-    de(i) = (E0_denominator(i) - hii)
-  end do
 
-  do i=1,N_states
-    val = 0D0
   do l_sd=1,n_minilist
     call i_h_j_s2(alpha,det_minilist(1,1,l_sd),N_int,hij, sij)
-      val += hij
+    a_h_i(l_sd, iproc) = hij
+    a_s2_i(l_sd, iproc) = sij
+    do i=1,Nstates
+      a_h_psi(i) += hij * psi_coef(minilist(l_sd), i)
+    end do
+  end do
+
+
+  do i=1,Nstates
+    de = E0_denominator(i) - haa
+    if(DABS(de) < 1D-5) cycle
+
+    c_alpha = a_h_psi(i) / de
+
+    do l_sd=1,n_minilist
+      hdress = c_alpha * a_h_i(l_sd, iproc)
+      sdress = c_alpha * a_s2_i(l_sd, iproc)
+      delta_ij_loc(i, minilist(l_sd), 1) += hdress
+      delta_ij_loc(i, minilist(l_sd), 2) += sdress
     end do
-    val = 2d0 * val
-    tmp = dsqrt(de(i)**2 + val**2)
-    if(de(i) < 0d0) tmp = -tmp
-    delta_ij_loc(i, minilist(l_sd), 1) += 0.5d0 * (tmp - de(i)) ! * psi_coef(minilist(l_sd), i)
   end do
 end subroutine
 
 
-
-
-
 BEGIN_PROVIDER [ logical, initialize_E0_denominator ]
     implicit none
     BEGIN_DOC
@@ -73,6 +87,7 @@ BEGIN_PROVIDER [ logical, initialize_E0_denominator ]
     initialize_E0_denominator = .True.
 END_PROVIDER
 
+
 BEGIN_PROVIDER [ double precision, E0_denominator, (N_states) ]
   implicit none
   BEGIN_DOC
@@ -88,3 +103,5 @@ BEGIN_PROVIDER [ double precision, E0_denominator, (N_states) ]
     E0_denominator = -huge(1.d0)
   endif
 END_PROVIDER
+
+