compiles

src/cipsi/NEED

@@ -3,3 +3,4 @@ zmq
 mpi
 davidson_undressed
 iterations
+two_body_rdm

src/cipsi/pert_rdm_providers.irp.f

@@ -0,0 +1,172 @@
+
+use bitmasks
+
+BEGIN_PROVIDER [logical , pert_2rdm ]
+ implicit none
+ pert_2rdm = .False.
+END_PROVIDER 
+
+BEGIN_PROVIDER [integer, n_orb_pert_rdm]
+ implicit none
+ n_orb_pert_rdm = n_act_orb
+END_PROVIDER 
+
+BEGIN_PROVIDER [integer, list_orb_reverse_pert_rdm, (mo_num)]
+ implicit none
+ list_orb_reverse_pert_rdm = list_act_reverse
+
+END_PROVIDER 
+
+BEGIN_PROVIDER [integer, list_orb_pert_rdm, (n_orb_pert_rdm)]
+ implicit none
+ list_orb_pert_rdm = list_act
+
+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)]
+ implicit none
+ pert_2rdm_provider = 0.d0
+
+END_PROVIDER
+
+subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf, psi_det_connection, psi_coef_connection, n_det_connection)
+  use bitmasks
+  use selection_types
+  implicit none
+  
+  integer, intent(in)           :: n_det_connection
+  double precision, intent(in)  :: psi_coef_connection(n_det_connection,N_states)
+  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
+  integer :: sze_buff
+  sze_buff = 5 * mo_num ** 2 
+  allocate(keys(4,sze_buff),values(sze_buff))
+  nkeys = 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
+
+  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
+
+      do istate=1,N_states
+        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)
+        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)
+
+        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_connection,psi_coef_connection,n_det_connection,nkeys,keys,values,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
+end
+
+

src/cipsi/run_selection_slave.irp.f

@@ -61,7 +61,6 @@ subroutine run_selection_slave(thread,iproc,energy)
         ! Only first time
         bsize = min(N, (elec_alpha_num * (mo_num-elec_alpha_num))**2)
         call create_selection_buffer(bsize, bsize*2, buf)
-!        call create_selection_buffer(N, N*2, buf2)
         buffer_ready = .True.
       else
         ASSERT (N == buf%N)

src/cipsi/selection.irp.f

@@ -1,3 +1,4 @@
+
 use bitmasks
 
 BEGIN_PROVIDER [ double precision, pt2_match_weight, (N_states) ]
@@ -248,6 +249,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
   integer,allocatable               :: tmp_array(:)
   integer(bit_kind), allocatable :: minilist(:, :, :), fullminilist(:, :, :)
   logical, allocatable           :: banned(:,:,:), bannedOrb(:,:)
+  double precision, allocatable  :: coef_fullminilist(:,:)
 
 
   double precision, allocatable   :: mat(:,:,:)
@@ -546,6 +548,12 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
       allocate (fullminilist (N_int, 2, fullinteresting(0)), &
                     minilist (N_int, 2,     interesting(0)) )
+      if(pert_2rdm)then
+       allocate(coef_fullminilist(fullinteresting(0),N_states)) 
+       do i=1,fullinteresting(0)
+         coef_fullminilist(i,:) = psi_coef_sorted(fullinteresting(i),:)
+       enddo
+      endif
       do i=1,fullinteresting(0)
         fullminilist(1:N_int,1:2,i) = psi_det_sorted(1:N_int,1:2,fullinteresting(i))
       enddo
@@ -597,12 +605,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
             
             call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
             
-            call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf)
+            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 
+             call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2, variance, norm, mat, buf,fullminilist, coef_fullminilist, fullinteresting(0))
+            endif
           end if
         enddo
         if(s1 /= s2) monoBdo = .false.
       enddo
       deallocate(fullminilist,minilist)
+      if(pert_2rdm)then
+       deallocate(coef_fullminilist)
+      endif
     enddo
   enddo
   deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)
@@ -633,6 +648,10 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
   double precision :: E_shift
 
   logical, external :: detEq
+  double precision, allocatable :: values(:)
+  integer, allocatable          :: keys(:,:)
+  integer                       :: nkeys
+  
 
   if(sp == 3) then
     s1 = 1
