Corrected MRCC bug

Conflicts:
	plugins/MRCC_Utils/mrcc_utils.irp.f

plugins/MRCC_Utils/mrcc_dress.irp.f

@@ -24,7 +24,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
   double precision, intent(inout) :: delta_ii_(Nstates,Ndet_ref)
 
   integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
-  integer                        :: i,j,k,l
+  integer                        :: i,j,k,l,m
   integer                        :: degree_alpha(psi_det_size)
   integer                        :: idx_alpha(0:psi_det_size)
   logical                        :: good, fullMatch
@@ -48,9 +48,14 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
   integer                        :: N_miniList, ni, leng
   double precision, allocatable  :: hij_cache(:)
   
+  integer(bit_kind), allocatable :: microlist(:,:,:), microlist_zero(:,:,:)
+  integer, allocatable           :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:)
+  integer :: mobiles(2), smallerlist
+  
+  
   
   leng = max(N_det_generators, N_det_non_ref)
-  allocate(miniList(Nint, 2, leng), idx_miniList(leng), hij_cache(N_det_non_ref))
+  allocate(miniList(Nint, 2, leng), idx_minilist(leng), hij_cache(N_det_non_ref))
   
   !create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint)
   call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
@@ -59,8 +64,21 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
     return
   end if
   
+  allocate(ptr_microlist(0:mo_tot_num*2+1),  &
+      N_microlist(0:mo_tot_num*2) )
+  allocate(   microlist(Nint,2,N_minilist*4),               &
+    idx_microlist(N_minilist*4))
   
+  if(key_mask(1,1) /= 0) then
+    call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint)
+    call find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
+  else
     call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
+  end if
+  
+  
+  
+  deallocate(microlist, idx_microlist)
   
   allocate (dIa_hla(Nstates,Ndet_non_ref))
   
@@ -69,17 +87,101 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
   ! |alpha>
   
   if(N_tq > 0) then
-    call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint)
+    call create_minilist(key_mask, psi_non_ref, miniList, idx_minilist, N_det_non_ref, N_minilist, Nint)
+    if(N_minilist == 0) return
+    
+    
+    if(key_mask(1,1) /= 0) then !!!!!!!!!!! PAS GENERAL !!!!!!!!!
+      allocate(microlist_zero(Nint,2,N_minilist), idx_microlist_zero(N_minilist))
+      
+      allocate(   microlist(Nint,2,N_minilist*4),               &
+        idx_microlist(N_minilist*4))
+      call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint)
+      
+      
+      do i=0,mo_tot_num*2
+        do k=ptr_microlist(i),ptr_microlist(i+1)-1
+          idx_microlist(k) = idx_minilist(idx_microlist(k))
+        end do
+      end do
+      
+      do l=1,N_microlist(0)
+        do k=1,Nint
+          microlist_zero(k,1,l) = microlist(k,1,l)
+          microlist_zero(k,2,l) = microlist(k,2,l)
+        enddo
+        idx_microlist_zero(l) = idx_microlist(l)
+      enddo
+    end if
   end if
   
   
-  do i_alpha=1,N_tq
-    !    call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
-    call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
       
+ 
+  do i_alpha=1,N_tq
+    if(key_mask(1,1) /= 0) then
+      call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint) 
+      
+      if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then
+        smallerlist = mobiles(1)
+      else
+        smallerlist = mobiles(2)
+      end if
+      
+    
+      do l=0,N_microlist(smallerlist)-1
+        microlist_zero(:,:,ptr_microlist(1) + l) = microlist(:,:,ptr_microlist(smallerlist) + l)
+        idx_microlist_zero(ptr_microlist(1) + l) = idx_microlist(ptr_microlist(smallerlist) + l)
+      end do
+      
+      call get_excitation_degree_vector(microlist_zero,tq(1,1,i_alpha),degree_alpha,Nint,N_microlist(smallerlist)+N_microlist(0),idx_alpha)
+      do j=1,idx_alpha(0)
+        idx_alpha(j) = idx_microlist_zero(idx_alpha(j))
+      end do
+      
+      
+!       i = 1
+!       j = 2
+!       do j = 2, idx_alpha_tmp(0)
+!         if(idx_alpha_tmp(j) < idx_alpha_tmp(j-1)) exit
+!       end do
+!       
+!       m = j
+!       
+!       idx_alpha(0) = idx_alpha_tmp(0)
+!           
+!       do l = 1, idx_alpha(0)      
+!         if(j > idx_alpha_tmp(0)) then
+!           k = i
+!           i += 1
+!         else if(i >= m) then
+!           k = j
+!           j += 1
+!         else if(idx_alpha_tmp(i) < idx_alpha_tmp(j)) then
+!           k = i
+!           i += 1
+!         else
+!           k = j
+!           j += 1
+!         end if
+! !  k=l       
+!         idx_alpha(l) = idx_alpha_tmp(k)
+!         degree_alpha(l) = degree_alpha_tmp(k)
+!       end do
+!       
+    else
+      call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
       do j=1,idx_alpha(0)
         idx_alpha(j) = idx_miniList(idx_alpha(j))
       end do
