Accelerated S2

src/Determinants/H_apply.irp.f

@@ -189,11 +189,15 @@ subroutine copy_H_apply_buffer_to_wf
   !$OMP BARRIER
   H_apply_buffer(j)%N_det = 0
   !$OMP END PARALLEL
-  call normalize(psi_coef,N_det)
   SOFT_TOUCH N_det psi_det psi_coef
   
   logical :: found_duplicates
   call remove_duplicates_in_psi_det(found_duplicates)
+  do k=1,N_states
+    call normalize(psi_coef(1,k),N_det)
+  enddo
+  SOFT_TOUCH N_det psi_det psi_coef
+  
 end
 
 subroutine remove_duplicates_in_psi_det(found_duplicates)
@@ -230,7 +234,7 @@ subroutine remove_duplicates_in_psi_det(found_duplicates)
     j = i+1
     do while (bit_tmp(j)==bit_tmp(i))
       if (duplicate(j)) then
-        j += 1
+        j = j+1
         if (j > N_det) then
           exit
         else
@@ -246,8 +250,10 @@ subroutine remove_duplicates_in_psi_det(found_duplicates)
         endif
       enddo
       if (dup) then
+        !$OMP CRITICAL
         duplicate(j) = .True.
         found_duplicates = .True.
+        !$OMP END CRITICAL
       endif
       j += 1
       if (j > N_det) then
@@ -265,17 +271,20 @@ subroutine remove_duplicates_in_psi_det(found_duplicates)
         k += 1
         psi_det(:,:,k) = psi_det_sorted_bit (:,:,i)
         psi_coef(k,:)  = psi_coef_sorted_bit(i,:)
-!      else
+      else
-!        call debug_det(psi_det_sorted_bit(1,1,i),N_int)
+        if (sum(abs(psi_coef_sorted_bit(i,:))) /= 0.d0 ) then
-!        stop 'duplicates in psi_det'
+          psi_coef(k,:)  = psi_coef_sorted_bit(i,:)
+        endif
       endif
     enddo
     N_det = k
-    call write_bool(6,found_duplicates,'Found duplicate determinants')
     psi_det_sorted_bit(:,:,1:N_det) = psi_det(:,:,1:N_det)
     psi_coef_sorted_bit(1:N_det,:) = psi_coef(1:N_det,:)
-    SOFT_TOUCH N_det psi_det psi_coef psi_det_sorted_bit psi_coef_sorted_bit
+    TOUCH N_det psi_det psi_coef psi_det_sorted_bit psi_coef_sorted_bit
   endif
+  psi_det = psi_det_sorted
+  psi_coef = psi_coef_sorted
+  TOUCH psi_det psi_coef psi_det_sorted_bit psi_coef_sorted_bit
   deallocate (duplicate,bit_tmp)
 end
 

src/Determinants/connected_to_ref.irp.f

@@ -7,11 +7,8 @@ integer*8 function det_search_key(det,Nint)
   integer, intent(in) :: Nint
   integer(bit_kind), intent(in) :: det(Nint,2)
   integer :: i
-  det_search_key = iand(det(1,1),det(1,2))
+  i = shiftr(elec_alpha_num, bit_kind_shift)+1
-  do i=2,Nint
+  det_search_key = shiftr(ior(det(i,1),det(i,2)),1)
-    det_search_key = ieor(det_search_key,iand(det(i,1),det(i,2)))
-  enddo
-  det_search_key = iand(huge(det(1,1)),det_search_key)
 end
 
 
@@ -24,11 +21,8 @@ integer*8 function occ_pattern_search_key(det,Nint)
   integer, intent(in) :: Nint
   integer(bit_kind), intent(in) :: det(Nint,2)
   integer :: i
