PT2 stoch fixed, new algo teeth

plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f

@@ -27,7 +27,7 @@ subroutine ZMQ_pt2(E, pt2,relative_error)
   double precision, external :: omp_get_wtime 
   double precision :: time
 
-  allocate(pt2_detail(N_states, N_det_generators), comb(N_det_generators/2), computed(N_det_generators), tbc(0:size_tbc))
+  allocate(pt2_detail(N_states, N_det_generators), comb(N_det_generators), computed(N_det_generators), tbc(0:size_tbc))
   sumabove = 0d0
   sum2above = 0d0
   Nabove = 0d0
@@ -54,7 +54,17 @@ subroutine ZMQ_pt2(E, pt2,relative_error)
     call create_selection_buffer(1, 1*2, b)
     
     Ncomb=size(comb)
-    call get_carlo_workbatch(computed, comb, Ncomb, tbc)
+!    i=N_det_generators
+!    do while (tbc(0) < i)
+      call get_carlo_workbatch(computed, comb, Ncomb, tbc)
+!      i=0
+!      do j=1,N_det_generators
+!        if (.not.computed(j)) then 
+!          i = i+1
+!        endif
+!      enddo
+!      i = i/2
+!    enddo
 
 
 
@@ -370,7 +380,6 @@ subroutine get_last_full_tooth(computed, last_tooth)
   last_tooth = 0
    combLoop : do i=comb_teeth, 1, -1
      missing = 1+ ishft(first_det_of_teeth(i+1)-first_det_of_teeth(i),-4) ! /16
-!     missing = 1+ ishft(first_det_of_teeth(i+1)-first_det_of_teeth(i),-14) ! /16384
      do j=first_det_of_teeth(i), first_det_of_teeth(i+1)-1
        if(.not.computed(j)) then
          missing -= 1
@@ -393,52 +402,34 @@ END_PROVIDER
 
 subroutine get_carlo_workbatch(computed, comb, Ncomb, tbc)
   implicit none
-  integer, intent(inout) :: Ncomb
+  integer, intent(inout)         :: Ncomb
-  double precision, intent(out) :: comb(Ncomb)
+  double precision, intent(out)  :: comb(Ncomb)
-  integer, intent(inout) :: tbc(0:size_tbc)
+  integer, intent(inout)         :: tbc(0:size_tbc)
-  logical, intent(inout) :: computed(N_det_generators)
+  logical, intent(inout)         :: computed(N_det_generators)
-  integer :: i, j, last_full, dets(comb_teeth), tbc_save
+  integer                        :: i, j, last_full, dets(comb_teeth), tbc_save
-  integer :: icount, n
+  integer                        :: icount, n
-!  n = tbc(0)
+  integer                        :: k, l
-!  icount = 1 
+  l=1
-!  call RANDOM_NUMBER(comb)
-!  do i=1,size(comb)
-!      comb(i) = comb(i) * comb_step
-!      tbc_save = tbc(0)
-!      !DIR$ FORCEINLINE
-!      call add_comb(comb(i), computed, tbc, size_tbc, comb_teeth)
-!      if (tbc(0) < size(tbc)) then
-!         Ncomb = i
-!      else
-!         tbc(0) = tbc_save
-!         return
-!      endif
-!      icount = icount + tbc(0) - tbc_save
-!      if ((i>1000).and.(icount > n)) then
-!        call get_filling_teeth(computed, tbc)
-!        icount = 0
-!        n = ishft(tbc_save,-4)
-!      endif
-!  enddo
-!  call get_filling_teeth(computed, tbc)
-
-  n = int(sqrt(dble(size(comb))))
- 
   call RANDOM_NUMBER(comb)
-  do j=1,size(comb),n
+  do i=1,size(comb)
-      do i=j,min(size(comb),j+n-1)
+    comb(i) = comb(i) * comb_step
-        comb(i) = comb(i) * comb_step
+    tbc_save = tbc(0)
-        tbc_save = tbc(0)
+    !DIR$ FORCEINLINE
-        !DIR$ FORCEINLINE
+    call add_comb(comb(i), computed, tbc, size_tbc, comb_teeth)
-        call add_comb(comb(i), computed, tbc, size_tbc, comb_teeth)
+    if ( (tbc(0) < size(tbc)-1).and.(l < first_det_of_teeth(comb_teeth)) ) then
-        if (tbc(0) < size(tbc)) then
+      Ncomb = i
-           Ncomb = i
+      do while (computed(l))
-        else
+        l=l+1
-           tbc(0) = tbc_save
+        if (l == size(computed)) exit
-           return
+      enddo
-        endif
+      k=tbc(0)+1
-      end do
+      tbc(k) = l
-      call get_filling_teeth(computed, tbc)
+      computed(l) = .True.
+      tbc(0) = k
+    else
+      tbc(0) = tbc_save
+      return
+    endif
   enddo
   
 end subroutine
@@ -563,7 +554,7 @@ end subroutine
   comb_step = 1d0/dfloat(comb_teeth)
   first_det_of_comb = 1
   do i=1,N_det_generators
-    if(pt2_weight(i)/norm_left < 0.5d0*comb_step) then
+    if(pt2_weight(i)/norm_left < .5d0*comb_step) then
       first_det_of_comb = i
       exit
     end if