fixed kpts cipsi

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

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
 

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)
 

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

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)

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

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

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
+

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

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