working on perturbative rdm

src/casscf/test_pert_2rdm.irp.f

@@ -2,6 +2,8 @@ program test_pert_2rdm
  implicit none
  read_wf = .True.
  touch read_wf
+ pert_2rdm = .True.
+ touch pert_2rdm
 
 !provide is_pert_2rdm_provided
  call test
@@ -22,21 +24,50 @@ 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
   enddo
  enddo
+ print*,'accu = ',accu
+ do ii = 1, n_orb_pert_rdm
+  i = list_orb_pert_rdm(ii) 
+  do kk = 1, n_orb_pert_rdm
+   k= list_orb_pert_rdm(kk) 
+   do jj = 1, n_core_inact_orb
+    j = list_core_inact(jj)
+    l = j
+    integral = get_two_e_integral(i,j,k,l,mo_integrals_map)
+    double precision :: exchange 
+    exchange = get_two_e_integral(i,j,l,k,mo_integrals_map)
+    accu    +=  pert_1rdm_provider(kk,ii) * ( 2.d0 * integral - exchange)
+   enddo
+  enddo
+ enddo
+ print*,'accu = ',accu
+ double precision :: accu_1
+ accu_1 = 0.d0
  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) 
-   accu +=  pert_1rdm_provider(j,i) * mo_one_e_integrals(j,i)
+   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)
+   endif
+   if(dabs(pert_1rdm_provider(jj,ii) * mo_one_e_integrals(j,i)).gt.1.d-10)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)
   enddo
  enddo
- print*,'accu = ',accu
+ print*,'accu_1 = ',accu_1
+ print*,'*********'
+ print*,'*********'
+ print*,'accu = ',accu + accu_1
+
+ print*,'pt2  = ',pt2_pert_2rdm
 end

src/cipsi/EZFIO.cfg

@@ -0,0 +1,5 @@
+[pert_2rdm]
+type: logical
+doc: If true, computes the one- and two-body rdms with perturbation theory
+interface: ezfio,provider,ocaml
+default: False

src/cipsi/pert_rdm_providers.irp.f

@@ -8,12 +8,19 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), pert_2rdm_lock]
   call omp_init_lock(pert_2rdm_lock)
 END_PROVIDER
 
-BEGIN_PROVIDER [logical , pert_2rdm ]
+BEGIN_PROVIDER [ integer(omp_lock_kind), pert_1rdm_lock]
+  use f77_zmq
   implicit none
- pert_2rdm = .True.
+  call omp_init_lock(pert_1rdm_lock)
 END_PROVIDER
 
-BEGIN_PROVIDER [logical, is_pert_2rdm_provided ]
+!BEGIN_PROVIDER [logical , pert_2rdm ]
+! implicit none
+! pert_2rdm = .False.
+!END_PROVIDER 
+
+ BEGIN_PROVIDER [logical, is_pert_2rdm_provided ]
+&BEGIN_PROVIDER [double precision, pt2_pert_2rdm, (N_states)]
  implicit none
  is_pert_2rdm_provided = .True.
  provide pert_2rdm_provider 
@@ -22,10 +29,12 @@ BEGIN_PROVIDER [logical, is_pert_2rdm_provided ]
   double precision :: pt2(N_states),relative_error, error(N_states),variance(N_states),norm(N_states)
   relative_error = 0.d0
   pert_2rdm_provider = 0.d0
-  call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, &                                         
+  pert_1rdm_provider = 0.d0
+  pert_1rdm_provider_bis = 0.d0
+  call ZMQ_pt2(psi_energy_with_nucl_rep, pt2_pert_2rdm,relative_error,error,variance, &                                         
           norm,0) ! Stochastic PT2
   print*,'is_pert_2rdm_provided = ',is_pert_2rdm_provided
-  print*,'pt2 = ',pt2
+  print*,'pt2 = ',pt2_pert_2rdm
  endif
 END_PROVIDER 
 
@@ -46,22 +55,8 @@ BEGIN_PROVIDER [integer, list_orb_pert_rdm, (n_orb_pert_rdm)]
 
 END_PROVIDER
 
-BEGIN_PROVIDER [double precision, pert_1rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm)]
- implicit none
- integer :: i,j,k,l
- if(is_pert_2rdm_provided)then
-  pert_1rdm_provider = 0.d0
-  do i = 1, n_orb_pert_rdm
-   do j = 1, n_orb_pert_rdm
-    do k = 1, n_orb_pert_rdm
-     pert_1rdm_provider(j,i) += 2.d0 * pert_2rdm_provider(i,k,j,k)
-    enddo
-   enddo
-  enddo
- endif
-END_PROVIDER
-
  BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm)]
+&BEGIN_PROVIDER [double precision, pert_1rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm)]
 &BEGIN_PROVIDER [double precision, pert_1rdm_provider_bis, (n_orb_pert_rdm,n_orb_pert_rdm)]
  implicit none
 
