updated NEEDED_MODULES in src repository

scripts/generate_h_apply.py

@@ -23,6 +23,7 @@ deinit_thread
 skip
 init_main
 filter_integrals
+filter2h2p
 """.split()
 
 class H_apply(object):
@@ -115,6 +116,13 @@ class H_apply(object):
       buffer = buffer.replace('$'+key, value)
     return buffer
 
+  def set_filter_2h_2p(self):
+    self["filter2h2p"] = """
+!    ! DIR$ FORCEINLINE
+     if(is_a_two_holes_two_particles(key))cycle
+    """
+
+
   def set_perturbation(self,pert):
     if self.perturbation is not None:
         raise

src/BiInts/map_integrals.irp.f

@@ -276,6 +276,7 @@ double precision function get_mo_bielec_integral(i,j,k,l,map)
   implicit none
   BEGIN_DOC
   ! Returns one integral <ij|kl> in the MO basis
+  ! i(1)j(1) 1/r12 k(2)l(2)
   END_DOC
   integer, intent(in)            :: i,j,k,l
   integer*8                      :: idx
@@ -292,6 +293,7 @@ double precision function mo_bielec_integral(i,j,k,l)
   implicit none
   BEGIN_DOC
   ! Returns one integral <ij|kl> in the MO basis
+  ! i(1)j(1) 1/r12 k(2)l(2)
   END_DOC
   integer, intent(in)            :: i,j,k,l
   double precision               :: get_mo_bielec_integral
@@ -306,6 +308,7 @@ subroutine get_mo_bielec_integrals(j,k,l,sze,out_val,map)
   BEGIN_DOC
   ! Returns multiple integrals <ij|kl> in the MO basis, all
   ! i for j,k,l fixed.
+  ! i(1)j(1) 1/r12 k(2)l(2)
   END_DOC
   integer, intent(in)            :: j,k,l, sze
   real(integral_kind), intent(out) :: out_val(sze)
@@ -327,6 +330,7 @@ subroutine get_mo_bielec_integrals_existing_ik(j,l,sze,out_array,map)
   implicit none
   BEGIN_DOC
   ! Returns multiple integrals <ij|kl> in the MO basis, all
+  ! i(1)j(1) 1/r12 k(2)l(2)
   ! i for j,k,l fixed.
   END_DOC
   integer, intent(in)            :: j,l, sze

src/Dets/H_apply_template.f

@@ -26,6 +26,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
   integer                        :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
   
   double precision               :: mo_bielec_integral
+  logical                        :: is_a_two_holes_two_particles
   integer, allocatable           :: ia_ja_pairs(:,:,:)
   integer, allocatable           :: ib_jb_pairs(:,:)
   double precision               :: diag_H_mat_elem
@@ -162,6 +163,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
           k = ishft(j_b-1,-bit_kind_shift)+1
           l = j_b-ishft(k-1,bit_kind_shift)-1
           key(k,other_spin) = ibset(key(k,other_spin),l)
+          $filter2h2p
           key_idx += 1
           do k=1,N_int
             keys_out(k,1,key_idx) = key(k,1)
@@ -210,6 +212,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
         k = ishft(j_b-1,-bit_kind_shift)+1
         l = j_b-ishft(k-1,bit_kind_shift)-1
         key(k,ispin) = ibset(key(k,ispin),l)
+        $filter2h2p
         key_idx += 1
         do k=1,N_int
           keys_out(k,1,key_idx) = key(k,1)
@@ -267,6 +270,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param
   integer                        :: kk,pp,other_spin,key_idx
   integer                        :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2)
   integer                        :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2)
+  logical                        :: is_a_two_holes_two_particles
   
   integer, allocatable           :: ia_ja_pairs(:,:,:)
   logical, allocatable           :: array_pairs(:,:)
@@ -333,6 +337,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param
       k_a = ishft(j_a-1,-bit_kind_shift)+1
       l_a = j_a-ishft(k_a-1,bit_kind_shift)-1
       hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a)
+      $filter2h2p
       key_idx += 1
       do k=1,N_int
         keys_out(k,1,key_idx) = hole(k,1)

src/Dets/program_beginer_determinants.irp.f

@@ -0,0 +1,138 @@
+program pouet
+ implicit none
+ print*,'HF energy = ',ref_bitmask_energy + nuclear_repulsion
+ call routine
+
+end
+subroutine routine 
+ use bitmasks
+ implicit none
+ integer :: i,j,k,l
+ double precision :: hij,get_mo_bielec_integral
+ double precision :: hmono,h_bi_ispin,h_bi_other_spin
+ integer(bit_kind),allocatable  :: key_tmp(:,:)
+ integer, allocatable :: occ(:,:)
+ integer       :: n_occ_alpha, n_occ_beta
+ ! First checks 
+ print*,'N_int            = ',N_int
+ print*,'mo_tot_num       = ',mo_tot_num
+ print*,'mo_tot_num / 64+1= ',mo_tot_num/64+1
+ ! We print the HF determinant
+ do i = 1, N_int
+  print*,'ref_bitmask(i,1) = ',ref_bitmask(i,1)
+  print*,'ref_bitmask(i,2) = ',ref_bitmask(i,2)
+ enddo
+ print*,''
+ print*,'Hartree Fock determinant ...'
+ call debug_det(ref_bitmask,N_int)
+ allocate(key_tmp(N_int,2))
+ ! We initialize key_tmp to the Hartree Fock one
+ key_tmp = ref_bitmask
+ integer :: i_hole,i_particle,ispin,i_ok,other_spin
+ ! We do a mono excitation on the top of the HF determinant
+ write(*,*)'Enter the (hole, particle) couple for the mono excitation ...'
+ read(5,*)i_hole,i_particle
+!!i_hole = 4
+!!i_particle = 20
+ write(*,*)'Enter the ispin variable ...'
+ write(*,*)'ispin = 1 ==> alpha '
+ write(*,*)'ispin = 2 ==> beta  '
+ read(5,*)ispin
+ if(ispin == 1)then
+  other_spin = 2
+ else if(ispin == 2)then 
+  other_spin = 1
+ else
+  print*,'PB !! '
+  print*,'ispin must be 1 or 2 !'
+  stop
+ endif
+!!ispin = 1
+ call do_mono_excitation(key_tmp,i_hole,i_particle,ispin,i_ok)
+ ! We check if it the excitation was possible with "i_ok"
+ if(i_ok == -1)then
+  print*,'i_ok = ',i_ok
+  print*,'You can not do this excitation because of Pauli principle ...'
+  print*,'check your hole particle couple, there must be something wrong ...'
+  stop
+
+ endif
+ print*,'New det = '
+ call debug_det(key_tmp,N_int)
+ call i_H_j(key_tmp,ref_bitmask,N_int,hij)
+ ! We calculate the H matrix element between the new determinant and HF
+ print*,'<D_i|H|D_j>    = ',hij
+ print*,''
+ print*,''
+ print*,'Recalculating it old school style ....'
+ print*,''
+ print*,''
+ ! We recalculate this old school style !!!
+ ! Mono electronic part
+ hmono = mo_mono_elec_integral(i_hole,i_particle)
+ print*,''
+ print*,'Mono electronic part '
+ print*,''
+ print*,'<D_i|h(1)|D_j> = ',hmono
+ h_bi_ispin = 0.d0
+ h_bi_other_spin = 0.d0
+ print*,''
+ print*,'Getting all the info for the calculation of the bi electronic part ...'
+ print*,''
+ allocate (occ(N_int*bit_kind_size,2))
+ ! We get the occupation of the alpha electrons in occ(:,1)
+ call bitstring_to_list(key_tmp(1,1), occ(1,1), n_occ_alpha, N_int)
+ print*,'n_occ_alpha    = ',n_occ_alpha
+ print*,'elec_alpha_num = ',elec_alpha_num
+ ! We get the occupation of the beta electrons in occ(:,2)
+ call bitstring_to_list(key_tmp(1,2), occ(1,2), n_occ_beta, N_int)
+ print*,'n_occ_beta    = ',n_occ_beta
+ print*,'elec_beta_num = ',elec_beta_num
+ ! We print the occupation of the alpha electrons 
+ print*,'Alpha electrons !'
+ do i = 1, n_occ_alpha
+  print*,'i = ',i
+  print*,'occ(i,1) = ',occ(i,1)
+ enddo
+ ! We print the occupation of the beta electrons 
+ print*,'Alpha electrons !'
+ do i = 1, n_occ_beta
+  print*,'i = ',i
+  print*,'occ(i,2) = ',occ(i,2)
+ enddo
+ integer                        :: exc(0:2,2,2),degree,h1,p1,h2,p2,s1,s2
+ double precision :: phase
+
+ call get_excitation_degree(key_tmp,ref_bitmask,degree,N_int)
+ print*,'degree = ',degree
+ call get_mono_excitation(ref_bitmask,key_tmp,exc,phase,N_int)
+ call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+ print*,'h1  = ',h1
+ print*,'p1  = ',p1
+ print*,'s1  = ',s1
+ print*,'phase = ',phase
+ do i = 1, elec_num_tab(ispin)
+  integer :: orb_occupied
+  orb_occupied = occ(i,ispin)
+  h_bi_ispin += get_mo_bielec_integral(i_hole,orb_occupied,i_particle,orb_occupied,mo_integrals_map) & 
+               -get_mo_bielec_integral(i_hole,i_particle,orb_occupied,orb_occupied,mo_integrals_map)
+ enddo
+ print*,'h_bi_ispin = ',h_bi_ispin
+
+ do i = 1, elec_num_tab(other_spin)
+  orb_occupied = occ(i,other_spin)
+  h_bi_other_spin += get_mo_bielec_integral(i_hole,orb_occupied,i_particle,orb_occupied,mo_integrals_map) 
+ enddo
+ print*,'h_bi_other_spin  = ',h_bi_other_spin
+ print*,'h_bi_ispin + h_bi_other_spin = ',h_bi_ispin + h_bi_other_spin 
+
+ print*,'Total matrix element =  ',phase*(h_bi_ispin + h_bi_other_spin + hmono)
+!i = 1
+!j = 1
+!k = 1
+!l = 1
+!hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
+!print*,'<ij|kl> = ',hij
+
+
+end

src/MonoInts/pot_mo_ints.irp.f

@@ -11,9 +11,9 @@ BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_to
  do i = 1, mo_tot_num
    do j = 1, mo_tot_num
     do i1 = 1,ao_num
-     c_i1 = mo_coef(i1,i)
+     c_i1 = mo_coef(i1,i)  ! <AO(i1)|MO(i)> 
      do j1 = 1,ao_num
-       c_j1 = c_i1*mo_coef(j1,j)
+       c_j1 = c_i1*mo_coef(j1,j)   ! <AO(j1)|MO(j)> 
        mo_nucl_elec_integral(j,i) = mo_nucl_elec_integral(j,i) + &
                              c_j1 * ao_nucl_elec_integral(j1,i1)
      enddo

src/NEEDED_MODULES

@@ -1 +1 @@
-AOs BiInts Bitmask Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS
+AOs BiInts Bitmask Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS Denstity_stuff CAS_SD_selected DDCI_selected Two_body_density_matrix