working but super slow mrcc_sto

plugins/Full_CI_ZMQ/selection.irp.f

@@ -71,7 +71,6 @@ subroutine select_connected(i_generator,E0,pt2,b,subset)
       hole_mask(k,2) = iand(generators_bitmask(k,2,s_hole,l), psi_det_generators(k,2,i_generator))
       particle_mask(k,1) = iand(generators_bitmask(k,1,s_part,l), not(psi_det_generators(k,1,i_generator)) )
       particle_mask(k,2) = iand(generators_bitmask(k,2,s_part,l), not(psi_det_generators(k,2,i_generator)) )
-
     enddo
     call select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2,b,subset)
   enddo
@@ -300,6 +299,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
     particle(k,2) = iand(not(psi_det_generators(k,2,i_generator)), particle_mask(k,2))
   enddo
   
+  
   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/alpha_factory.irp.f

@@ -31,6 +31,7 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
   
   integer                         :: h1,h2,s1,s2,s3,i1,i2,ib,sp,k,i,j,nt,ii,n
   integer(bit_kind)               :: hole(N_int,2), particle(N_int,2), mask(N_int, 2), pmask(N_int, 2)
+  integer(bit_kind)  :: mmask(N_int, 2)
   logical                         :: fullMatch, ok
   
   integer(bit_kind) :: mobMask(N_int, 2), negMask(N_int, 2)
@@ -58,12 +59,25 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
   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,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,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)
@@ -335,10 +349,10 @@ subroutine generate_singles_and_doubles(delta_ij_loc, i_generator, bitmask_index
             if(siz > size(abuf)) stop "buffer too small in alpha_factory"
             abuf = 0
             call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, indexes_end, abuf, interesting)
-            indexes_end(:,:) -= 1
+            !indexes_end(:,:) -= 1
-            do i=1,siz
+            !do i=1,siz
-              if(abuf(i) < 1 .or. abuf(i) > N_det) stop "foireous abuf"
+            !  if(abuf(i) < 1 .or. abuf(i) > N_det) stop "foireous abuf"
-            end do
+            !end do
             !print *, "IND1", indexes(1,:)
             !print *, "IND2", indexes_end(1,:)
             !stop
@@ -374,7 +388,7 @@ subroutine alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexe
   allocate(labuf(N_det), putten(N_det))
   putten = .false.
 
-  st1 = indexes_end(0,0)
+  st1 = indexes_end(0,0)-1 !!
   if(st1 > 0) labuf(:st1) = abuf(:st1)
   st1 += 1
 
@@ -382,19 +396,19 @@ subroutine alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexe
     s1 = 1
     s2 = 2
     lindex(:, 1) = indexes(1:,0)
-    lindex_end(:,1) = indexes_end(1:,0)
+    lindex_end(:,1) = indexes_end(1:,0)-1
     lindex(:, 2) = indexes(0, 1:)
-    lindex_end(:, 2) = indexes_end(0, 1:)
+    lindex_end(:, 2) = indexes_end(0, 1:)-1
   else if(sp == 2) then
     s1 = 2
     s2 = 2
     lindex(:, 2) = indexes(0, 1:)
-    lindex_end(:, 2) = indexes_end(0, 1:)
+    lindex_end(:, 2) = indexes_end(0, 1:)-1
   else if(sp == 1) then
     s1 = 1
     s2 = 1
     lindex(:, 1) = indexes(1:, 0)
-    lindex_end(:,1) = indexes_end(1:, 0)
+    lindex_end(:,1) = indexes_end(1:, 0)-1
   end if
 
     do i=1,mo_tot_num
@@ -432,8 +446,8 @@ subroutine alpha_callback_mask(delta_ij_loc, sp, mask, bannedOrb, banned, indexe
         end if
 
         if(indexes(i,j) /= 0) then
-          st4 = st3 + 1 + indexes_end(i,j)-indexes(i,j)
+          st4 = st3 + 1 + indexes_end(i,j)-indexes(i,j) -1!!
-          labuf(st3:st4-1) = abuf(indexes(i,j):indexes_end(i,j))
+          labuf(st3:st4-1) = abuf(indexes(i,j):indexes_end(i,j)-1) !!
         else
           st4 = st3
         end if

plugins/dress_zmq/dressing_vector.irp.f

@@ -1,3 +1,9 @@
+ BEGIN_PROVIDER [ integer, nalp ]
+&BEGIN_PROVIDER [ integer, ninc ]
+   nalp = 0
+   ninc = 0
+ END_PROVIDER
+
  BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ]
 &BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ]
  implicit none
@@ -11,9 +17,6 @@
  double precision :: f, tmp
  double precision, external :: u_dot_v
  
- print *, "DELTA_IJ", delta_ij(1,:10,1)
- print *, "DELTA_IJ div", delta_ij(1,:10,1) / psi_coef(dressed_column_idx(1),1)
- !stop
  do k=1,N_states
    l = dressed_column_idx(k)
    f = 1.d0/psi_coef(l,k)
