perturbative rdm gives an error of 1 percent with respect to the pt2

src/casscf/test_pert_2rdm.irp.f

@@ -24,10 +24,10 @@ program test_pert_2rdm
     do ll = 1, n_orb_pert_rdm
      l = list_orb_pert_rdm(ll) 
      integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
-     if(dabs(pert_2rdm_provider(ii,jj,kk,ll) * integral).gt.1.d-12)then
+    !if(dabs(pert_2rdm_provider(ii,jj,kk,ll) * integral).gt.1.d-12)then
-      print*,i,j,k,l
+    ! print*,i,j,k,l
-      print*,pert_2rdm_provider(ii,jj,kk,ll) , integral, pert_2rdm_provider(ii,jj,kk,ll)* integral
+    ! print*,pert_2rdm_provider(ii,jj,kk,ll) , integral, pert_2rdm_provider(ii,jj,kk,ll)* integral
-     endif
+    !endif
      accu    += pert_2rdm_provider(ii,jj,kk,ll) * integral
     enddo
    enddo
@@ -51,14 +51,22 @@ program test_pert_2rdm
  print*,'accu = ',accu
  double precision :: accu_1
  accu_1 = 0.d0
+!print*,'pert_1rdm_provider(1,2)',pert_1rdm_provider(1,2)
+!print*,'pert_1rdm_provider(2,1)',pert_1rdm_provider(2,1)
+!print*,'mo_one_e_integrals(2,1)',mo_one_e_integrals(2,1)
+!print*,'mo_one_e_integrals(1,2)',mo_one_e_integrals(1,2)
+!print*,'F(1,2)                 ',fock_operator_closed_shell_ref_bitmask(1,2)
+!print*,'(12|66) = <16|26>      ',get_two_e_integral(1,6,2,6,mo_integrals_map)
+!print*,'(16|26) = <16|62>      ',get_two_e_integral(1,6,6,2,mo_integrals_map)
+!print*,'(12|22) = <12|22>      ',get_two_e_integral(1,2,2,2,mo_integrals_map)
  do ii = 1, n_orb_pert_rdm
   i = list_orb_pert_rdm(ii) 
   do jj = 1, n_orb_pert_rdm
    j = list_orb_pert_rdm(jj) 
-   if(dabs(pert_1rdm_provider(jj,ii) - pert_1rdm_provider_bis(jj,ii)).gt.1.d-10)then
+  !if(dabs(pert_1rdm_provider(jj,ii) - pert_1rdm_provider_bis(jj,ii)).gt.1.d-10)then
-    print*,ii,jj,pert_1rdm_provider(jj,ii),pert_1rdm_provider_bis(jj,ii)
+  ! print*,ii,jj,pert_1rdm_provider(jj,ii),pert_1rdm_provider_bis(jj,ii)
-   endif
+  !endif
-   if(dabs(pert_1rdm_provider(jj,ii) * mo_one_e_integrals(j,i)).gt.1.d-10)then
+   if(dabs(pert_1rdm_provider(jj,ii) * mo_one_e_integrals(j,i)).gt.1.d-12)then
     print*,j,i,pert_1rdm_provider(jj,ii) , mo_one_e_integrals(j,i), pert_1rdm_provider(jj,ii) * mo_one_e_integrals(j,i)
    endif
    accu_1 +=  pert_1rdm_provider(jj,ii) * mo_one_e_integrals(j,i)

src/cipsi/pert_rdm_providers.irp.f

@@ -113,6 +113,9 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo
     E_shift = psi_det_Hii(i_generator) - psi_occ_pattern_Hii(j)
   endif
 
+  double precision, allocatable  :: hij(:)
+  allocate(hij(N_det_generators))
+
   do p1=1,mo_num
     if(bannedOrb(p1, s1)) cycle
     ib = 1
@@ -173,23 +176,8 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo
       Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
 
       sum_e_pert = 0d0
