mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-22 03:23:29 +01:00
normal ordering with DGEMM: OK
This commit is contained in:
Abdallah Ammar
parent
aafca191f1
commit
374a88bc62
@ -11,16 +11,15 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
|
|||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
|
|
||||||
integer :: i, h1, p1, h2, p2
|
integer :: i, ii, h1, p1, h2, p2, ipoint
|
||||||
integer :: hh1, hh2, pp1, pp2
|
integer :: hh1, hh2, pp1, pp2
|
||||||
integer :: Ne(2)
|
integer :: Ne(2)
|
||||||
double precision :: hthree_aaa, hthree_aab
|
double precision :: wall0, wall1, walli, wallf
|
||||||
double precision :: wall0, wall1
|
|
||||||
integer, allocatable :: occ(:,:)
|
integer, allocatable :: occ(:,:)
|
||||||
integer(bit_kind), allocatable :: key_i_core(:,:)
|
integer(bit_kind), allocatable :: key_i_core(:,:)
|
||||||
|
|
||||||
print*,' Providing normal_two_body_bi_orth ...'
|
print*,' Providing normal_two_body_bi_orth ...'
|
||||||
call wall_time(wall0)
|
call wall_time(walli)
|
||||||
|
|
||||||
if(read_tc_norm_ord) then
|
if(read_tc_norm_ord) then
|
||||||
|
|
||||||
@ -30,6 +29,11 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
|
|||||||
|
|
||||||
else
|
else
|
||||||
|
|
||||||
|
double precision, allocatable :: tmp_2d(:,:), tmp_3d(:,:,:)
|
||||||
|
double precision, allocatable :: tmp1(:,:,:), tmp2(:,:), tmp3(:,:,:)
|
||||||
|
double precision, allocatable :: tmpval_1(:), tmpval_2(:), tmpvec_1(:,:), tmpvec_2(:,:), tmpvec_3(:,:)
|
||||||
|
double precision, allocatable :: tmp(:,:,:,:)
|
||||||
|
|
||||||
PROVIDE N_int
|
PROVIDE N_int
|
||||||
|
|
||||||
allocate( occ(N_int*bit_kind_size,2) )
|
allocate( occ(N_int*bit_kind_size,2) )
|
||||||
@ -45,87 +49,15 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
|
|||||||
call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
|
call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
|
||||||
endif
|
endif
|
||||||
|
|
||||||
PROVIDE no_aba_contraction
|
allocate(tmp(mo_num,mo_num,mo_num,mo_num))
|
||||||
PROVIDE no_aab_contraction
|
|
||||||
PROVIDE no_aaa_contraction
|
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
! ---
|
||||||
!$OMP DEFAULT (NONE) &
|
! aba contraction
|
||||||
!$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, hthree_aab, hthree_aaa) &
|
|
||||||
!$OMP SHARED (N_int, n_act_orb, list_act, Ne, occ, normal_two_body_bi_orth, &
|
|
||||||
!$OMP no_aba_contraction, no_aab_contraction, no_aaa_contraction)
|
|
||||||
!$OMP DO SCHEDULE (static)
|
|
||||||
do hh1 = 1, n_act_orb
|
|
||||||
h1 = list_act(hh1)
|
|
||||||
|
|
||||||
do pp1 = 1, n_act_orb
|
print*,' Providing aba_contraction ...'
|
||||||
p1 = list_act(pp1)
|
|
||||||
|
|
||||||
do hh2 = 1, n_act_orb
|
|
||||||
h2 = list_act(hh2)
|
|
||||||
|
|
||||||
do pp2 = 1, n_act_orb
|
|
||||||
p2 = list_act(pp2)
|
|
||||||
|
|
||||||
normal_two_body_bi_orth(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) + no_aab_contraction(p2,h2,p1,h1) + no_aaa_contraction(p2,h2,p1,h1)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
!$OMP END DO
|
|
||||||
!$OMP END PARALLEL
|
|
||||||
|
|
||||||
deallocate( occ )
|
|
||||||
deallocate( key_i_core )
|
|
||||||
endif
|
|
||||||
|
|
||||||
if(write_tc_norm_ord.and.mpi_master) then
|
|
||||||
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/normal_two_body_bi_orth', action="write")
|
|
||||||
call ezfio_set_work_empty(.False.)
|
|
||||||
write(11) normal_two_body_bi_orth
|
|
||||||
close(11)
|
|
||||||
call ezfio_set_tc_keywords_io_tc_integ('Read')
|
|
||||||
endif
|
|
||||||
|
|
||||||
call wall_time(wall1)
|
|
||||||
print*,' Wall time for normal_two_body_bi_orth ', wall1-wall0
|
|
||||||
|
|
||||||
END_PROVIDER
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_num)]
|
|
||||||
|
|
||||||
use bitmasks ! you need to include the bitmasks_module.f90 features
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, ii, h1, p1, h2, p2, ipoint
|
|
||||||
integer :: Ne(2)
|
|
||||||
double precision :: wall0, wall1
|
|
||||||
integer, allocatable :: occ(:,:)
|
|
||||||
integer(bit_kind), allocatable :: key_i_core(:,:)
|
|
||||||
double precision, allocatable :: tmp_3d(:,:,:)
|
|
||||||
double precision, allocatable :: tmp1(:,:,:), tmp2(:,:)
|
|
||||||
double precision, allocatable :: tmpval_1(:), tmpval_2(:), tmpvec_1(:,:), tmpvec_2(:,:)
|
|
||||||
double precision, allocatable :: tmp_2d(:,:)
|
|
||||||
|
|
||||||
print*,' Providing no_aba_contraction ...'
|
|
||||||
call wall_time(wall0)
|
call wall_time(wall0)
|
||||||
|
|
||||||
PROVIDE N_int
|
tmp = 0.d0
|
||||||
|
|
||||||
allocate(occ(N_int*bit_kind_size,2))
|
|
||||||
allocate(key_i_core(N_int,2))
|
|
||||||
|
|
||||||
if(core_tc_op) then
|
|
||||||
do i = 1, N_int
|
|
||||||
key_i_core(i,1) = xor(ref_bitmask(i,1), core_bitmask(i,1))
|
|
||||||
key_i_core(i,2) = xor(ref_bitmask(i,2), core_bitmask(i,2))
|
|
||||||
enddo
|
|
||||||
call bitstring_to_list_ab(key_i_core, occ, Ne, N_int)
|
|
||||||
else
|
|
||||||
call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
|
|
||||||
endif
|
|
||||||
|
|
||||||
allocate(tmp_3d(mo_num,mo_num,mo_num))
|
allocate(tmp_3d(mo_num,mo_num,mo_num))
|
||||||
allocate(tmp1(n_points_final_grid,3,mo_num))
|
allocate(tmp1(n_points_final_grid,3,mo_num))
|
||||||
@ -136,7 +68,6 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
allocate(tmpvec_2(n_points_final_grid,3))
|
allocate(tmpvec_2(n_points_final_grid,3))
|
||||||
allocate(tmp_2d(mo_num,mo_num))
|
allocate(tmp_2d(mo_num,mo_num))
|
||||||
|
|
||||||
|
|
||||||
! purely closed shell part
|
! purely closed shell part
|
||||||
do ii = 1, Ne(2)
|
do ii = 1, Ne(2)
|
||||||
i = occ(ii,2)
|
i = occ(ii,2)
|
||||||
@ -195,7 +126,7 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
do p1 = 1, mo_num
|
do p1 = 1, mo_num
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aba_contraction(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -246,7 +177,7 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aba_contraction(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END PARALLEL DO
|
!$OMP END PARALLEL DO
|
||||||
@ -255,7 +186,6 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
enddo ! h1
|
enddo ! h1
|
||||||
enddo ! i
|
enddo ! i
|
||||||
|
|
||||||
|
|
||||||
! purely open-shell part
|
! purely open-shell part
|
||||||
if(Ne(2) < Ne(1)) then
|
if(Ne(2) < Ne(1)) then
|
||||||
do ii = Ne(2) + 1, Ne(1)
|
do ii = Ne(2) + 1, Ne(1)
|
||||||
@ -313,7 +243,7 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
do p1 = 1, mo_num
|
do p1 = 1, mo_num
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aba_contraction(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -363,7 +293,7 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aba_contraction(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END PARALLEL DO
|
!$OMP END PARALLEL DO
|
||||||
@ -373,53 +303,30 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
enddo !i
|
enddo !i
|
||||||
endif
|
endif
|
||||||
|
|
||||||
deallocate(tmp_2d, tmp_3d)
|
deallocate(tmp_3d)
|
||||||
deallocate(tmp1, tmp2)
|
deallocate(tmp1)
|
||||||
deallocate(tmpval_1, tmpval_2)
|
deallocate(tmp2)
|
||||||
deallocate(tmpvec_1, tmpvec_2)
|
deallocate(tmpval_1)
|
||||||
|
deallocate(tmpval_2)
|
||||||
|
deallocate(tmpvec_1)
|
||||||
|
deallocate(tmpvec_2)
|
||||||
|
deallocate(tmp_2d)
|
||||||
|
|
||||||
no_aba_contraction = -0.5d0 * no_aba_contraction
|
tmp = -0.5d0 * tmp
|
||||||
call sum_A_At(no_aba_contraction(1,1,1,1), mo_num*mo_num)
|
call sum_A_At(tmp(1,1,1,1), mo_num*mo_num)
|
||||||
|
|
||||||
call wall_time(wall1)
|
call wall_time(wall1)
|
||||||
print*,' Wall time for no_aba_contraction', wall1-wall0
|
print*,' Wall time for aba_contraction', wall1-wall0
|
||||||
|
|
||||||
END_PROVIDER
|
normal_two_body_bi_orth = tmp
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
! aab contraction
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, no_aab_contraction, (mo_num,mo_num,mo_num,mo_num)]
|
print*,' Providing aab_contraction ...'
|
||||||
|
|
||||||
use bitmasks ! you need to include the bitmasks_module.f90 features
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, ii, h1, p1, h2, p2, ipoint
|
|
||||||
integer :: Ne(2)
|
|
||||||
double precision :: wall0, wall1
|
|
||||||
integer, allocatable :: occ(:,:)
|
|
||||||
integer(bit_kind), allocatable :: key_i_core(:,:)
|
|
||||||
double precision, allocatable :: tmp_3d(:,:,:)
|
|
||||||
double precision, allocatable :: tmp1(:,:,:), tmp2(:,:)
|
|
||||||
double precision, allocatable :: tmpval_1(:), tmpvec_1(:,:)
|
|
||||||
double precision, allocatable :: tmp_2d(:,:)
|
|
||||||
|
|
||||||
print*,' Providing no_aab_contraction ...'
|
|
||||||
call wall_time(wall0)
|
call wall_time(wall0)
|
||||||
|
|
||||||
PROVIDE N_int
|
tmp = 0.d0
|
||||||
|
|
||||||
allocate(occ(N_int*bit_kind_size,2))
|
|
||||||
allocate(key_i_core(N_int,2))
|
|
||||||
|
|
||||||
if(core_tc_op) then
|
|
||||||
do i = 1, N_int
|
|
||||||
key_i_core(i,1) = xor(ref_bitmask(i,1), core_bitmask(i,1))
|
|
||||||
key_i_core(i,2) = xor(ref_bitmask(i,2), core_bitmask(i,2))
|
|
||||||
enddo
|
|
||||||
call bitstring_to_list_ab(key_i_core, occ, Ne, N_int)
|
|
||||||
else
|
|
||||||
call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
|
|
||||||
endif
|
|
||||||
|
|
||||||
allocate(tmp_2d(mo_num,mo_num))
|
allocate(tmp_2d(mo_num,mo_num))
|
||||||
allocate(tmp_3d(mo_num,mo_num,mo_num))
|
allocate(tmp_3d(mo_num,mo_num,mo_num))
|
||||||
@ -428,7 +335,6 @@ BEGIN_PROVIDER [ double precision, no_aab_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
allocate(tmpval_1(n_points_final_grid))
|
allocate(tmpval_1(n_points_final_grid))
|
||||||
allocate(tmpvec_1(n_points_final_grid,3))
|
allocate(tmpvec_1(n_points_final_grid,3))
|
||||||
|
|
||||||
|
|
||||||
! purely closed shell part
|
! purely closed shell part
|
||||||
do ii = 1, Ne(2)
|
do ii = 1, Ne(2)
|
||||||
i = occ(ii,2)
|
i = occ(ii,2)
|
||||||
@ -480,7 +386,7 @@ BEGIN_PROVIDER [ double precision, no_aab_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
do p1 = 1, mo_num
|
do p1 = 1, mo_num
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aab_contraction(p2,h2,p1,h1) = no_aab_contraction(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -528,7 +434,7 @@ BEGIN_PROVIDER [ double precision, no_aab_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aab_contraction(p2,h2,p1,h1) = no_aab_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END PARALLEL DO
|
!$OMP END PARALLEL DO
|
||||||
@ -537,24 +443,26 @@ BEGIN_PROVIDER [ double precision, no_aab_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
enddo ! h1
|
enddo ! h1
|
||||||
enddo ! i
|
enddo ! i
|
||||||
|
|
||||||
|
deallocate(tmp_2d)
|
||||||
deallocate(tmp_3d)
|
deallocate(tmp_3d)
|
||||||
deallocate(tmp1, tmp2)
|
deallocate(tmp1)
|
||||||
|
deallocate(tmp2)
|
||||||
deallocate(tmpval_1)
|
deallocate(tmpval_1)
|
||||||
deallocate(tmpvec_1)
|
deallocate(tmpvec_1)
|
||||||
|
|
||||||
no_aab_contraction = -0.5d0 * no_aab_contraction
|
tmp = -0.5d0 * tmp
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP DEFAULT (NONE) &
|
||||||
!$OMP PRIVATE (h1, h2, p1, p2) &
|
!$OMP PRIVATE (h1, h2, p1, p2) &
|
||||||
!$OMP SHARED (no_aab_contraction, mo_num)
|
!$OMP SHARED (tmp, mo_num)
|
||||||
|
|
||||||
!$OMP DO
|
!$OMP DO
|
||||||
do h1 = 1, mo_num
|
do h1 = 1, mo_num
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p1 = 1, mo_num
|
do p1 = 1, mo_num
|
||||||
do p2 = p1, mo_num
|
do p2 = p1, mo_num
|
||||||
no_aab_contraction(p2,h2,p1,h1) -= no_aab_contraction(p1,h2,p2,h1)
|
tmp(p2,h2,p1,h1) -= tmp(p1,h2,p2,h1)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -566,7 +474,7 @@ BEGIN_PROVIDER [ double precision, no_aab_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p1 = 2, mo_num
|
do p1 = 2, mo_num
|
||||||
do p2 = 1, p1-1
|
do p2 = 1, p1-1
|
||||||
no_aab_contraction(p2,h2,p1,h1) = -no_aab_contraction(p1,h2,p2,h1)
|
tmp(p2,h2,p1,h1) = -tmp(p1,h2,p2,h1)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -578,7 +486,7 @@ BEGIN_PROVIDER [ double precision, no_aab_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
do h2 = h1+1, mo_num
|
do h2 = h1+1, mo_num
|
||||||
do p1 = 2, mo_num
|
do p1 = 2, mo_num
|
||||||
do p2 = 1, p1-1
|
do p2 = 1, p1-1
|
||||||
no_aab_contraction(p2,h2,p1,h1) *= -1.d0
|
tmp(p2,h2,p1,h1) *= -1.d0
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -586,71 +494,19 @@ BEGIN_PROVIDER [ double precision, no_aab_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
call wall_time(wall1)
|
call wall_time(wall1)
|
||||||
print*,' Wall time for no_aab_contraction', wall1-wall0
|
print*,' Wall time for aab_contraction', wall1-wall0
|
||||||
|
|
||||||
END_PROVIDER
|
normal_two_body_bi_orth += tmp
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
! aaa contraction
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_num)]
|
if(Ne(2) .ge. 3) then
|
||||||
|
|
||||||
BEGIN_DOC
|
print*,' Providing aaa_contraction ...'
|
||||||
!
|
|
||||||
! if:
|
|
||||||
! h1 < h2
|
|
||||||
! p1 > p2
|
|
||||||
!
|
|
||||||
! no_aaa_contraction(p2,h2.p1,h1) = 0.5 [Ialpha(p2,h1,p1,h2) + Ibeta(p2,h1,p1,h2)]
|
|
||||||
! = -0.5 [Ialpha(p2,h2,p1,h1) + Ibeta(p2,h2,p1,h1)]
|
|
||||||
!
|
|
||||||
! else:
|
|
||||||
!
|
|
||||||
! no_aaa_contraction(p2,h2.p1,h1) = 0.5 [Ialpha(p2,h2,p1,h1) + Ibeta(p2,h2,p1,h1)]
|
|
||||||
!
|
|
||||||
!
|
|
||||||
! I(p2,h2,p1,h1) = J(p2,h2,p1,h1) - J(p1,h2,p2,h1)
|
|
||||||
! J(p2,h2,p1,h1) = \sum_i [ < i p2 p1 | i h2 h1 >
|
|
||||||
! + < p2 p1 i | i h2 h1 >
|
|
||||||
! + < p1 i p2 | i h2 h1 > ]
|
|
||||||
!
|
|
||||||
!
|
|
||||||
END_DOC
|
|
||||||
|
|
||||||
use bitmasks ! you need to include the bitmasks_module.f90 features
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
integer :: i, ii, h1, p1, h2, p2, ipoint
|
|
||||||
integer :: Ne(2)
|
|
||||||
double precision :: wall0, wall1
|
|
||||||
integer, allocatable :: occ(:,:)
|
|
||||||
integer(bit_kind), allocatable :: key_i_core(:,:)
|
|
||||||
double precision, allocatable :: tmp_2d(:,:), tmp_3d(:,:,:)
|
|
||||||
double precision, allocatable :: tmp1(:,:,:), tmp2(:,:), tmp3(:,:,:)
|
|
||||||
double precision, allocatable :: tmpval_1(:), tmpval_2(:), tmpvec_1(:,:), tmpvec_2(:,:), tmpvec_3(:,:)
|
|
||||||
|
|
||||||
print*,' Providing no_aaa_contraction ...'
|
|
||||||
call wall_time(wall0)
|
call wall_time(wall0)
|
||||||
|
|
||||||
PROVIDE N_int
|
tmp = 0.d0
|
||||||
|
|
||||||
allocate(occ(N_int*bit_kind_size,2))
|
|
||||||
allocate(key_i_core(N_int,2))
|
|
||||||
|
|
||||||
if(core_tc_op) then
|
|
||||||
do i = 1, N_int
|
|
||||||
key_i_core(i,1) = xor(ref_bitmask(i,1), core_bitmask(i,1))
|
|
||||||
key_i_core(i,2) = xor(ref_bitmask(i,2), core_bitmask(i,2))
|
|
||||||
enddo
|
|
||||||
call bitstring_to_list_ab(key_i_core, occ, Ne, N_int)
|
|
||||||
else
|
|
||||||
call bitstring_to_list_ab(ref_bitmask, occ, Ne, N_int)
|
|
||||||
endif
|
|
||||||
|
|
||||||
if(Ne(2) .lt. 3) then
|
|
||||||
|
|
||||||
no_aaa_contraction = 0.d0
|
|
||||||
|
|
||||||
else
|
|
||||||
|
|
||||||
allocate(tmp_2d(mo_num,mo_num))
|
allocate(tmp_2d(mo_num,mo_num))
|
||||||
allocate(tmp_3d(mo_num,mo_num,mo_num))
|
allocate(tmp_3d(mo_num,mo_num,mo_num))
|
||||||
@ -722,7 +578,7 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
do p1 = 1, mo_num
|
do p1 = 1, mo_num
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -754,7 +610,7 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
do p1 = 1, mo_num
|
do p1 = 1, mo_num
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_3d(p2,p1,h2)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,p1,h2)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -828,7 +684,7 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END PARALLEL DO
|
!$OMP END PARALLEL DO
|
||||||
@ -868,7 +724,7 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END PARALLEL DO
|
!$OMP END PARALLEL DO
|
||||||
@ -877,15 +733,12 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
enddo ! h1
|
enddo ! h1
|
||||||
enddo ! i
|
enddo ! i
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
! purely open-shell part
|
! purely open-shell part
|
||||||
if(Ne(2) < Ne(1)) then
|
if(Ne(2) < Ne(1)) then
|
||||||
|
|
||||||
do ii = Ne(2) + 1, Ne(1)
|
do ii = Ne(2) + 1, Ne(1)
|
||||||
i = occ(ii,1)
|
i = occ(ii,1)
|
||||||
|
|
||||||
|
|
||||||
! to avoid tmp(N^4)
|
! to avoid tmp(N^4)
|
||||||
do h1 = 1, mo_num
|
do h1 = 1, mo_num
|
||||||
|
|
||||||
@ -941,7 +794,7 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
do p1 = 1, mo_num
|
do p1 = 1, mo_num
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,h2,p1)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -973,7 +826,7 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
do p1 = 1, mo_num
|
do p1 = 1, mo_num
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_3d(p2,p1,h2)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_3d(p2,p1,h2)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -1047,7 +900,7 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END PARALLEL DO
|
!$OMP END PARALLEL DO
|
||||||
@ -1087,7 +940,7 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p2 = 1, mo_num
|
do p2 = 1, mo_num
|
||||||
no_aaa_contraction(p2,h2,p1,h1) = no_aaa_contraction(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
tmp(p2,h2,p1,h1) = tmp(p2,h2,p1,h1) + tmp_2d(p2,h2)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END PARALLEL DO
|
!$OMP END PARALLEL DO
|
||||||
@ -1097,24 +950,30 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
enddo !i
|
enddo !i
|
||||||
endif
|
endif
|
||||||
|
|
||||||
deallocate(tmp_2d, tmp_3d)
|
deallocate(tmp_2d)
|
||||||
deallocate(tmp1, tmp2, tmp3)
|
deallocate(tmp_3d)
|
||||||
deallocate(tmpval_1, tmpval_2)
|
deallocate(tmp1)
|
||||||
deallocate(tmpvec_1, tmpvec_2, tmpvec_3)
|
deallocate(tmp2)
|
||||||
|
deallocate(tmp3)
|
||||||
|
deallocate(tmpval_1)
|
||||||
|
deallocate(tmpval_2)
|
||||||
|
deallocate(tmpvec_1)
|
||||||
|
deallocate(tmpvec_2)
|
||||||
|
deallocate(tmpvec_3)
|
||||||
|
|
||||||
no_aaa_contraction = -0.5d0 * no_aaa_contraction
|
tmp = -0.5d0 * tmp
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP DEFAULT (NONE) &
|
||||||
!$OMP PRIVATE (h1, h2, p1, p2) &
|
!$OMP PRIVATE (h1, h2, p1, p2) &
|
||||||
!$OMP SHARED (no_aaa_contraction, mo_num)
|
!$OMP SHARED (tmp, mo_num)
|
||||||
|
|
||||||
!$OMP DO
|
!$OMP DO
|
||||||
do h1 = 1, mo_num
|
do h1 = 1, mo_num
|
||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p1 = 1, mo_num
|
do p1 = 1, mo_num
|
||||||
do p2 = p1, mo_num
|
do p2 = p1, mo_num
|
||||||
no_aaa_contraction(p2,h2,p1,h1) -= no_aaa_contraction(p1,h2,p2,h1)
|
tmp(p2,h2,p1,h1) -= tmp(p1,h2,p2,h1)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -1126,7 +985,7 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
do h2 = 1, mo_num
|
do h2 = 1, mo_num
|
||||||
do p1 = 2, mo_num
|
do p1 = 2, mo_num
|
||||||
do p2 = 1, p1-1
|
do p2 = 1, p1-1
|
||||||
no_aaa_contraction(p2,h2,p1,h1) = -no_aaa_contraction(p1,h2,p2,h1)
|
tmp(p2,h2,p1,h1) = -tmp(p1,h2,p2,h1)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -1138,18 +997,35 @@ BEGIN_PROVIDER [ double precision, no_aaa_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
do h2 = h1+1, mo_num
|
do h2 = h1+1, mo_num
|
||||||
do p1 = 2, mo_num
|
do p1 = 2, mo_num
|
||||||
do p2 = 1, p1-1
|
do p2 = 1, p1-1
|
||||||
no_aaa_contraction(p2,h2,p1,h1) *= -1.d0
|
tmp(p2,h2,p1,h1) *= -1.d0
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
call wall_time(wallf)
|
||||||
|
print*,' Wall time for aaa_contraction', wall1-wall0
|
||||||
|
|
||||||
|
normal_two_body_bi_orth += tmp
|
||||||
|
endif ! Ne(2) .ge. 3
|
||||||
|
|
||||||
|
deallocate(tmp)
|
||||||
|
|
||||||
|
endif ! read_tc_norm_ord
|
||||||
|
|
||||||
|
if(write_tc_norm_ord.and.mpi_master) then
|
||||||
|
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/normal_two_body_bi_orth', action="write")
|
||||||
|
call ezfio_set_work_empty(.False.)
|
||||||
|
write(11) normal_two_body_bi_orth
|
||||||
|
close(11)
|
||||||
|
call ezfio_set_tc_keywords_io_tc_integ('Read')
|
||||||
endif
|
endif
|
||||||
|
|
||||||
call wall_time(wall1)
|
call wall_time(wallf)
|
||||||
print*,' Wall time for no_aaa_contraction', wall1-wall0
|
print*,' Wall time for normal_two_body_bi_orth ', wallf-walli
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
|
|||||||
1062
src/tc_bi_ortho/normal_ordered_contractions.irp.f
Normal file
1062
src/tc_bi_ortho/normal_ordered_contractions.irp.f
Normal file
File diff suppressed because it is too large
Load Diff
Loading…
Reference in New Issue
Block a user