@@ -27,6 +30,8 @@
    tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det)
    dressing_column_s(l,k) -= tmp * f
  enddo
- !stop
+ print *, "NALP", nalp
+ print *, "NINC", ninc
+ print *, "DELTA_IJ", dressing_column_h(:10,1)
 END_PROVIDER
 

plugins/mrcc_sto/mrcc_sto.irp.f

@@ -4,18 +4,11 @@ program mrcc_sto
   BEGIN_DOC
 ! TODO
   END_DOC
-  print *, "========================"
-  print *, "========================"
-  print *, "========================"
-  print *, "MRCC_STO not implemented - acts as a unittest for dress_zmq"
-  print *, "========================"
-  print *, "========================"
-  print *, "========================"
   call dress_zmq()
 end
 
- BEGIN_PROVIDER [ integer, idx_non_ref_from_sorted, (N_det) ]
+ BEGIN_PROVIDER [ integer, psi_from_sorted, (N_det) ]
-&BEGIN_PROVIDER [ integer, psi_from_sorted, (N_det) ]
+&BEGIN_PROVIDER [ integer, idx_non_ref_from_sorted, (N_det) ]
   implicit none
   integer :: i,inpsisor
   
@@ -30,10 +23,17 @@ end
   end do
 END_PROVIDER
 
+
 subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha)
  use bitmasks
  implicit none
-
+  BEGIN_DOC
+  !delta_ij_loc(:,:,1) : dressing column for H
+  !delta_ij_loc(:,:,2) : dressing column for S2
+  !minilist : indices of determinants connected to alpha ( in psi_det_sorted )
+  !n_minilist : size of minilist
+  !alpha : alpha determinant
+  END_DOC
   integer(bit_kind), intent(in)   :: alpha(N_int,2)
   integer,intent(in)              :: minilist(n_minilist), n_minilist
   double precision, intent(inout) :: delta_ij_loc(N_states,N_det,2)
@@ -49,22 +49,36 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha)
   double precision               :: ci_inv(N_states)
   integer                        :: exc(0:2,2,2)
   integer                        :: h1,h2,p1,p2,s1,s2
-  integer(bit_kind)              :: tmp_det(N_int,2)
+  integer(bit_kind)              :: tmp_det(N_int,2), ctrl
   integer                        :: i_state, k_sd, l_sd, m_sd, ll_sd, i_I
   double precision, allocatable  :: hij_cache(:), sij_cache(:)
   double precision :: Delta_E_inv(N_states)
   double precision :: sdress, hdress
   double precision :: c0(N_states)
-  logical :: ok
+  logical :: ok, ok2
+  integer :: old_ninc
+  double precision :: shdress
   
 
+  if(n_minilist == 1) return
+
+  shdress = 0d0
+  old_ninc = ninc
+
   if (perturbative_triples) then
     PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
   endif
+  
+  do i_I=1,N_det_ref
+    call get_excitation_degree(alpha,psi_ref(1,1,i_I),degree1,N_int)
+    if(degree1 <= 2) return
+  end do
+
   allocate(hij_cache(N_det), sij_cache(N_det)) 
   allocate (dIa_hla(N_states,N_det), dIa_sla(N_states,N_det))
   allocate (idx_alpha(0:n_minilist))
   
+
   do i_state=1,N_states
     c0(i_state) = 1.d0/psi_coef(dressed_column_idx(i_state),i_state)
   enddo
@@ -82,17 +96,30 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha)
       end if
     end do
 
-    if( xor(ok, idx_non_ref_from_sorted(k_sd) > 0)) stop "BUGUE"
+    !if(ok) then
+    !  call get_excitation(psi_det_sorted(1,1,k_sd),alpha,exc,degree1,phase,N_int)
+    !  if(degree1 == 0 .or. degree1 > 2) stop "minilist error"
+    !  call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2)
+    !  
+    !  if(h1 > 10 .or. p1 < 7 .or. p1 == 8 .or. p1 == 9) ok = .false.
+    !  if(ok .and. degree1 == 2) then
+    !    if(h2 > 10 .or. p2 < 7 .or. p2 == 8 .or. p2 == 9) ok = .false.
+    !  end if
+    !  !if(degree1 == 0 .or. degree1 > 2) stop "minilist error"
+    !    !iand(xor(psi_det_sorted(i,2,k_sd), alpha(i,2)), alpha(i,2))
+    !end if
+    
+    !if( xor(ok, idx_non_ref_from_sorted(k_sd) > 0)) stop "BUGUE"
     if(ok) then
       ll_sd += 1
       idx_alpha(ll_sd) = k_sd
 !     call i_h_j(alpha,psi_non_ref(1,1,idx_alpha(l_sd)),N_int,hij_cache(k_sd))
 !     call get_s2(alpha,psi_non_ref(1,1,idx_alpha(l_sd)),N_int,sij_cache(k_sd))