-      integer                        :: exc(0:2,2,2)
-      integer                        :: degree
-      double precision               :: phase,hij,hij2
-      integer                        :: hh1,pp1,hh2,pp2,ss1,ss2
-      call get_excitation_degree(HF_bitmask,det,degree,N_int)
-      if(degree.ne.2)cycle
-      call get_excitation_degree(psi_det_generators(1,1,2),det,degree,N_int)
-      if(degree.ne.1)cycle
-      call get_excitation(HF_bitmask,det,exc,degree,phase,N_int) 
-      call decode_exc(exc,degree,hh1,pp1,hh2,pp2,ss1,ss2)
-      if(hh1 .ne. 1)cycle
-      if(pp1 .ne. 6)cycle
-      if(ss1 .ne. 1)cycle
-!     if (exc(0,1,1) .ne. 1) cycle !only double alpha/beta
 !!!!!!!!!!!!!!!!!! LOOP OVER STATES
-!     do istate=1,N_states
+      do istate=1,N_states
-        istate=1
         delta_E = E0(istate) - Hii + E_shift
         alpha_h_psi = mat(istate, p1, p2)
         val = alpha_h_psi + alpha_h_psi
@@ -198,44 +186,21 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo
             tmp = -tmp
         endif
         e_pert = 0.5d0 * (tmp - delta_E)
-        if(dabs(e_pert).lt.1.d-07)cycle
-!       if(dabs(e_pert).gt.1.d-06)cycle
-        write(*,*),'----'
-      print*,'hh1',hh1,'pp1',pp1,'ss1',ss1
-      print*,'hh2',hh2,'pp2',pp2,'ss2',ss2
-      call get_excitation(psi_det_generators(1,1,2),det,exc,degree,phase,N_int) 
-      call decode_exc(exc,degree,hh1,pp1,hh2,pp2,ss1,ss2)
-      print*,'hh1',hh1,'pp1',pp1,'ss1',ss1
-      print*,'hh2',hh2,'pp2',pp2,'ss2',ss2
 
         coef(istate) = e_pert / alpha_h_psi
         pt2(istate) = pt2(istate) + e_pert
         variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi
         norm(istate) = norm(istate) + coef(istate) * coef(istate)
 
-!       call i_H_j(HF_bitmask,det,N_int,hij)
+      end do
-!       call i_H_j(psi_det_generators(1,1,2),det,N_int,hij2)
+      call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,nkeys_1e,keys_1e,values_1e,sze_buff)
-        write(*,'(100(F16.13,X))'),hij,hij2
-        write(*,'(100(F16.13,X))'),phase,coef,alpha_h_psi,hij,e_pert
 
-        if (weight_selection /= 5) then
-          ! Energy selection
-          sum_e_pert = sum_e_pert + e_pert * selection_weight(istate)
-
-        else
-          ! Variance selection
-          sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
-        endif
-!     end do
-      call give_2rdm_pert_contrib(det,coef,psi_det_sorted,psi_coef_sorted_reverse,n_det,nkeys,keys,values,nkeys_1e,keys_1e,values_1e,sze_buff)
-
-      if(sum_e_pert <= buf%mini) then
-        call add_to_selection_buffer(buf, det, sum_e_pert)
-      end if
     end do
   end do
   call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
+  nkeys = 0
   call update_keys_values_1e(keys_1e,values_1e,nkeys_1e,n_orb_pert_rdm,pert_1rdm_provider,pert_1rdm_lock)
+  nkeys_1e = 0
 end
 
 

src/cipsi/routines_debug_pert_rdm.irp.f