-  occ_pattern_search_key = ieor(det(1,1),det(1,2))
+  i = shiftr(elec_alpha_num, bit_kind_shift)+1
-  do i=2,Nint
+  occ_pattern_search_key = shiftr(ior(det(i,1),det(i,2)),1)
-    occ_pattern_search_key = ieor(occ_pattern_search_key,iand(det(i,1),det(i,2)))
-  enddo
-  occ_pattern_search_key = iand(huge(det(1,1)),occ_pattern_search_key)
 end
 
 

src/Determinants/determinants.irp.f

@@ -336,19 +336,19 @@ end subroutine
    ! function.
    END_DOC
    
-   call sort_dets_by_det_search_key(N_det, psi_det, psi_coef,        &
+   call sort_dets_by_det_search_key(N_det, psi_det, psi_coef, size(psi_coef,1),       &
        psi_det_sorted_bit, psi_coef_sorted_bit, N_states)
 
 END_PROVIDER
  
-subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, det_out, coef_out, N_st)
+subroutine sort_dets_by_det_search_key(Ndet, det_in, coef_in, sze, det_out, coef_out, N_st)
    use bitmasks
    implicit none
-   integer, intent(in)            :: Ndet, N_st
+   integer, intent(in)            :: Ndet, N_st, sze
-   integer(bit_kind), intent(in)  :: det_in  (N_int,2,psi_det_size)
+   integer(bit_kind), intent(in)  :: det_in  (N_int,2,sze)
-   double precision , intent(in)  :: coef_in(psi_det_size,N_st)
+   double precision , intent(in)  :: coef_in(sze,N_st)
-   integer(bit_kind), intent(out) :: det_out (N_int,2,psi_det_size)
+   integer(bit_kind), intent(out) :: det_out (N_int,2,sze)
-   double precision , intent(out) :: coef_out(psi_det_size,N_st)
+   double precision , intent(out) :: coef_out(sze,N_st)
    BEGIN_DOC
    ! Determinants are sorted are sorted according to their det_search_key.
    ! Useful to accelerate the search of a random determinant in the wave

src/Determinants/occ_pattern.irp.f

@@ -3,7 +3,12 @@ subroutine det_to_occ_pattern(d,o,Nint)
   use bitmasks
   implicit none
   BEGIN_DOC
-  ! Transform a determinant to an occupation pattern
+  ! Transforms a determinant to an occupation pattern
+  !
+  ! occ(:,1) : Single occupations
+  ! 
+  ! occ(:,2) : Double occupations
+  !
   END_DOC
   integer          ,intent(in)   :: Nint
   integer(bit_kind),intent(in)   :: d(Nint,2)
@@ -17,6 +22,7 @@ subroutine det_to_occ_pattern(d,o,Nint)
   enddo
 end
 
+
 subroutine occ_pattern_to_dets_size(o,sze,n_alpha,Nint)
   use bitmasks
   implicit none
@@ -29,120 +35,77 @@ subroutine occ_pattern_to_dets_size(o,sze,n_alpha,Nint)
   integer                        :: amax,bmax,k
   double precision, external     :: binom_func
 
-  amax = n_alpha
   bmax = 0
+  amax = n_alpha
   do k=1,Nint
     bmax += popcnt( o(k,1) )
     amax -= popcnt( o(k,2) )
   enddo
-  sze = int( min(binom_func(bmax, amax), 1.d8) )
+  sze = binom_int(bmax, amax)
-  sze = 2*sze*sze + 64
-
 end
 
+
 subroutine occ_pattern_to_dets(o,d,sze,n_alpha,Nint)
   use bitmasks
   implicit none
   BEGIN_DOC
   ! Generate all possible determinants for a give occ_pattern
   END_DOC
-  integer          ,intent(in)   :: Nint, n_alpha
+  integer          ,intent(in)   :: Nint
-  integer         ,intent(inout) :: sze
+  integer          ,intent(in)   :: n_alpha        ! Number of alpha electrons
-  integer(bit_kind),intent(in)   :: o(Nint,2)
+  integer         ,intent(inout) :: sze            ! Dimension of the output dets
-  integer(bit_kind),intent(out)  :: d(Nint,2,sze)
+  integer(bit_kind),intent(in)   :: o(Nint,2)      ! Occ patters
+  integer(bit_kind),intent(out)  :: d(Nint,2,sze)  ! Output determinants
   