+    end if
+    
+    
+!     call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
+!     do j=1,idx_alpha(0)
+!       idx_alpha(j) = idx_miniList(idx_alpha(j))
+!     end do
+    !print *, idx_alpha(:idx_alpha(0))
     
     do l_sd=1,idx_alpha(0)
       k_sd = idx_alpha(l_sd)
@@ -133,32 +235,17 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
           tmp_det(k,1) = psi_ref(k,1,i_I)
           tmp_det(k,2) = psi_ref(k,2,i_I)
         enddo
-        ! Hole (see list_to_bitstring)
-        iint = ishft(h1-1,-bit_kind_shift) + 1
-        ipos = h1-ishft((iint-1),bit_kind_shift)-1
-        tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos)
         
-        ! Particle
+        logical :: ok
-        iint = ishft(p1-1,-bit_kind_shift) + 1
+        call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint)
-        ipos = p1-ishft((iint-1),bit_kind_shift)-1
+        if(.not. ok) cycle
-        tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
-        if (degree_alpha(k_sd) == 2) then
-          ! Hole (see list_to_bitstring)
-          iint = ishft(h2-1,-bit_kind_shift) + 1
-          ipos = h2-ishft((iint-1),bit_kind_shift)-1
-          tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos)
-          
-          ! Particle
-          iint = ishft(p2-1,-bit_kind_shift) + 1
-          ipos = p2-ishft((iint-1),bit_kind_shift)-1
-          tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
-        endif
         
         ! <I| \l/ |alpha>
         do i_state=1,Nstates
           dka(i_state) = 0.d0
         enddo
         do l_sd=k_sd+1,idx_alpha(0)
+
           call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
           if (degree == 0) then
             
@@ -172,11 +259,12 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
             if (.not.loop) then
               call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
               hIl = hij_mrcc(idx_alpha(l_sd),i_I)
-              !               call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
+!                             call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
               do i_state=1,Nstates
                 dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
               enddo
             endif
+
             exit
           endif
         enddo
@@ -211,21 +299,13 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
       call omp_unset_lock( psi_ref_lock(i_I) )
     enddo
   enddo
-  deallocate (dIa_hla,hij_cache)
+  !deallocate (dIa_hla,hij_cache)
-  deallocate(miniList, idx_miniList)
+  !deallocate(miniList, idx_miniList)
 end
 
 
 
- BEGIN_PROVIDER [ integer(bit_kind), gen_det_sorted,  (N_int,2,N_det_generators,2) ]
+
-&BEGIN_PROVIDER [ integer, gen_det_shortcut, (0:N_det_generators,2) ]
-&BEGIN_PROVIDER [ integer, gen_det_version, (N_int, N_det_generators,2) ]
-&BEGIN_PROVIDER [ integer, gen_det_idx, (N_det_generators,2) ]
-  gen_det_sorted(:,:,:,1) = psi_det_generators(:,:,:N_det_generators)
-  gen_det_sorted(:,:,:,2) = psi_det_generators(:,:,:N_det_generators)
-  call sort_dets_ab_v(gen_det_sorted(:,:,:,1), gen_det_idx(:,1), gen_det_shortcut(0:,1), gen_det_version(:,:,1), N_det_generators, N_int)
-  call sort_dets_ba_v(gen_det_sorted(:,:,:,2), gen_det_idx(:,2), gen_det_shortcut(0:,2), gen_det_version(:,:,2), N_det_generators, N_int)
-END_PROVIDER
 
 
 subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList)
@@ -258,6 +338,7 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
   N_tq = 0
   
   
+  
   i_loop : do i=1,N_selected
     if(is_connected_to(det_buffer(1,1,i), miniList, Nint, N_miniList)) then
       cycle
@@ -287,8 +368,84 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
 end
 
 