@@ -745,7 +764,6 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
   end do
 end
 
-
 subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting)
   use bitmasks
   implicit none

src/cipsi/update_2rdm.irp.f

@@ -0,0 +1,80 @@
+use bitmasks
+
+subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection,n_det_connection,nkeys,keys,values,sze_buff)
+ implicit none
+ integer, intent(in) :: n_det_connection,sze_buff
+ double precision, intent(in) :: coef(N_states)
+ integer(bit_kind), intent(in) :: det(N_int,2)
+ integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) 
+ double precision, intent(in)  :: psi_coef_connection(n_det_connection, N_states)
+ integer,           intent(inout) :: keys(4,sze_buff),nkeys
+ double precision,  intent(inout) :: values(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) 
+  if(degree.gt.2)cycle
+  contrib = 0.d0
+  do j = 1, N_states
+   contrib += state_average_weight(j) * psi_coef_connection(i,j) * phase * coef(j)
+  enddo
+  ! case of single excitations 
+  if(degree == 1)then
+   if (nkeys+ 2 * elec_alpha_num .ge. sze_buff)then
+    call update_rdms(nkeys,keys,values,sze_buff)
+    nkeys = 0
+   endif
+   call update_buffer_single_exc_rdm(det,psi_det_connection(1,1,i),exc,phase,contrib,nkeys,keys,values,sze_buff)
+  else 
+  ! case of double excitations 
+   if (nkeys+ 4  .ge. sze_buff)then
+    call update_rdms(nkeys,keys,values,sze_buff)
+    nkeys = 0
+   endif
+   call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff)
+  endif
+ enddo
+
+end
+
+subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,sze_buff)
+ implicit none
+ integer, intent(in) :: sze_buff
+ integer(bit_kind), intent(in) :: det1(N_int,2)
+ integer(bit_kind), intent(in) :: det2(N_int,2)
+ integer,intent(in)             :: exc(0:2,2,2)
+ double precision,intent(in)    :: phase, contrib
+ integer, intent(inout)         :: nkeys, keys(4,sze_buff)
+ double precision, intent(inout):: values(sze_buff)
+ 
+
+
+end
+
+subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff)
+ implicit none
+ integer, intent(in) :: sze_buff
+ integer,intent(in)             :: exc(0:2,2,2)
+ double precision,intent(in)    :: phase, contrib
+ integer, intent(inout)         :: nkeys, keys(4,sze_buff)
+ double precision, intent(inout):: values(sze_buff)
+
+
+end
+
+
+subroutine update_rdms(nkeys,keys,values,sze_buff)
+ implicit none
+ integer, intent(in) :: nkeys, keys(4,sze_buff),sze_buff
+ double precision, intent(in) :: values(sze_buff)
+ integer :: i,h1,h2,p1,p2
+ do i = 1, nkeys
+  h1 = keys(1,i)
+  h2 = keys(2,i)
+  p1 = keys(3,i)
+  p2 = keys(4,i)
+  pert_2rdm_provider(h1,h2,p1,p2) += values(i)
+ enddo
+end

src/tools/print_wf.irp.f

@@ -51,7 +51,7 @@ subroutine routine
   if(degree == 0)then
    print*,'Reference determinant '
    call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,h00)
-  else
+  else if(degree .le. 2)then
    call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hii)
    call i_H_j(psi_det(1,1,1),psi_det(1,1,i),N_int,hij)
    delta_e = hii - h00

src/two_body_rdm/orb_range_2_rdm_openmp.irp.f

@@ -1,6 +1,4 @@
 
-
-
  BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_alpha_alpha_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
  implicit none
  double precision, allocatable :: state_weights(:) 

src/two_body_rdm/orb_range_routines_openmp.irp.f

@@ -76,6 +76,7 @@ end
  BEGIN_TEMPLATE
 subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
    use bitmasks
+   use omp_lib
    implicit none
    BEGIN_DOC
    ! Computes the two rdm for the N_st vectors |u_t>
@@ -92,6 +93,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
    integer, intent(in)            :: dim1,norb,list_orb(norb),ispin
    double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
    