-  integer                        :: i, l, k, nt, na, nd, amax
+  integer                        :: i, k, n, ispin
-  integer                        :: list_todo(2*n_alpha)
-  integer                        :: list_a(2*n_alpha)
-  integer                        :: ishift
 
-  amax = n_alpha
+  ! Extract list of singly occupied MOs as (int,pos) pairs
-  do k=1,Nint
+  ! ------------------------------------------------------
-    amax -= popcnt( o(k,2) )
+
+  integer           :: iint(2*n_alpha), ipos(2*n_alpha)
+  integer(bit_kind) :: v, t, tt
+  integer           :: n_alpha_in_single
+
+  n=0
+  n_alpha_in_single = n_alpha
+  do i=1,Nint
+    v = o(i,1)
+    do while(v /= 0_bit_kind)
+      n = n+1
+      iint(n) = i
+      ipos(n) = trailz(v)
+      v = iand(v,v-1)
+    enddo
+    n_alpha_in_single = n_alpha_in_single - popcnt( o(i,2) )
   enddo
 
-  call bitstring_to_list(o(1,1), list_todo, nt, Nint)
+  v = shiftl(1,n_alpha_in_single) - 1
+  sze = binom_int(n,n_alpha_in_single)
+  do i=1,sze
+    ! Initialize with doubly occupied MOs
+    d(:,1,i) = o(:,2)
+    d(:,2,i) = o(:,2)
 
-  na = 0
+    do k=1,n
-  nd = 0
+      if (btest(v,k-1)) then
-  d = 0
+        ispin = 1
-  call rec_occ_pattern_to_dets(list_todo,nt,list_a,na,d,nd,sze,amax,Nint)
+      else
-
+        ispin = 2
-  sze = nd
-  
-  integer :: ne(2)
-  l=0
-  do i=1,nd
-    ne(1) = 0
-    ne(2) = 0
-    l=l+1
-    ! Doubly occupied orbitals
-    do k=1,Nint
-      d(k,1,l) = ior(d(k,1,i),o(k,2))
-      d(k,2,l) = ior(d(k,2,i),o(k,2))
-      ne(1) += popcnt(d(k,1,l))
-      ne(2) += popcnt(d(k,2,l))
-    enddo
-    if ( (ne(1) /= elec_alpha_num).or.(ne(2) /= elec_beta_num) ) then
-      l = l-1
-    endif
-  enddo
-  sze = l
-
-end
-
-
-recursive subroutine  rec_occ_pattern_to_dets(list_todo,nt,list_a,na,d,nd,sze,amax,Nint)
-  use bitmasks
-  implicit none
-
-  integer, intent(in)            :: nt, sze, amax, Nint,na
-  integer,intent(inout)          :: list_todo(nt)
-  integer, intent(inout)         :: list_a(na+1),nd
-  integer(bit_kind),intent(inout) :: d(Nint,2,sze)
-  integer :: iint, ipos, i,j,k
-
-  if (na == amax) then
-    nd += 1
-    if (nd > sze) then
-      print *,  irp_here, ': nd  = ', nd
-      print *,  irp_here, ': sze = ', sze
-      stop 'bug in rec_occ_pattern_to_dets'
-    endif
-    do i=1,na
-      iint = ishft(list_a(i)-1,-bit_kind_shift) + 1
-      ipos = list_a(i)-ishft((iint-1),bit_kind_shift)-1
-      d(iint,1,nd) = ibset( d(iint,1,nd), ipos )
-    enddo
-    do i=1,nt
-      iint = ishft(list_todo(i)-1,-bit_kind_shift) + 1
-      ipos = list_todo(i)-ishft((iint-1),bit_kind_shift)-1
-      d(iint,2,nd) = ibset( d(iint,2,nd), ipos )
-    enddo
-  else
-    integer, allocatable :: list_todo_tmp(:)
-    allocate (list_todo_tmp(nt))
-    do i=1,nt
-      if (na > 0) then
-        if (list_todo(i) < list_a(na)) then
-          cycle
-        endif
       endif
