i_h_j_s2 - optimization in mrcc_sto

plugins/dress_zmq/dressing_vector.irp.f

@@ -1,8 +1,3 @@
- BEGIN_PROVIDER [ integer, nalp ]
-&BEGIN_PROVIDER [ integer, ninc ]
-   nalp = 0
-   ninc = 0
- END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ]
 &BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ]
@@ -30,8 +25,5 @@
    tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det)
    dressing_column_s(l,k) -= tmp * f
  enddo
- print *, "NALP", nalp
- print *, "NINC", ninc
- print *, "DELTA_IJ", dressing_column_h(:10,1)
 END_PROVIDER
 

plugins/mrcc_sto/mrcc_sto.irp.f

@@ -12,6 +12,8 @@ end
 &BEGIN_PROVIDER [ double precision, sij_cache_, (N_det,Nproc) ]
 &BEGIN_PROVIDER [ double precision, dIa_hla_, (N_states,N_det,Nproc) ]
 &BEGIN_PROVIDER [ double precision, dIa_sla_, (N_states,N_det,Nproc) ]
+&BEGIN_PROVIDER [ integer, excs_ , (0:2,2,2,N_det,Nproc) ]
+&BEGIN_PROVIDER [ double precision, phases_, (N_det, Nproc) ]
 BEGIN_DOC
  ! temporay arrays for dress_with_alpha_buffer. Avoids realocation.
 END_DOC
@@ -46,8 +48,7 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil
   double precision :: Delta_E_inv(N_states)
   double precision :: sdress, hdress
   logical :: ok, ok2
-  integer :: old_ninc
+  integer :: canbediamond
-  double precision :: shdress
   PROVIDE mo_class
 
   
@@ -57,28 +58,33 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil
     if(idx_non_ref_rev(minilist(i)) == 0) return
   end do
 
-  shdress = 0d0
-  old_ninc = ninc
-
   if (perturbative_triples) then
     PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
   endif
-  
+  canbediamond = 0
   do l_sd=1,n_minilist
-   !call get_excitation(det_minilist(1,1,l_sd),alpha,exc,degree1,phase,N_int)
+   call get_excitation(det_minilist(1,1,l_sd),alpha,exc,degree1,phase,N_int)
-   !if(degree1 == 0 .or. degree1 > 2) stop "minilist error"
+   call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2)
-   !call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2)
-   !
-   !ok = (mo_class(h1)(1:1) == 'A' .or. mo_class(h1)(1:1) == 'I') .and. &
-   !     (mo_class(p1)(1:1) == 'A' .or. mo_class(p1)(1:1) == 'V') 
-   !if(ok .and. degree1 == 2) then
-   !        ok = (mo_class(h2)(1:1) == 'A' .or. mo_class(h2)(1:1) == 'I') .and. &
-   !             (mo_class(p2)(1:1) == 'A' .or. mo_class(p2)(1:1) == 'V') 
-   !end if
-    call i_h_j(alpha,det_minilist(1,1,l_sd),N_int,hij_cache_(l_sd,iproc))
-    call get_s2(alpha,det_minilist(1,1,l_sd),N_int,sij_cache_(l_sd,iproc))
-  enddo
    
+   ok = (mo_class(h1)(1:1) == 'A' .or. mo_class(h1)(1:1) == 'I') .and. &
+        (mo_class(p1)(1:1) == 'A' .or. mo_class(p1)(1:1) == 'V') 
+   if(ok .and. degree1 == 2) then
+           ok = (mo_class(h2)(1:1) == 'A' .or. mo_class(h2)(1:1) == 'I') .and. &
+                (mo_class(p2)(1:1) == 'A' .or. mo_class(p2)(1:1) == 'V') 
+   end if
+   
+   if(ok) then
+     canbediamond += 1
+     excs_(:,:,:,l_sd,iproc) = exc(:,:,:)
+     phases_(l_sd, iproc) = phase
+   else
+     phases_(l_sd, iproc) = 0d0
+   end if
+   !call i_h_j(alpha,det_minilist(1,1,l_sd),N_int,hij_cache_(l_sd,iproc))
+   !call get_s2(alpha,det_minilist(1,1,l_sd),N_int,sij_cache_(l_sd,iproc))
+   call i_h_j_s2(alpha,det_minilist(1,1,l_sd),N_int,hij_cache_(l_sd,iproc), sij_cache_(l_sd,iproc))
+  enddo
+  if(canbediamond <= 1) return
 
   do i_I=1,N_det_ref
     call get_excitation_degree(alpha,psi_ref(1,1,i_I),degree1,N_int)