@@ -0,0 +1,228 @@
+
+subroutine fill_buffer_double_rdm_debug(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf, psi_det_connection, psi_coef_connection_reverse, n_det_connection)
+  use bitmasks
+  use selection_types
+  implicit none
+  
+  integer, intent(in)           :: n_det_connection
+  double precision, intent(in)  :: psi_coef_connection_reverse(N_states,n_det_connection)
+  integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection)
+  integer, intent(in) :: i_generator, sp, h1, h2
+  double precision, intent(in) :: mat(N_states, mo_num, mo_num)
+  logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num)
+  double precision, intent(in)           :: fock_diag_tmp(mo_num)
+  double precision, intent(in)    :: E0(N_states)
+  double precision, intent(inout) :: pt2(N_states)
+  double precision, intent(inout) :: variance(N_states)
+  double precision, intent(inout) :: norm(N_states)
+  type(selection_buffer), intent(inout) :: buf
+  logical :: ok
+  integer :: s1, s2, p1, p2, ib, j, istate
+  integer(bit_kind) :: mask(N_int, 2), det(N_int, 2)
+  double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp, alpha_h_psi, coef(N_states)
+  double precision, external :: diag_H_mat_elem_fock
+  double precision :: E_shift
+
+  logical, external :: detEq
+  double precision, allocatable :: values(:)
+  integer, allocatable          :: keys(:,:)
+  integer                       :: nkeys
+  double precision, allocatable :: values_1e(:)
+  integer, allocatable          :: keys_1e(:,:)
+  integer                       :: nkeys_1e
+  integer :: sze_buff
+  sze_buff = 5 * mo_num ** 2 
+  allocate(keys(4,sze_buff),values(sze_buff),keys_1e(2,sze_buff),values_1e(sze_buff))
+  nkeys = 0
+  nkeys_1e = 0
+  if(sp == 3) then
+    s1 = 1
+    s2 = 2
+  else
+    s1 = sp
+    s2 = sp
+  end if
+  call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int)
+  E_shift = 0.d0
+
+  if (h0_type == 'SOP') then
+    j = det_to_occ_pattern(i_generator)
+    E_shift = psi_det_Hii(i_generator) - psi_occ_pattern_Hii(j)
+  endif
+
+  double precision, allocatable  :: hij(:)
+  allocate(hij(N_det_generators))
+
+  do p1=1,mo_num
+    if(bannedOrb(p1, s1)) cycle
+    ib = 1
+    if(sp /= 3) ib = p1+1
+
+    do p2=ib,mo_num
+
+! -----
+! /!\ Generating only single excited determinants doesn't work because a
+! determinant generated by a single excitation may be doubly excited wrt
+! to a determinant of the future. In that case, the determinant will be
+! detected as already generated when generating in the future with a
+! double excitation.
+!     
+!      if (.not.do_singles) then
+!        if ((h1 == p1) .or. (h2 == p2)) then
+!          cycle
+!        endif
+!      endif
+!
+!      if (.not.do_doubles) then
+!        if ((h1 /= p1).and.(h2 /= p2)) then
+!          cycle
+!        endif
+!      endif
+! -----
+
+      if(bannedOrb(p2, s2)) cycle
+      if(banned(p1,p2)) cycle
+
+
+      if( sum(abs(mat(1:N_states, p1, p2))) == 0d0) cycle
+      call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
+
+      if (do_only_cas) then
+        integer, external :: number_of_holes, number_of_particles
+        if (number_of_particles(det)>0) then
+          cycle
+        endif
+        if (number_of_holes(det)>0) then
+          cycle
+        endif
+      endif
+
+      if (do_ddci) then
+        logical, external  :: is_a_two_holes_two_particles
+        if (is_a_two_holes_two_particles(det)) then
+          cycle
+        endif
+      endif
+
+      if (do_only_1h1p) then
+        logical, external :: is_a_1h1p
+        if (.not.is_a_1h1p(det)) cycle
+      endif
+
+
+      Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
+
+      sum_e_pert = 0d0
+      integer                        :: exc(0:2,2,2)
+      integer                        :: degree
+      double precision               :: phase,hij1,hij2
+      integer                        :: hh1,pp1,hh2,pp2,ss1,ss2
+!     call get_excitation_degree(HF_bitmask,det,degree,N_int)
+!     if(degree.gt.2)cycle
+!     call get_excitation_degree(psi_det_generators(1,1,2),det,degree,N_int)
+!     if(degree.gt.2)cycle
+!     call get_excitation(HF_bitmask,det,exc,degree,phase,N_int) 
+!     call decode_exc(exc,degree,hh1,pp1,hh2,pp2,ss1,ss2)
+!     if(hh1 .ne. 1)cycle
+!     if(pp1 .ne. 6)cycle
+!     if(ss1 .ne. 1)cycle
+!     if (exc(0,1,1) .ne. 1) cycle !only double alpha/beta
+!!!!!!!!!!!!!!!!!! LOOP OVER STATES
+!     do istate=1,N_states
+        istate=1
+        delta_E = E0(istate) - Hii + E_shift
+        alpha_h_psi = mat(istate, p1, p2)
+        val = alpha_h_psi + alpha_h_psi
+        tmp = dsqrt(delta_E * delta_E + val * val)
+        if (delta_E < 0.d0) then
+            tmp = -tmp
+        endif
+        e_pert = 0.5d0 * (tmp - delta_E)
+!       if(dabs(e_pert).lt.1.d-07)cycle
+!       if(dabs(e_pert).gt.1.d-06)cycle
+!       write(*,*),'----'
+!     print*,'hh1',hh1,'pp1',pp1,'ss1',ss1
+!     print*,'hh2',hh2,'pp2',pp2,'ss2',ss2
+!     call get_excitation(psi_det_generators(1,1,2),det,exc,degree,phase,N_int) 
+!     call decode_exc(exc,degree,hh1,pp1,hh2,pp2,ss1,ss2)
+
+        coef(istate) = e_pert / alpha_h_psi
+        pt2(istate) = pt2(istate) + e_pert
+        variance(istate) = variance(istate) + alpha_h_psi * alpha_h_psi
+        norm(istate) = norm(istate) + coef(istate) * coef(istate)
+
+        integer :: igen_tmp
+        double precision :: alphahpsi
+        alphahpsi = 0.d0
+        hij = 0.d0
+        do igen_tmp = 1, N_det_generators
+         call get_excitation_degree(psi_det_generators(1,1,igen_tmp),det,degree,N_int)
+         if(degree.gt.2)cycle
+         call i_H_j(psi_det_generators(1,1,igen_tmp),det,N_int,hij(igen_tmp))
+         alphahpsi += psi_coef_generators(igen_tmp,istate) * hij(igen_tmp)
+        enddo
+        if(dabs(alphahpsi - alpha_h_psi).gt.1.d-12)then
+         print*,''
+         print*,''
+         print*,'alphhpsi   = ',alphahpsi,alpha_h_psi
+         print*,'<igen|H|k> = ',psi_coef_generators(i_generator,istate) * hij(i_generator)
+         call debug_det(det,N_int)
+         stop
+!        call get_excitation_degree(psi_det_generators(1,1,1),psi_det_generators(1,1,2),degree,N_int)
+!        if(degree.gt.2)cycle
+!        call get_excitation(psi_det_generators(1,1,1),det,exc,degree,phase,N_int) 
+!        call decode_exc(exc,degree,hh1,pp1,hh2,pp2,ss1,ss2)
+!        print*,'excitation between generators '
+!        print*,'degree = ',degree
+!        print*,'hh1',hh1,'pp1',pp1,'ss1',ss1
+!        if(hh2.ne.0)then
+!        print*,'hh2',hh2,'pp2',pp2,'ss2',ss2
+!        print*,''
+!        endif
+!        print*,''
+!        call get_excitation_degree(psi_det_generators(1,1,1),det,degree,N_int)
+!        if(degree.gt.2)cycle
+!        call get_excitation(psi_det_generators(1,1,1),det,exc,degree,phase,N_int) 
+!        call decode_exc(exc,degree,hh1,pp1,hh2,pp2,ss1,ss2)
+!        print*,'degree = ',degree
+!        print*,'hh1',hh1,'pp1',pp1,'ss1',ss1
+!        print*,'hh2',hh2,'pp2',pp2,'ss2',ss2
+!        print*,'coef, hij = ',psi_coef_generators(1,istate),hij1
+
+!        call get_excitation_degree(psi_det_generators(1,1,2),det,degree,N_int)
+!        if(degree.gt.2)cycle
+!        call get_excitation(psi_det_generators(1,1,2),det,exc,degree,phase,N_int) 
+!        call decode_exc(exc,degree,hh1,pp1,hh2,pp2,ss1,ss2)
+!        print*,'degree = ',degree
+!        print*,'hh1',hh1,'pp1',pp1,'ss1',ss1
+!        print*,'hh2',hh2,'pp2',pp2,'ss2',ss2
+!        print*,'coef, hij = ',psi_coef_generators(2,istate),hij2
+
+!        print*,'delta 1 = ',hij1 * psi_coef_generators(1,istate) * coef
+!        print*,'delta 2 = ',hij2* psi_coef_generators(2,istate) * coef
+         write(*,'(100(F16.13,X))')coef,alpha_h_psi,e_pert,coef * alpha_h_psi
+        endif
+
+        if (weight_selection /= 5) then
+          ! Energy selection
+          sum_e_pert = sum_e_pert + e_pert * selection_weight(istate)
+
+        else
+          ! Variance selection
+          sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
+        endif
+!     end do
+!     call give_2rdm_pert_contrib(det,coef,psi_det_generators,psi_coef_generators_reverse,n_det,nkeys,keys,values,nkeys_1e,keys_1e,values_1e,sze_buff)
+!     call give_2rdm_pert_contrib(det,coef,psi_det_sorted,psi_coef_sorted_reverse,n_det,nkeys,keys,values,nkeys_1e,keys_1e,values_1e,sze_buff)
+      call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,nkeys_1e,keys_1e,values_1e,sze_buff)
+
+     !if(sum_e_pert <= buf%mini) then
+     !  call add_to_selection_buffer(buf, det, sum_e_pert)
+     !end if
+    end do
+  end do
+  call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
+  call update_keys_values_1e(keys_1e,values_1e,nkeys_1e,n_orb_pert_rdm,pert_1rdm_provider,pert_1rdm_lock)
+end
+
+

