Debugged Davidson for large Ndet

src/Davidson/u0Hu0.irp.f

@@ -665,7 +665,7 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8)
   integer, allocatable           :: singles_b(:,:)
   integer, allocatable           :: idx(:), idx0(:)
   logical, allocatable           :: is_single_a(:)
-  integer                        :: maxab, n_singles_max, kcol_prev
+  integer                        :: maxab, n_singles_max, kcol_prev, nmax
   double precision, allocatable  :: u_t(:,:), v_t(:,:), s_t(:,:)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: v_t, s_t, u_t
 
@@ -702,7 +702,7 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8)
       !$OMP          singles_alpha_size, sze_8,                      &
       !$OMP          idx0, u_t, maxab, v_0, s_0)                     &
       !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,     &
-      !$OMP          lcol, lrow, is_single_a,l_a, l_b,               &
+      !$OMP          lcol, lrow, is_single_a,l_a, l_b, nmax,         &
       !$OMP          buffer, singles, doubles, n_singles, n_doubles, &
       !$OMP          tmp_det2, hij, sij, idx, l, kcol_prev, v_t, s_t)
   
@@ -756,11 +756,12 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8)
       tmp_det2(1:N_int,2) = psi_det_beta_unique(1:N_int, lcol)
 
       l_a = psi_bilinear_matrix_columns_loc(lcol)
-      do while (l_a <= k_a)
-        l_a += 1
-      enddo
 
       n_doubles=1
+
+      ! Loop over alpha singles
+      ! -----------------------
+
       do while ( l_a < psi_bilinear_matrix_columns_loc(lcol+1) )
         lrow = psi_bilinear_matrix_rows(l_a)
         if (.not.is_single_a(lrow)) then
@@ -778,8 +779,8 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8)
               call get_s2(tmp_det,tmp_det2,N_int,sij)
               do l=1,N_st
                 v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
-                s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,l_a)
                 v_t(l,l_a) = v_t(l,l_a) + hij * u_t(l,k_a)
+                s_t(l,k_a) = s_t(l,k_a) + sij * u_t(l,l_a)
                 s_t(l,l_a) = s_t(l,l_a) + sij * u_t(l,k_a)
               enddo
             enddo
@@ -811,6 +812,10 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8)
   enddo
   !$OMP END DO NOWAIT
 
+
+  ! Single and double alpha excitations
+  ! ===================================
+    
   !$OMP DO SCHEDULE(static,1)
   do k_a=1,N_det
     
@@ -828,23 +833,18 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8)
     
     k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
 
-    ! Get all single and double alpha excitations
-    ! ===========================================
-    
     spindet(1:N_int) = tmp_det(1:N_int,1)
     
     ! Loop inside the beta column to gather all the connected alphas
-    i=1
+    l_a = k_a+1
-    l_a = k_a
+    nmax = min(N_det_alpha_unique, N_det - l_a)
-    lcol = psi_bilinear_matrix_columns(l_a)
+    do i=1,nmax
-    do while (lcol == kcol)
+      lcol = psi_bilinear_matrix_columns(l_a)
+      if (lcol /= kcol) exit
       lrow = psi_bilinear_matrix_rows(l_a)
       buffer(1:N_int,i) = psi_det_alpha_unique(1:N_int, lrow)
       idx(i) = l_a
-      i   = i  +1
       l_a = l_a+1
-      if (l_a > N_det) exit
-      lcol = psi_bilinear_matrix_columns(l_a)
     enddo
     i = i-1
     
@@ -852,7 +852,6 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8)
         buffer, idx, spindet, N_int, i,                              &
         singles, doubles, n_singles, n_doubles )
 
-
     ! Compute Hij for all alpha singles
     ! ----------------------------------
 
@@ -884,26 +883,38 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8)
       enddo
     enddo
     
