diff --git a/src/casscf/test_pert_2rdm.irp.f b/src/casscf/test_pert_2rdm.irp.f index 917d44ba..4894ec89 100644 --- a/src/casscf/test_pert_2rdm.irp.f +++ b/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 - print*,i,j,k,l - print*,pert_2rdm_provider(ii,jj,kk,ll) , integral, pert_2rdm_provider(ii,jj,kk,ll)* integral - endif + !if(dabs(pert_2rdm_provider(ii,jj,kk,ll) * integral).gt.1.d-12)then + ! print*,i,j,k,l + ! print*,pert_2rdm_provider(ii,jj,kk,ll) , integral, pert_2rdm_provider(ii,jj,kk,ll)* integral + !endif accu += pert_2rdm_provider(ii,jj,kk,ll) * integral enddo enddo @@ -51,15 +51,23 @@ 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 - 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) + !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-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) enddo diff --git a/src/cipsi/pert_rdm_providers.irp.f b/src/cipsi/pert_rdm_providers.irp.f index 6a8d7f63..e901b002 100644 --- a/src/cipsi/pert_rdm_providers.irp.f +++ b/src/cipsi/pert_rdm_providers.irp.f @@ -31,7 +31,7 @@ END_PROVIDER pert_2rdm_provider = 0.d0 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, & + 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_pert_2rdm @@ -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 - istate=1 + 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 @@ -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) -! 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 + end do + 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 (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 diff --git a/src/cipsi/routines_debug_pert_rdm.irp.f b/src/cipsi/routines_debug_pert_rdm.irp.f new file mode 100644 index 00000000..930fed5e --- /dev/null +++ b/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*,' = ',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 + + diff --git a/src/cipsi/update_2rdm.irp.f b/src/cipsi/update_2rdm.irp.f index a7a22dc1..c9385573 100644 --- a/src/cipsi/update_2rdm.irp.f +++ b/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(nkeys_1e + 2 .ge. sze_buff)then + if(degree == 1)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(degree == 1)then ! check for the length of buffer for the TWO BODY DM - if (nkeys + 6 * elec_alpha_num .ge. sze_buff)then + if (nkeys + 6 * elec_alpha_num .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 @@ -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 diff --git a/src/generators_cas/generators.irp.f b/src/generators_cas/generators.irp.f index b2f58202..7dc0d9fe 100644 --- a/src/generators_cas/generators.irp.f +++ b/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 diff --git a/src/two_body_rdm/orb_range_routines_openmp.irp.f b/src/two_body_rdm/orb_range_routines_openmp.irp.f index d36b4cc6..e15751ba 100644 --- a/src/two_body_rdm/orb_range_routines_openmp.irp.f +++ b/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)