+subroutine find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
+
+ use bitmasks
+ implicit none
+
+  integer, intent(in)            :: i_generator,n_selected, Nint
+
+  integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
+  integer                        :: i,j,k,m
+  logical                        :: is_in_wavefunction
+  integer                        :: degree(psi_det_size)
+  integer                        :: idx(0:psi_det_size)
+  logical                        :: good
+
+  integer(bit_kind), intent(out) :: tq(Nint,2,n_selected)
+  integer, intent(out)           :: N_tq
+  
+  
+  integer                        :: nt,ni
+  logical, external              :: is_connected_to
+  
+  
+  integer(bit_kind),intent(in)  :: microlist(Nint,2,*)
+  integer,intent(in)  :: ptr_microlist(0:*)
+  integer,intent(in)            :: N_microlist(0:*)
+  integer(bit_kind),intent(in)   :: key_mask(Nint, 2)
+  
+  integer :: mobiles(2), smallerlist
+  
+  N_tq = 0
+  
+  
+  
+  i_loop : do i=1,N_selected
+    call getMobiles(det_buffer(1,1,i), key_mask, mobiles, Nint) 
+    if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then
+      smallerlist = mobiles(1)
+    else
+      smallerlist = mobiles(2)
+    end if
+    
+    if(N_microlist(smallerlist) > 0) then
+      if(is_connected_to(det_buffer(1,1,i), microlist(1,1,ptr_microlist(smallerlist)), Nint, N_microlist(smallerlist))) then
+        cycle
+      end if
+    end if
+    
+    if(N_microlist(0) > 0) then
+      if(is_connected_to(det_buffer(1,1,i), microlist, Nint, N_microlist(0))) then
+        cycle
+      end if
+    end if
+
+    ! Select determinants that are triple or quadruple excitations
+    ! from the ref
+    good = .True.
+    call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) 
+    !good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector
+    do k=1,idx(0)
+      if (degree(k) < 3) then
+        good = .False.
+        exit
+      endif
+    enddo
+    if (good) then
+      if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then
+        N_tq += 1
+        do k=1,N_int
+          tq(k,1,N_tq) = det_buffer(k,1,i)
+          tq(k,2,N_tq) = det_buffer(k,2,i)
+        enddo
+      endif
+    endif
+  enddo i_loop
+end
+
+
 
 
 
 
-

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -6,27 +6,21 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ]
   integer                        :: i,k
   double precision               :: ihpsi_current(N_states)
   integer                        :: i_pert_count
-  double precision :: hii, lambda_pert
   
   i_pert_count = 0
   lambda_mrcc = 0.d0
   
   do i=1,N_det_non_ref
-     call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref, &
+    call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,&
         size(psi_ref_coef,1), N_states,ihpsi_current)
-     call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
     do k=1,N_states
       if (ihpsi_current(k) == 0.d0) then
         ihpsi_current(k) = 1.d-32
       endif
       lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k)
-      if ( dabs(psi_non_ref_coef(i,k)*ihpsi_current(k)) < 1.d-8 ) then
+      if ( dabs(psi_non_ref_coef(i,k)*ihpsi_current(k)) < 1.d-6 ) then
         i_pert_count += 1
         lambda_mrcc(k,i) = 0.d0
-!          lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
-!          if((ihpsi_current(k) * lambda_pert) < 0.5d0 * psi_non_ref_coef_restart(i,k) ) then
-!             lambda_mrcc(k,i) = 0.d0
-!          endif
       endif
       double precision, parameter    :: x = 2.d0
       if (lambda_mrcc(k,i) > x) then
@@ -45,8 +39,6 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ]
 END_PROVIDER
 
 
-
-
 !BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ]
 !implicit none
 !BEGIN_DOC

src/Determinants/slater_rules.irp.f

@@ -1728,3 +1728,55 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
   deallocate (shortcut, sort_idx, sorted, version)
 end
 
+
+subroutine apply_excitation(det, exc, res, ok, Nint)
+  use bitmasks
+  implicit none
+  
+  integer, intent(in) :: Nint
+  integer, intent(in) :: exc(0:2,2,2)
+  integer(bit_kind),intent(in) :: det(Nint, 2)
+  integer(bit_kind),intent(out) :: res(Nint, 2)
+  logical, intent(out) :: ok
+  integer :: h1,p1,h2,p2,s1,s2,degree
+  integer :: ii, pos 
+  
+  
+  ok = .false.
+  degree = exc(0,1,1) + exc(0,1,2)
+  
+  if(.not. (degree > 0 .and. degree <= 2)) then
+    print *, degree
+    print *, "apply ex"
+    STOP
+  endif
+  
+  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+  res = det 
+  
+  ii = (h1-1)/bit_kind_size + 1 
+  pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64
+  if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return
+  res(ii, s1) = ibclr(res(ii, s1), pos)
+  
+  ii = (p1-1)/bit_kind_size + 1 
+  pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1)
+  if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return
+  res(ii, s1) = ibset(res(ii, s1), pos)
+  
+  if(degree == 2) then
+    ii = (h2-1)/bit_kind_size + 1 
+    pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1)
+    if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return
+    res(ii, s2) = ibclr(res(ii, s2), pos)
+    
+    ii = (p2-1)/bit_kind_size + 1 
+    pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1)
+    if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return
+    res(ii, s2) = ibset(res(ii, s2), pos)
+  endif
+
+  ok = .true.
+end subroutine
+
+