LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-26 22:42:06 +02:00
Code Releases Activity

two rdm seems to work with buffer, ready for openmp

Browse Source
This commit is contained in:
Emmanuel Giner 2019-07-04 17:34:56 +02:00
parent 59aaf3806d
commit 887afe97b4
5 changed files with 440 additions and 335 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
src/bitmask/core_inact_act_virt.irp.f
View File

@ -322,6 +322,7 @@ END_PROVIDER
enddo enddo
print *, 'Active MOs:' print *, 'Active MOs:'
print *, list_act(1:n_act_orb) print *, list_act(1:n_act_orb)
print*, list_act_reverse(1:n_act_orb)
END_PROVIDER END_PROVIDER

2
src/two_body_rdm/orb_range_routines.irp.f
View File

@ -30,6 +30,7 @@ subroutine orb_range_two_rdm_state_av(big_array,dim1,norb,list_orb,list_orb_reve
u_t, & u_t, &
size(u_t, 1), & size(u_t, 1), &
N_det, N_st) 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) 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) 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 stop
endif endif
PROVIDE N_int PROVIDE N_int
call list_to_bitstring( orb_bitmask, list_orb, norb, N_int) call list_to_bitstring( orb_bitmask, list_orb, norb, N_int)

35
src/two_body_rdm/orb_range_routines_openmp.irp.f
View File

@ -145,7 +145,7 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
PROVIDE N_int PROVIDE N_int
call list_to_bitstring( orb_bitmask, list_orb, norb, 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 list_orb_reverse = -1000
do i = 1, norb do i = 1, norb
list_orb_reverse(list_orb(i)) = i 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 enddo
if(alpha_beta.or.spin_trace.or.alpha_alpha)then if(alpha_beta.or.spin_trace.or.alpha_alpha)then
! increment the alpha/beta part for single excitations ! increment the alpha/beta part for single excitations
if (nkeys+norb .ge. size(values)) then if (nkeys+ 2 * norb .ge. size(values)) then
call update_keys_values(keys,values,size(values),nkeys,dim1,big_array) call update_keys_values(keys,values,size(values),nkeys,dim1,big_array)
nkeys = 0 nkeys = 0
endif 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,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
! increment the alpha/alpha part for single excitations ! 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 endif
enddo 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_2(l) = u_t(l,k_a)
c_average += c_1(l) * c_2(l) * state_weights(l) c_average += c_1(l) * c_2(l) * state_weights(l)
enddo 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 enddo
endif endif
@ -455,7 +463,11 @@ subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,lis
endif 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,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
! increment the beta /beta part for single excitations ! 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 endif
enddo 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_2(l) = u_t(l,k_a)
c_average += c_1(l) * c_2(l) * state_weights(l) c_average += c_1(l) * c_2(l) * state_weights(l)
enddo 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) ASSERT (l_a <= N_det)
enddo enddo

38
src/two_body_rdm/routines_compute_2rdm_orb_range.irp.f
View File

