diff --git a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f index 784ca8fa..eb64fc2f 100644 --- a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f +++ b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f @@ -258,7 +258,7 @@ subroutine pt2_collector(E, b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, time = omp_get_wtime() - if(time - timeLast > 1d1 .or. more /= 1) then + if(time - timeLast > 3d0 .or. more /= 1) then timeLast = time do i=1, first_det_of_teeth(1)-1 if(.not.(actually_computed(i))) then @@ -331,7 +331,7 @@ end function BEGIN_PROVIDER [ integer, comb_teeth ] implicit none - comb_teeth = 50 + comb_teeth = 100 END_PROVIDER @@ -369,7 +369,7 @@ 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),-12) ! /4096 + 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 @@ -543,7 +543,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 < comb_step*.5d0) then + if(pt2_weight(i)/norm_left < comb_step*.25d0) then first_det_of_comb = i exit end if diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f index c277cf58..fae1e644 100644 --- a/plugins/Full_CI_ZMQ/selection.irp.f +++ b/plugins/Full_CI_ZMQ/selection.irp.f @@ -509,7 +509,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d logical :: ok integer :: s1, s2, p1, p2, ib, j, istate integer(bit_kind) :: mask(N_int, 2), det(N_int, 2) - double precision :: e_pert, delta_E, val, Hii, max_e_pert,tmp + double precision :: e_pert, delta_E, val, Hii, min_e_pert,tmp double precision, external :: diag_H_mat_elem_fock logical, external :: detEq @@ -536,7 +536,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d 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) - max_e_pert = 0d0 + min_e_pert = 0d0 do istate=1,N_states delta_E = E0(istate) - Hii @@ -545,14 +545,14 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d if (delta_E < 0.d0) then tmp = -tmp endif - e_pert = 0.5d0 * ( tmp - delta_E) + e_pert = 0.5d0 * (tmp - delta_E) pt2(istate) = pt2(istate) + e_pert - max_e_pert = min(e_pert,max_e_pert) + min_e_pert = min(e_pert,min_e_pert) ! ci(istate) = e_pert / mat(istate, p1, p2) end do - if(dabs(max_e_pert) > buf%mini) then - call add_to_selection_buffer(buf, det, max_e_pert) + if(min_e_pert <= buf%mini) then + call add_to_selection_buffer(buf, det, min_e_pert) end if end do end do diff --git a/plugins/Full_CI_ZMQ/selection_buffer.irp.f b/plugins/Full_CI_ZMQ/selection_buffer.irp.f index 902e2af7..d0bc05dd 100644 --- a/plugins/Full_CI_ZMQ/selection_buffer.irp.f +++ b/plugins/Full_CI_ZMQ/selection_buffer.irp.f @@ -25,7 +25,7 @@ subroutine add_to_selection_buffer(b, det, val) double precision, intent(in) :: val integer :: i - if(dabs(val) >= b%mini) then + if(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 @@ -41,33 +41,25 @@ subroutine sort_selection_buffer(b) implicit none type(selection_buffer), intent(inout) :: b - double precision, allocatable:: absval(:) integer, allocatable :: iorder(:) - double precision, pointer :: vals(:) integer(bit_kind), pointer :: detmp(:,:,:) integer :: i, nmwen logical, external :: detEq nmwen = min(b%N, b%cur) - - allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3)), absval(b%cur), vals(size(b%val))) - absval = b%val(:b%cur) + allocate(iorder(b%cur), detmp(N_int, 2, size(b%det,3))) do i=1,b%cur iorder(i) = i end do - call dsort(absval, iorder, b%cur) + call dsort(b%val, iorder, b%cur) do i=1, nmwen detmp(1:N_int,1,i) = b%det(1:N_int,1,iorder(i)) detmp(1:N_int,2,i) = b%det(1:N_int,2,iorder(i)) - vals(i) = b%val(iorder(i)) end do - do i=nmwen+1, size(vals) - vals(i) = 0.d0 - enddo - deallocate(b%det, b%val) + deallocate(b%det) b%det => detmp - b%val => vals - b%mini = max(b%mini,dabs(b%val(b%N))) + b%mini = min(b%mini,b%val(b%N)) b%cur = nmwen + deallocate(iorder) end subroutine diff --git a/src/Utils/sort.irp.f b/src/Utils/sort.irp.f index 1ebf3b17..a9594d6c 100644 --- a/src/Utils/sort.irp.f +++ b/src/Utils/sort.irp.f @@ -27,6 +27,56 @@ BEGIN_TEMPLATE enddo end subroutine insertion_$Xsort + subroutine quick_$Xsort(x, iorder, isize) + implicit none + BEGIN_DOC + ! Sort array x(isize) using the quicksort algorithm. + ! iorder in input should be (1,2,3,...,isize), and in output + ! contains the new order of the elements. + END_DOC + integer,intent(in) :: isize + $type,intent(inout) :: x(isize) + integer,intent(inout) :: iorder(isize) + call rec_$X_quicksort(x,iorder,isize,1,isize) + end + + recursive subroutine rec_$X_quicksort(x, iorder, isize, first, last) + implicit none + integer, intent(in) :: isize, first, last + integer,intent(inout) :: iorder(isize) + $type, intent(inout) :: x(isize) + $type :: c, tmp + integer :: itmp + integer :: i, j + + c = x( ishft(first+last,-1) ) + i = first + j = last + do + do while (x(i) < c) + i=i+1 + end do + do while (c < x(j)) + j=j-1 + end do + if (i >= j) exit + tmp = x(i) + x(i) = x(j) + x(j) = tmp + itmp = iorder(i) + iorder(i) = iorder(j) + iorder(j) = itmp + i=i+1 + j=j-1 + enddo + if (first < i-1) then + call rec_$X_quicksort(x, iorder, isize, first, i-1) + endif + if (j+1 < last) then + call rec_$X_quicksort(x, iorder, isize, j+1, last) + endif + end + subroutine heap_$Xsort(x,iorder,isize) implicit none BEGIN_DOC @@ -206,10 +256,11 @@ BEGIN_TEMPLATE ! if (isize == n) then ! return ! endif - if ( isize < 16) then + if ( isize < 32) then call insertion_$Xsort(x,iorder,isize) else - call heap_$Xsort(x,iorder,isize) +! call heap_$Xsort(x,iorder,isize) + call quick_$Xsort(x,iorder,isize) endif end subroutine $Xsort