1003 lines
29 KiB
Fortran
1003 lines
29 KiB
Fortran
subroutine orb_range_diag_to_all_states_2_rdm_trans_buffer(det_1,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
|
|
use bitmasks
|
|
BEGIN_DOC
|
|
! routine that update the DIAGONAL PART of the two body trans_rdms in a specific range of orbitals for a given determinant det_1
|
|
!
|
|
! c_1 is the array of the contributions to the trans_rdm for all states
|
|
!
|
|
! 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-trans_rdm you will update
|
|
!
|
|
! ispin == 1 :: alpha/ alpha
|
|
! ispin == 2 :: beta / beta
|
|
! ispin == 3 :: alpha/ beta
|
|
! ispin == 4 :: spin traced <=> total two-trans_rdm
|
|
END_DOC
|
|
implicit none
|
|
integer, intent(in) :: ispin,sze_buff,N_st
|
|
integer, intent(in) :: list_orb_reverse(mo_num)
|
|
integer(bit_kind), intent(in) :: det_1(N_int,2)
|
|
integer(bit_kind), intent(in) :: orb_bitmask(N_int)
|
|
double precision, intent(in) :: c_1(N_st,N_st)
|
|
double precision, intent(out) :: values(N_st,N_st,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
|
|
integer(bit_kind) :: det_1_act(N_int,2)
|
|
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
|
|
do i = 1, N_int
|
|
det_1_act(i,1) = iand(det_1(i,1),orb_bitmask(i))
|
|
det_1_act(i,2) = iand(det_1(i,2),orb_bitmask(i))
|
|
enddo
|
|
|
|
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_act, occ, n_occ_ab, N_int)
|
|
logical :: is_integer_in_string
|
|
integer :: i1,i2,istate
|
|
integer :: jstate
|
|
if(alpha_beta)then
|
|
do i = 1, n_occ_ab(1)
|
|
i1 = occ(i,1)
|
|
do j = 1, n_occ_ab(2)
|
|
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
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h1
|
|
keys(4,nkeys) = h2
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = h1
|
|
enddo
|
|
enddo
|
|
|
|
else if (alpha_alpha)then
|
|
do i = 1, n_occ_ab(1)
|
|
i1 = occ(i,1)
|
|
do j = 1, n_occ_ab(1)
|
|
i2 = occ(j,1)
|
|
h1 = list_orb_reverse(i1)
|
|
h2 = list_orb_reverse(i2)
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h1
|
|
keys(4,nkeys) = h2
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = h1
|
|
enddo
|
|
enddo
|
|
else if (beta_beta)then
|
|
do i = 1, n_occ_ab(2)
|
|
i1 = occ(i,2)
|
|
do j = 1, n_occ_ab(2)
|
|
i2 = occ(j,2)
|
|
h1 = list_orb_reverse(i1)
|
|
h2 = list_orb_reverse(i2)
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h1
|
|
keys(4,nkeys) = h2
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = h1
|
|
enddo
|
|
enddo
|
|
else if(spin_trace)then
|
|
! 0.5 * (alpha beta + beta alpha)
|
|
do i = 1, n_occ_ab(1)
|
|
i1 = occ(i,1)
|
|
do j = 1, n_occ_ab(2)
|
|
i2 = occ(j,2)
|
|
h1 = list_orb_reverse(i1)
|
|
h2 = list_orb_reverse(i2)
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h1
|
|
keys(4,nkeys) = h2
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = h1
|
|
enddo
|
|
enddo
|
|
do i = 1, n_occ_ab(1)
|
|
i1 = occ(i,1)
|
|
do j = 1, n_occ_ab(1)
|
|
i2 = occ(j,1)
|
|
h1 = list_orb_reverse(i1)
|
|
h2 = list_orb_reverse(i2)
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h1
|
|
keys(4,nkeys) = h2
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = h1
|
|
enddo
|
|
enddo
|
|
do i = 1, n_occ_ab(2)
|
|
i1 = occ(i,2)
|
|
do j = 1, n_occ_ab(2)
|
|
i2 = occ(j,2)
|
|
h1 = list_orb_reverse(i1)
|
|
h2 = list_orb_reverse(i2)
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h1
|
|
keys(4,nkeys) = h2
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = -0.5d0 * c_1(istate,jstate)
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = h1
|
|
enddo
|
|
enddo
|
|
endif
|
|
end
|
|
|
|
|
|
subroutine orb_range_off_diag_double_to_all_states_ab_trans_rdm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
|
|
use bitmasks
|
|
BEGIN_DOC
|
|
! routine that update the OFF DIAGONAL PART of the two body trans_rdms in a specific range of orbitals for
|
|
!
|
|
! a given couple of determinant det_1, det_2 being a alpha/beta DOUBLE excitation with respect to one another
|
|
!
|
|
! c_1 is the array of the contributions to the trans_rdm for all states
|
|
!
|
|
! 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-trans_rdm you will update
|
|
!
|
|
! ispin == 1 :: alpha/ alpha
|
|
! ispin == 2 :: beta / beta
|
|
! ispin == 3 :: alpha/ beta
|
|
! ispin == 4 :: spin traced <=> total two-trans_rdm
|
|
!
|
|
! here, only ispin == 3 or 4 will do something
|
|
END_DOC
|
|
implicit none
|
|
integer, intent(in) :: ispin,sze_buff,N_st
|
|
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(N_st,N_st)
|
|
double precision, intent(out) :: values(N_st,N_st,sze_buff)
|
|
integer , intent(out) :: keys(4,sze_buff)
|
|
integer , intent(inout):: nkeys
|
|
integer :: i,j,h1,h2,p1,p2,istate
|
|
integer :: exc(0:2,2,2)
|
|
double precision :: phase
|
|
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
|
|
logical :: is_integer_in_string
|
|
integer :: jstate
|
|
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(det_1,det_2,exc,phase,N_int)
|
|
h1 = exc(1,1,1)
|
|
if(list_orb_reverse(h1).lt.0)return
|
|
h1 = list_orb_reverse(h1)
|
|
h2 = exc(1,1,2)
|
|
if(list_orb_reverse(h2).lt.0)return
|
|
h2 = list_orb_reverse(h2)
|
|
p1 = exc(1,2,1)
|
|
if(list_orb_reverse(p1).lt.0)return
|
|
p1 = list_orb_reverse(p1)
|
|
p2 = exc(1,2,2)
|
|
if(list_orb_reverse(p2).lt.0)return
|
|
p2 = list_orb_reverse(p2)
|
|
if(alpha_beta)then
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = p2
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = p2
|
|
keys(4,nkeys) = p1
|
|
else if(spin_trace)then
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = p2
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = p2
|
|
keys(4,nkeys) = p1
|
|
endif
|
|
end
|
|
|
|
subroutine orb_range_off_diag_single_to_all_states_ab_trans_rdm_buffer(det_1,det_2,c_1,N_st,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 trans_rdms in a specific range of orbitals for
|
|
!
|
|
! a given couple of determinant det_1, det_2 being a 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 trans_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-trans_rdm you will update
|
|
!
|
|
! ispin == 1 :: alpha/ alpha
|
|
! ispin == 2 :: beta / beta
|
|
! ispin == 3 :: alpha/ beta
|
|
! ispin == 4 :: spin traced <=> total two-trans_rdm
|
|
!
|
|
! here, only ispin == 3 or 4 will do something
|
|
END_DOC
|
|
implicit none
|
|
integer, intent(in) :: ispin,sze_buff,N_st
|
|
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(N_st,N_st)
|
|
double precision, intent(out) :: values(N_st,N_st,sze_buff)
|
|
integer , intent(out) :: keys(4,sze_buff)
|
|
integer , intent(inout):: nkeys
|
|
|
|
integer :: jstate
|
|
integer :: occ(N_int*bit_kind_size,2)
|
|
integer :: n_occ_ab(2)
|
|
integer :: i,j,h1,h2,p1,istate
|
|
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_beta)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(2)
|
|
h2 = occ(i,2)
|
|
if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
|
|
h2 = list_orb_reverse(h2)
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = h2
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = p1
|
|
enddo
|
|
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 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)
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = h2
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = p1
|
|
enddo
|
|
endif
|
|
else if(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(2)
|
|
h2 = occ(i,2)
|
|
if(.not.is_integer_in_string(h2,orb_bitmask,N_int))cycle
|
|
h2 = list_orb_reverse(h2)
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = h2
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = p1
|
|
enddo
|
|
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 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)
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = h2
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = p1
|
|
enddo
|
|
endif
|
|
endif
|
|
end
|
|
|
|
subroutine orb_range_off_diag_single_to_all_states_aa_trans_rdm_buffer(det_1,det_2,c_1,N_st,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
|
|
BEGIN_DOC
|
|
! routine that update the OFF DIAGONAL PART of the two body trans_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 trans_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-trans_rdm you will update
|
|
!
|
|
! ispin == 1 :: alpha/ alpha
|
|
! ispin == 2 :: beta / beta
|
|
! ispin == 3 :: alpha/ beta
|
|
! ispin == 4 :: spin traced <=> total two-trans_rdm
|
|
!
|
|
! here, only ispin == 1 or 4 will do something
|
|
END_DOC
|
|
use bitmasks
|
|
implicit none
|
|
integer, intent(in) :: ispin,sze_buff,N_st
|
|
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(N_st,N_st)
|
|
double precision, intent(out) :: values(N_st,N_st,sze_buff)
|
|
integer , intent(out) :: keys(4,sze_buff)
|
|
integer , intent(inout):: nkeys
|
|
|
|
integer :: jstate
|
|
integer :: occ(N_int*bit_kind_size,2)
|
|
integer :: n_occ_ab(2)
|
|
integer :: i,j,h1,h2,p1,istate
|
|
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)
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = h2
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = p1
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = p1
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
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_all_states_bb_trans_rdm_buffer(det_1,det_2,c_1,N_st,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 trans_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 trans_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-trans_rdm you will update
|
|
!
|
|
! ispin == 1 :: alpha/ alpha
|
|
! ispin == 2 :: beta / beta
|
|
! ispin == 3 :: alpha/ beta
|
|
! ispin == 4 :: spin traced <=> total two-trans_rdm
|
|
!
|
|
! here, only ispin == 2 or 4 will do something
|
|
END_DOC
|
|
implicit none
|
|
integer, intent(in) :: ispin,sze_buff,N_st
|
|
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(N_st,N_st)
|
|
double precision, intent(out) :: values(N_st,N_st,sze_buff)
|
|
integer , intent(out) :: keys(4,sze_buff)
|
|
integer , intent(inout):: nkeys
|
|
|
|
integer :: jstate
|
|
integer :: occ(N_int*bit_kind_size,2)
|
|
integer :: n_occ_ab(2)
|
|
integer :: i,j,h1,h2,p1,istate
|
|
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 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)
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = h2
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = p1
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = h2
|
|
keys(4,nkeys) = p1
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = h2
|
|
enddo
|
|
endif
|
|
endif
|
|
end
|
|
|
|
|
|
subroutine orb_range_off_diag_double_to_all_states_aa_trans_rdm_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
|
|
use bitmasks
|
|
BEGIN_DOC
|
|
! routine that update the OFF DIAGONAL PART of the two body trans_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 trans_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-trans_rdm you will update
|
|
!
|
|
! ispin == 1 :: alpha/ alpha
|
|
! ispin == 2 :: beta / beta
|
|
! ispin == 3 :: alpha/ beta
|
|
! ispin == 4 :: spin traced <=> total two-trans_rdm
|
|
!
|
|
! here, only ispin == 1 or 4 will do something
|
|
END_DOC
|
|
implicit none
|
|
integer, intent(in) :: ispin,sze_buff,N_st
|
|
integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int)
|
|
integer, intent(in) :: list_orb_reverse(mo_num)
|
|
double precision, intent(in) :: c_1(N_st,N_st)
|
|
double precision, intent(out) :: values(N_st,N_st,sze_buff)
|
|
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 :: jstate
|
|
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
|
|
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = p2
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p2
|
|
keys(4,nkeys) = p1
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = p2
|
|
keys(4,nkeys) = p1
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
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_all_states_trans_rdm_bb_buffer(det_1,det_2,c_1,N_st,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
|
|
use bitmasks
|
|
BEGIN_DOC
|
|
! routine that update the OFF DIAGONAL PART of the two body trans_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 trans_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-trans_rdm you will update
|
|
!
|
|
! ispin == 1 :: alpha/ alpha
|
|
! ispin == 2 :: beta / beta
|
|
! ispin == 3 :: alpha/ beta
|
|
! ispin == 4 :: spin traced <=> total two-trans_rdm
|
|
!
|
|
! here, only ispin == 2 or 4 will do something
|
|
END_DOC
|
|
implicit none
|
|
|
|
integer, intent(in) :: ispin,sze_buff,N_st
|
|
integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int)
|
|
integer, intent(in) :: list_orb_reverse(mo_num)
|
|
double precision, intent(in) :: c_1(N_st,N_st)
|
|
double precision, intent(out) :: values(N_st,N_st,sze_buff)
|
|
integer , intent(out) :: keys(4,sze_buff)
|
|
integer , intent(inout):: nkeys
|
|
|
|
integer :: jstate
|
|
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(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
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = p2
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h1
|
|
keys(2,nkeys) = h2
|
|
keys(3,nkeys) = p2
|
|
keys(4,nkeys) = p1
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = p2
|
|
keys(4,nkeys) = p1
|
|
|
|
nkeys += 1
|
|
do jstate = 1, N_st
|
|
do istate = 1, N_st
|
|
values(istate,jstate,nkeys) = - 0.5d0 * c_1(istate,jstate) * phase
|
|
enddo
|
|
enddo
|
|
keys(1,nkeys) = h2
|
|
keys(2,nkeys) = h1
|
|
keys(3,nkeys) = p1
|
|
keys(4,nkeys) = p2
|
|
endif
|
|
end
|
|
|