+  end do
+  !$OMP END DO NOWAIT
+
+
+  ! Single and double beta excitations
+  ! ==================================
+
+  !$OMP DO SCHEDULE(static,1)
+  do k_b=1,N_det
+
+    ! Initial determinant is at k_b in beta-major representation
+    ! -----------------------------------------------------------------------
     
+    krow = psi_bilinear_matrix_transp_rows(k_b)
+    kcol = psi_bilinear_matrix_transp_columns(k_b)
     
-    ! Get all single and double beta excitations
+    tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int, krow)
-    ! ===========================================
+    tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int, kcol)
     
     spindet(1:N_int) = tmp_det(1:N_int,2)
+    k_a = psi_bilinear_matrix_transp_order(k_b)
     
     ! Loop inside the alpha row to gather all the connected betas
-    i=1
+    l_b = k_b+1
-    l_b = k_b
+    nmax = min(N_det_beta_unique, N_det - l_b)
-
+    do i=1,nmax
-    lrow = psi_bilinear_matrix_transp_rows(l_b)
+      lrow = psi_bilinear_matrix_transp_rows(l_b)
-    do while (lrow == krow)
+      if (lrow /= krow) exit
       lcol = psi_bilinear_matrix_transp_columns(l_b)
       buffer(1:N_int,i) = psi_det_beta_unique(1:N_int, lcol)
       idx(i) = l_b
-      i   = i  +1
       l_b = l_b+1
-      if (l_b > N_det) exit
-      lrow = psi_bilinear_matrix_transp_rows(l_b)
     enddo
     i = i-1
   
@@ -943,9 +954,25 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8)
       enddo
     enddo
 
-    ! Diagonal contribution
+  end do
-    ! ---------------------
+  !$OMP END DO NOWAIT
 
+
+  ! Diagonal contribution
+  ! =====================
+
+  !$OMP DO SCHEDULE(static,1)
+  do k_a=1,N_det
+    
+    ! Initial determinant is at k_a in alpha-major representation
+    ! -----------------------------------------------------------------------
+    
+    krow = psi_bilinear_matrix_rows(k_a)
+    kcol = psi_bilinear_matrix_columns(k_a)
+    
+    tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int, krow)
+    tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int, kcol)
+    
     double precision, external :: diag_H_mat_elem, diag_S_mat_elem
   
     hij = diag_H_mat_elem(tmp_det,N_int) 
@@ -956,7 +983,7 @@ subroutine H_S2_u_0_nstates_bilinear_order(v_0,s_0,u_0,N_st,sze_8)
     enddo
 
   end do
-  !$OMP END DO
+  !$OMP END DO NOWAIT
 
   !$OMP CRITICAL
   do l=1,N_st