@ -13,7 +13,7 @@
double precision, intent(in) :: c_1 double precision, intent(in) :: c_1
integer :: occ(N_int*bit_kind_size,2) integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(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) call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
do i = 1, n_occ_ab(1) do i = 1, n_occ_ab(1)
h1 = occ(i,1) h1 = occ(i,1)
@ -53,7 +53,7 @@
integer :: occ(N_int*bit_kind_size,2) integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(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) integer(bit_kind) :: det_1_act(N_int,2)
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
do i = 1, N_int do i = 1, N_int
@ -193,7 +193,7 @@
integer(bit_kind), intent(in) :: orb_bitmask(N_int) integer(bit_kind), intent(in) :: orb_bitmask(N_int)
integer, intent(in) :: list_orb_reverse(mo_num) integer, intent(in) :: list_orb_reverse(mo_num)
double precision, intent(in) :: c_1 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) integer :: exc(0:2,2,2)
double precision :: phase double precision :: phase
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
@ -278,7 +278,7 @@
integer :: occ(N_int*bit_kind_size,2) integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(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) integer :: exc(0:2,2,2)
double precision :: phase double precision :: phase
@ -397,7 +397,7 @@
integer :: occ(N_int*bit_kind_size,2) integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(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) integer :: exc(0:2,2,2)
double precision :: phase double precision :: phase
@ -477,7 +477,7 @@
integer :: occ(N_int*bit_kind_size,2) integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(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) integer :: exc(0:2,2,2)
double precision :: phase double precision :: phase
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
@ -510,18 +510,16 @@
p1 = exc(1,2,2) p1 = exc(1,2,2)
if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
p1 = list_orb_reverse(p1) p1 = list_orb_reverse(p1)
do istate = 1, N_states do i = 1, n_occ_ab(2)
do i = 1, n_occ_ab(2) h2 = occ(i,2)
h2 = occ(i,2) if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle h2 = list_orb_reverse(h2)
h2 = list_orb_reverse(h2) big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase
big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase big_array(h1,h2,h2,p1) -= 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,h2,p1) += 0.5d0 * c_1 * phase
big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase
enddo enddo
enddo
endif endif
endif endif
end end
@ -557,7 +555,7 @@
integer, intent(in) :: list_orb_reverse(mo_num) integer, intent(in) :: list_orb_reverse(mo_num)
double precision, intent(in) :: c_1 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) integer :: exc(0:2,2)
double precision :: phase double precision :: phase
@ -590,13 +588,11 @@
if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return
p2 = list_orb_reverse(p2) p2 = list_orb_reverse(p2)
if(alpha_alpha.or.spin_trace)then 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,p1,p2) += 0.5d0 * c_1 * phase
big_array(h1,h2,p2,p1) -= 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,p2,p1) += 0.5d0 * c_1 * phase
big_array(h2,h1,p1,p2) -= 0.5d0 * c_1 * phase big_array(h2,h1,p1,p2) -= 0.5d0 * c_1 * phase
enddo
endif endif
end end
@ -631,7 +627,7 @@
integer, intent(in) :: list_orb_reverse(mo_num) integer, intent(in) :: list_orb_reverse(mo_num)
double precision, intent(in) :: c_1 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) integer :: exc(0:2,2)
double precision :: phase double precision :: phase
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace

699
src/two_body_rdm/routines_compute_2rdm_orb_range_openmp.irp.f
View File