-      list_a(na+1) = list_todo(i)
+      d(iint(k), ispin, i) = ibset( d(iint(k), ispin, i), ipos(k) )
-      k=1
-      do j=1,nt
-        if (i/=j) then
-          list_todo_tmp(k) = list_todo(j)
-          k += 1
-        endif
-      enddo
-      call rec_occ_pattern_to_dets(list_todo_tmp,nt-1,list_a,na+1,d,nd,sze,amax,Nint)
     enddo
-    deallocate(list_todo_tmp)
+
-  endif
+    ! Generate next permutation with Anderson's algorithm
+    t = ior(v,v-1)
+    tt = t+1
+    v = ior(tt, shiftr( and(not(t),tt) - 1, trailz(v)+1) )
+  enddo
+
 
 end
 
+
+
  BEGIN_PROVIDER [ integer(bit_kind), psi_occ_pattern, (N_int,2,psi_det_size) ]
 &BEGIN_PROVIDER [ integer, N_occ_pattern ]
  implicit none
@@ -348,13 +311,14 @@ subroutine make_s2_eigenfunction
     endif
     call occ_pattern_to_dets(psi_occ_pattern(1,1,i),d,s,elec_alpha_num,N_int)
     do j=1,s
-      if (.not. is_in_wavefunction(d(1,1,j), N_int) ) then
+      if ( is_in_wavefunction(d(1,1,j), N_int) ) then
-        N_det_new += 1
+        cycle
-        det_buffer(:,:,N_det_new) = d(:,:,j)
+      endif
-        if (N_det_new == bufsze) then
+      N_det_new += 1
-          call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,ithread)
+      det_buffer(:,:,N_det_new) = d(:,:,j)
-          N_det_new = 0
+      if (N_det_new == bufsze) then
-        endif
+        call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,ithread)
+        N_det_new = 0
       endif
     enddo
   enddo

src/Determinants/psi_cas.irp.f

@@ -53,7 +53,7 @@ END_PROVIDER
  ! CAS determinants sorted to accelerate the search of a random determinant in the wave
  ! function.
  END_DOC