src/cipsi/update_2rdm.irp.f

@@ -18,33 +18,24 @@ subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connectio
  do i = 1, n_det_connection
   call get_excitation_degree(det,psi_det_connection(1,1,i),degree,N_int)
   if(degree.gt.2)cycle
+  if(degree.eq.0)then
+   print*,'PB !! there is a perturbative determinant already in the WF !!'
+   call debug_det(det,N_int)
+   call debug_det(psi_det_connection(1,1,i),N_int)
+  endif
   call get_excitation(det,psi_det_connection(1,1,i),exc,degree,phase,N_int) 
-  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-  integer :: h1,p1,h2,p2,s1,s2
-! print*,''
-! print*,''
-! print*,''
-! print*,''
-! print*,'n_det_connection = ',n_det_connection
-
-!  call debug_det(det,N_int)
-!  call debug_det(psi_det_connection(1,1,i),N_int)
-!  print*,'degree = ',degree
-!  print*,'h1,p1,s1',h1,p1,s1
-!  print*,'h2,p2,s2',h2,p2,s2
 
   contrib = 0.d0
   do j = 1, N_states
    contrib += state_average_weight(j) * psi_coef_connection_reverse(j,i)  * coef(j)
   enddo
   ! case of single excitations 
-  if(degree == 1)then  ! check for the length of buffer for the ONE BODY DM
+  if(degree == 1)then 
-   if(nkeys_1e + 2 .ge. sze_buff)then
+   if(nkeys_1e + 2 .ge. sze_buff)then ! check for the length of buffer for the ONE BODY DM
     call update_keys_values_1e(keys_1e,values_1e,nkeys_1e,n_orb_pert_rdm,pert_1rdm_provider,pert_1rdm_lock)