@@ -986,7 +1013,7 @@ subroutine H_S2_u_0_nstates_test(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze
 
   double precision, allocatable :: vt(:,:)
   integer, allocatable :: idx(:)
-  integer                        :: i,j, jj
+  integer                        :: i,j, jj, l
   double precision               :: hij
 
   do i=1,n
@@ -995,6 +1022,7 @@ subroutine H_S2_u_0_nstates_test(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze
 
   allocate(idx(0:n), vt(N_st,n))
   Vt = 0.d0
+  !$OMP PARALLEL DO DEFAULT(shared) PRIVATE(i,idx,jj,j,degree,exc,phase,hij,l) SCHEDULE(static,1)
   do i=2,n
     idx(0) = i
     call filter_connected(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
@@ -1012,11 +1040,19 @@ subroutine H_S2_u_0_nstates_test(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze
 !       if ((degree == 2).and.(exc(0,1,1)==1)) cycle
 !      if ((degree > 1)) cycle
 !      if (exc(0,1,2) /= 0) cycle
+!      if (exc(0,1,1) == 2) cycle
+!      if (exc(0,1,2) == 2) cycle
+!      if ((degree==1).and.(exc(0,1,2) == 1)) cycle
       call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
-      vt (:,i) = vt (:,i) + hij*u_0(j,:)
+      do l=1,N_st
-      vt (:,j) = vt (:,j) + hij*u_0(i,:)
+        !$OMP ATOMIC
+        vt (l,i) = vt (l,i) + hij*u_0(j,l)
+        !$OMP ATOMIC
+        vt (l,j) = vt (l,j) + hij*u_0(i,l)
+      enddo
     enddo
   enddo
+  !$OMP END PARALLEL DO
   do i=1,n
     v_0(i,:) = v_0(i,:) + vt(:,i)
   enddo

src/Determinants/spindeterminants.irp.f

@@ -389,8 +389,6 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_rows   , (N_det) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns, (N_det) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order  , (N_det) ]
-&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_reverse , (N_det) ]
-&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns_loc, (N_det_beta_unique+1) ]
   use bitmasks
   implicit none
   BEGIN_DOC
@@ -408,7 +406,7 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
 
   PROVIDE psi_coef_sorted_bit
 
-  integer, allocatable :: to_sort(:)
+  integer*8, allocatable :: to_sort(:)
   integer, external :: get_index_in_psi_det_alpha_unique
   integer, external :: get_index_in_psi_det_beta_unique
   allocate(to_sort(N_det))
@@ -421,16 +419,47 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
     enddo
     psi_bilinear_matrix_rows(k) = i
     psi_bilinear_matrix_columns(k) = j
-    to_sort(k) = N_det_alpha_unique * (j-1) + i
+    to_sort(k) = int(N_det_alpha_unique,8) * int(j-1,8) + int(i,8)
     psi_bilinear_matrix_order(k) = k
   enddo
-  call isort(to_sort, psi_bilinear_matrix_order, N_det)
+  call i8sort(to_sort, psi_bilinear_matrix_order, N_det)
   call iset_order(psi_bilinear_matrix_rows,psi_bilinear_matrix_order,N_det)
   call iset_order(psi_bilinear_matrix_columns,psi_bilinear_matrix_order,N_det)
   do l=1,N_states
     call dset_order(psi_bilinear_matrix_values(1,l),psi_bilinear_matrix_order,N_det)
   enddo
-  psi_bilinear_matrix_columns_loc(1) = 1
+  deallocate(to_sort)
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_reverse , (N_det) ]
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+! Order which allors to go from psi_bilinear_matrix to psi_det
+  END_DOC
+  integer                        :: k
+  do k=1,N_det
+    psi_bilinear_matrix_order_reverse(psi_bilinear_matrix_order(k)) = k
+  enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns_loc, (N_det_beta_unique+1) ]
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
+!  D_a^t C D_b
+!
+! Rows are alpha determinants and columns are beta.
+!
+! Order refers to psi_det
+  END_DOC
+  integer                        :: i,j,k, l
+
+  l = psi_bilinear_matrix_columns(1)
+  psi_bilinear_matrix_columns_loc(l) = 1
   do k=2,N_det
     if (psi_bilinear_matrix_columns(k) == psi_bilinear_matrix_columns(k-1)) then
       cycle
@@ -440,35 +469,27 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_values, (N_det,N_states)
     endif
   enddo
   psi_bilinear_matrix_columns_loc(N_det_beta_unique+1) = N_det+1
-  do k=1,N_det
-    psi_bilinear_matrix_order_reverse(psi_bilinear_matrix_order(k)) = k
-  enddo
-  deallocate(to_sort)
 END_PROVIDER
 
-
 BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_states) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows   , (N_det) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_columns, (N_det) ]
 &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_order  , (N_det) ]
-&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_transp_reverse , (N_det) ]
   use bitmasks
   implicit none
   BEGIN_DOC
-! Sparse coefficient matrix if the wave function is expressed in a bilinear form :
+! Transpose of psi_bilinear_matrix
-!  D_a^t C D_b
+!  D_b^t C^t D_a
 !
 ! Rows are Alpha determinants and columns are beta, but the matrix is stored in row major
 ! format