@ -26,7 +26,7 @@
integer :: occ(N_int*bit_kind_size,2) integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(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) integer(bit_kind) :: det_1_act(N_int,2)
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
do i = 1, N_int do i = 1, N_int
@ -201,7 +201,7 @@
double precision, intent(out) :: values(sze_buff) double precision, intent(out) :: values(sze_buff)
integer , intent(out) :: keys(4,sze_buff) integer , intent(out) :: keys(4,sze_buff)
integer , intent(inout):: nkeys 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) integer :: exc(0:2,2,2)
double precision :: phase double precision :: phase
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
@ -288,7 +288,7 @@
integer :: occ(N_int*bit_kind_size,2) integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(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) integer :: exc(0:2,2,2)
double precision :: phase double precision :: phase
@ -409,310 +409,399 @@
endif endif
end end
! subroutine orb_range_off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) 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 BEGIN_DOC
!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for ! 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 ! 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 ! 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 ! 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 ! 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 determines which spin-spin component of the two-rdm you will update
!! !
!! ispin == 1 :: alpha/ alpha ! ispin == 1 :: alpha/ alpha
!! ispin == 2 :: beta / beta ! ispin == 2 :: beta / beta
!! ispin == 3 :: alpha/ beta ! ispin == 3 :: alpha/ beta
!! ispin == 4 :: spin traced <=> total two-rdm ! ispin == 4 :: spin traced <=> total two-rdm
!! !
!! here, only ispin == 1 or 4 will do something ! here, only ispin == 1 or 4 will do something
! END_DOC END_DOC
! use bitmasks use bitmasks
! implicit none implicit none
! integer, intent(in) :: dim1,ispin integer, intent(in) :: ispin,sze_buff
! 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) :: 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
! integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(out) :: values(sze_buff)
! double precision, intent(in) :: c_1 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 :: occ(N_int*bit_kind_size,2)
! integer :: i,j,h1,h2,istate,p1 integer :: n_occ_ab(2)
! integer :: exc(0:2,2,2) integer :: i,j,h1,h2,p1
! double precision :: phase integer :: exc(0:2,2,2)
! double precision :: phase
! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
! logical :: is_integer_in_string logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
! alpha_alpha = .False. logical :: is_integer_in_string
! beta_beta = .False. alpha_alpha = .False.
! alpha_beta = .False. beta_beta = .False.
! spin_trace = .False. alpha_beta = .False.
! if( ispin == 1)then spin_trace = .False.
! alpha_alpha = .True. if( ispin == 1)then
! else if(ispin == 2)then alpha_alpha = .True.
! beta_beta = .True. else if(ispin == 2)then
! else if(ispin == 3)then beta_beta = .True.
! alpha_beta = .True. else if(ispin == 3)then
! else if(ispin == 4)then alpha_beta = .True.
! spin_trace = .True. else if(ispin == 4)then
! endif 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) call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
! if(alpha_alpha.or.spin_trace)then call get_single_excitation(det_1,det_2,exc,phase,N_int)
! if (exc(0,1,1) == 1) then if(alpha_alpha.or.spin_trace)then
! ! Mono alpha if (exc(0,1,1) == 1) then
! h1 = exc(1,1,1) ! Mono alpha
! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return h1 = exc(1,1,1)
! h1 = list_orb_reverse(h1) if(list_orb_reverse(h1).lt.0)return
! p1 = exc(1,2,1) h1 = list_orb_reverse(h1)
! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return p1 = exc(1,2,1)
! p1 = list_orb_reverse(p1) if(list_orb_reverse(p1).lt.0)return
! do i = 1, n_occ_ab(1) p1 = list_orb_reverse(p1)
! h2 = occ(i,1) do i = 1, n_occ_ab(1)
! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle h2 = occ(i,1)
! h2 = list_orb_reverse(h2) if(list_orb_reverse(h2).lt.0)return
! big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase h2 = list_orb_reverse(h2)
! 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_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) nkeys += 1
! use bitmasks values(nkeys) = 0.5d0 * c_1 * phase
! BEGIN_DOC keys(1,nkeys) = h1
!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for keys(2,nkeys) = h2
!! keys(3,nkeys) = p1
!! a given couple of determinant det_1, det_2 being a BETA SINGLE excitation with respect to one another keys(4,nkeys) = h2
!!
!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 nkeys += 1
!! values(nkeys) = - 0.5d0 * c_1 * phase
!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation keys(1,nkeys) = h1
!! keys(2,nkeys) = h2
!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals keys(3,nkeys) = h2
!! keys(4,nkeys) = p1
!! ispin determines which spin-spin component of the two-rdm you will update
!! nkeys += 1
!! ispin == 1 :: alpha/ alpha values(nkeys) = 0.5d0 * c_1 * phase
!! ispin == 2 :: beta / beta keys(1,nkeys) = h2
!! ispin == 3 :: alpha/ beta keys(2,nkeys) = h1
!! ispin == 4 :: spin traced <=> total two-rdm keys(3,nkeys) = h2
!! keys(4,nkeys) = p1
!! here, only ispin == 2 or 4 will do something
! END_DOC nkeys += 1
! implicit none values(nkeys) = - 0.5d0 * c_1 * phase
! integer, intent(in) :: dim1,ispin keys(1,nkeys) = h2
! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) keys(2,nkeys) = h1
! integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2) keys(3,nkeys) = p1
! integer(bit_kind), intent(in) :: orb_bitmask(N_int) keys(4,nkeys) = h2
! integer, intent(in) :: list_orb_reverse(mo_num) enddo
! double precision, intent(in) :: c_1 else
! return
! endif
! integer :: occ(N_int*bit_kind_size,2) endif
! integer :: n_occ_ab(2) end
! integer :: i,j,h1,h2,istate,p1
! integer :: exc(0:2,2,2) 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)
! double precision :: phase use bitmasks
! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace BEGIN_DOC
! logical :: is_integer_in_string ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for
! alpha_alpha = .False. !
! beta_beta = .False. ! a given couple of determinant det_1, det_2 being a BETA SINGLE excitation with respect to one another
! alpha_beta = .False. !
! spin_trace = .False. ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
! if( ispin == 1)then !
! alpha_alpha = .True. ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation
! else if(ispin == 2)then !
! beta_beta = .True. ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
! else if(ispin == 3)then !
! alpha_beta = .True. ! ispin determines which spin-spin component of the two-rdm you will update
! else if(ispin == 4)then !
! spin_trace = .True. ! ispin == 1 :: alpha/ alpha
! endif ! ispin == 2 :: beta / beta
! ! ispin == 3 :: alpha/ beta
! ! ispin == 4 :: spin traced <=> total two-rdm
! 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) ! here, only ispin == 2 or 4 will do something
! if(beta_beta.or.spin_trace)then END_DOC
! if (exc(0,1,1) == 1) then implicit none
! return integer, intent(in) :: ispin,sze_buff
! else integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
! ! Mono beta integer, intent(in) :: list_orb_reverse(mo_num)
! h1 = exc(1,1,2) double precision, intent(in) :: c_1
! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return double precision, intent(out) :: values(sze_buff)
! h1 = list_orb_reverse(h1) integer , intent(out) :: keys(4,sze_buff)
! p1 = exc(1,2,2) integer , intent(inout):: nkeys
! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return
! p1 = list_orb_reverse(p1) integer :: occ(N_int*bit_kind_size,2)
! do istate = 1, N_states integer :: n_occ_ab(2)
! do i = 1, n_occ_ab(2) integer :: i,j,h1,h2,p1
! h2 = occ(i,2) integer :: exc(0:2,2,2)
! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle double precision :: phase
! h2 = list_orb_reverse(h2) logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
! big_array(h1,h2,p1,h2) += 0.5d0 * c_1 * phase logical :: is_integer_in_string
! big_array(h1,h2,h2,p1) -= 0.5d0 * c_1 * phase alpha_alpha = .False.
! beta_beta = .False.
! big_array(h2,h1,h2,p1) += 0.5d0 * c_1 * phase alpha_beta = .False.
! big_array(h2,h1,p1,h2) -= 0.5d0 * c_1 * phase spin_trace = .False.
! enddo if( ispin == 1)then
! enddo alpha_alpha = .True.
! endif else if(ispin == 2)then
! endif beta_beta = .True.
! end 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) 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 use bitmasks
! BEGIN_DOC BEGIN_DOC
!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for ! 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 ! 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 ! 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 ! 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 ! 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 determines which spin-spin component of the two-rdm you will update
!! !
!! ispin == 1 :: alpha/ alpha ! ispin == 1 :: alpha/ alpha
!! ispin == 2 :: beta / beta ! ispin == 2 :: beta / beta
!! ispin == 3 :: alpha/ beta ! ispin == 3 :: alpha/ beta
!! ispin == 4 :: spin traced <=> total two-rdm ! ispin == 4 :: spin traced <=> total two-rdm
!! !
!! here, only ispin == 1 or 4 will do something ! here, only ispin == 1 or 4 will do something
! END_DOC END_DOC
! implicit none implicit none
! integer, intent(in) :: dim1,ispin integer, intent(in) :: ispin,sze_buff
! 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) :: det_1(N_int),det_2(N_int) integer, intent(in) :: list_orb_reverse(mo_num)
! integer(bit_kind), intent(in) :: orb_bitmask(N_int) double precision, intent(in) :: c_1
! integer, intent(in) :: list_orb_reverse(mo_num) double precision, intent(out) :: values(sze_buff)
! double precision, intent(in) :: c_1 integer , intent(out) :: keys(4,sze_buff)
! integer , intent(inout):: nkeys
! integer :: i,j,h1,h2,p1,p2,istate
! integer :: exc(0:2,2)
! double precision :: phase integer :: i,j,h1,h2,p1,p2
! integer :: exc(0:2,2)
! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace double precision :: phase
! logical :: is_integer_in_string
! alpha_alpha = .False. logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
! beta_beta = .False. logical :: is_integer_in_string
! alpha_beta = .False. alpha_alpha = .False.
! spin_trace = .False. beta_beta = .False.
! if( ispin == 1)then alpha_beta = .False.
! alpha_alpha = .True. spin_trace = .False.
! else if(ispin == 2)then if( ispin == 1)then
! beta_beta = .True. alpha_alpha = .True.
! else if(ispin == 3)then else if(ispin == 2)then
! alpha_beta = .True. beta_beta = .True.
! else if(ispin == 4)then else if(ispin == 3)then
! spin_trace = .True. alpha_beta = .True.
! endif else if(ispin == 4)then
! call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) spin_trace = .True.
! h1 =exc(1,1) endif
! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return call get_double_excitation_spin(det_1,det_2,exc,phase,N_int)
! h1 = list_orb_reverse(h1) h1 =exc(1,1)
! h2 =exc(2,1) if(list_orb_reverse(h1).lt.0)return
! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return h1 = list_orb_reverse(h1)
! h2 = list_orb_reverse(h2) h2 =exc(2,1)
! p1 =exc(1,2) if(list_orb_reverse(h2).lt.0)return
! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return h2 = list_orb_reverse(h2)
! p1 = list_orb_reverse(p1) p1 =exc(1,2)
! p2 =exc(2,2) if(list_orb_reverse(p1).lt.0)return
! if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return p1 = list_orb_reverse(p1)
! p2 = list_orb_reverse(p2) p2 =exc(2,2)
! if(alpha_alpha.or.spin_trace)then if(list_orb_reverse(p2).lt.0)return
! do istate = 1, N_states p2 = list_orb_reverse(p2)
! big_array(h1,h2,p1,p2) += 0.5d0 * c_1 * phase if(alpha_alpha.or.spin_trace)then
! big_array(h1,h2,p2,p1) -= 0.5d0 * c_1 * phase nkeys += 1
! values(nkeys) = 0.5d0 * c_1 * phase
! big_array(h2,h1,p2,p1) += 0.5d0 * c_1 * phase keys(1,nkeys) = h1
! big_array(h2,h1,p1,p2) -= 0.5d0 * c_1 * phase keys(2,nkeys) = h2
! enddo keys(3,nkeys) = p1
! endif keys(4,nkeys) = p2
! end
! subroutine orb_range_off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,gorb_bitmask,list_orb_reverse,ispin) nkeys += 1
! use bitmasks values(nkeys) = - 0.5d0 * c_1 * phase
! BEGIN_DOC keys(1,nkeys) = h1
!! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for keys(2,nkeys) = h2
!! keys(3,nkeys) = p2
!! a given couple of determinant det_1, det_2 being a BETA /BETA DOUBLE excitation with respect to one another keys(4,nkeys) = p1
!!
!! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1 nkeys += 1
!! values(nkeys) = 0.5d0 * c_1 * phase
!! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation keys(1,nkeys) = h2
!! keys(2,nkeys) = h1
!! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals keys(3,nkeys) = p2
!! keys(4,nkeys) = p1
!! ispin determines which spin-spin component of the two-rdm you will update
!! nkeys += 1
!! ispin == 1 :: alpha/ alpha values(nkeys) = - 0.5d0 * c_1 * phase
!! ispin == 2 :: beta / beta keys(1,nkeys) = h2
!! ispin == 3 :: alpha/ beta keys(2,nkeys) = h1
!! ispin == 4 :: spin traced <=> total two-rdm keys(3,nkeys) = p1
!! keys(4,nkeys) = p2
!! here, only ispin == 2 or 4 will do something endif
! END_DOC end
! implicit none
! 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)
! integer, intent(in) :: dim1,ispin use bitmasks
! double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1) BEGIN_DOC
! integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int) ! routine that update the OFF DIAGONAL PART of the two body rdms in a specific range of orbitals for
! integer(bit_kind), intent(in) :: orb_bitmask(N_int) !
! integer, intent(in) :: list_orb_reverse(mo_num) ! a given couple of determinant det_1, det_2 being a BETA /BETA DOUBLE excitation with respect to one another
! double precision, intent(in) :: c_1 !
! ! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
! integer :: i,j,h1,h2,p1,p2,istate !
! integer :: exc(0:2,2) ! big_array(dim1,dim1,dim1,dim1) is the two-body rdm to be updated in physicist notation
! double precision :: phase !
! logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace ! orb_bitmask(N_int) is the bitmask for the orbital range, list_orb_reverse(mo_num) is the inverse range of orbitals
! logical :: is_integer_in_string !
! alpha_alpha = .False. ! ispin determines which spin-spin component of the two-rdm you will update
! beta_beta = .False. !
! alpha_beta = .False. ! ispin == 1 :: alpha/ alpha
! spin_trace = .False. ! ispin == 2 :: beta / beta
! if( ispin == 1)then ! ispin == 3 :: alpha/ beta
! alpha_alpha = .True. ! ispin == 4 :: spin traced <=> total two-rdm
! else if(ispin == 2)then !
! beta_beta = .True. ! here, only ispin == 2 or 4 will do something
! else if(ispin == 3)then END_DOC
! alpha_beta = .True. implicit none
! else if(ispin == 4)then
! spin_trace = .True. integer, intent(in) :: ispin,sze_buff
! endif integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
! integer, intent(in) :: list_orb_reverse(mo_num)
! call get_double_excitation_spin(det_1,det_2,exc,phase,N_int) double precision, intent(in) :: c_1
! h1 =exc(1,1) double precision, intent(out) :: values(sze_buff)
! if(.not.is_integer_in_string(h1,orb_bitmask,N_int))return integer , intent(out) :: keys(4,sze_buff)
! h1 = list_orb_reverse(h1) integer , intent(inout):: nkeys
! h2 =exc(2,1)
! if(.not.is_integer_in_string(h2,orb_bitmask,N_int))return integer :: i,j,h1,h2,p1,p2
! h2 = list_orb_reverse(h2) integer :: exc(0:2,2)
! p1 =exc(1,2) double precision :: phase
! if(.not.is_integer_in_string(p1,orb_bitmask,N_int))return logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
! p1 = list_orb_reverse(p1) logical :: is_integer_in_string
! p2 =exc(2,2) alpha_alpha = .False.
! if(.not.is_integer_in_string(p2,orb_bitmask,N_int))return beta_beta = .False.
! p2 = list_orb_reverse(p2) alpha_beta = .False.
! if(beta_beta.or.spin_trace)then spin_trace = .False.
! big_array(h1,h2,p1,p2) += 0.5d0 * c_1* phase if( ispin == 1)then
! big_array(h1,h2,p2,p1) -= 0.5d0 * c_1* phase alpha_alpha = .True.
! else if(ispin == 2)then
! big_array(h2,h1,p2,p1) += 0.5d0 * c_1* phase beta_beta = .True.
! big_array(h2,h1,p1,p2) -= 0.5d0 * c_1* phase else if(ispin == 3)then
! endif alpha_beta = .True.
! end 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