+    nkeys_1e = 0
    endif
-  endif
+   if (nkeys + 6 * elec_alpha_num .ge. sze_buff)then ! check for the length of buffer for the TWO BODY DM
-  if(degree == 1)then ! check for the length of buffer for the TWO BODY DM
-   if (nkeys + 6 * elec_alpha_num .ge. sze_buff)then
     call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
     nkeys = 0
    endif
@@ -98,7 +89,7 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v
  if(list_orb_reverse_pert_rdm(h1).le.0)return
  h1 = list_orb_reverse_pert_rdm(h1)
  if(list_orb_reverse_pert_rdm(p1).le.0)return
-!write(*,'(100(I3,X))')list_orb_reverse_pert_rdm(:)
+
  p1 = list_orb_reverse_pert_rdm(p1)
  pert_1rdm_provider_bis(h1,p1) += 0.5d0 * contrib * phase
  pert_1rdm_provider_bis(p1,h1) += 0.5d0 * contrib * phase
@@ -110,6 +101,7 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v
  values_1e(nkeys_1e) = 0.5d0 * contrib * phase
  keys_1e(1,nkeys_1e) = p1
  keys_1e(2,nkeys_1e) = h1
+
  !update the alpha/beta part
  do i = 1, n_occ_ab(other_spin)
   h2 = occ(i,other_spin)
