From 75dbda613aac542465550a4da4678e3c8d4cdae0 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Mon, 13 Jul 2020 10:52:09 -0500 Subject: [PATCH] cleaner k-point orbital indexing --- src/ao_basis/aos_cplx.irp.f | 16 ++++++++ src/cipsi/selection.irp.f | 39 ++++++------------ src/determinants/density_matrix_cplx.irp.f | 12 ++---- src/determinants/single_excitations.irp.f | 9 ++-- src/determinants/slater_rules.irp.f | 23 +++++------ src/mo_basis/mos.irp.f | 16 ++++++++ src/scf_utils/fock_matrix_cplx.irp.f | 48 ++++++++-------------- 7 files changed, 81 insertions(+), 82 deletions(-) diff --git a/src/ao_basis/aos_cplx.irp.f b/src/ao_basis/aos_cplx.irp.f index f571b28d..da1adb94 100644 --- a/src/ao_basis/aos_cplx.irp.f +++ b/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 diff --git a/src/cipsi/selection.irp.f b/src/cipsi/selection.irp.f index a3703a62..8bbf41c7 100644 --- a/src/cipsi/selection.irp.f +++ b/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 - khfix = (hfix-1)/mo_num_per_kpt + 1 - kp1 = (p1-1)/mo_num_per_kpt + 1 - kp2 = (p2-1)/mo_num_per_kpt + 1 - 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_kpt_idx_mo(puti,kputi,iputi) + call get_kpt_idx_mo(hfix,khfix,ihfix) + call get_kpt_idx_mo(p1,kp1,ip1) + call get_kpt_idx_mo(p2,kp2,ip2) 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 - ipfix = mod(pfix-1,mo_num_per_kpt) + 1 + call get_kpt_idx_mo(pfix,kpfix,ipfix) 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 - khfix = (hfix-1)/mo_num_per_kpt + 1 - kp1 = (p1-1)/mo_num_per_kpt + 1 - kp2 = (p2-1)/mo_num_per_kpt + 1 - 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_kpt_idx_mo(puti,kputi,iputi) + call get_kpt_idx_mo(hfix,khfix,ihfix) + call get_kpt_idx_mo(p1,kp1,ip1) + call get_kpt_idx_mo(p2,kp2,ip2) 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 - khfix = (hfix-1)/mo_num_per_kpt + 1 - kp1 = (p1-1)/mo_num_per_kpt + 1 - kp2 = (p2-1)/mo_num_per_kpt + 1 - 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 + call get_kpt_idx_mo(pfix,kpfix,ipfix) + call get_kpt_idx_mo(hfix,khfix,ihfix) + call get_kpt_idx_mo(p1,kp1,ip1) + call get_kpt_idx_mo(p2,kp2,ip2) tmp_row = (0.d0,0.d0) tmp_row2 = (0.d0,0.d0) !tmp_row_kpts = (0.d0,0.d0) diff --git a/src/determinants/density_matrix_cplx.irp.f b/src/determinants/density_matrix_cplx.irp.f index 882b73ee..e5d74347 100644 --- a/src/determinants/density_matrix_cplx.irp.f +++ b/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 - ip1 = mod(p1-1,mo_num_per_kpt)+1 - kh1 = (h1-1)/mo_num_per_kpt + 1 - kp1 = (p1-1)/mo_num_per_kpt + 1 + call get_kpt_idx_mo(h1,kh1,ih1) + call get_kpt_idx_mo(p1,kp1,ip1) 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 - ip1 = mod(p1-1,mo_num_per_kpt)+1 - kh1 = (h1-1)/mo_num_per_kpt + 1 - kp1 = (p1-1)/mo_num_per_kpt + 1 + call get_kpt_idx_mo(h1,kh1,ih1) + call get_kpt_idx_mo(p1,kp1,ip1) if (kh1.ne.kp1) then print *,'problem in: ',irp_here,'b' print *,' h1 = ',h1 diff --git a/src/determinants/single_excitations.irp.f b/src/determinants/single_excitations.irp.f index 044c7d06..192681b3 100644 --- a/src/determinants/single_excitations.irp.f +++ b/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 - h = mod(ih-1,mo_num_per_kpt)+1 + + call get_kpt_idx_mo(ip,khp,p) + 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 diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index e9b9aca9..fb77fb45 100644 --- a/src/determinants/slater_rules.irp.f +++ b/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 - ih2 = mod(exc(1,1,2)-1,mo_num_per_kpt)+1 - kh1 = (exc(1,1,1)-1)/mo_num_per_kpt+1 - kh2 = (exc(1,1,2)-1)/mo_num_per_kpt+1 - ip1 = mod(exc(1,2,1)-1,mo_num_per_kpt)+1 - kp1 = (exc(1,2,1)-1)/mo_num_per_kpt+1 + !h1(a) = p2(b) + call get_kpt_idx_mo(exc(1,1,1),kh1,ih1) + call get_kpt_idx_mo(exc(1,1,2),kh2,ih2) + call get_kpt_idx_mo(exc(1,2,1),kp1,ip1) + 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 - kh1 = (exc(1,1,1)-1)/mo_num_per_kpt+1 - ip1 = mod(exc(1,2,1)-1,mo_num_per_kpt)+1 - kp1 = (exc(1,2,1)-1)/mo_num_per_kpt+1 - ip2 = mod(exc(1,2,2)-1,mo_num_per_kpt)+1 - kp2 = (exc(1,2,2)-1)/mo_num_per_kpt+1 + !p1(a)==h2(b) + call get_kpt_idx_mo(exc(1,1,1),kh1,ih1) + call get_kpt_idx_mo(exc(1,2,1),kp1,ip1) + call get_kpt_idx_mo(exc(1,2,2),kp2,ip2) if(kp2.ne.kh1) then print*,'problem with hij kpts: ',irp_here print*,is_allowed diff --git a/src/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f index 440d1703..f5310696 100644 --- a/src/mo_basis/mos.irp.f +++ b/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 diff --git a/src/scf_utils/fock_matrix_cplx.irp.f b/src/scf_utils/fock_matrix_cplx.irp.f index cc0dc4af..b59465f9 100644 --- a/src/scf_utils/fock_matrix_cplx.irp.f +++ b/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 - kpt_j = (j-1)/ao_num_per_kpt +1 - kpt_k = (k-1)/ao_num_per_kpt +1 - kpt_l = (l-1)/ao_num_per_kpt +1 - 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 + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) integral = i_sign(k2)*values(k1) !for klij and lkji, take complex conjugate !G_a(i,k) += D_{ab}(l,j)*() @@ -636,14 +632,10 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) - kpt_i = (i-1)/ao_num_per_kpt +1 - kpt_j = (j-1)/ao_num_per_kpt +1 - kpt_k = (k-1)/ao_num_per_kpt +1 - kpt_l = (l-1)/ao_num_per_kpt +1 - 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 + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) 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 - kpt_j = (j-1)/ao_num_per_kpt +1 - kpt_k = (k-1)/ao_num_per_kpt +1 - kpt_l = (l-1)/ao_num_per_kpt +1 - 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 + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) integral = i_sign(k2)*values(k1) ! for klij and lkji, take conjugate !G_a(i,k) += D_{ab}(l,j)*() @@ -757,14 +745,10 @@ END_PROVIDER j = jj(k2) k = kk(k2) l = ll(k2) - kpt_i = (i-1)/ao_num_per_kpt +1 - kpt_j = (j-1)/ao_num_per_kpt +1 - kpt_k = (k-1)/ao_num_per_kpt +1 - kpt_l = (l-1)/ao_num_per_kpt +1 - 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 + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) integral = values(k1) if (kpt_l.eq.kpt_j) then