Accelerated MRSC2

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -652,14 +652,12 @@ END_PROVIDER
   allocate(rho_mrcc_init(N_det_non_ref))
   allocate(x_new(hh_nex))
   allocate(x(hh_nex), AtB(hh_nex))
-  x = 0d0
-        
 
   do s=1,N_states
 
     AtB(:) = 0.d0
     !$OMP PARALLEL default(none) shared(k, psi_non_ref_coef, active_excitation_to_determinants_idx,&
-        !$OMP   active_excitation_to_determinants_val, x, N_det_ref, hh_nex, N_det_non_ref)          &
+        !$OMP   active_excitation_to_determinants_val, N_det_ref, hh_nex, N_det_non_ref)          &
         !$OMP private(at_row, a_col, i, j, r1, r2, wk, A_ind_mwen, A_val_mwen, a_coll, at_roww)&
         !$OMP shared(N_states,mrcc_col_shortcut, mrcc_N_col, AtB, mrcc_AtA_val, mrcc_AtA_ind, s, n_exc_active, active_pp_idx)
     
@@ -721,14 +719,14 @@ END_PROVIDER
     factor = 1.d0
     resold = huge(1.d0)
 
-    do k=0,hh_nex*hh_nex
+    do k=0,10*hh_nex
       !$OMP PARALLEL default(shared) private(cx, i, a_col, a_coll)
       
       !$OMP DO
       do i=1,N_det_non_ref
         rho_mrcc(i,s) = rho_mrcc_init(i) 
       enddo
-      !$OMP END DO NOWAIT
+      !$OMP END DO
       
       !$OMP DO
       do a_coll = 1, n_exc_active

src/Determinants/determinants.irp.f

@@ -764,37 +764,85 @@ subroutine apply_excitation(det, exc, res, ok, Nint)
   ok = .false.
   degree = exc(0,1,1) + exc(0,1,2)
   
-  if(.not. (degree > 0 .and. degree <= 2)) then
+!  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+! INLINE
+  select case(degree)
+    case(2)
+      if (exc(0,1,1) == 2) then
+        h1 = exc(1,1,1)
+        h2 = exc(2,1,1)
+        p1 = exc(1,2,1)
+        p2 = exc(2,2,1)
+        s1 = 1
+        s2 = 1
+      else if (exc(0,1,2) == 2) then
+        h1 = exc(1,1,2)
+        h2 = exc(2,1,2)
+        p1 = exc(1,2,2)
+        p2 = exc(2,2,2)
+        s1 = 2
+        s2 = 2
+      else
+        h1 = exc(1,1,1)
+        h2 = exc(1,1,2)
+        p1 = exc(1,2,1)
+        p2 = exc(1,2,2)
+        s1 = 1
+        s2 = 2
+      endif
+    case(1)
+      if (exc(0,1,1) == 1) then
+        h1 = exc(1,1,1)
+        h2 = 0
+        p1 = exc(1,2,1)
+        p2 = 0
+        s1 = 1
+        s2 = 0
+      else
+        h1 = exc(1,1,2)
+        h2 = 0
+        p1 = exc(1,2,2)
+        p2 = 0
+        s1 = 2
+        s2 = 0
+      endif
+    case(0)
+      h1 = 0
+      p1 = 0
+      h2 = 0
+      p2 = 0
+      s1 = 0
+      s2 = 0
+    case default
       print *, degree
       print *, "apply ex"
       STOP
-  endif
+  end select
+! END INLINE
 
-  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
   res = det 
   
-  ii = (h1-1)/bit_kind_size + 1 
+  ii = ishft(h1-1,-bit_kind_shift) + 1 
-  pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64
+  pos = h1-1-ishft(ii-1,bit_kind_shift)
-  if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return
+  if(iand(det(ii, s1), ibset(0_bit_kind, pos)) == 0_8) return
   res(ii, s1) = ibclr(res(ii, s1), pos)
   
-  ii = (p1-1)/bit_kind_size + 1 
+  ii = ishft(p1-1,-bit_kind_shift) + 1 
-  pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1)
+  pos = p1-1-ishft(ii-1,bit_kind_shift)
   if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return
   res(ii, s1) = ibset(res(ii, s1), pos)
   
   if(degree == 2) then
-    ii = (h2-1)/bit_kind_size + 1 
+    ii = ishft(h2-1,-bit_kind_shift) + 1 
-    pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1)
+    pos = h2-1-ishft(ii-1,bit_kind_shift)
     if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return
     res(ii, s2) = ibclr(res(ii, s2), pos)
     
-    ii = (p2-1)/bit_kind_size + 1 
+    ii = ishft(p2-1,-bit_kind_shift) + 1 
-    pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1)
+    pos = p2-1-ishft(ii-1,bit_kind_shift)
     if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return
     res(ii, s2) = ibset(res(ii, s2), pos)
   endif
-
   ok = .true.
 end subroutine
 
@@ -814,13 +862,13 @@ subroutine apply_particles(det, s1, p1, s2, p2, res, ok, Nint)
   
   if(p1 /= 0) then
   ii = (p1-1)/bit_kind_size + 1 
-  pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1)
+  pos = mod(p1-1, bit_kind_size)!iand(p1-1,bit_kind_size-1)
   if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return
   res(ii, s1) = ibset(res(ii, s1), pos)
   end if
 
   ii = (p2-1)/bit_kind_size + 1 
-  pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1)
+  pos = mod(p2-1, bit_kind_size)!iand(p2-1,bit_kind_size-1)
   if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return
   res(ii, s2) = ibset(res(ii, s2), pos)
 
@@ -843,13 +891,13 @@ subroutine apply_holes(det, s1, h1, s2, h2, res, ok, Nint)
   
   if(h1 /= 0) then
   ii = (h1-1)/bit_kind_size + 1 
-  pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1)
+  pos = mod(h1-1, bit_kind_size)!iand(h1-1,bit_kind_size-1)
   if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return
   res(ii, s1) = ibclr(res(ii, s1), pos)
   end if
 
   ii = (h2-1)/bit_kind_size + 1 
-  pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1)
+  pos = mod(h2-1, bit_kind_size)!iand(h2-1,bit_kind_size-1)
   if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return
   res(ii, s2) = ibclr(res(ii, s2), pos)
 
@@ -870,7 +918,7 @@ subroutine apply_particle(det, s1, p1, res, ok, Nint)
   res = det 
   
   ii = (p1-1)/bit_kind_size + 1 
-  pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1)
+  pos = mod(p1-1, bit_kind_size)!iand(p1-1,bit_kind_size-1)
   if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return
   res(ii, s1) = ibset(res(ii, s1), pos)
 
@@ -892,7 +940,7 @@ subroutine apply_hole(det, s1, h1, res, ok, Nint)
   res = det 
   
   ii = (h1-1)/bit_kind_size + 1 
-  pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1)
+  pos = mod(h1-1, bit_kind_size)!iand(h1-1,bit_kind_size-1)
   if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return
   res(ii, s1) = ibclr(res(ii, s1), pos)