@@ -95,10 +90,14 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo
   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))
+  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
@@ -176,12 +175,21 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo
       sum_e_pert = 0d0
       integer                        :: exc(0:2,2,2)
       integer                        :: degree
-      double precision               :: phase
+      double precision               :: phase,hij,hij2
-      call get_excitation(det,HF_bitmask,exc,degree,phase,N_int) 
+      integer                        :: hh1,pp1,hh2,pp2,ss1,ss2
-      call get_excitation_degree(det,HF_bitmask,degree,N_int)
+      call get_excitation_degree(HF_bitmask,det,degree,N_int)
-      if(degree == 1)cycle
+      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
-      do istate=1,N_states
+!!!!!!!!!!!!!!!!!! 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
@@ -190,12 +198,26 @@ 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
-!       print*,coef,alpha_h_psi,e_pert
         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)
+!       call i_H_j(psi_det_generators(1,1,2),det,N_int,hij2)
+        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)
@@ -204,8 +226,8 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo
           ! Variance selection
           sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
         endif
-      end do
+!     end do
-      call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,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)
 
       if(sum_e_pert <= buf%mini) then
         call add_to_selection_buffer(buf, det, sum_e_pert)
@@ -213,6 +235,7 @@ subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fo
     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/selection.irp.f

@@ -616,6 +616,10 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
             if(.not.pert_2rdm)then
              call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf)
             else 
+!            print*,'fullinteresting(0) just before ',fullinteresting(0)
+!            print*,'i_generator = ',i_generator
+!            call debug_det(psi_det_generators(1,1,i_generator),N_int)
+!            pause
              call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0))
             endif
           end if

src/cipsi/update_2rdm.irp.f

@@ -1,6 +1,6 @@
 use bitmasks
 
-subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff)
+subroutine 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)
  implicit none
  integer, intent(in) :: n_det_connection,sze_buff
  double precision, intent(in) :: coef(N_states)
@@ -9,27 +9,49 @@ subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connectio
  double precision, intent(in)  :: psi_coef_connection_reverse(N_states,n_det_connection)
  integer,           intent(inout) :: keys(4,sze_buff),nkeys
  double precision,  intent(inout) :: values(sze_buff)
+ integer,           intent(inout) :: keys_1e(2,sze_buff),nkeys_1e
+ double precision,  intent(inout) :: values_1e(sze_buff)
  integer :: i,j
  integer                        :: exc(0:2,2,2)
  integer                        :: degree
  double precision               :: phase, contrib
  do i = 1, n_det_connection
-  call get_excitation(det,psi_det_connection(1,1,i),exc,degree,phase,N_int) 
+  call get_excitation_degree(det,psi_det_connection(1,1,i),degree,N_int)
   if(degree.gt.2)cycle
+  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
+  if(degree == 1)then  ! check for the length of buffer for the ONE BODY DM
+   if(nkeys_1e + 2 .ge. sze_buff)then
+    call update_keys_values_1e(keys_1e,values_1e,nkeys_1e,n_orb_pert_rdm,pert_1rdm_provider,pert_1rdm_lock)
+   endif
+  endif
+  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
-   call update_buffer_single_exc_rdm(det,psi_det_connection(1,1,i),exc,phase,contrib,nkeys,keys,values,sze_buff)
+   call update_buffer_single_exc_rdm(det,psi_det_connection(1,1,i),exc,phase,contrib,nkeys,keys,values,keys_1e,values_1e,nkeys_1e,sze_buff)
   else 
   ! case of double excitations 
-   if (nkeys + 4 .ge. sze_buff)then
+   if (nkeys + 4 .ge. sze_buff)then! check for the length of buffer for the TWO BODY DM
     call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
     nkeys = 0
    endif
@@ -38,10 +60,12 @@ subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connectio
  enddo
  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
 
-subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,sze_buff)
+subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,keys_1e,values_1e,nkeys_1e,sze_buff)
  implicit none
  integer, intent(in) :: sze_buff
  integer(bit_kind), intent(in) :: det1(N_int,2)
@@ -50,6 +74,8 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v
  double precision,intent(in)    :: phase, contrib
  integer, intent(inout)         :: nkeys, keys(4,sze_buff)
  double precision, intent(inout):: values(sze_buff)
+ integer, intent(inout)         :: nkeys_1e, keys_1e(2,sze_buff)
+ double precision, intent(inout):: values_1e(sze_buff)
  
  integer                        :: occ(N_int*bit_kind_size,2)
  integer                        :: n_occ_ab(2),ispin,other_spin
@@ -69,15 +95,25 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v
   ispin = 2
   other_spin = 1
  endif
-!print*,'h1,p1,s1',h1,p1,ispin
+ if(list_orb_reverse_pert_rdm(h1).le.0)return
- if(list_orb_reverse_pert_rdm(h1).lt.0)return
  h1 = list_orb_reverse_pert_rdm(h1)
