LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-21 11:03:29 +01:00
Code Releases Activity

beginning the cleaning of two_body_rdm

Browse Source
This commit is contained in:
Emmanuel Giner 2020-03-18 15:13:49 +01:00
parent be3aa1af71
commit 3e0ada9538
16 changed files with 50 additions and 2662 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/density_for_dft/EZFIO.cfg
View File

@ -11,10 +11,10 @@ interface: ezfio,provider,ocaml
default: 0.5
[no_core_density]
type: character*(32)
doc: Type of density. If [no_core_dm] then all elements of the density matrix involving at least one orbital set as core are set to zero
type: logical
doc: If [no_core_density] then all elements of the density matrix involving at least one orbital set as core are set to zero. The default is False in order to take all the density.
interface: ezfio, provider, ocaml
default: full_density
default: False
[normalize_dm]
type: logical

4
src/density_for_dft/density_for_dft.irp.f
View File

@ -22,7 +22,7 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_alpha_for_dft, (mo_num,mo_num, N_s
one_e_dm_mo_alpha_for_dft(:,:,1) = one_e_dm_mo_alpha_average(:,:)
endif
if(no_core_density .EQ. "no_core_dm")then
if(no_core_density)then
integer :: ii,i,j
do ii = 1, n_core_orb
i = list_core(ii)
@ -73,7 +73,7 @@ BEGIN_PROVIDER [double precision, one_e_dm_mo_beta_for_dft, (mo_num,mo_num, N_st
one_e_dm_mo_beta_for_dft(:,:,1) = one_e_dm_mo_beta_average(:,:)
endif
if(no_core_density .EQ. "no_core_dm")then
if(no_core_density)then
integer :: ii,i,j
do ii = 1, n_core_orb
i = list_core(ii)

1
src/two_body_rdm/README.rst
View File

@ -6,3 +6,4 @@ Contains the two rdms $\alpha\alpha$, $\beta\beta$ and $\alpha\beta$ stored as
arrays, with pysicists notation, consistent with the two-electron integrals in the
MO basis.

402
src/two_body_rdm/ab_only_routines.irp.f
View File

@ -1,402 +0,0 @@
subroutine two_rdm_ab_nstates(big_array,dim1,dim2,dim3,dim4,u_0,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes the alpha/beta part of the two-body density matrix IN CHEMIST NOTATIONS
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
integer, intent(in) :: N_st,sze
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states)
double precision, intent(inout) :: u_0(sze,N_st)
integer :: k
double precision, allocatable :: u_t(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
allocate(u_t(N_st,N_det))
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
enddo
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_st)
call two_rdm_ab_nstates_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1)
deallocate(u_t)
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
enddo
end
subroutine two_rdm_ab_nstates_work(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes the alpha/beta part of the two-body density matrix
!
! Default should be 1,N_det,0,1
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states)
double precision, intent(in) :: u_t(N_st,N_det)
PROVIDE N_int
select case (N_int)
case (1)
call two_rdm_ab_nstates_work_1(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
case (2)
call two_rdm_ab_nstates_work_2(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
case (3)
call two_rdm_ab_nstates_work_3(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
case (4)
call two_rdm_ab_nstates_work_4(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
case default
call two_rdm_ab_nstates_work_N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
end select
end
BEGIN_TEMPLATE
subroutine two_rdm_ab_nstates_work_$N_int(big_array,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states)
double precision, intent(in) :: u_t(N_st,N_det)
double precision :: hij, sij
integer :: i,j,k,l
integer :: k_a, k_b, l_a, l_b, m_a, m_b
integer :: istate
integer :: krow, kcol, krow_b, kcol_b
integer :: lrow, lcol
integer :: mrow, mcol
integer(bit_kind) :: spindet($N_int)
integer(bit_kind) :: tmp_det($N_int,2)
integer(bit_kind) :: tmp_det2($N_int,2)
integer(bit_kind) :: tmp_det3($N_int,2)
integer(bit_kind), allocatable :: buffer(:,:)
integer :: n_doubles
integer, allocatable :: doubles(:)
integer, allocatable :: singles_a(:)
integer, allocatable :: singles_b(:)
integer, allocatable :: idx(:), idx0(:)
integer :: maxab, n_singles_a, n_singles_b, kcol_prev, nmax
integer*8 :: k8
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
allocate(idx0(maxab))
do i=1,maxab
idx0(i) = i
enddo
! Prepare the array of all alpha single excitations
! -------------------------------------------------
PROVIDE N_int nthreads_davidson
! Alpha/Beta double excitations
! =============================
allocate( buffer($N_int,maxab), &
singles_a(maxab), &
singles_b(maxab), &
doubles(maxab), &
idx(maxab))
kcol_prev=-1
ASSERT (iend <= N_det)
ASSERT (istart > 0)
ASSERT (istep > 0)
do k_a=istart+ishift,iend,istep
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
if (kcol /= kcol_prev) then
call get_all_spin_singles_$N_int( &
psi_det_beta_unique, idx0, &
tmp_det(1,2), N_det_beta_unique, &
singles_b, n_singles_b)
endif
kcol_prev = kcol
! Loop over singly excited beta columns
! -------------------------------------
do i=1,n_singles_b
lcol = singles_b(i)
tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
l_a = psi_bilinear_matrix_columns_loc(lcol)
ASSERT (l_a <= N_det)
do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)
ASSERT (l_a <= N_det)
idx(j) = l_a
l_a = l_a+1
enddo
j = j-1
call get_all_spin_singles_$N_int( &
buffer, idx, tmp_det(1,1), j, &
singles_a, n_singles_a )
! Loop over alpha singles
! -----------------------
do k = 1,n_singles_a
l_a = singles_a(k)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
!!!!!!!!!!!!!!!!!! ALPHA BETA
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
enddo
call off_diagonal_double_to_two_rdm_ab_dm(tmp_det,tmp_det2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
enddo
enddo
enddo
do k_a=istart+ishift,iend,istep
! Single and double alpha excitations
! ===================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
! Initial determinant is at k_b in beta-major representation
! ----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
spindet(1:$N_int) = tmp_det(1:$N_int,1)
! Loop inside the beta column to gather all the connected alphas
lcol = psi_bilinear_matrix_columns(k_a)
l_a = psi_bilinear_matrix_columns_loc(lcol)
do i=1,N_det_alpha_unique
if (l_a > N_det) exit
lcol = psi_bilinear_matrix_columns(l_a)
if (lcol /= kcol) exit
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
idx(i) = l_a
l_a = l_a+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_a, doubles, n_singles_a, n_doubles )
! Compute Hij for all alpha singles
! ----------------------------------
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
do i=1,n_singles_a
l_a = singles_a(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
!!!! MONO SPIN
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
enddo
call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
enddo
!! Compute Hij for all alpha doubles
!! ----------------------------------
!
!do i=1,n_doubles
! l_a = doubles(i)
! ASSERT (l_a <= N_det)
! lrow = psi_bilinear_matrix_rows(l_a)
! ASSERT (lrow <= N_det_alpha_unique)
! call i_H_j_double_spin_erf( tmp_det(1,1), psi_det_alpha_unique(1, lrow), $N_int, hij)
! do l=1,N_st
! v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
! ! same spin => sij = 0
! enddo
!enddo
! Single and double beta excitations
! ==================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
kcol = psi_bilinear_matrix_columns(k_a)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
spindet(1:$N_int) = tmp_det(1:$N_int,2)
! Initial determinant is at k_b in beta-major representation
! -----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
! Loop inside the alpha row to gather all the connected betas
lrow = psi_bilinear_matrix_transp_rows(k_b)
l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
do i=1,N_det_beta_unique
if (l_b > N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l_b)
if (lrow /= krow) exit
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
idx(i) = l_b
l_b = l_b+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_b, doubles, n_singles_b, n_doubles )
! Compute Hij for all beta singles
! ----------------------------------
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
do i=1,n_singles_b
l_b = singles_b(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
l_a = psi_bilinear_matrix_transp_order(l_b)
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
enddo
call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
ASSERT (l_a <= N_det)
enddo
!
!! Compute Hij for all beta doubles
!! ----------------------------------
!
!do i=1,n_doubles
! l_b = doubles(i)
! ASSERT (l_b <= N_det)
! lcol = psi_bilinear_matrix_transp_columns(l_b)
! ASSERT (lcol <= N_det_beta_unique)
! call i_H_j_double_spin_erf( tmp_det(1,2), psi_det_beta_unique(1, lcol), $N_int, hij)
! l_a = psi_bilinear_matrix_transp_order(l_b)
! ASSERT (l_a <= N_det)
! do l=1,N_st
! v_t(l,k_a) = v_t(l,k_a) + hij * u_t(l,l_a)
! ! same spin => sij = 0
! enddo
!enddo
! Diagonal contribution
! =====================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
double precision, external :: diag_H_mat_elem_erf, diag_S_mat_elem
double precision :: c_1(N_states),c_2(N_states)
do l = 1, N_states
c_1(l) = u_t(l,k_a)
enddo
call diagonal_contrib_to_two_rdm_ab_dm(tmp_det,c_1,big_array,dim1,dim2,dim3,dim4)
end do
deallocate(buffer, singles_a, singles_b, doubles, idx)
end
SUBST [ N_int ]
1;;
2;;
3;;
4;;
N_int;;
END_TEMPLATE

442
src/two_body_rdm/all_2rdm_routines.irp.f
View File

@ -1,442 +0,0 @@
subroutine all_two_rdm_dm_nstates(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_0,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes the alpha/alpha, beta/beta and alpha/beta part of the two-body density matrix IN CHEMIST NOTATIONS
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
integer, intent(in) :: N_st,sze
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states)
double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states)
double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states)
double precision, intent(inout) :: u_0(sze,N_st)
integer :: k
double precision, allocatable :: u_t(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
allocate(u_t(N_st,N_det))
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
enddo
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_st)
call all_two_rdm_dm_nstates_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,1,N_det,0,1)
deallocate(u_t)
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
enddo
end
subroutine all_two_rdm_dm_nstates_work(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes two-rdm
!
! Default should be 1,N_det,0,1
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states)
double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states)
double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states)
double precision, intent(in) :: u_t(N_st,N_det)
PROVIDE N_int
select case (N_int)
case (1)
call all_two_rdm_dm_nstates_work_1(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
case (2)
call all_two_rdm_dm_nstates_work_2(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
case (3)
call all_two_rdm_dm_nstates_work_3(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
case (4)
call all_two_rdm_dm_nstates_work_4(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
case default
call all_two_rdm_dm_nstates_work_N_int(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
end select
end
BEGIN_TEMPLATE
subroutine all_two_rdm_dm_nstates_work_$N_int(big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes $v_t = H | u_t \\rangle$ and $s_t = S^2 | u_t \\rangle$
!
! Default should be 1,N_det,0,1
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det)
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states)
double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states)
double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states)
integer :: i,j,k,l
integer :: k_a, k_b, l_a, l_b, m_a, m_b
integer :: istate
integer :: krow, kcol, krow_b, kcol_b
integer :: lrow, lcol
integer :: mrow, mcol
integer(bit_kind) :: spindet($N_int)
integer(bit_kind) :: tmp_det($N_int,2)
integer(bit_kind) :: tmp_det2($N_int,2)
integer(bit_kind) :: tmp_det3($N_int,2)
integer(bit_kind), allocatable :: buffer(:,:)
integer :: n_doubles
integer, allocatable :: doubles(:)
integer, allocatable :: singles_a(:)
integer, allocatable :: singles_b(:)
integer, allocatable :: idx(:), idx0(:)
integer :: maxab, n_singles_a, n_singles_b, kcol_prev
integer*8 :: k8
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
allocate(idx0(maxab))
do i=1,maxab
idx0(i) = i
enddo
! Prepare the array of all alpha single excitations
! -------------------------------------------------
PROVIDE N_int nthreads_davidson
!!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) &
! !$OMP SHARED(psi_bilinear_matrix_rows, N_det, &
! !$OMP psi_bilinear_matrix_columns, &
! !$OMP psi_det_alpha_unique, psi_det_beta_unique,&
! !$OMP n_det_alpha_unique, n_det_beta_unique, N_int,&
! !$OMP psi_bilinear_matrix_transp_rows, &
! !$OMP psi_bilinear_matrix_transp_columns, &
! !$OMP psi_bilinear_matrix_transp_order, N_st, &
! !$OMP psi_bilinear_matrix_order_transp_reverse, &
! !$OMP psi_bilinear_matrix_columns_loc, &
! !$OMP psi_bilinear_matrix_transp_rows_loc, &
! !$OMP istart, iend, istep, irp_here, v_t, s_t, &
! !$OMP ishift, idx0, u_t, maxab) &
! !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,&
! !$OMP lcol, lrow, l_a, l_b, &
! !$OMP buffer, doubles, n_doubles, &
! !$OMP tmp_det2, idx, l, kcol_prev, &
! !$OMP singles_a, n_singles_a, singles_b, &
! !$OMP n_singles_b, k8)
! Alpha/Beta double excitations
! =============================
allocate( buffer($N_int,maxab), &
singles_a(maxab), &
singles_b(maxab), &
doubles(maxab), &
idx(maxab))
kcol_prev=-1
ASSERT (iend <= N_det)
ASSERT (istart > 0)
ASSERT (istep > 0)
!!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
if (kcol /= kcol_prev) then
call get_all_spin_singles_$N_int( &
psi_det_beta_unique, idx0, &
tmp_det(1,2), N_det_beta_unique, &
singles_b, n_singles_b)
endif
kcol_prev = kcol
! Loop over singly excited beta columns
! -------------------------------------
do i=1,n_singles_b
lcol = singles_b(i)
tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
l_a = psi_bilinear_matrix_columns_loc(lcol)
ASSERT (l_a <= N_det)
do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)
ASSERT (l_a <= N_det)
idx(j) = l_a
l_a = l_a+1
enddo
j = j-1
call get_all_spin_singles_$N_int( &
buffer, idx, tmp_det(1,1), j, &
singles_a, n_singles_a )
! Loop over alpha singles
! -----------------------
do k = 1,n_singles_a
l_a = singles_a(k)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
!call i_H_j_double_alpha_beta(tmp_det,tmp_det2,$N_int,hij)
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
enddo
call off_diagonal_double_to_two_rdm_ab_dm(tmp_det,tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4)
enddo
enddo
enddo
! !$OMP END DO
! !$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
! Single and double alpha exitations
! ===================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
! Initial determinant is at k_b in beta-major representation
! ----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
spindet(1:$N_int) = tmp_det(1:$N_int,1)
! Loop inside the beta column to gather all the connected alphas
lcol = psi_bilinear_matrix_columns(k_a)
l_a = psi_bilinear_matrix_columns_loc(lcol)
do i=1,N_det_alpha_unique
if (l_a > N_det) exit
lcol = psi_bilinear_matrix_columns(l_a)
if (lcol /= kcol) exit
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
idx(i) = l_a
l_a = l_a+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_a, doubles, n_singles_a, n_doubles )
! Compute Hij for all alpha singles
! ----------------------------------
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
do i=1,n_singles_a
l_a = singles_a(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
enddo
! increment the alpha/beta part for single excitations
call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4)
! increment the alpha/alpha part for single excitations
call off_diagonal_single_to_two_rdm_aa_dm(tmp_det,tmp_det2,c_1,c_2,big_array_aa,dim1,dim2,dim3,dim4)
enddo
! Compute Hij for all alpha doubles
! ----------------------------------
do i=1,n_doubles
l_a = doubles(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
enddo
call off_diagonal_double_to_two_rdm_aa_dm(tmp_det(1,1),psi_det_alpha_unique(1, lrow),c_1,c_2,big_array_aa,dim1,dim2,dim3,dim4)
enddo
! Single and double beta excitations
! ==================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
kcol = psi_bilinear_matrix_columns(k_a)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
spindet(1:$N_int) = tmp_det(1:$N_int,2)
! Initial determinant is at k_b in beta-major representation
! -----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
! Loop inside the alpha row to gather all the connected betas
lrow = psi_bilinear_matrix_transp_rows(k_b)
l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
do i=1,N_det_beta_unique
if (l_b > N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l_b)
if (lrow /= krow) exit
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
idx(i) = l_b
l_b = l_b+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_b, doubles, n_singles_b, n_doubles )
! Compute Hij for all beta singles
! ----------------------------------
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
do i=1,n_singles_b
l_b = singles_b(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
l_a = psi_bilinear_matrix_transp_order(l_b)
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
enddo
! increment the alpha/beta part for single excitations
call off_diagonal_single_to_two_rdm_ab_dm(tmp_det, tmp_det2,c_1,c_2,big_array_ab,dim1,dim2,dim3,dim4)
! increment the beta /beta part for single excitations
call off_diagonal_single_to_two_rdm_bb_dm(tmp_det, tmp_det2,c_1,c_2,big_array_bb,dim1,dim2,dim3,dim4)
enddo
! Compute Hij for all beta doubles
! ----------------------------------
do i=1,n_doubles
l_b = doubles(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
l_a = psi_bilinear_matrix_transp_order(l_b)
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
enddo
call off_diagonal_double_to_two_rdm_bb_dm(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_1,c_2,big_array_bb,dim1,dim2,dim3,dim4)
ASSERT (l_a <= N_det)
enddo
! Diagonal contribution
! =====================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
double precision, external :: diag_wee_mat_elem, diag_S_mat_elem
double precision :: c_1(N_states),c_2(N_states)
do l = 1, N_states
c_1(l) = u_t(l,k_a)
enddo
call diagonal_contrib_to_all_two_rdm_dm(tmp_det,c_1,big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4)
end do
!!$OMP END DO
deallocate(buffer, singles_a, singles_b, doubles, idx)
!!$OMP END PARALLEL
end
SUBST [ N_int ]
1;;
2;;
3;;
4;;
N_int;;
END_TEMPLATE

0
src/two_body_rdm/all_states_routines.irp.f → src/two_body_rdm/all_states_act_2_rdm_dav_routines.irp.f
View File

34
src/two_body_rdm/all_states_2_rdm.irp.f → src/two_body_rdm/all_states_act_2_rdm_prov.irp.f
View File

@ -5,8 +5,11 @@
implicit none
double precision, allocatable :: state_weights(:)
BEGIN_DOC
! all_states_act_two_rdm_alpha_alpha_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-alpha electron pairs
! = <Psi| a^{\dagger}_i a^{\dagger}_j a_l a_k |Psi>
! all_states_act_two_rdm_alpha_alpha_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha electrons
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \alpha} |Psi>
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
allocate(state_weights(N_states))
state_weights = 1.d0/dble(N_states)
@ -20,11 +23,14 @@
BEGIN_PROVIDER [double precision, all_states_act_two_rdm_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)]
implicit none
double precision, allocatable :: state_weights(:)
BEGIN_DOC
! all_states_act_two_rdm_beta_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for beta-beta electron pairs
! = <Psi| a^{\dagger}_i a^{\dagger}_j a_l a_k |Psi>
! all_states_act_two_rdm_beta_beta_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of beta electrons
!
! <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi>
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
double precision, allocatable :: state_weights(:)
allocate(state_weights(N_states))
state_weights = 1.d0/dble(N_states)
integer :: ispin
@ -39,8 +45,11 @@
implicit none
double precision, allocatable :: state_weights(:)
BEGIN_DOC
! all_states_act_two_rdm_alpha_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs
! = <Psi| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi>
! all_states_act_two_rdm_alpha_beta_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM of alpha/beta electrons
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
allocate(state_weights(N_states))
state_weights = 1.d0/dble(N_states)
@ -61,10 +70,15 @@
BEGIN_PROVIDER [double precision, all_states_act_two_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb,N_states)]
implicit none
BEGIN_DOC
! all_states_act_two_rdm_spin_trace_mo(i,j,k,l) = state average physicist spin trace two-body rdm restricted to the ACTIVE indices
! The active part of the two-electron energy can be computed as:
! all_states_act_two_rdm_spin_trace_mo(i,j,k,l,istate) = STATE SPECIFIC physicist notation for 2RDM
!
! \sum_{i,j,k,l = 1, n_act_orb} all_states_act_two_rdm_spin_trace_mo(i,j,k,l) * < ii jj | kk ll >
! \sum_{\sigma, \sigma'} <Psi| a^{\dagger}_{i \sigma} a^{\dagger}_{j \sigma'} a_{l \sigma'} a_{k \sigma} |Psi>
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
!
! The active part of the two-electron energy for the state istate can be computed as:
!
! \sum_{i,j,k,l = 1, n_act_orb} all_states_act_two_rdm_spin_trace_mo(i,j,k,l,istate) * < ii jj | kk ll >
!
! with ii = list_act(i), jj = list_act(j), kk = list_act(k), ll = list_act(l)
END_DOC

0
src/two_body_rdm/compute_all_states.irp.f → src/two_body_rdm/all_states_act_2_rdm_update_routines.irp.f
View File

269
src/two_body_rdm/compute.irp.f
View File

@ -1,269 +0,0 @@
subroutine diagonal_contrib_to_two_rdm_ab_dm(det_1,c_1,big_array,dim1,dim2,dim3,dim4)
use bitmasks
BEGIN_DOC
! routine that update the DIAGONAL PART of the alpha/beta two body rdm IN CHEMIST NOTATIONS
END_DOC
implicit none
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states)
integer(bit_kind), intent(in) :: det_1(N_int,2)
double precision, intent(in) :: c_1(N_states)
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2)
integer :: i,j,h1,h2,istate
double precision :: c_1_bis
call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
do istate = 1, N_states
c_1_bis = c_1(istate) * c_1(istate)
do i = 1, n_occ_ab(1)
h1 = occ(i,1)
do j = 1, n_occ_ab(2)
h2 = occ(j,2)
big_array(h1,h1,h2,h2,istate) += c_1_bis
enddo
enddo
enddo
end
subroutine diagonal_contrib_to_all_two_rdm_dm(det_1,c_1,big_array_aa,big_array_bb,big_array_ab,dim1,dim2,dim3,dim4)
use bitmasks
BEGIN_DOC
! routine that update the DIAGONAL PART of ALL THREE two body rdm IN CHEMIST NOTATIONS
END_DOC
implicit none
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array_ab(dim1,dim2,dim3,dim4,N_states)
double precision, intent(inout) :: big_array_aa(dim1,dim2,dim3,dim4,N_states)
double precision, intent(inout) :: big_array_bb(dim1,dim2,dim3,dim4,N_states)
integer(bit_kind), intent(in) :: det_1(N_int,2)
double precision, intent(in) :: c_1(N_states)
integer :: occ(N_int*bit_kind_size,2)
integer :: n_occ_ab(2)
integer :: i,j,h1,h2,istate
double precision :: c_1_bis
BEGIN_DOC
! no factor 1/2 have to be taken into account as the permutations are already taken into account
END_DOC
call bitstring_to_list_ab(det_1, occ, n_occ_ab, N_int)
do istate = 1, N_states
c_1_bis = c_1(istate) * c_1(istate)
do i = 1, n_occ_ab(1)
h1 = occ(i,1)
do j = 1, n_occ_ab(2)
h2 = occ(j,2)
big_array_ab(h1,h1,h2,h2,istate) += c_1_bis
enddo
do j = 1, n_occ_ab(1)
h2 = occ(j,1)
big_array_aa(h1,h1,h2,h2,istate) += 0.5d0 * c_1_bis
big_array_aa(h1,h2,h2,h1,istate) -= 0.5d0 * c_1_bis
enddo
enddo
do i = 1, n_occ_ab(2)
h1 = occ(i,2)
do j = 1, n_occ_ab(2)
h2 = occ(j,2)
big_array_bb(h1,h1,h2,h2,istate) += 0.5d0 * c_1_bis
big_array_bb(h1,h2,h2,h1,istate) -= 0.5d0 * c_1_bis
enddo
enddo
enddo
end
subroutine off_diagonal_double_to_two_rdm_ab_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
use bitmasks
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for DOUBLE EXCITATIONS IN CHEMIST NOTATIONS
END_DOC
implicit none
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states)
integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
double precision, intent(in) :: c_1(N_states),c_2(N_states)
integer :: i,j,h1,h2,p1,p2,istate
integer :: exc(0:2,2,2)
double precision :: phase
call get_double_excitation(det_1,det_2,exc,phase,N_int)
h1 = exc(1,1,1)
h2 = exc(1,1,2)
p1 = exc(1,2,1)
p2 = exc(1,2,2)
do istate = 1, N_states
big_array(h1,p1,h2,p2,istate) += c_1(istate) * phase * c_2(istate)
! big_array(p1,h1,p2,h2,istate) += c_1(istate) * phase * c_2(istate)
enddo
end
subroutine off_diagonal_single_to_two_rdm_ab_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
use bitmasks
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the alpha/beta 2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS
END_DOC
implicit none
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states)
integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
double precision, intent(in) :: c_1(N_states),c_2(N_states)
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
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 (exc(0,1,1) == 1) then
! Mono alpha
h1 = exc(1,1,1)
p1 = exc(1,2,1)
do istate = 1, N_states
do i = 1, n_occ_ab(2)
h2 = occ(i,2)
big_array(h1,p1,h2,h2,istate) += 1.d0 * c_1(istate) * c_2(istate) * phase
enddo
enddo
else
! Mono beta
h1 = exc(1,1,2)
p1 = exc(1,2,2)
do istate = 1, N_states
do i = 1, n_occ_ab(1)
h2 = occ(i,1)
big_array(h2,h2,h1,p1,istate) += 1.d0 * c_1(istate) * c_2(istate) * phase
enddo
enddo
endif
end
subroutine off_diagonal_single_to_two_rdm_aa_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS
END_DOC
use bitmasks
implicit none
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states)
integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
double precision, intent(in) :: c_1(N_states),c_2(N_states)
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
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 (exc(0,1,1) == 1) then
! Mono alpha
h1 = exc(1,1,1)
p1 = exc(1,2,1)
do istate = 1, N_states
do i = 1, n_occ_ab(1)
h2 = occ(i,1)
big_array(h1,p1,h2,h2,istate) += 0.5d0 * c_1(istate) * c_2(istate) * phase
big_array(h1,h2,h2,p1,istate) -= 0.5d0 * c_1(istate) * c_2(istate) * phase
big_array(h2,h2,h1,p1,istate) += 0.5d0 * c_1(istate) * c_2(istate) * phase
big_array(h2,p1,h1,h2,istate) -= 0.5d0 * c_1(istate) * c_2(istate) * phase
enddo
enddo
else
return
endif
end
subroutine off_diagonal_single_to_two_rdm_bb_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
use bitmasks
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for SINGLE EXCITATIONS IN CHEMIST NOTATIONS
END_DOC
implicit none
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states)
integer(bit_kind), intent(in) :: det_1(N_int,2),det_2(N_int,2)
double precision, intent(in) :: c_1(N_states),c_2(N_states)
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
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 (exc(0,1,1) == 1) then
return
else
! Mono beta
h1 = exc(1,1,2)
p1 = exc(1,2,2)
do istate = 1, N_states
do i = 1, n_occ_ab(2)
h2 = occ(i,2)
big_array(h1,p1,h2,h2,istate) += 0.5d0 * c_1(istate) * c_2(istate) * phase
big_array(h1,h2,h2,p1,istate) -= 0.5d0 * c_1(istate) * c_2(istate) * phase
big_array(h2,h2,h1,p1,istate) += 0.5d0 * c_1(istate) * c_2(istate) * phase
big_array(h2,p1,h1,h2,istate) -= 0.5d0 * c_1(istate) * c_2(istate) * phase
enddo
enddo
endif
end
subroutine off_diagonal_double_to_two_rdm_aa_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
use bitmasks
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the alpha/alpha 2RDM only for DOUBLE EXCITATIONS IN CHEMIST NOTATIONS
END_DOC
implicit none
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states)
integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int)
double precision, intent(in) :: c_1(N_states),c_2(N_states)
integer :: i,j,h1,h2,p1,p2,istate
integer :: exc(0:2,2)
double precision :: phase
call get_double_excitation_spin(det_1,det_2,exc,phase,N_int)
h1 =exc(1,1)
h2 =exc(2,1)
p1 =exc(1,2)
p2 =exc(2,2)
!print*,'h1,p1,h2,p2',h1,p1,h2,p2,c_1(istate) * phase * c_2(istate)
do istate = 1, N_states
big_array(h1,p1,h2,p2,istate) += 0.5d0 * c_1(istate) * phase * c_2(istate)
big_array(h1,p2,h2,p1,istate) -= 0.5d0 * c_1(istate) * phase * c_2(istate)
big_array(h2,p2,h1,p1,istate) += 0.5d0 * c_1(istate) * phase * c_2(istate)
big_array(h2,p1,h1,p2,istate) -= 0.5d0 * c_1(istate) * phase * c_2(istate)
enddo
end
subroutine off_diagonal_double_to_two_rdm_bb_dm(det_1,det_2,c_1,c_2,big_array,dim1,dim2,dim3,dim4)
use bitmasks
BEGIN_DOC
! routine that update the OFF DIAGONAL PART of the beta /beta 2RDM only for DOUBLE EXCITATIONS
END_DOC
implicit none
integer, intent(in) :: dim1,dim2,dim3,dim4
double precision, intent(inout) :: big_array(dim1,dim2,dim3,dim4,N_states)
integer(bit_kind), intent(in) :: det_1(N_int),det_2(N_int)
double precision, intent(in) :: c_1(N_states),c_2(N_states)
integer :: i,j,h1,h2,p1,p2,istate
integer :: exc(0:2,2)
double precision :: phase
call get_double_excitation_spin(det_1,det_2,exc,phase,N_int)
h1 =exc(1,1)
h2 =exc(2,1)
p1 =exc(1,2)
p2 =exc(2,2)
!print*,'h1,p1,h2,p2',h1,p1,h2,p2,c_1(istate) * phase * c_2(istate)
do istate = 1, N_states
big_array(h1,p1,h2,p2,istate) += 0.5d0 * c_1(istate) * phase * c_2(istate)
big_array(h1,p2,h2,p1,istate) -= 0.5d0 * c_1(istate) * phase * c_2(istate)
big_array(h2,p2,h1,p1,istate) += 0.5d0 * c_1(istate) * phase * c_2(istate)
big_array(h2,p1,h1,p2,istate) -= 0.5d0 * c_1(istate) * phase * c_2(istate)
enddo
end

