diff --git a/src/bitmask/core_inact_act_virt.irp.f b/src/bitmask/core_inact_act_virt.irp.f index ff7ee2de..b016f1fd 100644 --- a/src/bitmask/core_inact_act_virt.irp.f +++ b/src/bitmask/core_inact_act_virt.irp.f @@ -322,6 +322,7 @@ END_PROVIDER enddo print *, 'Active MOs:' print *, list_act(1:n_act_orb) + print*, list_act_reverse(1:n_act_orb) END_PROVIDER diff --git a/src/two_body_rdm/orb_range_routines.irp.f b/src/two_body_rdm/orb_range_routines.irp.f index d5bd7d1c..d63a0390 100644 --- a/src/two_body_rdm/orb_range_routines.irp.f +++ b/src/two_body_rdm/orb_range_routines.irp.f @@ -30,6 +30,7 @@ subroutine orb_range_two_rdm_state_av(big_array,dim1,norb,list_orb,list_orb_reve u_t, & size(u_t, 1), & N_det, N_st) + call orb_range_two_rdm_state_av_work(big_array,dim1,norb,list_orb,list_orb_reverse,state_weights,ispin,u_t,N_st,sze,1,N_det,0,1) deallocate(u_t) @@ -135,6 +136,7 @@ subroutine orb_range_two_rdm_state_av_work_$N_int(big_array,dim1,norb,list_orb,l stop endif + PROVIDE N_int call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) 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 82649a23..b4ff7405 100644 --- a/src/two_body_rdm/orb_range_routines_openmp.irp.f +++ b/src/two_body_rdm/orb_range_routines_openmp.irp.f @@ -145,7 +145,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis PROVIDE N_int call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) - sze_buff = norb ** 3 + sze_buff = norb ** 3 + 6 * norb list_orb_reverse = -1000 do i = 1, norb list_orb_reverse(list_orb(i)) = i @@ -353,13 +353,17 @@ 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+norb .ge. size(values)) then - call update_keys_values(keys,values,size(values),nkeys,dim1,big_array) - nkeys = 0 - endif + if (nkeys+ 2 * norb .ge. size(values)) then + call update_keys_values(keys,values,size(values),nkeys,dim1,big_array) + 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 -!!!! call orb_range_off_diagonal_single_to_two_rdm_aa_dm(tmp_det,tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) + if (nkeys+2 * norb .ge. size(values)) then + call update_keys_values(keys,values,size(values),nkeys,dim1,big_array) + 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) endif enddo @@ -382,7 +386,11 @@ 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 -!!!! call orb_range_off_diagonal_double_to_two_rdm_aa_dm(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) + if (nkeys+4 .ge. size(values)) then + call update_keys_values(keys,values,size(values),nkeys,dim1,big_array) + 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) enddo endif @@ -455,7 +463,11 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis 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 -!!!! call orb_range_off_diagonal_single_to_two_rdm_bb_dm(tmp_det, tmp_det2,c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) + if (nkeys+norb .ge. size(values)) then + call update_keys_values(keys,values,size(values),nkeys,dim1,big_array) + 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) endif enddo @@ -477,7 +489,12 @@ 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 -!!!! call orb_range_off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) +! call orb_range_off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_average,big_array,dim1,orb_bitmask,list_orb_reverse,ispin) + if (nkeys+4 .ge. size(values)) then + call update_keys_values(keys,values,size(values),nkeys,dim1,big_array) + nkeys = 0 + endif + call orb_range_off_diag_double_to_two_rdm_bb_dm_buffer(tmp_det(1,2),psi_det_alpha_unique(1, lcol),c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values) ASSERT (l_a <= N_det) enddo diff --git a/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f index a3c7a76d..52cccbf3 100644 --- a/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f +++ b/src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f @@ -13,7 +13,7 @@ double precision, intent(in) :: c_1 integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) - integer :: i,j,h1,h2,istate + integer :: i,j,h1,h2 call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) do i = 1, n_occ_ab(1) h1 = occ(i,1) @@ -53,7 +53,7 @@ integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) - integer :: i,j,h1,h2,istate + integer :: i,j,h1,h2 integer(bit_kind) :: det_1_act(N_int,2) logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace do i = 1, N_int @@ -193,7 +193,7 @@ integer(bit_kind), intent(in) :: orb_bitmask(N_int) integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(in) :: c_1 - integer :: i,j,h1,h2,p1,p2,istate + integer :: i,j,h1,h2,p1,p2 integer :: exc(0:2,2,2) double precision :: phase logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace @@ -278,7 +278,7 @@ integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) - integer :: i,j,h1,h2,istate,p1 + integer :: i,j,h1,h2,p1 integer :: exc(0:2,2,2) double precision :: phase @@ -397,7 +397,7 @@ integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) - integer :: i,j,h1,h2,istate,p1 + integer :: i,j,h1,h2,p1 integer :: exc(0:2,2,2) double precision :: phase @@ -477,7 +477,7 @@ integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) - integer :: i,j,h1,h2,istate,p1 + integer :: i,j,h1,h2,p1 integer :: exc(0:2,2,2) double precision :: phase logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace @@ -510,18 +510,16 @@ p1 = exc(1,2,2) if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return p1 = list_orb_reverse(p1) - do istate = 1, N_states - do i = 1, n_occ_ab(2) - h2 = occ(i,2) - if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle - h2 = list_orb_reverse(h2) - big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase - big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase + do i = 1, n_occ_ab(2) + h2 = occ(i,2) + if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle + h2 = list_orb_reverse(h2) + big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase + big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase - big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase - big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase - enddo - enddo + big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase + big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase + enddo endif endif end @@ -557,7 +555,7 @@ integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(in) :: c_1 - integer :: i,j,h1,h2,p1,p2,istate + integer :: i,j,h1,h2,p1,p2 integer :: exc(0:2,2) double precision :: phase @@ -590,13 +588,11 @@ if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return p2 = list_orb_reverse(p2) if(alpha_alpha.or.spin_trace)then - do istate = 1, N_states big_array(h1,h2,p1,p2) += 0.5d0 * c_1 * phase big_array(h1,h2,p2,p1) -= 0.5d0 * c_1 * phase big_array(h2,h1,p2,p1) += 0.5d0 * c_1 * phase big_array(h2,h1,p1,p2) -= 0.5d0 * c_1 * phase - enddo endif end @@ -631,7 +627,7 @@ integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(in) :: c_1 - integer :: i,j,h1,h2,p1,p2,istate + integer :: i,j,h1,h2,p1,p2 integer :: exc(0:2,2) double precision :: phase logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace diff --git a/src/two_body_rdm/routines_compute_2rdm_orb_range_openmp.irp.f b/src/two_body_rdm/routines_compute_2rdm_orb_range_openmp.irp.f index b5bc66b4..ffbb2711 100644 --- a/src/two_body_rdm/routines_compute_2rdm_orb_range_openmp.irp.f +++ b/src/two_body_rdm/routines_compute_2rdm_orb_range_openmp.irp.f @@ -26,7 +26,7 @@ integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) - integer :: i,j,h1,h2,istate + integer :: i,j,h1,h2 integer(bit_kind) :: det_1_act(N_int,2) logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace do i = 1, N_int @@ -201,7 +201,7 @@ double precision, intent(out) :: values(sze_buff) integer , intent(out) :: keys(4,sze_buff) integer , intent(inout):: nkeys - integer :: i,j,h1,h2,p1,p2,istate + integer :: i,j,h1,h2,p1,p2 integer :: exc(0:2,2,2) double precision :: phase logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace @@ -288,7 +288,7 @@ integer :: occ(N_int*bit_kind_size,2) integer :: n_occ_ab(2) - integer :: i,j,h1,h2,istate,p1 + integer :: i,j,h1,h2,p1 integer :: exc(0:2,2,2) double precision :: phase @@ -409,310 +409,399 @@ endif end -! subroutine orb_range_off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) -! BEGIN_DOC -!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for -!! -!! a given couple of determinant det_1, det_2 being a ALPHA SINGLE excitation with respect to one another -!! -!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 -!! -!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation -!! -!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals -!! -!! ispin determines which spin-spin component of the two-rdm you will update -!! -!! ispin == 1 :: alpha/ alpha -!! ispin == 2 :: beta / beta -!! ispin == 3 :: alpha/ beta -!! ispin == 4 :: spin traced <=> total two-rdm -!! -!! here, only ispin == 1 or 4 will do something -! END_DOC -! use bitmasks -! implicit none -! integer, intent(in) :: dim1,ispin -! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) -! integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) -! integer(bit_kind), intent(in) :: orb_bitmask(N_int) -! integer, intent(in) :: list_orb_reverse(mo_num) -! double precision, intent(in) :: c_1 -! -! integer :: occ(N_int*bit_kind_size,2) -! integer :: n_occ_ab(2) -! integer :: i,j,h1,h2,istate,p1 -! integer :: exc(0:2,2,2) -! double precision :: phase -! -! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace -! logical :: is_integer_in_string -! alpha_alpha = .False. -! beta_beta = .False. -! alpha_beta = .False. -! spin_trace = .False. -! if( ispin == 1)then -! alpha_alpha = .True. -! else if(ispin == 2)then -! beta_beta = .True. -! else if(ispin == 3)then -! alpha_beta = .True. -! else if(ispin == 4)then -! spin_trace = .True. -! endif -! -! call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) -! call get_single_excitation(det_1,det_2,exc,phase,N_int) -! if(alpha_alpha.or.spin_trace)then -! if (exc(0,1,1) == 1) then -! ! Mono alpha -! h1 = exc(1,1,1) -! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return -! h1 = list_orb_reverse(h1) -! p1 = exc(1,2,1) -! 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(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle -! h2 = list_orb_reverse(h2) -! big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase -! big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase -! -! big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase -! big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase -! enddo -! else -! return -! endif -! 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) + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a ALPHA SINGLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm + ! + ! here, only ispin == 1 or 4 will do something + END_DOC + use bitmasks + implicit none + 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) + double precision, intent(in) :: c_1 + double precision, intent(out) :: values(sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,p1 + integer :: exc(0:2,2,2) + double precision :: phase + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + + call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) + call get_single_excitation(det_1,det_2,exc,phase,N_int) + if(alpha_alpha.or.spin_trace)then + if (exc(0,1,1) == 1) then + ! Mono alpha + h1 = exc(1,1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,1) + if(list_orb_reverse(p1).lt.0)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 + h2 = list_orb_reverse(h2) -! subroutine orb_range_off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) -! use bitmasks -! BEGIN_DOC -!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for -!! -!! a given couple of determinant det_1, det_2 being a BETA SINGLE excitation with respect to one another -!! -!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 -!! -!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation -!! -!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals -!! -!! ispin determines which spin-spin component of the two-rdm you will update -!! -!! ispin == 1 :: alpha/ alpha -!! ispin == 2 :: beta / beta -!! ispin == 3 :: alpha/ beta -!! ispin == 4 :: spin traced <=> total two-rdm -!! -!! here, only ispin == 2 or 4 will do something -! END_DOC -! implicit none -! integer, intent(in) :: dim1,ispin -! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) -! integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) -! integer(bit_kind), intent(in) :: orb_bitmask(N_int) -! integer, intent(in) :: list_orb_reverse(mo_num) -! double precision, intent(in) :: c_1 -! -! -! integer :: occ(N_int*bit_kind_size,2) -! integer :: n_occ_ab(2) -! integer :: i,j,h1,h2,istate,p1 -! integer :: exc(0:2,2,2) -! double precision :: phase -! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace -! logical :: is_integer_in_string -! alpha_alpha = .False. -! beta_beta = .False. -! alpha_beta = .False. -! spin_trace = .False. -! if( ispin == 1)then -! alpha_alpha = .True. -! else if(ispin == 2)then -! beta_beta = .True. -! else if(ispin == 3)then -! alpha_beta = .True. -! else if(ispin == 4)then -! spin_trace = .True. -! endif -! -! -! call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) -! call get_single_excitation(det_1,det_2,exc,phase,N_int) -! if(beta_beta.or.spin_trace)then -! if (exc(0,1,1) == 1) then -! return -! else -! ! Mono beta -! h1 = exc(1,1,2) -! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return -! h1 = list_orb_reverse(h1) -! p1 = exc(1,2,2) -! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return -! p1 = list_orb_reverse(p1) -! do istate = 1, N_states -! do i = 1, n_occ_ab(2) -! h2 = occ(i,2) -! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle -! h2 = list_orb_reverse(h2) -! big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase -! big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase -! -! big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase -! big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase -! enddo -! enddo -! endif -! endif -! end + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + enddo + else + return + endif + 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) + use bitmasks + BEGIN_DOC + ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for + ! + ! a given couple of determinant det_1, det_2 being a BETA SINGLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm + ! + ! here, only ispin == 2 or 4 will do something + END_DOC + implicit none + 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) + double precision, intent(in) :: c_1 + double precision, intent(out) :: values(sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab(2) + integer :: i,j,h1,h2,p1 + integer :: exc(0:2,2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + + + call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int) + call get_single_excitation(det_1,det_2,exc,phase,N_int) + if(beta_beta.or.spin_trace)then + if (exc(0,1,1) == 1) then + return + else + ! Mono beta + h1 = exc(1,1,2) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + p1 = exc(1,2,2) + if(list_orb_reverse(p1).lt.0)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 + h2 = list_orb_reverse(h2) + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = h2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = h2 + enddo + endif + endif + end -! subroutine orb_range_off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) -! use bitmasks -! BEGIN_DOC -!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for -!! -!! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another -!! -!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 -!! -!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation -!! -!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals -!! -!! ispin determines which spin-spin component of the two-rdm you will update -!! -!! ispin == 1 :: alpha/ alpha -!! ispin == 2 :: beta / beta -!! ispin == 3 :: alpha/ beta -!! ispin == 4 :: spin traced <=> total two-rdm -!! -!! here, only ispin == 1 or 4 will do something -! END_DOC -! implicit none -! integer, intent(in) :: dim1,ispin -! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) -! integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) -! integer(bit_kind), intent(in) :: orb_bitmask(N_int) -! integer, intent(in) :: list_orb_reverse(mo_num) -! double precision, intent(in) :: c_1 -! -! integer :: i,j,h1,h2,p1,p2,istate -! integer :: exc(0:2,2) -! double precision :: phase -! -! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace -! logical :: is_integer_in_string -! alpha_alpha = .False. -! beta_beta = .False. -! alpha_beta = .False. -! spin_trace = .False. -! if( ispin == 1)then -! alpha_alpha = .True. -! else if(ispin == 2)then -! beta_beta = .True. -! else if(ispin == 3)then -! alpha_beta = .True. -! else if(ispin == 4)then -! spin_trace = .True. -! endif -! call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) -! h1 =exc(1,1) -! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return -! h1 = list_orb_reverse(h1) -! h2 =exc(2,1) -! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return -! h2 = list_orb_reverse(h2) -! p1 =exc(1,2) -! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return -! p1 = list_orb_reverse(p1) -! p2 =exc(2,2) -! if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return -! p2 = list_orb_reverse(p2) -! if(alpha_alpha.or.spin_trace)then -! do istate = 1, N_states -! big_array(h1,h2,p1,p2) += 0.5d0 * c_1 * phase -! big_array(h1,h2,p2,p1) -= 0.5d0 * c_1 * phase -! -! big_array(h2,h1,p2,p1) += 0.5d0 * c_1 * phase -! big_array(h2,h1,p1,p2) -= 0.5d0 * c_1 * phase -! enddo -! endif -! end + subroutine orb_range_off_diag_double_to_two_rdm_aa_dm_buffer(det_1,det_2,c_1,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 + ! + ! a given couple of determinant det_1, det_2 being a ALPHA/ALPHA DOUBLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm + ! + ! here, only ispin == 1 or 4 will do something + END_DOC + implicit none + 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) + double precision, intent(in) :: c_1 + double precision, intent(out) :: values(sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + + integer :: i,j,h1,h2,p1,p2 + integer :: exc(0:2,2) + double precision :: phase + + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) + h1 =exc(1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + h2 =exc(2,1) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + p1 =exc(1,2) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + p2 =exc(2,2) + if(list_orb_reverse(p2).lt.0)return + p2 = list_orb_reverse(p2) + if(alpha_alpha.or.spin_trace)then + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 -! subroutine orb_range_off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) -! use bitmasks -! BEGIN_DOC -!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for -!! -!! a given couple of determinant det_1, det_2 being a BETA /BETA DOUBLE excitation with respect to one another -!! -!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 -!! -!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation -!! -!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals -!! -!! ispin determines which spin-spin component of the two-rdm you will update -!! -!! ispin == 1 :: alpha/ alpha -!! ispin == 2 :: beta / beta -!! ispin == 3 :: alpha/ beta -!! ispin == 4 :: spin traced <=> total two-rdm -!! -!! here, only ispin == 2 or 4 will do something -! END_DOC -! implicit none -! -! integer, intent(in) :: dim1,ispin -! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) -! integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) -! integer(bit_kind), intent(in) :: orb_bitmask(N_int) -! integer, intent(in) :: list_orb_reverse(mo_num) -! double precision, intent(in) :: c_1 -! -! integer :: i,j,h1,h2,p1,p2,istate -! integer :: exc(0:2,2) -! double precision :: phase -! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace -! logical :: is_integer_in_string -! alpha_alpha = .False. -! beta_beta = .False. -! alpha_beta = .False. -! spin_trace = .False. -! if( ispin == 1)then -! alpha_alpha = .True. -! else if(ispin == 2)then -! beta_beta = .True. -! else if(ispin == 3)then -! alpha_beta = .True. -! else if(ispin == 4)then -! spin_trace = .True. -! endif -! -! call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) -! h1 =exc(1,1) -! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return -! h1 = list_orb_reverse(h1) -! h2 =exc(2,1) -! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return -! h2 = list_orb_reverse(h2) -! p1 =exc(1,2) -! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return -! p1 = list_orb_reverse(p1) -! p2 =exc(2,2) -! if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return -! p2 = list_orb_reverse(p2) -! if(beta_beta.or.spin_trace)then -! big_array(h1,h2,p1,p2) += 0.5d0 * c_1* phase -! big_array(h1,h2,p2,p1) -= 0.5d0 * c_1* phase -! -! big_array(h2,h1,p2,p1) += 0.5d0 * c_1* phase -! big_array(h2,h1,p1,p2) -= 0.5d0 * c_1* phase -! endif -! end + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + endif + end + + subroutine orb_range_off_diag_double_to_two_rdm_bb_dm_buffer(det_1,det_2,c_1,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 + ! + ! a given couple of determinant det_1, det_2 being a BETA /BETA DOUBLE excitation with respect to one another + ! + ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 + ! + ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation + ! + ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals + ! + ! ispin determines which spin-spin component of the two-rdm you will update + ! + ! ispin == 1 :: alpha/ alpha + ! ispin == 2 :: beta / beta + ! ispin == 3 :: alpha/ beta + ! ispin == 4 :: spin traced <=> total two-rdm + ! + ! here, only ispin == 2 or 4 will do something + END_DOC + implicit none + + 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) + double precision, intent(in) :: c_1 + double precision, intent(out) :: values(sze_buff) + integer , intent(out) :: keys(4,sze_buff) + integer , intent(inout):: nkeys + + integer :: i,j,h1,h2,p1,p2 + integer :: exc(0:2,2) + double precision :: phase + logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace + logical :: is_integer_in_string + alpha_alpha = .False. + beta_beta = .False. + alpha_beta = .False. + spin_trace = .False. + if( ispin == 1)then + alpha_alpha = .True. + else if(ispin == 2)then + beta_beta = .True. + else if(ispin == 3)then + alpha_beta = .True. + else if(ispin == 4)then + spin_trace = .True. + endif + + call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) + h1 =exc(1,1) + if(list_orb_reverse(h1).lt.0)return + h1 = list_orb_reverse(h1) + h2 =exc(2,1) + if(list_orb_reverse(h2).lt.0)return + h2 = list_orb_reverse(h2) + p1 =exc(1,2) + if(list_orb_reverse(p1).lt.0)return + p1 = list_orb_reverse(p1) + p2 =exc(2,2) + if(list_orb_reverse(p2).lt.0)return + p2 = list_orb_reverse(p2) + if(beta_beta.or.spin_trace)then + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h1 + keys(2,nkeys) = h2 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p2 + keys(4,nkeys) = p1 + + nkeys += 1 + values(nkeys) = - 0.5d0 * c_1 * phase + keys(1,nkeys) = h2 + keys(2,nkeys) = h1 + keys(3,nkeys) = p1 + keys(4,nkeys) = p2 + endif + end