- if(list_orb_reverse_pert_rdm(p1).lt.0)return
+ 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
+ nkeys_1e += 1
+ values_1e(nkeys_1e) = 0.5d0 * contrib * phase
+ keys_1e(1,nkeys_1e) = h1
+ keys_1e(2,nkeys_1e) = p1
+ nkeys_1e += 1
+ 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)
-  if(list_orb_reverse_pert_rdm(h2).lt.0)return
+  if(list_orb_reverse_pert_rdm(h2).le.0)cycle
   h2 = list_orb_reverse_pert_rdm(h2)
 
   nkeys += 1
@@ -96,7 +132,7 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v
  !update the same spin part 
  do i = 1, n_occ_ab(ispin)
   h2 = occ(i,ispin)
-  if(list_orb_reverse_pert_rdm(h2).lt.0)return
+  if(list_orb_reverse_pert_rdm(h2).le.0)cycle
   h2 = list_orb_reverse_pert_rdm(h2)
 
   nkeys += 1
@@ -129,6 +165,8 @@ subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,v
  enddo 
  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
 
@@ -148,13 +186,13 @@ subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_
   p1 = exc(1,2,1)
   p2 = exc(1,2,2)
   ! check if the orbitals involved are within the orbital range
-  if(list_orb_reverse_pert_rdm(h1).lt.0)return
+  if(list_orb_reverse_pert_rdm(h1).le.0)return
   h1 = list_orb_reverse_pert_rdm(h1)
-  if(list_orb_reverse_pert_rdm(h2).lt.0)return
+  if(list_orb_reverse_pert_rdm(h2).le.0)return
   h2 = list_orb_reverse_pert_rdm(h2)
-  if(list_orb_reverse_pert_rdm(p1).lt.0)return
+  if(list_orb_reverse_pert_rdm(p1).le.0)return
   p1 = list_orb_reverse_pert_rdm(p1)
-  if(list_orb_reverse_pert_rdm(p2).lt.0)return
+  if(list_orb_reverse_pert_rdm(p2).le.0)return
   p2 = list_orb_reverse_pert_rdm(p2)
   nkeys += 1
   values(nkeys) = 0.5d0 * contrib * phase
@@ -186,13 +224,14 @@ 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
-  if(list_orb_reverse_pert_rdm(h1).lt.0)return
+ !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).lt.0)return
+  if(list_orb_reverse_pert_rdm(h2).le.0)return
   h2 = list_orb_reverse_pert_rdm(h2)
-  if(list_orb_reverse_pert_rdm(p1).lt.0)return
+  if(list_orb_reverse_pert_rdm(p1).le.0)return
   p1 = list_orb_reverse_pert_rdm(p1)
-  if(list_orb_reverse_pert_rdm(p2).lt.0)return
+  if(list_orb_reverse_pert_rdm(p2).le.0)return
   p2 = list_orb_reverse_pert_rdm(p2)
   nkeys += 1
   values(nkeys) = 0.5d0 * contrib * phase

src/determinants/determinants.irp.f

@@ -306,6 +306,16 @@ END_PROVIDER
 
 END_PROVIDER
 
+ BEGIN_PROVIDER [double precision, psi_coef_sorted_reverse, (N_states,psi_det_size)
+ implicit none
+ integer :: i,j
+ do i = 1, N_states
+  do j = 1, N_det
+   psi_coef_sorted_reverse(i,j) = psi_coef_sorted(j,i)
+  enddo
+ enddo
+ END_PROVIDER 
+
  BEGIN_PROVIDER [ integer(bit_kind), psi_det_sorted_bit, (N_int,2,psi_det_size) ]
 &BEGIN_PROVIDER [ double precision, psi_coef_sorted_bit, (psi_det_size,N_states) ]
    implicit none

src/two_body_rdm/orb_range_routines_openmp.irp.f

@@ -566,3 +566,24 @@ subroutine update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
 
 end
 
+
+subroutine update_keys_values_1e(keys_1e,values_1e,nkeys_1e,dim1,big_array,lock_1rdm)
+ use omp_lib
+ implicit none
+ integer, intent(in) :: nkeys_1e,dim1
+ integer, intent(in) :: keys_1e(2,nkeys_1e)
+ double precision, intent(in) :: values_1e(nkeys_1e)
+ double precision, intent(inout) :: big_array(dim1,dim1)
+ 
+ integer(omp_lock_kind),intent(inout):: lock_1rdm
+ integer :: i,h1,p1
+ call omp_set_lock(lock_1rdm)
+ do i = 1, nkeys_1e
+  h1 = keys_1e(1,i)
+  p1 = keys_1e(2,i)
+  big_array(h1,p1) += values_1e(i)
+ enddo
+ call omp_unset_lock(lock_1rdm)
+
+end
+