+   integer(omp_lock_kind)          :: lock_2rdm
    integer                        :: i,j,k,l
    integer                        :: k_a, k_b, l_a, l_b
    integer                        :: krow, kcol
@@ -148,30 +150,31 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
    do i=1,maxab
      idx0(i) = i
    enddo
+   call omp_init_lock(lock_2rdm)
    
    ! Prepare the array of all alpha single excitations
    ! -------------------------------------------------
    
-   PROVIDE N_int nthreads_davidson
-   !!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson)      &
-       !    !$OMP   SHARED(psi_bilinear_matrix_rows, N_det,          &
-       !    !$OMP          psi_bilinear_matrix_columns,              &
-       !    !$OMP          psi_det_alpha_unique, psi_det_beta_unique,&
-       !    !$OMP          n_det_alpha_unique, n_det_beta_unique, N_int,&
-       !    !$OMP          psi_bilinear_matrix_transp_rows,          &
-       !    !$OMP          psi_bilinear_matrix_transp_columns,       &
-       !    !$OMP          psi_bilinear_matrix_transp_order, N_st,   &
-       !    !$OMP          psi_bilinear_matrix_order_transp_reverse, &
-       !    !$OMP          psi_bilinear_matrix_columns_loc,          &
-       !    !$OMP          psi_bilinear_matrix_transp_rows_loc,      &
-       !    !$OMP          istart, iend, istep, irp_here, v_t, s_t,  &
-       !    !$OMP          ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin)     &
-       !    !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,&
-       !    !$OMP          lcol, lrow, l_a, l_b,                     &
-       !    !$OMP          buffer, doubles, n_doubles,               &
-       !    !$OMP          tmp_det2, idx, l, kcol_prev,              &
-       !    !$OMP          singles_a, n_singles_a, singles_b,        &
-       !    !$OMP          n_singles_b, nkeys, keys, valus, c_average)
+   PROVIDE N_int nthreads_davidson elec_alpha_num 
+    !$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson)      &
+            !$OMP   SHARED(psi_bilinear_matrix_rows, N_det,lock_2rdm,&
+            !$OMP          psi_bilinear_matrix_columns,              &
+            !$OMP          psi_det_alpha_unique, psi_det_beta_unique,&
+            !$OMP          n_det_alpha_unique, n_det_beta_unique, N_int,&
+            !$OMP          psi_bilinear_matrix_transp_rows,          &
+            !$OMP          psi_bilinear_matrix_transp_columns,       &
+            !$OMP          psi_bilinear_matrix_transp_order, N_st,   &
+            !$OMP          psi_bilinear_matrix_order_transp_reverse, &
+            !$OMP          psi_bilinear_matrix_columns_loc,          &
+            !$OMP          psi_bilinear_matrix_transp_rows_loc,elec_alpha_num, &
+            !$OMP          istart, iend, istep, irp_here,list_orb_reverse, n_states, state_weights, dim1, & 
+            !$OMP          ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask)     &
+            !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1, c_2, &
+            !$OMP          lcol, lrow, l_a, l_b,                     &
+            !$OMP          buffer, doubles, n_doubles,               &
+            !$OMP          tmp_det2, idx, l, kcol_prev,              &
+            !$OMP          singles_a, n_singles_a, singles_b,        &
+            !$OMP          n_singles_b, nkeys, keys, values, c_average)
    
    ! Alpha/Beta double excitations
    ! =============================
@@ -189,7 +192,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
    ASSERT (istart > 0)
    ASSERT (istep  > 0)
    
-   !!$OMP DO SCHEDULE(dynamic,64)
+    !$OMP DO SCHEDULE(dynamic,64)
    do k_a=istart+ishift,iend,istep
      
      krow = psi_bilinear_matrix_rows(k_a)
@@ -257,13 +260,13 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
          if(alpha_beta)then
           ! only ONE contribution 
           if (nkeys+1 .ge. size(values)) then
-            call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
+            call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
             nkeys = 0
           endif
          else if (spin_trace)then
           ! TWO contributions 
           if (nkeys+2 .ge. size(values)) then
-            call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
+            call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
             nkeys = 0
           endif
          endif