-     call i_h_j(alpha,psi_det_sorted(1,1,k_sd),N_int,hij_cache(k_sd))
+      call i_h_j(alpha,psi_det_sorted(1,1,k_sd),N_int,hij_cache(k_sd))
-     call get_s2(alpha,psi_det_sorted(1,1,k_sd),N_int,sij_cache(k_sd))
+      call get_s2(alpha,psi_det_sorted(1,1,k_sd),N_int,sij_cache(k_sd))
     end if
   enddo
-  
+  if(ll_sd <= 1) return 
   idx_alpha(0) = ll_sd
   
 
@@ -126,9 +153,12 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha)
       enddo
       call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, N_int)
       
+      ok2 = .false.
       do i_state=1,N_states
         dIK(i_state) = dij(i_I, idx_non_ref_from_sorted(idx_alpha(k_sd)), i_state)
+        if(dIK(i_state) /= 0d0) ok2 = .true.
       enddo
+      if(.not. ok2) cycle
       
       ! <I| \l/ |alpha>
       do i_state=1,N_states
@@ -174,6 +204,12 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha)
       enddo
     enddo
     
+    ok2 = .false.
+    do i_state=1,N_states
+      if(dIa(i_state) /= 0d0) ok2 = .true.
+    enddo
+    if(.not. ok2) cycle
+
     do i_state=1,N_states
       ci_inv(i_state) = psi_ref_coef_inv(i_I,i_state)
     enddo
@@ -192,9 +228,15 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha)
         !print *, i_state, idx_alpha(l_sd)
         k_sd = idx_alpha(l_sd)
         m_sd = psi_from_sorted(k_sd)
-        if(psi_det(1,1,m_sd) /= psi_det_sorted(1,1,k_sd)) stop "psi_from_sorted foireous"
+        hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state)! * c0(i_state)
-        hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state)
+        sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state)! * c0(i_state)
-        sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) * c0(i_state)
+        !!$OMP ATOMIC
+        !shdress += 1d0
+        nalp += 1
+        if(hdress /= 0d0) then
+          ninc = ninc + 1
+          !print *, "grepme2", hdress, shdress
+        end if
         !$OMP ATOMIC
         delta_ij_loc(i_state,m_sd,1) += hdress
         !$OMP ATOMIC
@@ -203,6 +245,11 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, n_minilist, alpha)
       enddo
     enddo
   enddo
+
+  !if(ninc /= old_ninc) then
+  !  nalp = nalp + 1
+  !  !print "(A8,I20,I20,E15.5)", "grepme",  alpha(1,1), alpha(1,2), shdress
+  !end if
 end subroutine
 
 

plugins/mrcepa0/dressing.irp.f

@@ -343,7 +343,10 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
   !double precision, external :: get_dij, get_dij_index
   double precision :: Delta_E_inv(N_states)
   double precision, intent(inout) :: contrib(N_states)
-  double precision :: sdress, hdress
+  double precision :: sdress, hdress, shdress
+  integer :: old_ninc
+
+  old_ninc = ninc
 
   if (perturbative_triples) then
     PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
@@ -408,6 +411,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
       
  
   do i_alpha=1,N_tq
+    old_ninc = ninc
+    shdress = 0d0
+    
     if(key_mask(1,1) /= 0) then
       call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint)
       
@@ -545,6 +551,9 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
           k_sd = idx_alpha(l_sd)
           hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) 
           sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) 
+          !!$OMP ATOMIC
+          if(hdress /= 0d0) ninc = ninc + 1
+          shdress += hdress
           !$OMP ATOMIC
           contrib(i_state) += hdress * psi_coef(dressed_column_idx(i_state), i_state) * psi_non_ref_coef(k_sd, i_state)
           !$OMP ATOMIC
@@ -554,7 +563,18 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
         enddo
       enddo
     enddo
+    if(ninc /= old_ninc) then
+      nalp = nalp + 1
+      !print "(A8,I20,I20,E15.5)", "grepme", tq(1,1,i_alpha), tq(1,2,i_alpha), shdress
+      !if(tq(1,1,i_alpha) == 1007 .and. tq(1,2,i_alpha) == 17301943) then
+      ! print *, "foinder", i_generator
+      ! call debug_det(psi_det_generators(1,1, i_generator), N_int)
+      ! call debug_det(tq(1,1,i_alpha), N_int)
+      ! stop
+     ! end if
+    end if
   enddo
+
   deallocate (dIa_hla,dIa_sla,hij_cache,sij_cache)
   deallocate(miniList, idx_miniList)
 end

plugins/mrcepa0/dressing_vector.irp.f

@@ -1,3 +1,10 @@
+
+ BEGIN_PROVIDER [ integer, nalp ]
+&BEGIN_PROVIDER [ integer, ninc ]
+  nalp = 0
+  ninc = 0
+END_PROVIDER
+
  BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ]
 &BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ]
  implicit none
@@ -24,6 +31,8 @@
    tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det)
    dressing_column_s(l,k) -= tmp * f
  enddo
+ print *, "NALP", nalp
+ print *, "NINC", ninc
  print *, "DRESS", dressing_column_h(:10,1)
 ! stop
 END_PROVIDER