@@ -91,12 +97,16 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil
     enddo
 
     do k_sd=1,n_minilist
+      if(phases_(k_sd,iproc) == 0d0) cycle
       call get_excitation_degree(psi_ref(1,1,i_I),det_minilist(1,1,k_sd),degree,N_int)
       if (degree > 2) then
         cycle
       endif
       
-      call get_excitation(det_minilist(1,1,k_sd),alpha,exc,degree2,phase,N_int)
+      !call get_excitation(det_minilist(1,1,k_sd),alpha,exc,degree2,phase,N_int)
+      phase = phases_(k_sd, iproc)
+      exc(:,:,:) = excs_(:,:,:,k_sd,iproc)
+      degree2 = exc(0,1,1) + exc(0,1,2)
       call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, N_int)
       
       if((.not. ok) .and. (.not. perturbative_triples)) cycle
@@ -118,6 +128,7 @@ subroutine dress_with_alpha_buffer(delta_ij_loc, minilist, det_minilist, n_minil
       if (ok) then
         phase2 = 0d0
         do l_sd=k_sd+1,n_minilist
+          if(phases_(l_sd, iproc) == 0d0) cycle
           call get_excitation_degree(tmp_det,det_minilist(1,1,l_sd),degree,N_int)
           if (degree == 0) then
             do i_state=1,N_states
@@ -204,11 +215,11 @@ subroutine test_minilist(minilist, n_minilist, alpha)
   refc = 0
   testc = 0
   do i=1,N_det
-    call get_excitation_degree(psi_det_sorted(1,1,i), alpha, deg, N_int) 
+    call get_excitation_degree(psi_det(1,1,i), alpha, deg, N_int) 
     if(deg <= 2) refc(i) = refc(i) + 1
   end do
   do i=1,n_minilist
-    call get_excitation_degree(psi_det_sorted(1,1,minilist(i)), alpha, deg, N_int) 
+    call get_excitation_degree(psi_det(1,1,minilist(i)), alpha, deg, N_int) 
     if(deg <= 2) then
       testc(minilist(i)) += 1
     else

plugins/mrcepa0/dressing.irp.f

@@ -343,10 +343,8 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
   !double precision, external :: get_dij, get_dij_index
   double precision :: Delta_E_inv(N_states)
   double precision, intent(inout) :: contrib(N_states)
-  double precision :: sdress, hdress, shdress
+  double precision :: sdress, hdress
-  integer :: old_ninc
 
-  old_ninc = ninc
 
   if (perturbative_triples) then
     PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
@@ -411,9 +409,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
       
  
   do i_alpha=1,N_tq
-    old_ninc = ninc
-    shdress = 0d0
-    
     if(key_mask(1,1) /= 0) then
       call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint)
       
@@ -552,8 +547,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
           hdress = dIa_hla(i_state,k_sd) * psi_ref_coef(i_I,i_state) 
           sdress = dIa_sla(i_state,k_sd) * psi_ref_coef(i_I,i_state) 
           !!$OMP ATOMIC
-          if(hdress /= 0d0) ninc = ninc + 1
-          shdress += hdress
           !$OMP ATOMIC
           contrib(i_state) += hdress * psi_coef(dressed_column_idx(i_state), i_state) * psi_non_ref_coef(k_sd, i_state)
           !$OMP ATOMIC
@@ -563,16 +556,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ij_s2_, i_generator,n_selected,det_b
         enddo
       enddo
     enddo
-    if(ninc /= old_ninc) then
-      nalp = nalp + 1
-      !print "(A8,I20,I20,E15.5)", "grepme", tq(1,1,i_alpha), tq(1,2,i_alpha), shdress
-      !if(tq(1,1,i_alpha) == 1007 .and. tq(1,2,i_alpha) == 17301943) then
-      ! print *, "foinder", i_generator
-      ! call debug_det(psi_det_generators(1,1, i_generator), N_int)
-      ! call debug_det(tq(1,1,i_alpha), N_int)
-      ! stop
-     ! end if
-    end if
   enddo
 
   deallocate (dIa_hla,dIa_sla,hij_cache,sij_cache)

plugins/mrcepa0/dressing_vector.irp.f

@@ -1,9 +1,4 @@
 
- BEGIN_PROVIDER [ integer, nalp ]
-&BEGIN_PROVIDER [ integer, ninc ]
-  nalp = 0
-  ninc = 0
-END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, dressing_column_h, (N_det,N_states) ]
 &BEGIN_PROVIDER [ double precision, dressing_column_s, (N_det,N_states) ]
