diff --git a/src/two_body_rdm/all_states_act_2_rdm_update_routines.irp.f b/src/two_body_rdm/all_states_act_2_rdm_update_routines.irp.f index a42f2d79..3e4a070c 100644 --- a/src/two_body_rdm/all_states_act_2_rdm_update_routines.irp.f +++ b/src/two_body_rdm/all_states_act_2_rdm_update_routines.irp.f @@ -73,7 +73,6 @@ else if(ispin == 4)then spin_trace = .True. endif -! call debug_det(det_1_act,N_int) call bitstring_to_list_ab(det_1_act, occ, n_occ_ab, N_int) logical :: is_integer_in_string integer :: i1,i2 diff --git a/src/two_body_rdm/compute_orb_range_omp.irp.f b/src/two_body_rdm/compute_orb_range_omp.irp.f index 0ba934d7..9d0f3fe8 100644 --- a/src/two_body_rdm/compute_orb_range_omp.irp.f +++ b/src/two_body_rdm/compute_orb_range_omp.irp.f @@ -57,6 +57,8 @@ i2 = occ(j,2) h1 = list_orb_reverse(i1) h2 = list_orb_reverse(i2) + ! If alpha/beta, electron 1 is alpha, electron 2 is beta + ! Therefore you don't necessayr have symmetry between electron 1 and 2 nkeys += 1 values(nkeys) = c_1 keys(1,nkeys) = h1 @@ -255,7 +257,7 @@ endif end - subroutine orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + subroutine orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(det_1,det_2,c_1,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) use bitmasks BEGIN_DOC ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for @@ -281,6 +283,7 @@ integer, intent(in) :: ispin,sze_buff integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) double precision, intent(in) :: c_1 double precision, intent(out) :: values(sze_buff) integer , intent(out) :: keys(4,sze_buff) @@ -314,14 +317,14 @@ if (exc(0,1,1) == 1) then ! Mono alpha h1 = exc(1,1,1) - if(list_orb_reverse(h1).lt.0)return + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return h1 = list_orb_reverse(h1) p1 = exc(1,2,1) - if(list_orb_reverse(p1).lt.0)return + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return p1 = list_orb_reverse(p1) do i = 1, n_occ_ab(2) h2 = occ(i,2) - if(list_orb_reverse(h2).lt.0)return + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle h2 = list_orb_reverse(h2) nkeys += 1 values(nkeys) = c_1 * phase @@ -333,14 +336,14 @@ else ! Mono beta h1 = exc(1,1,2) - if(list_orb_reverse(h1).lt.0)return + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return h1 = list_orb_reverse(h1) p1 = exc(1,2,2) - if(list_orb_reverse(p1).lt.0)return + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return p1 = list_orb_reverse(p1) do i = 1, n_occ_ab(1) h2 = occ(i,1) - if(list_orb_reverse(h2).lt.0)return + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle h2 = list_orb_reverse(h2) nkeys += 1 values(nkeys) = c_1 * phase @@ -354,14 +357,14 @@ if (exc(0,1,1) == 1) then ! Mono alpha h1 = exc(1,1,1) - if(list_orb_reverse(h1).lt.0)return + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return h1 = list_orb_reverse(h1) p1 = exc(1,2,1) - if(list_orb_reverse(p1).lt.0)return + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return p1 = list_orb_reverse(p1) do i = 1, n_occ_ab(2) h2 = occ(i,2) - if(list_orb_reverse(h2).lt.0)return + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle h2 = list_orb_reverse(h2) nkeys += 1 values(nkeys) = 0.5d0 * c_1 * phase @@ -379,19 +382,15 @@ else ! Mono beta h1 = exc(1,1,2) - if(list_orb_reverse(h1).lt.0)return + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return h1 = list_orb_reverse(h1) p1 = exc(1,2,2) - if(list_orb_reverse(p1).lt.0)return + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return p1 = list_orb_reverse(p1) - !print*,'****************' - !print*,'****************' - !print*,'h1,p1',h1,p1 do i = 1, n_occ_ab(1) h2 = occ(i,1) - if(list_orb_reverse(h2).lt.0)return + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle h2 = list_orb_reverse(h2) - ! print*,'h2 = ',h2 nkeys += 1 values(nkeys) = 0.5d0 * c_1 * phase keys(1,nkeys) = h1 @@ -409,7 +408,7 @@ endif end - subroutine orb_range_off_diag_single_to_two_rdm_aa_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + subroutine orb_range_off_diag_single_to_two_rdm_aa_dm_buffer(det_1,det_2,c_1,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) BEGIN_DOC ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for ! @@ -435,6 +434,7 @@ integer, intent(in) :: ispin,sze_buff integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) double precision, intent(in) :: c_1 double precision, intent(out) :: values(sze_buff) integer , intent(out) :: keys(4,sze_buff) @@ -468,14 +468,14 @@ if (exc(0,1,1) == 1) then ! Mono alpha h1 = exc(1,1,1) - if(list_orb_reverse(h1).lt.0)return + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return h1 = list_orb_reverse(h1) p1 = exc(1,2,1) - if(list_orb_reverse(p1).lt.0)return + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return p1 = list_orb_reverse(p1) do i = 1, n_occ_ab(1) h2 = occ(i,1) - if(list_orb_reverse(h2).lt.0)return + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle h2 = list_orb_reverse(h2) nkeys += 1 @@ -512,7 +512,7 @@ endif end - subroutine orb_range_off_diag_single_to_two_rdm_bb_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) + subroutine orb_range_off_diag_single_to_two_rdm_bb_dm_buffer(det_1,det_2,c_1,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) use bitmasks BEGIN_DOC ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for @@ -538,6 +538,7 @@ integer, intent(in) :: ispin,sze_buff integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) integer, intent(in) :: list_orb_reverse(mo_num) + integer(bit_kind), intent(in) :: orb_bitmask(N_int) double precision, intent(in) :: c_1 double precision, intent(out) :: values(sze_buff) integer , intent(out) :: keys(4,sze_buff) @@ -573,14 +574,14 @@ else ! Mono beta h1 = exc(1,1,2) - if(list_orb_reverse(h1).lt.0)return + if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return h1 = list_orb_reverse(h1) p1 = exc(1,2,2) - if(list_orb_reverse(p1).lt.0)return + if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return p1 = list_orb_reverse(p1) do i = 1, n_occ_ab(2) h2 = occ(i,2) - if(list_orb_reverse(h2).lt.0)return + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle h2 = list_orb_reverse(h2) nkeys += 1 values(nkeys) = 0.5d0 * c_1 * phase diff --git a/src/two_body_rdm/orb_range_routines_omp.irp.f b/src/two_body_rdm/orb_range_routines_omp.irp.f index b6e59540..bb195454 100644 --- a/src/two_body_rdm/orb_range_routines_omp.irp.f +++ b/src/two_body_rdm/orb_range_routines_omp.irp.f @@ -271,6 +271,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis endif endif call orb_range_off_diag_double_to_two_rdm_ab_dm_buffer(tmp_det,tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) +! print*,'todo orb_range_off_diag_double_to_two_rdm_ab_dm_buffer' enddo endif @@ -352,13 +353,13 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis 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) + call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) ! increment the alpha/alpha part for single excitations 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) + call orb_range_off_diag_single_to_two_rdm_aa_dm_buffer(tmp_det,tmp_det2,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) endif enddo @@ -456,13 +457,13 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis 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) + call orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(tmp_det, tmp_det2,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) ! increment the beta /beta part for single excitations 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) + call orb_range_off_diag_single_to_two_rdm_bb_dm_buffer(tmp_det, tmp_det2,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) endif enddo @@ -489,6 +490,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis 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) +! print*,'to do orb_range_off_diag_double_to_two_rdm_bb_dm_buffer' ASSERT (l_a <= N_det) enddo diff --git a/src/two_body_rdm/test_2_rdm.irp.f b/src/two_body_rdm/test_2_rdm.irp.f index e993da24..27a25024 100644 --- a/src/two_body_rdm/test_2_rdm.irp.f +++ b/src/two_body_rdm/test_2_rdm.irp.f @@ -2,7 +2,7 @@ program test_2_rdm implicit none read_wf = .True. touch read_wf - call routine_full_mos +! call routine_full_mos call routine_active_only end @@ -16,6 +16,12 @@ subroutine routine_active_only double precision :: vijkl,rdmaa,get_two_e_integral,rdmab,rdmbb,rdmtot double precision :: accu_aa(N_states),accu_bb(N_states),accu_ab(N_states),accu_tot(N_states) + double precision :: accu_ab_omp,rdmab_omp + double precision :: accu_bb_omp,rdmbb_omp + double precision :: accu_aa_omp,rdmaa_omp + accu_ab_omp = 0.d0 + accu_bb_omp = 0.d0 + accu_aa_omp = 0.d0 accu_aa = 0.d0 accu_ab = 0.d0 accu_bb = 0.d0 @@ -34,11 +40,17 @@ subroutine routine_active_only vijkl = get_two_e_integral(lorb,korb,jorb,iorb,mo_integrals_map) + rdmab_omp = state_av_act_two_rdm_openmp_alpha_beta_mo(l,k,j,i) + rdmbb_omp = state_av_act_two_rdm_openmp_beta_beta_mo(l,k,j,i) + rdmaa_omp = state_av_act_two_rdm_openmp_alpha_alpha_mo(l,k,j,i) rdmaa = all_states_act_two_rdm_alpha_alpha_mo(l,k,j,i,istate) rdmbb = all_states_act_two_rdm_beta_beta_mo(l,k,j,i,istate) rdmab = all_states_act_two_rdm_alpha_beta_mo(l,k,j,i,istate) rdmtot = all_states_act_two_rdm_spin_trace_mo(l,k,j,i,istate) + accu_ab_omp += vijkl * rdmab_omp + accu_bb_omp += vijkl * rdmbb_omp + accu_aa_omp += vijkl * rdmaa_omp accu_ab(istate) += vijkl * rdmab accu_aa(istate) += vijkl * rdmaa accu_bb(istate) += vijkl * rdmbb @@ -50,8 +62,11 @@ subroutine routine_active_only print*,'' print*,'Active space only energy ' print*,'accu_aa(istate) = ',accu_aa(istate) + print*,'accu_aa_omp = ',accu_aa_omp print*,'accu_bb(istate) = ',accu_bb(istate) + print*,'accu_bb_omp = ',accu_bb_omp print*,'accu_ab(istate) = ',accu_ab(istate) + print*,'accu_ab_omp = ',accu_ab_omp print*,'' print*,'sum (istate) = ',accu_aa(istate) + accu_bb(istate) + accu_ab(istate) print*,'accu_tot(istate) = ',accu_tot(istate)