cleaner k-point orbital indexing

src/ao_basis/aos_cplx.irp.f

@@ -5,3 +5,19 @@
 ! END_DOC
 ! ao_num_per_kpt = ao_num/kpt_num
 !END_PROVIDER
+
+subroutine get_kpt_idx_ao(idx_full,k,i)
+  implicit none
+  BEGIN_DOC
+  ! idx_full is ao index in full range (up to ao_num)
+  ! k is index of the k-point for this ao
+  ! i is index of this ao within k-point k
+  ! this assumes that all kpts have the same number of aos
+  END_DOC
+
+  integer, intent(in) :: idx_full
+  integer, intent(out) :: i,k
+  i = mod(idx_full-1,ao_num_per_kpt)+1
+  k = (idx_full-1)/ao_num_per_kpt+1
+  ASSERT (k <= kpt_num)
+end

src/cipsi/selection.irp.f

@@ -2918,14 +2918,10 @@ subroutine get_d1_kpts(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, c
     hfix = h(1,ma)
     p1 = p(1,ma)
     p2 = p(2,ma)
-    kputi = (puti-1)/mo_num_per_kpt + 1
+    call get_kpt_idx_mo(puti,kputi,iputi)
-    khfix = (hfix-1)/mo_num_per_kpt + 1
+    call get_kpt_idx_mo(hfix,khfix,ihfix)
-    kp1   =   (p1-1)/mo_num_per_kpt + 1
+    call get_kpt_idx_mo(p1,kp1,ip1)
-    kp2   =   (p2-1)/mo_num_per_kpt + 1
+    call get_kpt_idx_mo(p2,kp2,ip2)
-    iputi = mod(puti-1,mo_num_per_kpt) + 1
-    ihfix = mod(hfix-1,mo_num_per_kpt) + 1
-    ip1   = mod(p1-1,  mo_num_per_kpt) + 1
-    ip2   = mod(p2-1,  mo_num_per_kpt) + 1
 
     if(.not. bannedOrb(puti, mi)) then
       !==================
@@ -3059,8 +3055,7 @@ subroutine get_d1_kpts(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, c
 
     !MOVE MI
     pfix = p(1,mi)
-    kpfix = (pfix-1)/mo_num_per_kpt + 1
+    call get_kpt_idx_mo(pfix,kpfix,ipfix)
-    ipfix = mod(pfix-1,mo_num_per_kpt) + 1
     tmp_row = (0.d0,0.d0)
     tmp_row2 = (0.d0,0.d0)
     !tmp_row_kpts = (0.d0,0.d0)
@@ -3270,14 +3265,10 @@ subroutine get_d1_kpts(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, c
         puti = p(i, ma)
         p1 = p(turn3(1,i), ma)
         p2 = p(turn3(2,i), ma)
-        kputi = (puti-1)/mo_num_per_kpt + 1
+        call get_kpt_idx_mo(puti,kputi,iputi)
-        khfix = (hfix-1)/mo_num_per_kpt + 1
+        call get_kpt_idx_mo(hfix,khfix,ihfix)
-        kp1   =   (p1-1)/mo_num_per_kpt + 1
+        call get_kpt_idx_mo(p1,kp1,ip1)
-        kp2   =   (p2-1)/mo_num_per_kpt + 1
+        call get_kpt_idx_mo(p2,kp2,ip2)
-        iputi = mod(puti-1,mo_num_per_kpt) + 1
-        ihfix = mod(hfix-1,mo_num_per_kpt) + 1
-        ip1   = mod(p1-1,  mo_num_per_kpt) + 1
-        ip2   = mod(p2-1,  mo_num_per_kpt) + 1
         call get_mo_two_e_integrals_complex(hfix,p1,p2,mo_num,hij_cache(1,1),mo_integrals_map,mo_integrals_map_2)
         call get_mo_two_e_integrals_complex(hfix,p2,p1,mo_num,hij_cache(1,2),mo_integrals_map,mo_integrals_map_2)
         call get_mo_two_e_integrals_kpts(hfix,ihfix,khfix,p1,ip1,kp1,p2,ip2,kp2,mo_num_per_kpt,hij_cache2(1,1),mo_integrals_map,mo_integrals_map_2)
@@ -3425,14 +3416,10 @@ subroutine get_d1_kpts(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, c
       pfix = p(1,mi)
       p1 = p(1,ma)
       p2 = p(2,ma)
-      kpfix = (pfix-1)/mo_num_per_kpt + 1
+      call get_kpt_idx_mo(pfix,kpfix,ipfix)
-      khfix = (hfix-1)/mo_num_per_kpt + 1
+      call get_kpt_idx_mo(hfix,khfix,ihfix)
-      kp1   =   (p1-1)/mo_num_per_kpt + 1
+      call get_kpt_idx_mo(p1,kp1,ip1)
-      kp2   =   (p2-1)/mo_num_per_kpt + 1
+      call get_kpt_idx_mo(p2,kp2,ip2)
-      ipfix = mod(pfix-1,mo_num_per_kpt) + 1
-      ihfix = mod(hfix-1,mo_num_per_kpt) + 1
-      ip1   = mod(p1-1,  mo_num_per_kpt) + 1
-      ip2   = mod(p2-1,  mo_num_per_kpt) + 1
       tmp_row = (0.d0,0.d0)
       tmp_row2 = (0.d0,0.d0)
       !tmp_row_kpts = (0.d0,0.d0)

src/determinants/density_matrix_cplx.irp.f

@@ -495,10 +495,8 @@ END_PROVIDER
         call decode_exc_spin(exc,h1,p1,h2,p2)
         ! h1 occ in k
         ! p1 occ in l
-        ih1 = mod(h1-1,mo_num_per_kpt)+1
+        call get_kpt_idx_mo(h1,kh1,ih1)
-        ip1 = mod(p1-1,mo_num_per_kpt)+1
+        call get_kpt_idx_mo(p1,kp1,ip1)
-        kh1 = (h1-1)/mo_num_per_kpt + 1
-        kp1 = (p1-1)/mo_num_per_kpt + 1
         if (kh1.ne.kp1) then
           print *,'problem in: ',irp_here,'a'
           print *,' h1 = ',h1
@@ -577,10 +575,8 @@ END_PROVIDER
         exc = 0
         call get_single_excitation_spin(tmp_det(1,2),tmp_det2,exc,phase,N_int)
         call decode_exc_spin(exc,h1,p1,h2,p2)
-        ih1 = mod(h1-1,mo_num_per_kpt)+1
+        call get_kpt_idx_mo(h1,kh1,ih1)
-        ip1 = mod(p1-1,mo_num_per_kpt)+1
+        call get_kpt_idx_mo(p1,kp1,ip1)
-        kh1 = (h1-1)/mo_num_per_kpt + 1
-        kp1 = (p1-1)/mo_num_per_kpt + 1
         if (kh1.ne.kp1) then
           print *,'problem in: ',irp_here,'b'
           print *,' h1 = ',h1

src/determinants/single_excitations.irp.f

@@ -449,11 +449,12 @@ subroutine get_single_excitation_from_fock_kpts(det_1,det_2,ih,ip,spin,phase,hij
   integer :: occ_partcl(N_int*bit_kind_size,2)
   integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2)
   integer :: i0,i,h,p
-  integer :: ki,khp
+  integer :: ki,khp,kh
   complex*16 :: buffer_c(mo_num_per_kpt),buffer_x(mo_num_per_kpt)
-  khp = (ip-1)/mo_num_per_kpt+1
+
-  p = mod(ip-1,mo_num_per_kpt)+1
+  call get_kpt_idx_mo(ip,khp,p)
-  h = mod(ih-1,mo_num_per_kpt)+1
+  call get_kpt_idx_mo(ih,kh,h)
+  ASSERT (kh==khp)
   !todo: omp kpts
   do ki=1,kpt_num
     do i=1, mo_num_per_kpt

src/determinants/slater_rules.irp.f

@@ -2443,18 +2443,19 @@ subroutine i_H_j_complex(key_i,key_j,Nint,hij)
       if (exc(0,1,1) == 1) then
         call double_allowed_mo_kpts(exc(1,1,1),exc(1,1,2),exc(1,2,1),exc(1,2,2),is_allowed)
         if (.not.is_allowed) then
+          ! excitation doesn't conserve momentum
           hij = (0.d0,0.d0)
           return
         endif
         ! Single alpha, single beta
         if(exc(1,1,1) == exc(1,2,2) )then
-          ih1 = mod(exc(1,1,1)-1,mo_num_per_kpt)+1
+          !h1(a) = p2(b)
-          ih2 = mod(exc(1,1,2)-1,mo_num_per_kpt)+1
+          call get_kpt_idx_mo(exc(1,1,1),kh1,ih1)
-          kh1 = (exc(1,1,1)-1)/mo_num_per_kpt+1
+          call get_kpt_idx_mo(exc(1,1,2),kh2,ih2)
-          kh2 = (exc(1,1,2)-1)/mo_num_per_kpt+1
+          call get_kpt_idx_mo(exc(1,2,1),kp1,ip1)
-          ip1 = mod(exc(1,2,1)-1,mo_num_per_kpt)+1
+
-          kp1 = (exc(1,2,1)-1)/mo_num_per_kpt+1
           if(kp1.ne.kh2) then
+            !if h1==p2 then kp1==kh2
             print*,'problem with hij kpts: ',irp_here
             print*,is_allowed
             print*,exc(1,1,1),exc(1,1,2),exc(1,2,1),exc(1,2,2)
@@ -2464,12 +2465,10 @@ subroutine i_H_j_complex(key_i,key_j,Nint,hij)
           hij = phase * big_array_exchange_integrals_kpts(ih1,kh1,ih2,ip1,kp1)
           !hij = phase * big_array_exchange_integrals_complex(exc(1,1,1),exc(1,1,2),exc(1,2,1))
         else if (exc(1,2,1) ==exc(1,1,2))then
-          ih1 = mod(exc(1,1,1)-1,mo_num_per_kpt)+1
+          !p1(a)==h2(b)
-          kh1 = (exc(1,1,1)-1)/mo_num_per_kpt+1
+          call get_kpt_idx_mo(exc(1,1,1),kh1,ih1)
-          ip1 = mod(exc(1,2,1)-1,mo_num_per_kpt)+1
+          call get_kpt_idx_mo(exc(1,2,1),kp1,ip1)
-          kp1 = (exc(1,2,1)-1)/mo_num_per_kpt+1
+          call get_kpt_idx_mo(exc(1,2,2),kp2,ip2)
-          ip2 = mod(exc(1,2,2)-1,mo_num_per_kpt)+1
-          kp2 = (exc(1,2,2)-1)/mo_num_per_kpt+1
           if(kp2.ne.kh1) then
             print*,'problem with hij kpts: ',irp_here
             print*,is_allowed

src/mo_basis/mos.irp.f

@@ -80,6 +80,22 @@ BEGIN_PROVIDER [ integer, mo_num_per_kpt ]
 
 END_PROVIDER
 
+subroutine get_kpt_idx_mo(idx_full,k,i)
+  implicit none
+  BEGIN_DOC
+  ! idx_full is mo index in full range (up to mo_num)
+  ! k is index of the k-point for this mo
+  ! i is index of this mo within k-point k
+  ! this assumes that all kpts have the same number of mos
+  END_DOC
+
+  integer, intent(in) :: idx_full
+  integer, intent(out) :: i,k
+  i = mod(idx_full-1,mo_num_per_kpt)+1
+  k = (idx_full-1)/mo_num_per_kpt+1
+  ASSERT (k <= kpt_num)
+end
+
 
 BEGIN_PROVIDER [ double precision, mo_coef, (ao_num,mo_num) ]
   implicit none

src/scf_utils/fock_matrix_cplx.irp.f

@@ -593,14 +593,10 @@ END_PROVIDER
           j = jj(k2)
           k = kk(k2)
           l = ll(k2)
-          kpt_i = (i-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(i,kpt_i,idx_i)
-          kpt_j = (j-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(j,kpt_j,idx_j)
-          kpt_k = (k-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(k,kpt_k,idx_k)
-          kpt_l = (l-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(l,kpt_l,idx_l)
-          idx_i = mod(i-1,ao_num_per_kpt)+1
-          idx_j = mod(j-1,ao_num_per_kpt)+1
-          idx_k = mod(k-1,ao_num_per_kpt)+1
-          idx_l = mod(l-1,ao_num_per_kpt)+1
           integral = i_sign(k2)*values(k1) !for klij and lkji, take complex conjugate
 
           !G_a(i,k) += D_{ab}(l,j)*(<ij|kl>)
@@ -636,14 +632,10 @@ END_PROVIDER
           j = jj(k2)
           k = kk(k2)
           l = ll(k2)
-          kpt_i = (i-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(i,kpt_i,idx_i)
-          kpt_j = (j-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(j,kpt_j,idx_j)
-          kpt_k = (k-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(k,kpt_k,idx_k)
-          kpt_l = (l-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(l,kpt_l,idx_l)
-          idx_i = mod(i-1,ao_num_per_kpt)+1
-          idx_j = mod(j-1,ao_num_per_kpt)+1
-          idx_k = mod(k-1,ao_num_per_kpt)+1
-          idx_l = mod(l-1,ao_num_per_kpt)+1
           integral = values(k1)
 
           if (kpt_l.eq.kpt_j) then
@@ -714,14 +706,10 @@ END_PROVIDER
           j = jj(k2)
           k = kk(k2)
           l = ll(k2)
-          kpt_i = (i-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(i,kpt_i,idx_i)
-          kpt_j = (j-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(j,kpt_j,idx_j)
-          kpt_k = (k-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(k,kpt_k,idx_k)
-          kpt_l = (l-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(l,kpt_l,idx_l)
-          idx_i = mod(i-1,ao_num_per_kpt)+1
-          idx_j = mod(j-1,ao_num_per_kpt)+1
-          idx_k = mod(k-1,ao_num_per_kpt)+1
-          idx_l = mod(l-1,ao_num_per_kpt)+1
           integral = i_sign(k2)*values(k1) ! for klij and lkji, take conjugate
 
           !G_a(i,k) += D_{ab}(l,j)*(<ij|kl>)
@@ -757,14 +745,10 @@ END_PROVIDER
           j = jj(k2)
           k = kk(k2)
           l = ll(k2)
-          kpt_i = (i-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(i,kpt_i,idx_i)
-          kpt_j = (j-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(j,kpt_j,idx_j)
-          kpt_k = (k-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(k,kpt_k,idx_k)
-          kpt_l = (l-1)/ao_num_per_kpt +1
+          call get_kpt_idx_ao(l,kpt_l,idx_l)
-          idx_i = mod(i-1,ao_num_per_kpt)+1
-          idx_j = mod(j-1,ao_num_per_kpt)+1
-          idx_k = mod(k-1,ao_num_per_kpt)+1
-          idx_l = mod(l-1,ao_num_per_kpt)+1
           integral = values(k1)
 
           if (kpt_l.eq.kpt_j) then