From 9fa523fe6621c7a789b159c67ca058f3826e6c4c Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Thu, 26 Mar 2020 11:30:15 -0500 Subject: [PATCH] fixed kpts cipsi --- src/bitmask/bitmasks_routines.irp.f | 31 +++++++++++ src/bitmask/core_inact_act_virt.irp.f | 4 ++ src/cipsi/stochastic_cipsi.irp.f | 1 + src/determinants/density_matrix_cplx.irp.f | 11 ++-- src/determinants/slater_rules.irp.f | 54 ++++++++++++++++++-- src/iterations/print_summary.irp.f | 12 +++++ src/mo_two_e_ints/mo_bi_integrals_cplx.irp.f | 12 +++++ src/nuclei/kconserv_cplx.irp.f | 9 ++++ src/utils_complex/dump_mo_2e_cplx.irp.f | 4 +- src/utils_complex/qp2-pbc-diff.txt | 9 ++++ 10 files changed, 139 insertions(+), 8 deletions(-) diff --git a/src/bitmask/bitmasks_routines.irp.f b/src/bitmask/bitmasks_routines.irp.f index 5c4bf347..8a374e94 100644 --- a/src/bitmask/bitmasks_routines.irp.f +++ b/src/bitmask/bitmasks_routines.irp.f @@ -214,6 +214,37 @@ subroutine print_spindet(string,Nint) end +subroutine debug_single_spindet(string,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Subroutine to print the content of a determinant in '+-' notation and + ! hexadecimal representation. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: string(Nint) + character*(2048) :: output(1) + call bitstring_to_hexa( output(1), string(1), Nint ) + print *, trim(output(1)) + call print_single_spindet(string,Nint) + +end + +subroutine print_single_spindet(string,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Subroutine to print the content of a determinant using the '+-' notation + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: string(Nint) + character*(2048) :: output(1) + + call bitstring_to_str( output(1), string(1), Nint ) + print *, trim(output(1)) + +end + logical function is_integer_in_string(bite,string,Nint) use bitmasks implicit none diff --git a/src/bitmask/core_inact_act_virt.irp.f b/src/bitmask/core_inact_act_virt.irp.f index 337e275b..5fea3418 100644 --- a/src/bitmask/core_inact_act_virt.irp.f +++ b/src/bitmask/core_inact_act_virt.irp.f @@ -426,12 +426,16 @@ BEGIN_PROVIDER [ integer(bit_kind), kpts_bitmask , (N_int,kpt_num) ] integer :: k,i,di integer :: tmp_mo_list(mo_num_per_kpt) kpts_bitmask = 0_bit_kind + print*,'kpts bitmask' do k=1,kpt_num di=(k-1)*mo_num_per_kpt do i=1,mo_num_per_kpt tmp_mo_list(i) = i+di enddo call list_to_bitstring( kpts_bitmask(1,k), tmp_mo_list, mo_num_per_kpt, N_int) + !debugging + print*,'k' + call debug_single_spindet(kpts_bitmask(1,k),N_int) enddo END_PROVIDER diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 7a07577a..7c16cd41 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -110,6 +110,7 @@ subroutine run_stochastic_cipsi call write_double(6,correlation_energy_ratio, 'Correlation ratio') call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2) + !call print_debug_fci() call save_energy(psi_energy_with_nucl_rep, rpt2) diff --git a/src/determinants/density_matrix_cplx.irp.f b/src/determinants/density_matrix_cplx.irp.f index c1b63cf2..d7281e76 100644 --- a/src/determinants/density_matrix_cplx.irp.f +++ b/src/determinants/density_matrix_cplx.irp.f @@ -508,9 +508,11 @@ END_PROVIDER print *,'kh1 = ',kh1 print *,'kp1 = ',kp1 !call debug_det(tmp_det,N_int) - !call debug_spindet(tmp_det2,N_int) + call debug_single_spindet(tmp_det(1,1),N_int) + call debug_single_spindet(tmp_det2,N_int) + call debug_single_spindet(tmp_det(1,2),N_int) !call print_spindet(tmp_det2,N_int) - !stop -2 + stop -2 endif do m=1,N_states ckl = dconjg(psi_bilinear_matrix_values_complex(k_a,m))*psi_bilinear_matrix_values_complex(l,m) * phase @@ -587,7 +589,10 @@ END_PROVIDER print *,'ip1 = ',ip1 print *,'kh1 = ',kh1 print *,'kp1 = ',kp1 - !stop -3 + call debug_single_spindet(tmp_det(1,2),N_int) + call debug_single_spindet(tmp_det2,N_int) + call debug_single_spindet(tmp_det(1,1),N_int) + stop -3 endif do m=1,N_states ckl = dconjg(psi_bilinear_matrix_transp_values_complex(k_b,m))*psi_bilinear_matrix_transp_values_complex(l,m) * phase diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index ddd469b1..a4fc267a 100644 --- a/src/determinants/slater_rules.irp.f +++ b/src/determinants/slater_rules.irp.f @@ -2418,9 +2418,11 @@ subroutine i_H_j_complex(key_i,key_j,Nint,hij) complex*16 :: get_two_e_integral_complex integer :: m,n,p,q integer :: i,j,k + integer :: ih1,ih2,ip1,ip2,kh1,kh2,kp1,kp2 integer :: occ(Nint*bit_kind_size,2) double precision :: diag_H_mat_elem, phase integer :: n_occ_ab(2) + logical :: is_allowed PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals_complex ASSERT (Nint > 0) @@ -2439,11 +2441,38 @@ subroutine i_H_j_complex(key_i,key_j,Nint,hij) case (2) call get_double_excitation(key_i,key_j,exc,phase,Nint) 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 + hij = (0.d0,0.d0) + return + endif ! Single alpha, single beta if(exc(1,1,1) == exc(1,2,2) )then - hij = phase * big_array_exchange_integrals_complex(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + 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 + if(kp1.ne.kh2) then + print*,'problem with hij kpts: ',irp_here + stop -4 + endif + 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 - hij = phase * big_array_exchange_integrals_complex(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + 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 + if(kp2.ne.kh1) then + print*,'problem with hij kpts: ',irp_here + stop -4 + endif + hij = phase * big_array_exchange_integrals_kpts(ip1,kp1,ih1,ip2,kp2) + !hij = phase * big_array_exchange_integrals_complex(exc(1,2,1),exc(1,1,1),exc(1,2,2)) else hij = phase*get_two_e_integral_complex( & exc(1,1,1), & @@ -2452,6 +2481,11 @@ subroutine i_H_j_complex(key_i,key_j,Nint,hij) exc(1,2,2) ,mo_integrals_map,mo_integrals_map_2) endif else if (exc(0,1,1) == 2) then + call double_allowed_mo_kpts(exc(1,1,1),exc(2,1,1),exc(1,2,1),exc(2,2,1),is_allowed) + if (.not.is_allowed) then + hij = (0.d0,0.d0) + return + endif ! Double alpha hij = phase*(get_two_e_integral_complex( & exc(1,1,1), & @@ -2464,6 +2498,11 @@ subroutine i_H_j_complex(key_i,key_j,Nint,hij) exc(2,2,1), & exc(1,2,1) ,mo_integrals_map,mo_integrals_map_2) ) else if (exc(0,1,2) == 2) then + call double_allowed_mo_kpts(exc(1,1,2),exc(2,1,2),exc(1,2,2),exc(2,2,2),is_allowed) + if (.not.is_allowed) then + hij = (0.d0,0.d0) + return + endif ! Double beta hij = phase*(get_two_e_integral_complex( & exc(1,1,2), & @@ -2491,6 +2530,11 @@ subroutine i_H_j_complex(key_i,key_j,Nint,hij) p = exc(1,2,2) spin = 2 endif + !if m,p not from same kpt, single not allowed + if (int((m-1)/mo_num_per_kpt + 1).ne.int((p-1)/mo_num_per_kpt + 1)) then + hij = (0.d0,0.d0) + return + endif !call get_single_excitation_from_fock_complex(key_i,key_j,m,p,spin,phase,hij) call get_single_excitation_from_fock_kpts(key_i,key_j,m,p,spin,phase,hij) @@ -2775,10 +2819,12 @@ subroutine i_H_j_single_spin_complex(key_i,key_j,Nint,spin,hij) integer :: exc(0:2,2) double precision :: phase - PROVIDE big_array_exchange_integrals_complex mo_two_e_integrals_in_map + !PROVIDE big_array_exchange_integrals_complex mo_two_e_integrals_in_map + PROVIDE big_array_exchange_integrals_kpts mo_two_e_integrals_in_map call get_single_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint) - call get_single_excitation_from_fock_complex(key_i,key_j,exc(1,1),exc(1,2),spin,phase,hij) + !call get_single_excitation_from_fock_complex(key_i,key_j,exc(1,1),exc(1,2),spin,phase,hij) + call get_single_excitation_from_fock_kpts(key_i,key_j,exc(1,1),exc(1,2),spin,phase,hij) end subroutine i_H_j_double_spin_complex(key_i,key_j,Nint,hij) diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index ad87bc8e..87bc28fa 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -102,3 +102,15 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_ end subroutine +subroutine print_debug_fci + implicit none + integer :: i + do i=1,n_det + print'(2((F25.15),2X))',psi_coef_complex(i,1) + call debug_det(psi_det(1,1,i),n_int) + enddo + print*,'hamiltonian' + do i=1,n_det + print '(1000(F25.15))',h_matrix_all_dets_complex(i,:) + enddo +end subroutine diff --git a/src/mo_two_e_ints/mo_bi_integrals_cplx.irp.f b/src/mo_two_e_ints/mo_bi_integrals_cplx.irp.f index 632ff591..8b96c498 100644 --- a/src/mo_two_e_ints/mo_bi_integrals_cplx.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals_cplx.irp.f @@ -1,4 +1,16 @@ +subroutine double_allowed_mo_kpts(h1,h2,p1,p2,is_allowed) + implicit none + integer, intent(in) :: h1,h2,p1,p2 + logical, intent(out) :: is_allowed + integer :: kh1,kh2,kp1,kp2 + + kh1 = (h1-1)/mo_num_per_kpt+1 + kh2 = (h2-1)/mo_num_per_kpt+1 + kp1 = (p1-1)/mo_num_per_kpt+1 + kp2 = (p2-1)/mo_num_per_kpt+1 + call double_allowed_kpts(kh1,kh2,kp1,kp2,is_allowed) +end subroutine subroutine add_integrals_to_map_complex(mask_ijkl) use map_module diff --git a/src/nuclei/kconserv_cplx.irp.f b/src/nuclei/kconserv_cplx.irp.f index 8978ed9b..616ba779 100644 --- a/src/nuclei/kconserv_cplx.irp.f +++ b/src/nuclei/kconserv_cplx.irp.f @@ -29,3 +29,12 @@ BEGIN_PROVIDER [integer, kconserv, (kpt_num,kpt_num,kpt_num)] print *, 'kconserv written to disk' endif END_PROVIDER + +subroutine double_allowed_kpts(kh1,kh2,kp1,kp2,is_allowed) + implicit none + integer, intent(in) :: kh1,kh2,kp1,kp2 + logical, intent(out) :: is_allowed + + is_allowed = (kconserv(kh1,kh2,kp1) == kp2) +end subroutine + diff --git a/src/utils_complex/dump_mo_2e_cplx.irp.f b/src/utils_complex/dump_mo_2e_cplx.irp.f index 80dba969..44bb706d 100644 --- a/src/utils_complex/dump_mo_2e_cplx.irp.f +++ b/src/utils_complex/dump_mo_2e_cplx.irp.f @@ -15,7 +15,9 @@ subroutine run do k=1,mo_num do l=1,mo_num tmp_cmplx = get_two_e_integral_complex(i,j,k,l,mo_integrals_map,mo_integrals_map_2) - print'(4(I4),2(E23.15))',i,j,k,l,tmp_cmplx + if (cdabs(tmp_cmplx).gt. 1d-12) then + print'(4(I4),2(E23.15))',i,j,k,l,tmp_cmplx + endif enddo enddo enddo diff --git a/src/utils_complex/qp2-pbc-diff.txt b/src/utils_complex/qp2-pbc-diff.txt index e7345240..8148248c 100644 --- a/src/utils_complex/qp2-pbc-diff.txt +++ b/src/utils_complex/qp2-pbc-diff.txt @@ -1,4 +1,13 @@ +kpts: + changed matrices to block diagonal (1-e ints, fock, 1rdm) + double_allowed_mo_kpts(h1,h2,p1,p2,is_allowed) + {h,p}{1,2} indices in total mo_num (not per kpt) + double_allowed_kpts(kh1,kh2,kp1,kp2,is_allowed) + k{h,p}{1,2} k-point indices + + only allow momentum-conserving excitations + todo: change everything to be blocked by kpt