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

indentation

Browse Source
This commit is contained in:
Anthony Scemama 2019-06-27 21:41:17 +02:00
parent 9bb66d5b3a
commit 3e38912dcb
4 changed files with 516 additions and 513 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
src/casscf/test_two_rdm.irp.f
View File

@ -8,11 +8,11 @@ program print_two_rdm
double precision :: accu,twodm double precision :: accu,twodm
accu = 0.d0 accu = 0.d0
do i=1,mo_num do i=1,n_act_orb
do j=1,mo_num do j=1,n_act_orb
do k=1,mo_num do k=1,n_act_orb
do l=1,mo_num do l=1,n_act_orb
twodm = coussin_peter_two_rdm_mo(i,j,k,l,1) twodm = coussin_peter_two_rdm_mo(list_act(i),list_act(j),list_act(k),list_act(l),1)
if(dabs(twodm - P0tuvx(i,j,k,l)).gt.thr)then if(dabs(twodm - P0tuvx(i,j,k,l)).gt.thr)then
print*,'' print*,''
print*,'sum' print*,'sum'

4
src/two_body_rdm/README.rst
View File

@ -2,5 +2,7 @@
two_body_rdm two_body_rdm
============ ============
Contains the two rdms (aa,bb,ab) stored as plain arrays Contains the two rdms $\alpha\alpha$, $\beta\beta$ and $\alpha\beta$ stored as
maps, with pysicists notation, consistent with the two-electron integrals in the
MO basis.

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

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

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

@ -1,84 +1,86 @@
BEGIN_PROVIDER [double precision, coussin_peter_two_rdm_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
BEGIN_PROVIDER [double precision, coussin_peter_two_rdm_mo, (mo_num,mo_num,mo_num,mo_num,N_states)] implicit none
implicit none BEGIN_DOC
BEGIN_DOC ! coussin_peter_two_rdm_mo(i,j,k,l) = the two rdm that peter wants for his CASSCF
! coussin_peter_two_rdm_mo(i,j,k,l) = the two rdm that peter wants for his CASSCF END_DOC
END_DOC integer :: i,j,k,l, istate
integer :: i,j,k,l do istate = 1,N_states
do l = 1, mo_num do l = 1, mo_num
do k = 1, mo_num do k = 1, mo_num
do j = 1, mo_num do j = 1, mo_num
do i = 1, mo_num do i = 1, mo_num
coussin_peter_two_rdm_mo(i,j,k,l,:) = 0.5d0 * (two_rdm_alpha_beta_mo(i,j,k,l,:) + two_rdm_alpha_beta_mo(i,j,k,l,:)) & coussin_peter_two_rdm_mo (i,j,k,l,istate) = &
+ two_rdm_alpha_alpha_mo(i,j,k,l,:) & two_rdm_alpha_beta_mo (i,j,k,l,istate) + &
+ two_rdm_beta_beta_mo(i,j,k,l,:) two_rdm_alpha_alpha_mo(i,j,k,l,istate) + &
two_rdm_beta_beta_mo (i,j,k,l,istate)
enddo
enddo
enddo
enddo enddo
enddo
enddo enddo
enddo
END_PROVIDER END_PROVIDER
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_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_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)] &BEGIN_PROVIDER [double precision, two_rdm_beta_beta_mo, (mo_num,mo_num,mo_num,mo_num,N_states)]
implicit none implicit none
BEGIN_DOC 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> ! 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 ! 1 1 2 2 = chemist notations
! note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry ! note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry
! !
END_DOC END_DOC
integer :: dim1,dim2,dim3,dim4 integer :: dim1,dim2,dim3,dim4
double precision :: cpu_0,cpu_1 double precision :: cpu_0,cpu_1
dim1 = mo_num dim1 = mo_num
dim2 = mo_num dim2 = mo_num
dim3 = mo_num dim3 = mo_num
dim4 = mo_num dim4 = mo_num
two_rdm_alpha_beta_mo = 0.d0 two_rdm_alpha_beta_mo = 0.d0
two_rdm_alpha_alpha_mo= 0.d0 two_rdm_alpha_alpha_mo= 0.d0
two_rdm_beta_beta_mo = 0.d0 two_rdm_beta_beta_mo = 0.d0
print*,'providing two_rdm_alpha_beta ...' print*,'providing two_rdm_alpha_beta ...'
call wall_time(cpu_0) call wall_time(cpu_0)
call all_two_rdm_dm_nstates_openmp(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 all_two_rdm_dm_nstates_openmp(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) call wall_time(cpu_1)
print*,'two_rdm_alpha_beta provided in',dabs(cpu_1-cpu_0) print*,'two_rdm_alpha_beta provided in',dabs(cpu_1-cpu_0)
END_PROVIDER 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_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_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)] &BEGIN_PROVIDER [double precision, two_rdm_beta_beta_mo_physicist, (mo_num,mo_num,mo_num,mo_num,N_states)]
implicit none implicit none
BEGIN_DOC 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> ! 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 ! 1 2 1 2 = physicist notations
! note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry ! note that no 1/2 factor is introduced in order to take into acccount for the spin symmetry
! !
END_DOC END_DOC
integer :: i,j,k,l,istate integer :: i,j,k,l,istate
double precision :: cpu_0,cpu_1 double precision :: cpu_0,cpu_1
two_rdm_alpha_beta_mo_physicist = 0.d0 two_rdm_alpha_beta_mo_physicist = 0.d0
print*,'providing two_rdm_alpha_beta_mo_physicist ...' print*,'providing two_rdm_alpha_beta_mo_physicist ...'
call wall_time(cpu_0) call wall_time(cpu_0)
do istate = 1, N_states do istate = 1, N_states
do i = 1, mo_num do i = 1, mo_num
do j = 1, mo_num do j = 1, mo_num
do k = 1, mo_num do k = 1, mo_num
do l = 1, mo_num do l = 1, mo_num
! 1 2 1 2 1 1 2 2 ! 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_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_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) 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 enddo
enddo
enddo enddo
enddo call wall_time(cpu_1)
call wall_time(cpu_1) print*,'two_rdm_alpha_beta_mo_physicist provided in',dabs(cpu_1-cpu_0)
print*,'two_rdm_alpha_beta_mo_physicist provided in',dabs(cpu_1-cpu_0)
END_PROVIDER
END_PROVIDER