- call sort_dets_by_det_search_key(N_det_cas, psi_cas, psi_cas_coef, &
+ call sort_dets_by_det_search_key(N_det_cas, psi_cas, psi_cas_coef, size(psi_cas_coef,1), &
      psi_cas_sorted_bit, psi_cas_coef_sorted_bit, N_states)
 
 END_PROVIDER
@@ -106,7 +106,7 @@ END_PROVIDER
  ! CAS determinants sorted to accelerate the search of a random determinant in the wave
  ! function.
  END_DOC
- call sort_dets_by_det_search_key(N_det_cas, psi_non_cas, psi_non_cas_coef, &
+ call sort_dets_by_det_search_key(N_det_cas, psi_non_cas, psi_non_cas_coef, size(psi_non_cas_coef,1), &
      psi_non_cas_sorted_bit, psi_non_cas_coef_sorted_bit, N_states)
 
 END_PROVIDER

src/FCI/selection_buffer.irp.f

@@ -137,3 +137,42 @@ subroutine sort_selection_buffer(b)
   b%cur = nmwen
 end subroutine
 
+subroutine make_selection_buffer_s2(b)
+  use selection_types
+  type(selection_buffer), intent(inout) :: b
+
+  integer(bit_kind), pointer :: o(:,:,:)
+  double precision, pointer  :: val(:)
+
+  integer :: n_d
+  integer :: i,k,sze,n_alpha,j,n
+ 
+  ! Create occupation patterns
+  n_d = 0
+  allocate(o(N_int,2,b%cur))
+  do i=1,b%cur
+    do k=1,N_int
+      o(k,1,i) = ieor(b%det(k,1,i), b%det(k,2,i))
+      o(k,2,i) = iand(b%det(k,1,i), b%det(k,2,i))
+    enddo
+    call occ_pattern_to_dets_size(o(1,1,i),sze,elec_alpha_num,N_int)
+    n_d = n_d + sze
+  enddo
+  deallocate(b%det)
+
+  allocate(b%det(N_int,2,n_d), val(n_d))
+  k=1
+  do i=1,b%cur
+    n=n_d
+    call occ_pattern_to_dets(o(1,1,i),b%det(1,1,k),n,elec_alpha_num,N_int)
+    val(k) = b%val(i)
+    do j=k+1,k+n-1
+      val(j) = 0.d0
+    enddo
+    k = k+n
+  enddo
+  deallocate(o,b%val)
+  b%val => val
+  b%N = k-1
+  b%cur = k-1
+end

src/FCI/zmq_selection.irp.f

@@ -120,11 +120,11 @@ subroutine ZMQ_selection(N_in, pt2, variance, norm)
     norm(i) = 0.d0
   enddo
   if (N_in > 0) then
+    if (s2_eig.or.(N_states > 1) ) then
+      call make_selection_buffer_s2(b)
+    endif
     call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) 
     call copy_H_apply_buffer_to_wf()
-    if (s2_eig.or.(N_states > 1) ) then
-      call make_s2_eigenfunction
-    endif
     call save_wavefunction
   endif
   call delete_selection_buffer(b)

src/Psiref_Utils/psi_ref_utils.irp.f

@@ -151,7 +151,7 @@ END_PROVIDER
  ! Reference determinants sorted to accelerate the search of a random determinant in the wave
  ! function.
  END_DOC
- call sort_dets_by_det_search_key(N_det_ref, psi_non_ref, psi_non_ref_coef, &
+ call sort_dets_by_det_search_key(N_det_ref, psi_non_ref, psi_non_ref_coef, size(psi_non_ref_coef,1), &
      psi_non_ref_sorted_bit, psi_non_ref_coef_sorted_bit, N_states)
 
 END_PROVIDER

src/Utils/util.irp.f

@@ -28,6 +28,27 @@ double precision function binom_func(i,j)
 end
 
 
+FUNCTION n_combinations(n,k) RESULT(r)
+    IMPLICIT NONE
+
+    INTEGER(4), INTENT(IN) :: n,k
+    INTEGER(8)             :: r
+    INTEGER(4)             :: d, n0
+
+    IF (k > n) THEN
+        r = 0
+        RETURN
+    ELSE
+        r = 1
+    ENDIF
+
+    n0 = n
+    DO d=1, k
+        r = (r*n0) / d
+        n0 = n0 - 1
+    ENDDO
+END FUNCTION n_combinations
+
  BEGIN_PROVIDER [ double precision, binom, (0:40,0:40) ]
 &BEGIN_PROVIDER [ double precision, binom_transp, (0:40,0:40) ]
   implicit none
@@ -45,6 +66,22 @@ end
 END_PROVIDER
 
 
+ BEGIN_PROVIDER [ integer*8, binom_int, (0:40,0:40) ]
+&BEGIN_PROVIDER [ integer*8, binom_int_transp, (0:40,0:40) ]
+  implicit none
+  BEGIN_DOC
+  ! Binomial coefficients, as integers*8
+  END_DOC
+  integer                        :: k,l
+  double precision               :: logfact
+  do l=0,40
+    do k=0,40
+      binom_int(k,l) = int(binom(k,l)+0.1d0,8)
+    enddo
+  enddo
+END_PROVIDER
+
+
 
 double precision function fact(n)
   implicit none