807
src/two_body_rdm/compute_orb_range_omp.irp.f
View File

@ -1,807 +0,0 @@
subroutine orb_range_diag_to_all_two_rdm_dm_buffer(det_1,c_1,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 rdms in a specific range of orbitals for a given determinant det_1
!
! c_1 is supposed to be a scalar quantity, such as state averaged coef of the determinant det_1
!
! 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
END_DOC
implicit none
integer, intent(in) :: ispin,sze_buff
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
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
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
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)
nkeys += 1
values(nkeys) = c_1
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = h1
keys(4,nkeys) = h2
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
values(nkeys) = 0.5d0 * c_1
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = h1
keys(4,nkeys) = h2
nkeys += 1
values(nkeys) = -0.5d0 * c_1
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
values(nkeys) = 0.5d0 * c_1
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = h1
keys(4,nkeys) = h2
nkeys += 1
values(nkeys) = -0.5d0 * c_1
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
values(nkeys) = 0.5d0 * c_1
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = h1
keys(4,nkeys) = h2
nkeys += 1
values(nkeys) = 0.5d0 * c_1
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
values(nkeys) = 0.5d0 * c_1
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = h1
keys(4,nkeys) = h2
nkeys += 1
values(nkeys) = -0.5d0 * c_1
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
values(nkeys) = 0.5d0 * c_1
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = h1
keys(4,nkeys) = h2
nkeys += 1
values(nkeys) = -0.5d0 * c_1
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_two_rdm_ab_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/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
!
! 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 == 3 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,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(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
values(nkeys) = c_1 * phase
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = p2
else if(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) = p1
keys(2,nkeys) = p2
keys(3,nkeys) = h1
keys(4,nkeys) = h2
endif
end
subroutine orb_range_off_diag_single_to_two_rdm_ab_dm_buffer(det_1,det_2,c_1,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
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 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 == 3 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(alpha_beta)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(2)
h2 = occ(i,2)
if(list_orb_reverse(h2).lt.0)return
h2 = list_orb_reverse(h2)
nkeys += 1
values(nkeys) = c_1 * phase
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = h2
enddo
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(1)
h2 = occ(i,1)
if(list_orb_reverse(h2).lt.0)return
h2 = list_orb_reverse(h2)
nkeys += 1
values(nkeys) = c_1 * phase
keys(1,nkeys) = h1
keys(2,nkeys) = h2
keys(3,nkeys) = p1
keys(4,nkeys) = h2
enddo
endif
else if(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(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) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo
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)
!print*,'****************'
!print*,'****************'
!print*,'h1,p1',h1,p1
do i = 1, n_occ_ab(1)
h2 = occ(i,1)
if(list_orb_reverse(h2).lt.0)return
h2 = list_orb_reverse(h2)
! print*,'h2 = ',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) = h2
keys(2,nkeys) = h1
keys(3,nkeys) = h2
keys(4,nkeys) = p1
enddo
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)
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_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),det_2(N_int)
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
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),det_2(N_int)
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