@@ -275,9 +278,9 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
      enddo
      
    enddo
-   ! !$OMP END DO
+     !$OMP END DO
    
-   ! !$OMP DO SCHEDULE(dynamic,64)
+     !$OMP DO SCHEDULE(dynamic,64)
    do k_a=istart+ishift,iend,istep
      
      
@@ -345,14 +348,14 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
        enddo
        if(alpha_beta.or.spin_trace.or.alpha_alpha)then
        ! increment the alpha/beta part for single excitations
-       if (nkeys+ 2 * norb .ge. size(values)) then
-         call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
+       if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then
+         call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
          nkeys = 0
        endif
        call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
        ! increment the alpha/alpha part for single excitations
-       if (nkeys+4 * norb .ge. size(values)) then
-         call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
+       if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then
+         call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
          nkeys = 0
        endif
        call orb_range_off_diag_single_to_two_rdm_aa_dm_buffer(tmp_det,tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
@@ -378,8 +381,8 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
           c_2(l) = u_t(l,k_a)
           c_average += c_1(l) * c_2(l) * state_weights(l)
         enddo
-        if (nkeys+4 .ge. size(values)) then
-          call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
+        if (nkeys+4 .ge. sze_buff) then
+          call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
           nkeys = 0
         endif
         call orb_range_off_diag_double_to_two_rdm_aa_dm_buffer(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
@@ -449,14 +452,14 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
        enddo
        if(alpha_beta.or.spin_trace.or.beta_beta)then
         ! increment the alpha/beta  part for single excitations
-        if (nkeys+2 * norb .ge. size(values)) then
-          call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
+        if (nkeys+2 * elec_alpha_num .ge. sze_buff ) then
+          call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
           nkeys = 0
         endif
         call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
         ! increment the beta /beta  part for single excitations
-        if (nkeys+4 * norb .ge. size(values)) then
-          call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
+        if (nkeys+4 * elec_alpha_num .ge. sze_buff) then
+          call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
           nkeys = 0
         endif
         call orb_range_off_diag_single_to_two_rdm_bb_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
@@ -481,8 +484,8 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
           c_2(l) = u_t(l,k_a)
           c_average += c_1(l) * c_2(l) * state_weights(l)
         enddo
-        if (nkeys+4 .ge. size(values)) then
-          call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
+        if (nkeys+4 .ge. sze_buff) then
+          call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
           nkeys = 0
         endif
         call orb_range_off_diag_double_to_two_rdm_bb_dm_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
@@ -517,16 +520,16 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
        c_average += c_1(l) * c_1(l) * state_weights(l) 
      enddo
      
-     call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
+     call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
      nkeys = 0
      call orb_range_diag_to_all_two_rdm_dm_buffer(tmp_det,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
-     call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
+     call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
      nkeys = 0
      
    end do
-   !!$OMP END DO
+    !$OMP END DO
    deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values)
-   !!$OMP END PARALLEL
+    !$OMP END PARALLEL
    
 end
  
@@ -541,14 +544,17 @@ end
  END_TEMPLATE
  
 
-subroutine update_keys_values(keys,values,size_buff,nkeys,dim1,big_array)
+subroutine update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
+ use omp_lib
  implicit none
- integer, intent(in) :: size_buff,nkeys,dim1
- integer, intent(in) :: keys(4,size_buff)
- double precision, intent(in) :: values(size_buff)
+ integer, intent(in) :: nkeys,dim1
+ integer, intent(in) :: keys(4,nkeys)
+ double precision, intent(in) :: values(nkeys)
  double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
  
+ integer(omp_lock_kind),intent(inout):: lock_2rdm
  integer :: i,h1,h2,p1,p2
+ call omp_set_lock(lock_2rdm)
  do i = 1, nkeys
   h1 = keys(1,i)
   h2 = keys(2,i)
@@ -556,5 +562,7 @@ subroutine update_keys_values(keys,values,size_buff,nkeys,dim1,big_array)
   p2 = keys(4,i)
   big_array(h1,h2,p1,p2) += values(i)
  enddo
+ call omp_unset_lock(lock_2rdm)
 
 end
+