-!
-! Order refers to psi_bilinear_matrix
   END_DOC
   integer                        :: i,j,k,l
 
 
   PROVIDE psi_coef_sorted_bit
 
-  integer, allocatable :: to_sort(:)
+  integer*8, allocatable :: to_sort(:)
   allocate(to_sort(N_det))
   do l=1,N_states
     do k=1,N_det
@@ -480,19 +501,50 @@ BEGIN_PROVIDER  [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_
     psi_bilinear_matrix_transp_rows   (k) = psi_bilinear_matrix_rows   (k)
     i = psi_bilinear_matrix_transp_columns(k) 
     j = psi_bilinear_matrix_transp_rows   (k)
-    to_sort(k) = N_det_beta_unique * (j-1) + i
+    to_sort(k) = int(N_det_beta_unique,8) * int(j-1,8) + int(i,8)
     psi_bilinear_matrix_transp_order(k) = k
   enddo
-  call isort(to_sort, psi_bilinear_matrix_transp_order, N_det)
+  call i8sort(to_sort, psi_bilinear_matrix_transp_order, N_det)
   call iset_order(psi_bilinear_matrix_transp_rows,psi_bilinear_matrix_transp_order,N_det)
   call iset_order(psi_bilinear_matrix_transp_columns,psi_bilinear_matrix_transp_order,N_det)
   do l=1,N_states
     call dset_order(psi_bilinear_matrix_transp_values(1,l),psi_bilinear_matrix_transp_order,N_det)
   enddo
+  deallocate(to_sort)
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows_loc, (N_det_alpha_unique+1) ]
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+! Location of the columns in the psi_bilinear_matrix
+  END_DOC
+  integer                        :: i,j,k, l
+
+  l = psi_bilinear_matrix_transp_rows(1)
+  psi_bilinear_matrix_transp_rows_loc(l) = 1
+  do k=2,N_det
+    if (psi_bilinear_matrix_transp_rows(k) == psi_bilinear_matrix_transp_rows(k-1)) then
+      cycle
+    else
+      l = psi_bilinear_matrix_transp_rows(k)
+      psi_bilinear_matrix_transp_rows_loc(l) = k
+    endif
+  enddo
+  psi_bilinear_matrix_transp_rows_loc(N_det_beta_unique+1) = N_det+1
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order_transp_reverse , (N_det) ]
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+! Order which allows to go from psi_bilinear_matrix_order_transp to psi_bilinear_matrix
+  END_DOC
+  integer                        :: k
+
   do k=1,N_det
     psi_bilinear_matrix_order_transp_reverse(psi_bilinear_matrix_transp_order(k)) = k
   enddo
-  deallocate(to_sort)
 END_PROVIDER
 
 
@@ -641,17 +693,17 @@ subroutine get_all_spin_singles_and_doubles(buffer, idx, spindet, Nint, size_buf
 
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xorvec, degree
 
-!  select case (Nint)
+  select case (Nint)
-!    case (1)
+    case (1)
-!      call get_all_spin_singles_and_doubles_1(buffer, idx, spindet(1), size_buffer, singles, doubles, n_singles, n_doubles)
+      call get_all_spin_singles_and_doubles_1(buffer, idx, spindet(1), size_buffer, singles, doubles, n_singles, n_doubles)
-!      return
+      return
-!    case (2)
+    case (2)
-!      call get_all_spin_singles_and_doubles_2(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
+      call get_all_spin_singles_and_doubles_2(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
-!      return
+      return
-!    case (3)
+    case (3)
-!      call get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
+      call get_all_spin_singles_and_doubles_3(buffer, idx, spindet, size_buffer, singles, doubles, n_singles, n_doubles)
-!      return
+      return
-!  end select
+  end select
 
 
   size_buffer_align = align_double(size_buffer)