@@ -200,14 +192,12 @@ subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_
   keys(2,nkeys) = h2
   keys(3,nkeys) = p1
   keys(4,nkeys) = p2
-! print*,contrib 
   nkeys += 1
   values(nkeys) = 0.5d0 * contrib * phase
   keys(1,nkeys) = p1
   keys(2,nkeys) = p2
   keys(3,nkeys) = h1
   keys(4,nkeys) = h2 
-! print*,contrib 
 
  else 
   if (exc(0,1,1) == 2) then
@@ -224,7 +214,6 @@ subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_
    p2 = exc(2,2,2)
   endif
   ! check if the orbitals involved are within the orbital range
- !print*,'h1 = ',h1
   if(list_orb_reverse_pert_rdm(h1).le.0)return
   h1 = list_orb_reverse_pert_rdm(h1)
   if(list_orb_reverse_pert_rdm(h2).le.0)return

src/generators_cas/generators.irp.f

@@ -82,3 +82,12 @@ BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ]
   select_max = huge(1.d0)
 END_PROVIDER
 
+ BEGIN_PROVIDER [double precision, psi_coef_generators_reverse, (N_states,psi_det_size)
+ implicit none
+ integer :: i,j
+ do i = 1, N_states
+  do j = 1, N_det_generators
+   psi_coef_generators_reverse(i,j) = psi_coef_generators(j,i)
+  enddo
+ enddo
+ END_PROVIDER 

src/two_body_rdm/orb_range_routines_openmp.irp.f

@@ -578,10 +578,12 @@ subroutine update_keys_values_1e(keys_1e,values_1e,nkeys_1e,dim1,big_array,lock_
  integer(omp_lock_kind),intent(inout):: lock_1rdm
  integer :: i,h1,p1
  call omp_set_lock(lock_1rdm)
+!print*,'coucoucoucocu'
  do i = 1, nkeys_1e
   h1 = keys_1e(1,i)
   p1 = keys_1e(2,i)
   big_array(h1,p1) += values_1e(i)
+! print*,h1,p1,big_array(h1,p1)
  enddo
  call omp_unset_lock(lock_1rdm)