@@ -31,9 +26,6 @@ END_PROVIDER
    tmp = u_dot_v(dressing_column_s(1,k), psi_coef(1,k), N_det)
    dressing_column_s(l,k) -= tmp * f
  enddo
- print *, "NALP", nalp
- print *, "NINC", ninc
- print *, "DRESS", dressing_column_h(:10,1)
 ! stop
 END_PROVIDER
 

src/Determinants/slater_rules.irp.f

@@ -474,6 +474,107 @@ subroutine bitstring_to_list_ab_old( string, list, n_elements, Nint)
 end
 
 
+subroutine i_H_j_s2(key_i,key_j,Nint,hij,s2)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  ! Returns <i|H|j> where i and j are determinants
+  END_DOC
+  integer, intent(in)            :: Nint
+  integer(bit_kind), intent(in)  :: key_i(Nint,2), key_j(Nint,2)
+  double precision, intent(out)  :: hij, s2
+  
+  integer                        :: exc(0:2,2,2)
+  integer                        :: degree
+  double precision               :: get_mo_bielec_integral
+  integer                        :: m,n,p,q
+  integer                        :: i,j,k
+  integer                        :: occ(Nint*bit_kind_size,2)
+  double precision               :: diag_H_mat_elem, phase,phase_2
+  integer                        :: n_occ_ab(2)
+  PROVIDE mo_bielec_integrals_in_map mo_integrals_map big_array_exchange_integrals
+  
+  ASSERT (Nint > 0)
+  ASSERT (Nint == N_int)
+  ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num)
+  ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num)
+  ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num)
+  ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num)
+  
+  hij = 0.d0
+  s2 = 0d0
+  !DIR$ FORCEINLINE
+  call get_excitation_degree(key_i,key_j,degree,Nint)
+  integer :: spin
+  select case (degree)
+    case (2)
+      call get_double_excitation(key_i,key_j,exc,phase,Nint)
+      
+      if (exc(0,1,1) == 1) then
+        ! Mono alpha, mono beta
+        if(exc(1,1,1) == exc(1,2,2) )then
+         if(exc(1,1,2) == exc(1,2,1)) s2 = -phase !!!!!
+         hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1))
+        else if (exc(1,2,1) ==exc(1,1,2))then
+         hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2))
+        else
+         hij = phase*get_mo_bielec_integral(                          &
+             exc(1,1,1),                                              &
+             exc(1,1,2),                                              &
+             exc(1,2,1),                                              &
+             exc(1,2,2) ,mo_integrals_map)
+        endif
+      else if (exc(0,1,1) == 2) then
+        ! Double alpha
+        hij = phase*(get_mo_bielec_integral(                         &
+            exc(1,1,1),                                              &
+            exc(2,1,1),                                              &
+            exc(1,2,1),                                              &
+            exc(2,2,1) ,mo_integrals_map) -                          &
+            get_mo_bielec_integral(                                  &
+            exc(1,1,1),                                              &
+            exc(2,1,1),                                              &
+            exc(2,2,1),                                              &
+            exc(1,2,1) ,mo_integrals_map) )
+      else if (exc(0,1,2) == 2) then
+        ! Double beta
+        hij = phase*(get_mo_bielec_integral(                         &
+            exc(1,1,2),                                              &
+            exc(2,1,2),                                              &
+            exc(1,2,2),                                              &
+            exc(2,2,2) ,mo_integrals_map) -                          &
+            get_mo_bielec_integral(                                  &
+            exc(1,1,2),                                              &
+            exc(2,1,2),                                              &
+            exc(2,2,2),                                              &
+            exc(1,2,2) ,mo_integrals_map) )
+      endif
+    case (1)
+      call get_mono_excitation(key_i,key_j,exc,phase,Nint)
+      !DIR$ FORCEINLINE
+      call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
+      if (exc(0,1,1) == 1) then
+        ! Mono alpha
+        m = exc(1,1,1)
+        p = exc(1,2,1)
+        spin = 1
+      else
+        ! Mono beta
+        m = exc(1,1,2)
+        p = exc(1,2,2)
+        spin = 2
+      endif
+      call get_mono_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij)
+      
+    case (0)
+      print *," ZERO"
+      double precision, external :: diag_S_mat_elem
+      s2 = diag_S_mat_elem(key_i,Nint)
+      hij = diag_H_mat_elem(key_i,Nint)
+  end select
+end
+
+
 
 subroutine i_H_j(key_i,key_j,Nint,hij)
   use bitmasks