85
src/two_body_rdm/orb_range_omp.irp.f
View File

@ -1,85 +0,0 @@
BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_alpha_alpha_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
implicit none
double precision, allocatable :: state_weights(:)
BEGIN_DOC
! state_av_act_two_rdm_openmp_alpha_alpha_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-alpha electron pairs
! = <Psi| a^{\dagger}_i a^{\dagger}_j a_l a_k |Psi>
END_DOC
allocate(state_weights(N_states))
state_weights = state_average_weight
integer :: ispin
! condition for alpha/beta spin
ispin = 1
state_av_act_two_rdm_openmp_alpha_alpha_mo = 0.D0
call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_alpha_alpha_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
END_PROVIDER
BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_beta_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
implicit none
double precision, allocatable :: state_weights(:)
BEGIN_DOC
! state_av_act_two_rdm_openmp_beta_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for beta-beta electron pairs
! = <Psi| a^{\dagger}_i a^{\dagger}_j a_l a_k |Psi>
END_DOC
allocate(state_weights(N_states))
state_weights = state_average_weight
integer :: ispin
! condition for alpha/beta spin
ispin = 2
state_av_act_two_rdm_openmp_beta_beta_mo = 0.d0
call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_beta_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
END_PROVIDER
BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_alpha_beta_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
implicit none
double precision, allocatable :: state_weights(:)
BEGIN_DOC
! state_av_act_two_rdm_openmp_alpha_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs
! = <Psi| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi>
END_DOC
allocate(state_weights(N_states))
state_weights = state_average_weight
integer :: ispin
! condition for alpha/beta spin
print*,''
print*,''
print*,''
print*,'providint state_av_act_two_rdm_openmp_alpha_beta_mo '
ispin = 3
print*,'ispin = ',ispin
state_av_act_two_rdm_openmp_alpha_beta_mo = 0.d0
call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_alpha_beta_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
END_PROVIDER
BEGIN_PROVIDER [double precision, state_av_act_two_rdm_openmp_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
implicit none
BEGIN_DOC
! state_av_act_two_rdm_openmp_spin_trace_mo(i,j,k,l) = state average physicist spin trace two-body rdm restricted to the ACTIVE indices
! The active part of the two-electron energy can be computed as:
!
! \sum_{i,j,k,l = 1, n_act_orb} state_av_act_two_rdm_openmp_spin_trace_mo(i,j,k,l) * < ii jj | kk ll >
!
! with ii = list_act(i), jj = list_act(j), kk = list_act(k), ll = list_act(l)
END_DOC
double precision, allocatable :: state_weights(:)
allocate(state_weights(N_states))
state_weights = state_average_weight
integer :: ispin
! condition for alpha/beta spin
ispin = 4
state_av_act_two_rdm_openmp_spin_trace_mo = 0.d0
integer :: i
double precision :: wall_0,wall_1
call wall_time(wall_0)
print*,'providing the state average TWO-RDM ...'
call orb_range_two_rdm_state_av_openmp(state_av_act_two_rdm_openmp_spin_trace_mo,n_act_orb,n_act_orb,list_act,state_weights,ispin,psi_coef,size(psi_coef,2),size(psi_coef,1))
call wall_time(wall_1)
print*,'Time to provide the state average TWO-RDM',wall_1 - wall_0
END_PROVIDER

568
src/two_body_rdm/orb_range_routines_omp.irp.f
View File

@ -1,568 +0,0 @@
subroutine orb_range_two_rdm_state_av_openmp(big_array,dim1,norb,list_orb,state_weights,ispin,u_0,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! if ispin == 1 :: alpha/alpha 2rdm
! == 2 :: beta /beta 2rdm
! == 3 :: alpha/beta 2rdm
! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba))
!
! Assumes that the determinants are in psi_det
!
! istart, iend, ishift, istep are used in ZMQ parallelization.
END_DOC
integer, intent(in) :: N_st,sze
integer, intent(in) :: dim1,norb,list_orb(norb),ispin
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
double precision, intent(in) :: u_0(sze,N_st),state_weights(N_st)
integer :: k
double precision, allocatable :: u_t(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: u_t
allocate(u_t(N_st,N_det))
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order,N_det)
enddo
call dtranspose( &
u_0, &
size(u_0, 1), &
u_t, &
size(u_t, 1), &
N_det, N_st)
call orb_range_two_rdm_state_av_openmp_work(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,1,N_det,0,1)
deallocate(u_t)
do k=1,N_st
call dset_order(u_0(1,k),psi_bilinear_matrix_order_reverse,N_det)
enddo
end
subroutine orb_range_two_rdm_state_av_openmp_work(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
implicit none
BEGIN_DOC
! Computes two-rdm
!
! Default should be 1,N_det,0,1
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
integer, intent(in) :: dim1,norb,list_orb(norb),ispin
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
double precision, intent(in) :: u_t(N_st,N_det),state_weights(N_st)
integer :: k
PROVIDE N_int
select case (N_int)
case (1)
call orb_range_two_rdm_state_av_openmp_work_1(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
case (2)
call orb_range_two_rdm_state_av_openmp_work_2(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
case (3)
call orb_range_two_rdm_state_av_openmp_work_3(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
case (4)
call orb_range_two_rdm_state_av_openmp_work_4(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
case default
call orb_range_two_rdm_state_av_openmp_work_N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
end select
end
BEGIN_TEMPLATE
subroutine orb_range_two_rdm_state_av_openmp_work_$N_int(big_array,dim1,norb,list_orb,state_weights,ispin,u_t,N_st,sze,istart,iend,ishift,istep)
use bitmasks
use omp_lib
implicit none
BEGIN_DOC
! Computes the two rdm for the N_st vectors |u_t>
! if ispin == 1 :: alpha/alpha 2rdm
! == 2 :: beta /beta 2rdm
! == 3 :: alpha/beta 2rdm
! == 4 :: spin traced 2rdm :: aa + bb + 0.5 (ab + ba))
! The 2rdm will be computed only on the list of orbitals list_orb, which contains norb
! In any cases, the state average weights will be used with an array state_weights
! Default should be 1,N_det,0,1 for istart,iend,ishift,istep
END_DOC
integer, intent(in) :: N_st,sze,istart,iend,ishift,istep
double precision, intent(in) :: u_t(N_st,N_det),state_weights(N_st)
integer, intent(in) :: dim1,norb,list_orb(norb),ispin
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
integer(omp_lock_kind) :: lock_2rdm
integer :: i,j,k,l
integer :: k_a, k_b, l_a, l_b
integer :: krow, kcol
integer :: lrow, lcol
integer(bit_kind) :: spindet($N_int)
integer(bit_kind) :: tmp_det($N_int,2)
integer(bit_kind) :: tmp_det2($N_int,2)
integer(bit_kind) :: tmp_det3($N_int,2)
integer(bit_kind), allocatable :: buffer(:,:)
integer :: n_doubles
integer, allocatable :: doubles(:)
integer, allocatable :: singles_a(:)
integer, allocatable :: singles_b(:)
integer, allocatable :: idx(:), idx0(:)
integer :: maxab, n_singles_a, n_singles_b, kcol_prev
double precision :: c_average
logical :: alpha_alpha,beta_beta,alpha_beta,spin_trace
integer(bit_kind) :: orb_bitmask($N_int)
integer :: list_orb_reverse(mo_num)
integer, allocatable :: keys(:,:)
double precision, allocatable :: values(:)
integer :: nkeys,sze_buff
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.
else
print*,'Wrong parameter for ispin in general_two_rdm_state_av_openmp_work'
print*,'ispin = ',ispin
stop
endif
PROVIDE N_int
call list_to_bitstring( orb_bitmask, list_orb, norb, N_int)
sze_buff = norb ** 3 + 6 * norb
list_orb_reverse = -1000
do i = 1, norb
list_orb_reverse(list_orb(i)) = i
enddo
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
allocate(idx0(maxab))
do i=1,maxab
idx0(i) = i
enddo
call omp_init_lock(lock_2rdm)
! Prepare the array of all alpha single excitations
! -------------------------------------------------
PROVIDE N_int nthreads_davidson elec_alpha_num
!$OMP PARALLEL DEFAULT(NONE) NUM_THREADS(nthreads_davidson) &
!$OMP SHARED(psi_bilinear_matrix_rows, N_det,lock_2rdm,&
!$OMP psi_bilinear_matrix_columns, &
!$OMP psi_det_alpha_unique, psi_det_beta_unique,&
!$OMP n_det_alpha_unique, n_det_beta_unique, N_int,&
!$OMP psi_bilinear_matrix_transp_rows, &
!$OMP psi_bilinear_matrix_transp_columns, &
!$OMP psi_bilinear_matrix_transp_order, N_st, &
!$OMP psi_bilinear_matrix_order_transp_reverse, &
!$OMP psi_bilinear_matrix_columns_loc, &
!$OMP psi_bilinear_matrix_transp_rows_loc,elec_alpha_num, &
!$OMP istart, iend, istep, irp_here,list_orb_reverse, n_states, state_weights, dim1, &
!$OMP ishift, idx0, u_t, maxab, alpha_alpha,beta_beta,alpha_beta,spin_trace,ispin,big_array,sze_buff,orb_bitmask) &
!$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,c_1, c_2, &
!$OMP lcol, lrow, l_a, l_b, &
!$OMP buffer, doubles, n_doubles, &
!$OMP tmp_det2, idx, l, kcol_prev, &
!$OMP singles_a, n_singles_a, singles_b, &
!$OMP n_singles_b, nkeys, keys, values, c_average)
! Alpha/Beta double excitations
! =============================
nkeys = 0
allocate( keys(4,sze_buff), values(sze_buff))
allocate( buffer($N_int,maxab), &
singles_a(maxab), &
singles_b(maxab), &
doubles(maxab), &
idx(maxab))
kcol_prev=-1
ASSERT (iend <= N_det)
ASSERT (istart > 0)
ASSERT (istep > 0)
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
if (kcol /= kcol_prev) then
call get_all_spin_singles_$N_int( &
psi_det_beta_unique, idx0, &
tmp_det(1,2), N_det_beta_unique, &
singles_b, n_singles_b)
endif
kcol_prev = kcol
! Loop over singly excited beta columns
! -------------------------------------
do i=1,n_singles_b
lcol = singles_b(i)
tmp_det2(1:$N_int,2) = psi_det_beta_unique(1:$N_int, lcol)
l_a = psi_bilinear_matrix_columns_loc(lcol)
ASSERT (l_a <= N_det)
do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - l_a
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)
ASSERT (l_a <= N_det)
idx(j) = l_a
l_a = l_a+1
enddo
j = j-1
call get_all_spin_singles_$N_int( &
buffer, idx, tmp_det(1,1), j, &
singles_a, n_singles_a )
! Loop over alpha singles
! -----------------------
if(alpha_beta.or.spin_trace)then
do k = 1,n_singles_a
l_a = singles_a(k)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
c_average = 0.d0
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
c_average += c_1(l) * c_2(l) * state_weights(l)
enddo
if(alpha_beta)then
! only ONE contribution
if (nkeys+1 .ge. size(values)) then
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
nkeys = 0
endif
else if (spin_trace)then
! TWO contributions
if (nkeys+2 .ge. size(values)) then
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
nkeys = 0
endif
endif
call orb_range_off_diag_double_to_two_rdm_ab_dm_buffer(tmp_det,tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
enddo
endif
enddo
enddo
!$OMP END DO
!$OMP DO SCHEDULE(dynamic,64)
do k_a=istart+ishift,iend,istep
! Single and double alpha exitations
! ===================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
! Initial determinant is at k_b in beta-major representation
! ----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
spindet(1:$N_int) = tmp_det(1:$N_int,1)
! Loop inside the beta column to gather all the connected alphas
lcol = psi_bilinear_matrix_columns(k_a)
l_a = psi_bilinear_matrix_columns_loc(lcol)
do i=1,N_det_alpha_unique
if (l_a > N_det) exit
lcol = psi_bilinear_matrix_columns(l_a)
if (lcol /= kcol) exit
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
buffer(1:$N_int,i) = psi_det_alpha_unique(1:$N_int, lrow)
idx(i) = l_a
l_a = l_a+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_a, doubles, n_singles_a, n_doubles )
! Compute Hij for all alpha singles
! ----------------------------------
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
do i=1,n_singles_a
l_a = singles_a(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, lrow)
c_average = 0.d0
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
c_average += c_1(l) * c_2(l) * state_weights(l)
enddo
if(alpha_beta.or.spin_trace.or.alpha_alpha)then
! increment the alpha/beta part for single excitations
if (nkeys+ 2 * elec_alpha_num .ge. sze_buff) then
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
nkeys = 0
endif
call orb_range_off_diag_single_to_two_rdm_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
if (nkeys+4 * elec_alpha_num .ge. sze_buff ) then
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
nkeys = 0
endif
call orb_range_off_diag_single_to_two_rdm_aa_dm_buffer(tmp_det,tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
endif
enddo
! Compute Hij for all alpha doubles
! ----------------------------------
if(alpha_alpha.or.spin_trace)then
do i=1,n_doubles
l_a = doubles(i)
ASSERT (l_a <= N_det)
lrow = psi_bilinear_matrix_rows(l_a)
ASSERT (lrow <= N_det_alpha_unique)
c_average = 0.d0
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
c_average += c_1(l) * c_2(l) * state_weights(l)
enddo
if (nkeys+4 .ge. sze_buff) then
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
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
! Single and double beta excitations
! ==================================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
kcol = psi_bilinear_matrix_columns(k_a)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
spindet(1:$N_int) = tmp_det(1:$N_int,2)
! Initial determinant is at k_b in beta-major representation
! -----------------------------------------------------------------------
k_b = psi_bilinear_matrix_order_transp_reverse(k_a)
ASSERT (k_b <= N_det)
! Loop inside the alpha row to gather all the connected betas
lrow = psi_bilinear_matrix_transp_rows(k_b)
l_b = psi_bilinear_matrix_transp_rows_loc(lrow)
do i=1,N_det_beta_unique
if (l_b > N_det) exit
lrow = psi_bilinear_matrix_transp_rows(l_b)
if (lrow /= krow) exit
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
buffer(1:$N_int,i) = psi_det_beta_unique(1:$N_int, lcol)
idx(i) = l_b
l_b = l_b+1
enddo
i = i-1
call get_all_spin_singles_and_doubles_$N_int( &
buffer, idx, spindet, i, &
singles_b, doubles, n_singles_b, n_doubles )
! Compute Hij for all beta singles
! ----------------------------------
tmp_det2(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
do i=1,n_singles_b
l_b = singles_b(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
tmp_det2(1:$N_int,2) = psi_det_beta_unique (1:$N_int, lcol)
l_a = psi_bilinear_matrix_transp_order(l_b)
c_average = 0.d0
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
c_average += c_1(l) * c_2(l) * state_weights(l)
enddo
if(alpha_beta.or.spin_trace.or.beta_beta)then
! increment the alpha/beta part for single excitations
if (nkeys+2 * elec_alpha_num .ge. sze_buff ) then
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
nkeys = 0
endif
call orb_range_off_diag_single_to_two_rdm_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
if (nkeys+4 * elec_alpha_num .ge. sze_buff) then
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
nkeys = 0
endif
call orb_range_off_diag_single_to_two_rdm_bb_dm_buffer(tmp_det, tmp_det2,c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
endif
enddo
! Compute Hij for all beta doubles
! ----------------------------------
if(beta_beta.or.spin_trace)then
do i=1,n_doubles
l_b = doubles(i)
ASSERT (l_b <= N_det)
lcol = psi_bilinear_matrix_transp_columns(l_b)
ASSERT (lcol <= N_det_beta_unique)
l_a = psi_bilinear_matrix_transp_order(l_b)
c_average = 0.d0
do l= 1, N_states
c_1(l) = u_t(l,l_a)
c_2(l) = u_t(l,k_a)
c_average += c_1(l) * c_2(l) * state_weights(l)
enddo
if (nkeys+4 .ge. sze_buff) then
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
nkeys = 0
endif
call orb_range_off_diag_double_to_two_rdm_bb_dm_buffer(tmp_det(1,2),psi_det_beta_unique(1, lcol),c_average,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
ASSERT (l_a <= N_det)
enddo
endif
! Diagonal contribution
! =====================
! Initial determinant is at k_a in alpha-major representation
! -----------------------------------------------------------------------
krow = psi_bilinear_matrix_rows(k_a)
ASSERT (krow <= N_det_alpha_unique)
kcol = psi_bilinear_matrix_columns(k_a)
ASSERT (kcol <= N_det_beta_unique)
tmp_det(1:$N_int,1) = psi_det_alpha_unique(1:$N_int, krow)
tmp_det(1:$N_int,2) = psi_det_beta_unique (1:$N_int, kcol)
double precision, external :: diag_wee_mat_elem, diag_S_mat_elem
double precision :: c_1(N_states),c_2(N_states)
c_average = 0.d0
do l = 1, N_states
c_1(l) = u_t(l,k_a)
c_average += c_1(l) * c_1(l) * state_weights(l)
enddo
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
nkeys = 0
call orb_range_diag_to_all_two_rdm_dm_buffer(tmp_det,c_average,orb_bitmask,list_orb_reverse,ispin,sze_buff,nkeys,keys,values)
call update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
nkeys = 0
end do
!$OMP END DO
deallocate(buffer, singles_a, singles_b, doubles, idx, keys, values)
!$OMP END PARALLEL
end
SUBST [ N_int ]
1;;
2;;
3;;
4;;
N_int;;
END_TEMPLATE
subroutine update_keys_values(keys,values,nkeys,dim1,big_array,lock_2rdm)
use omp_lib
implicit none
integer, intent(in) :: nkeys,dim1
integer, intent(in) :: keys(4,nkeys)
double precision, intent(in) :: values(nkeys)
double precision, intent(inout) :: big_array(dim1,dim1,dim1,dim1)
integer(omp_lock_kind),intent(inout):: lock_2rdm
integer :: i,h1,h2,p1,p2
call omp_set_lock(lock_2rdm)
do i = 1, nkeys
h1 = keys(1,i)
h2 = keys(2,i)
p1 = keys(3,i)
p2 = keys(4,i)
big_array(h1,h2,p1,p2) += values(i)
enddo
call omp_unset_lock(lock_2rdm)
end

0
src/two_body_rdm/orb_range_routines.irp.f → src/two_body_rdm/state_av_act_2_rdm_dav_routines.irp.f
View File

28
src/two_body_rdm/orb_range.irp.f → src/two_body_rdm/state_av_act_2_rdm_prov.irp.f
View File

@ -5,8 +5,11 @@
implicit none
double precision, allocatable :: state_weights(:)
BEGIN_DOC
! state_av_act_two_rdm_alpha_alpha_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-alpha electron pairs
! = <Psi| a^{\dagger}_i a^{\dagger}_j a_l a_k |Psi>
! state_av_act_two_rdm_alpha_alpha_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of alpha electrons
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \alpha} a_{l \alpha} a_{k \alpha} |Psi>
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
allocate(state_weights(N_states))
state_weights = state_average_weight
@ -22,8 +25,11 @@
implicit none
double precision, allocatable :: state_weights(:)
BEGIN_DOC
! state_av_act_two_rdm_beta_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for beta-beta electron pairs
! = <Psi| a^{\dagger}_i a^{\dagger}_j a_l a_k |Psi>
! state_av_act_two_rdm_beta_beta_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of beta electrons
!
! <Psi| a^{\dagger}_{i \beta} a^{\dagger}_{j \beta} a_{l \beta} a_{k \beta} |Psi>
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
allocate(state_weights(N_states))
state_weights = state_average_weight
@ -39,8 +45,11 @@
implicit none
double precision, allocatable :: state_weights(:)
BEGIN_DOC
! state_av_act_two_rdm_alpha_beta_mo(i,j,k,l) = state average physicist two-body rdm restricted to the ACTIVE indices for alpha-beta electron pairs
! = <Psi| a^{\dagger}_{i,alpha} a^{\dagger}_{j,beta} a_{l,beta} a_{k,alpha} |Psi>
! state_av_act_two_rdm_alpha_beta_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM of alpha/beta electrons
!
! <Psi| a^{\dagger}_{i \alpha} a^{\dagger}_{j \beta} a_{l \beta} a_{k \alpha} |Psi>
!
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
allocate(state_weights(N_states))
state_weights = state_average_weight
@ -61,12 +70,11 @@
BEGIN_PROVIDER [double precision, state_av_act_two_rdm_spin_trace_mo, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)]
implicit none
BEGIN_DOC
! state_av_act_two_rdm_spin_trace_mo(i,j,k,l) = state average physicist spin trace two-body rdm restricted to the ACTIVE indices
! The active part of the two-electron energy can be computed as:
! state_av_act_two_rdm_spin_trace_mo(i,j,k,l) = STATE AVERAGE physicist notation for 2RDM
!
! \sum_{i,j,k,l = 1, n_act_orb} state_av_act_two_rdm_spin_trace_mo(i,j,k,l) * < ii jj | kk ll >
! \sum_{\sigma, \sigma'} <Psi| a^{\dagger}_{i \sigma} a^{\dagger}_{j \sigma'} a_{l \sigma'} a_{k \sigma} |Psi>
!
! with ii = list_act(i), jj = list_act(j), kk = list_act(k), ll = list_act(l)
! !!!!! WARNING !!!!! ALL SLATER DETERMINANTS IN PSI_DET MUST BELONG TO AN ACTIVE SPACE DEFINED BY "list_act"
END_DOC
double precision, allocatable :: state_weights(:)
allocate(state_weights(N_states))

0
src/two_body_rdm/compute_orb_range.irp.f → src/two_body_rdm/state_av_act_2_rdm_update_routines.irp.f
View File

62
src/two_body_rdm/two_rdm.irp.f
View File

@ -1,62 +0,0 @@
BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
&BEGIN_PROVIDER [double precision, two_rdm_alpha_alpha_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
&BEGIN_PROVIDER [double precision, two_rdm_beta_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
implicit none
BEGIN_DOC
! two_rdm_alpha_beta(i,j,k,l) = <Psi| a^{dagger}_{j,alpha} a^{dagger}_{l,beta} a_{k,beta} a_{i,alpha} | Psi>
! 1 1 2 2 = chemist notations
! note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry
!
END_DOC
integer :: dim1,dim2,dim3,dim4
double precision :: cpu_0,cpu_1
dim1 = mo_num
dim2 = mo_num
dim3 = mo_num
dim4 = mo_num
two_rdm_alpha_beta_mo = 0.d0
two_rdm_alpha_alpha_mo= 0.d0
two_rdm_beta_beta_mo = 0.d0
print*,'providing two_rdm_alpha_beta ...'
call wall_time(cpu_0)
call all_two_rdm_dm_nstates(two_rdm_alpha_alpha_mo,two_rdm_beta_beta_mo,two_rdm_alpha_beta_mo,dim1,dim2,dim3,dim4,psi_coef,size(psi_coef,2),size(psi_coef,1))
call wall_time(cpu_1)
print*,'two_rdm_alpha_beta provided in',dabs(cpu_1-cpu_0)
END_PROVIDER
BEGIN_PROVIDER [double precision, two_rdm_alpha_beta_mo_physicist, (mo_num,mo_num,mo_num,mo_num,N_states)]
&BEGIN_PROVIDER [double precision, two_rdm_alpha_alpha_mo_physicist, (mo_num,mo_num,mo_num,mo_num,N_states)]
&BEGIN_PROVIDER [double precision, two_rdm_beta_beta_mo_physicist, (mo_num,mo_num,mo_num,mo_num,N_states)]
implicit none
BEGIN_DOC
! two_rdm_alpha_beta_mo_physicist,(i,j,k,l) = <Psi| a^{dagger}_{k,alpha} a^{dagger}_{l,beta} a_{j,beta} a_{i,alpha} | Psi>
! 1 2 1 2 = physicist notations
! note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry
!
END_DOC
integer :: i,j,k,l,istate
double precision :: cpu_0,cpu_1
two_rdm_alpha_beta_mo_physicist = 0.d0
print*,'providing two_rdm_alpha_beta_mo_physicist ...'
call wall_time(cpu_0)
do istate = 1, N_states
do i = 1, mo_num
do j = 1, mo_num
do k = 1, mo_num
do l = 1, mo_num
! 1 2 1 2 1 1 2 2
two_rdm_alpha_beta_mo_physicist(l,k,i,j,istate) = two_rdm_alpha_beta_mo(i,l,j,k,istate)
two_rdm_alpha_alpha_mo_physicist(l,k,i,j,istate) = two_rdm_alpha_alpha_mo(i,l,j,k,istate)
two_rdm_beta_beta_mo_physicist(l,k,i,j,istate) = two_rdm_beta_beta_mo(i,l,j,k,istate)
enddo
enddo
enddo
enddo
enddo
call wall_time(cpu_1)
print*,'two_rdm_alpha_beta_mo_physicist provided in',dabs(cpu_1-cpu_0)
END_PROVIDER