mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-22 03:23:29 +01:00
normal ordering: aba-DGEMM OK
This commit is contained in:
Abdallah Ammar
parent
34ea28d462
commit
4712836349
@ -22,8 +22,6 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
|
|||||||
print*,' Providing normal_two_body_bi_orth ...'
|
print*,' Providing normal_two_body_bi_orth ...'
|
||||||
call wall_time(wall0)
|
call wall_time(wall0)
|
||||||
|
|
||||||
PROVIDE N_int
|
|
||||||
|
|
||||||
if(read_tc_norm_ord) then
|
if(read_tc_norm_ord) then
|
||||||
|
|
||||||
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/normal_two_body_bi_orth', action="read")
|
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/normal_two_body_bi_orth', action="read")
|
||||||
@ -48,12 +46,13 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
|
|||||||
endif
|
endif
|
||||||
|
|
||||||
! opposite spin double excitations : s1 /= s2
|
! opposite spin double excitations : s1 /= s2
|
||||||
normal_two_body_bi_orth(:,:,:,:) = no_aba_contraction(:,:,:,:)
|
PROVIDE no_aba_contraction
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP DEFAULT (NONE) &
|
||||||
!$OMP PRIVATE (hh1, h1, hh2, h2, pp1, p1, pp2, p2, hthree_aab, hthree_aaa) &
|
!$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 SHARED (N_int, n_act_orb, list_act, Ne, occ, normal_two_body_bi_orth, &
|
||||||
|
!$OMP no_aba_contraction)
|
||||||
!$OMP DO SCHEDULE (static)
|
!$OMP DO SCHEDULE (static)
|
||||||
do hh1 = 1, n_act_orb
|
do hh1 = 1, n_act_orb
|
||||||
h1 = list_act(hh1)
|
h1 = list_act(hh1)
|
||||||
@ -97,7 +96,7 @@ BEGIN_PROVIDER [ double precision, normal_two_body_bi_orth, (mo_num, mo_num, mo_
|
|||||||
|
|
||||||
endif
|
endif
|
||||||
|
|
||||||
normal_two_body_bi_orth(p2,h2,p1,h1) = 0.5d0*(hthree_aab + hthree_aaa)
|
normal_two_body_bi_orth(p2,h2,p1,h1) = no_aba_contraction(p2,h2,p1,h1) + 0.5d0*(hthree_aab + hthree_aaa)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
@ -124,103 +123,103 @@ END_PROVIDER
|
|||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
|
!subroutine give_aaa_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
|
||||||
|
!
|
||||||
BEGIN_DOC
|
! BEGIN_DOC
|
||||||
! pure same spin contribution to same spin double excitation s1=h1,p1, s2=h2,p2, with s1==s2
|
! ! pure same spin contribution to same spin double excitation s1=h1,p1, s2=h2,p2, with s1==s2
|
||||||
END_DOC
|
! END_DOC
|
||||||
|
!
|
||||||
use bitmasks ! you need to include the bitmasks_module.f90 features
|
! use bitmasks ! you need to include the bitmasks_module.f90 features
|
||||||
|
!
|
||||||
implicit none
|
! implicit none
|
||||||
integer, intent(in) :: Nint, h1, h2, p1, p2
|
! integer, intent(in) :: Nint, h1, h2, p1, p2
|
||||||
integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
|
! integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
|
||||||
double precision, intent(out) :: hthree
|
! double precision, intent(out) :: hthree
|
||||||
integer :: ii,i
|
! integer :: ii,i
|
||||||
double precision :: int_direct,int_exc_12,int_exc_13,int_exc_23
|
! double precision :: int_direct,int_exc_12,int_exc_13,int_exc_23
|
||||||
double precision :: integral,int_exc_l,int_exc_ll
|
! double precision :: integral,int_exc_l,int_exc_ll
|
||||||
|
!
|
||||||
hthree = 0.d0
|
! hthree = 0.d0
|
||||||
do ii = 1, Ne(2) ! purely closed shell part
|
! do ii = 1, Ne(2) ! purely closed shell part
|
||||||
i = occ(ii,2)
|
! i = occ(ii,2)
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
|
||||||
int_direct = -1.d0 * integral
|
! int_direct = -1.d0 * integral
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(p2, p1, i, i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(p2, p1, i, i, h2, h1, integral)
|
||||||
int_exc_l = -1.d0 * integral
|
! int_exc_l = -1.d0 * integral
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral)
|
||||||
int_exc_ll= -1.d0 * integral
|
! int_exc_ll= -1.d0 * integral
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
|
||||||
int_exc_12= -1.d0 * integral
|
! int_exc_12= -1.d0 * integral
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
|
||||||
int_exc_13= -1.d0 * integral
|
! int_exc_13= -1.d0 * integral
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral)
|
||||||
int_exc_23= -1.d0 * integral
|
! int_exc_23= -1.d0 * integral
|
||||||
|
!
|
||||||
hthree += 1.d0 * int_direct + int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23)
|
! hthree += 1.d0 * int_direct + int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23)
|
||||||
enddo
|
! enddo
|
||||||
|
!
|
||||||
do ii = Ne(2)+1,Ne(1) ! purely open-shell part
|
! do ii = Ne(2)+1,Ne(1) ! purely open-shell part
|
||||||
i = occ(ii,1)
|
! i = occ(ii,1)
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
|
||||||
int_direct = -1.d0 * integral
|
! int_direct = -1.d0 * integral
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(p2, p1, i , i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(p2, p1, i , i, h2, h1, integral)
|
||||||
int_exc_l = -1.d0 * integral
|
! int_exc_l = -1.d0 * integral
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(p1, i, p2, i, h2, h1, integral)
|
||||||
int_exc_ll = -1.d0 * integral
|
! int_exc_ll = -1.d0 * integral
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
|
||||||
int_exc_12 = -1.d0 * integral
|
! int_exc_12 = -1.d0 * integral
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
|
||||||
int_exc_13 = -1.d0 * integral
|
! int_exc_13 = -1.d0 * integral
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral)
|
! call give_integrals_3_body_bi_ort(i, p1, p2, i, h2, h1, integral)
|
||||||
int_exc_23 = -1.d0 * integral
|
! int_exc_23 = -1.d0 * integral
|
||||||
|
!
|
||||||
hthree += 1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23))
|
! hthree += 1.d0 * int_direct + 0.5d0 * (int_exc_l + int_exc_ll - (int_exc_12 + int_exc_13 + int_exc_23))
|
||||||
enddo
|
! enddo
|
||||||
|
!
|
||||||
return
|
! return
|
||||||
end
|
!end
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
subroutine give_aab_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
|
!subroutine give_aab_contraction(Nint, h1, h2, p1, p2, Ne, occ, hthree)
|
||||||
|
!
|
||||||
use bitmasks ! you need to include the bitmasks_module.f90 features
|
! use bitmasks ! you need to include the bitmasks_module.f90 features
|
||||||
|
!
|
||||||
implicit none
|
! implicit none
|
||||||
integer, intent(in) :: Nint, h1, h2, p1, p2
|
! integer, intent(in) :: Nint, h1, h2, p1, p2
|
||||||
integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
|
! integer, intent(in) :: Ne(2), occ(Nint*bit_kind_size,2)
|
||||||
double precision, intent(out) :: hthree
|
! double precision, intent(out) :: hthree
|
||||||
integer :: ii, i
|
! integer :: ii, i
|
||||||
double precision :: int_direct, int_exc_12, int_exc_13, int_exc_23
|
! double precision :: int_direct, int_exc_12, int_exc_13, int_exc_23
|
||||||
double precision :: integral, int_exc_l, int_exc_ll
|
! double precision :: integral, int_exc_l, int_exc_ll
|
||||||
|
!
|
||||||
hthree = 0.d0
|
! hthree = 0.d0
|
||||||
do ii = 1, Ne(2) ! purely closed shell part
|
! do ii = 1, Ne(2) ! purely closed shell part
|
||||||
i = occ(ii,2)
|
! i = occ(ii,2)
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(p2, p1, i, h2, h1, i, integral)
|
! call give_integrals_3_body_bi_ort(p2, p1, i, h2, h1, i, integral)
|
||||||
int_direct = -1.d0 * integral
|
! int_direct = -1.d0 * integral
|
||||||
|
!
|
||||||
call give_integrals_3_body_bi_ort(p1, p2, i, h2, h1, i, integral)
|
! call give_integrals_3_body_bi_ort(p1, p2, i, h2, h1, i, integral)
|
||||||
int_exc_23= -1.d0 * integral
|
! int_exc_23= -1.d0 * integral
|
||||||
|
!
|
||||||
hthree += 1.d0 * int_direct - int_exc_23
|
! hthree += 1.d0 * int_direct - int_exc_23
|
||||||
enddo
|
! enddo
|
||||||
|
!
|
||||||
return
|
! return
|
||||||
end
|
!end
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
@ -264,6 +263,10 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
allocate(tmpvec_1(n_points_final_grid,3))
|
allocate(tmpvec_1(n_points_final_grid,3))
|
||||||
allocate(tmpvec_2(n_points_final_grid,3))
|
allocate(tmpvec_2(n_points_final_grid,3))
|
||||||
|
|
||||||
|
double precision, allocatable :: tmp_2d(:,:)
|
||||||
|
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)
|
||||||
@ -313,9 +316,10 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 &
|
call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 1.d0 &
|
||||||
, int2_grad1_u12_bimo_t, 3*n_points_final_grid, tmp1, 3*n_points_final_grid &
|
, int2_grad1_u12_bimo_t(1,1,1,1), 3*n_points_final_grid &
|
||||||
, 0.d0, tmp_3d, mo_num)
|
, tmp1(1,1,1), 3*n_points_final_grid &
|
||||||
|
, 0.d0, tmp_3d(1,1,1), mo_num*mo_num)
|
||||||
|
|
||||||
!$OMP PARALLEL DO PRIVATE(p1,h2,p2)
|
!$OMP PARALLEL DO PRIVATE(p1,h2,p2)
|
||||||
do p1 = 1, mo_num
|
do p1 = 1, mo_num
|
||||||
@ -364,38 +368,163 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 &
|
call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 1.d0 &
|
||||||
, mos_l_in_r_array_transp, n_points_final_grid, tmp2, n_points_final_grid &
|
, mos_l_in_r_array_transp(1,1), n_points_final_grid &
|
||||||
, 1.d0, no_aba_contraction(p2,h2,1,1), mo_num*mo_num)
|
, tmp2(1,1), n_points_final_grid &
|
||||||
|
, 0.d0, tmp_2d(1,1), mo_num)
|
||||||
|
|
||||||
|
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
||||||
|
do h2 = 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)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END PARALLEL DO
|
||||||
|
|
||||||
enddo ! p1
|
enddo ! p1
|
||||||
enddo ! h1
|
enddo ! h1
|
||||||
enddo ! i
|
enddo ! i
|
||||||
|
|
||||||
|
|
||||||
double precision :: integral, int_direct, int_exc_13, int_exc_12
|
|
||||||
|
|
||||||
! TODO
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
! 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)
|
||||||
|
|
||||||
call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
|
do h1 = 1, mo_num
|
||||||
int_direct = -1.d0 * integral
|
|
||||||
|
|
||||||
call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
|
!$OMP PARALLEL &
|
||||||
int_exc_13 = -1.d0 * integral
|
!$OMP DEFAULT (NONE) &
|
||||||
|
!$OMP PRIVATE (ipoint) &
|
||||||
|
!$OMP SHARED (n_points_final_grid, i, h1, &
|
||||||
|
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||||
|
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
|
||||||
|
!$OMP tmpval_1, tmpval_2, tmpvec_1, tmpvec_2)
|
||||||
|
!$OMP DO
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmpval_1(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint, i)
|
||||||
|
tmpval_2(ipoint) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,h1)
|
||||||
|
tmpvec_1(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i, i) * mos_r_in_r_array_transp(ipoint,h1)
|
||||||
|
tmpvec_1(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i, i) * mos_r_in_r_array_transp(ipoint,h1)
|
||||||
|
tmpvec_1(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i, i) * mos_r_in_r_array_transp(ipoint,h1)
|
||||||
|
tmpvec_2(ipoint,1) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,i,h1) * mos_r_in_r_array_transp(ipoint, i)
|
||||||
|
tmpvec_2(ipoint,2) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,i,h1) * mos_r_in_r_array_transp(ipoint, i)
|
||||||
|
tmpvec_2(ipoint,3) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,i,h1) * mos_r_in_r_array_transp(ipoint, i)
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
|
!$OMP PARALLEL &
|
||||||
int_exc_12 = -1.d0 * integral
|
!$OMP DEFAULT (NONE) &
|
||||||
|
!$OMP PRIVATE (p1, ipoint) &
|
||||||
|
!$OMP SHARED (mo_num, n_points_final_grid, h1, i, &
|
||||||
|
!$OMP mos_l_in_r_array_transp, int2_grad1_u12_bimo_t, &
|
||||||
|
!$OMP tmpval_1, tmpval_2, tmpvec_1, tmpvec_2, tmp1)
|
||||||
|
!$OMP DO
|
||||||
|
do p1 = 1, mo_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmp1(ipoint,1,p1) = mos_l_in_r_array_transp(ipoint,p1) * (tmpvec_1(ipoint,1) - tmpvec_2(ipoint,1)) &
|
||||||
|
+ tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,h1) - tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,p1,i)
|
||||||
|
tmp1(ipoint,2,p1) = mos_l_in_r_array_transp(ipoint,p1) * (tmpvec_1(ipoint,2) - tmpvec_2(ipoint,2)) &
|
||||||
|
+ tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,h1) - tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,p1,i)
|
||||||
|
tmp1(ipoint,3,p1) = mos_l_in_r_array_transp(ipoint,p1) * (tmpvec_1(ipoint,3) - tmpvec_2(ipoint,3)) &
|
||||||
|
+ tmpval_1(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,h1) - tmpval_2(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,p1,i)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
no_aba_contraction(p2,h2,p1,h1) += 1.d0 * int_direct - 0.5d0 * (int_exc_13 + int_exc_12)
|
call dgemm( 'T', 'N', mo_num*mo_num, mo_num, 3*n_points_final_grid, 0.5d0 &
|
||||||
enddo
|
, int2_grad1_u12_bimo_t(1,1,1,1), 3*n_points_final_grid &
|
||||||
|
, tmp1(1,1,1), 3*n_points_final_grid &
|
||||||
|
, 0.d0, tmp_3d(1,1,1), mo_num*mo_num)
|
||||||
|
|
||||||
|
!$OMP PARALLEL DO PRIVATE(p1,h2,p2)
|
||||||
|
do p1 = 1, mo_num
|
||||||
|
do h2 = 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)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END PARALLEL DO
|
||||||
|
|
||||||
|
do p1 = 1, mo_num
|
||||||
|
|
||||||
|
! to minimize the number of operations
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT (NONE) &
|
||||||
|
!$OMP PRIVATE (ipoint) &
|
||||||
|
!$OMP SHARED (n_points_final_grid, i, h1, p1, &
|
||||||
|
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
|
||||||
|
!$OMP tmpval_1)
|
||||||
|
!$OMP DO
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmpval_1(ipoint) = final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1, i,i) * int2_grad1_u12_bimo_t(ipoint,1,p1,h1) &
|
||||||
|
+ int2_grad1_u12_bimo_t(ipoint,2, i,i) * int2_grad1_u12_bimo_t(ipoint,2,p1,h1) &
|
||||||
|
+ int2_grad1_u12_bimo_t(ipoint,3, i,i) * int2_grad1_u12_bimo_t(ipoint,3,p1,h1) &
|
||||||
|
- int2_grad1_u12_bimo_t(ipoint,1,p1,i) * int2_grad1_u12_bimo_t(ipoint,1, i,h1) &
|
||||||
|
- int2_grad1_u12_bimo_t(ipoint,2,p1,i) * int2_grad1_u12_bimo_t(ipoint,2, i,h1) &
|
||||||
|
- int2_grad1_u12_bimo_t(ipoint,3,p1,i) * int2_grad1_u12_bimo_t(ipoint,3, i,h1) )
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT (NONE) &
|
||||||
|
!$OMP PRIVATE (h2, ipoint) &
|
||||||
|
!$OMP SHARED (mo_num, n_points_final_grid, &
|
||||||
|
!$OMP mos_r_in_r_array_transp, &
|
||||||
|
!$OMP tmpval_1, tmp2)
|
||||||
|
!$OMP DO
|
||||||
|
do h2 = 1, mo_num
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
tmp2(ipoint,h2) = mos_r_in_r_array_transp(ipoint,h2) * tmpval_1(ipoint)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
call dgemm( 'T', 'N', mo_num, mo_num, n_points_final_grid, 0.5d0 &
|
||||||
|
, mos_l_in_r_array_transp(1,1), n_points_final_grid &
|
||||||
|
, tmp2(1,1), n_points_final_grid &
|
||||||
|
, 0.d0, tmp_2d(1,1), mo_num)
|
||||||
|
|
||||||
|
!$OMP PARALLEL DO PRIVATE(h2,p2)
|
||||||
|
do h2 = 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)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END PARALLEL DO
|
||||||
|
|
||||||
|
enddo ! p1
|
||||||
|
enddo ! h1
|
||||||
|
enddo !i
|
||||||
endif
|
endif
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
deallocate(tmp_3d)
|
deallocate(tmp_3d)
|
||||||
deallocate(tmp1, tmp2)
|
deallocate(tmp1, tmp2)
|
||||||
@ -403,17 +532,121 @@ BEGIN_PROVIDER [ double precision, no_aba_contraction, (mo_num,mo_num,mo_num,mo_
|
|||||||
deallocate(tmpvec_1, tmpvec_2)
|
deallocate(tmpvec_1, tmpvec_2)
|
||||||
|
|
||||||
|
|
||||||
!$OMP PARALLEL DO PRIVATE(h1,h2,p1,p2)
|
|
||||||
do h1 = 1, mo_num
|
|
||||||
do p1 = 1, mo_num
|
|
||||||
do h2 = 1, mo_num
|
|
||||||
do p2 = 1, mo_num
|
|
||||||
no_aba_contraction(p2,h2,p1,h1) = -0.5d0 * (no_aba_contraction(p2,h2,p1,h1) + no_aba_contraction(p1,h1,p2,h2))
|
|
||||||
enddo
|
|
||||||
enddo
|
no_aba_contraction = -0.5d0 * no_aba_contraction
|
||||||
enddo
|
call sum_A_At(no_aba_contraction(1,1,1,1), mo_num*mo_num)
|
||||||
enddo
|
|
||||||
!$OMP END PARALLEL DO
|
! do h1 = 1, mo_num
|
||||||
|
! do p1 = 1, mo_num
|
||||||
|
! do h2 = 1, mo_num
|
||||||
|
! do p2 = 1, mo_num
|
||||||
|
! no_aba_contraction(p2,h2,p1,h1) = -0.5d0 * (tmp_4d(p2,h2,p1,h1) + tmp_4d(p1,h1,p2,h2))
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
double precision :: integral, int_direct, int_exc_13, int_exc_12
|
||||||
|
|
||||||
|
! no_aba_contraction = 0.d0
|
||||||
|
!
|
||||||
|
! ! purely closed shell part
|
||||||
|
! do ii = 1, Ne(2)
|
||||||
|
! i = occ(ii,1)
|
||||||
|
!
|
||||||
|
! !$OMP PARALLEL &
|
||||||
|
! !$OMP DEFAULT (NONE) &
|
||||||
|
! !$OMP PRIVATE (h1, h2, p1, p2, int_direct, int_exc_13, int_exc_12, integral) &
|
||||||
|
! !$OMP SHARED (mo_num, i, no_aba_contraction)
|
||||||
|
! !$OMP DO SCHEDULE (static)
|
||||||
|
! do h1 = 1, mo_num
|
||||||
|
! do p1 = 1, mo_num
|
||||||
|
! do h2 = 1, mo_num
|
||||||
|
! do p2 = 1, mo_num
|
||||||
|
!
|
||||||
|
! call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
|
||||||
|
! int_direct = -1.d0 * integral
|
||||||
|
!
|
||||||
|
! call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
|
||||||
|
! int_exc_13 = -1.d0 * integral
|
||||||
|
!
|
||||||
|
! call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
|
||||||
|
! int_exc_12 = -1.d0 * integral
|
||||||
|
!
|
||||||
|
! !no_aba_contraction(p2,h2,p1,h1) += 1.d0 * int_direct - 0.5d0 * (int_exc_13 + int_exc_12)
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! !$OMP END DO
|
||||||
|
! !$OMP END PARALLEL
|
||||||
|
! enddo
|
||||||
|
|
||||||
|
! ! purely open-shell part
|
||||||
|
! if(Ne(2) < Ne(1)) then
|
||||||
|
!
|
||||||
|
! do ii = Ne(2) + 1, Ne(1)
|
||||||
|
! i = occ(ii,1)
|
||||||
|
!
|
||||||
|
! !$OMP PARALLEL &
|
||||||
|
! !$OMP DEFAULT (NONE) &
|
||||||
|
! !$OMP PRIVATE (h1, h2, p1, p2, int_direct, int_exc_13, int_exc_12, integral) &
|
||||||
|
! !$OMP SHARED (mo_num, i, no_aba_contraction)
|
||||||
|
! !$OMP DO SCHEDULE (static)
|
||||||
|
! do h1 = 1, mo_num
|
||||||
|
! do p1 = 1, mo_num
|
||||||
|
! do h2 = 1, mo_num
|
||||||
|
! do p2 = 1, mo_num
|
||||||
|
!
|
||||||
|
! call give_integrals_3_body_bi_ort(i, p2, p1, i, h2, h1, integral)
|
||||||
|
! int_direct = -1.d0 * integral
|
||||||
|
!
|
||||||
|
! call give_integrals_3_body_bi_ort(p1, p2, i, i, h2, h1, integral)
|
||||||
|
! int_exc_13 = -1.d0 * integral
|
||||||
|
!
|
||||||
|
! call give_integrals_3_body_bi_ort(p2, i, p1, i, h2, h1, integral)
|
||||||
|
! int_exc_12 = -1.d0 * integral
|
||||||
|
!
|
||||||
|
! no_aba_contraction(p2,h2,p1,h1) += 0.5d0 * int_direct - 0.25d0 * (int_exc_13 + int_exc_12)
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! !$OMP END DO
|
||||||
|
! !$OMP END PARALLEL
|
||||||
|
! enddo
|
||||||
|
! endif
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
! !$OMP PARALLEL &
|
||||||
|
! !$OMP DEFAULT (NONE) &
|
||||||
|
! !$OMP PRIVATE (h1, h2, p1, p2, integral) &
|
||||||
|
! !$OMP SHARED (mo_num, N_int,Ne, occ, no_aba_contraction)
|
||||||
|
! !$OMP DO SCHEDULE (static)
|
||||||
|
! do h1 = 1, mo_num
|
||||||
|
! do p1 = 1, mo_num
|
||||||
|
! do h2 = 1, mo_num
|
||||||
|
! do p2 = 1, mo_num
|
||||||
|
! call give_aba_contraction(N_int, h1, h2, p1, p2, Ne, occ, integral)
|
||||||
|
! no_aba_contraction(p2,h2,p1,h1) = 0.5d0 * integral
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! enddo
|
||||||
|
! !$OMP END DO
|
||||||
|
! !$OMP END PARALLEL
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|||||||
@ -11,12 +11,14 @@ program tc_bi_ortho
|
|||||||
touch read_wf
|
touch read_wf
|
||||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||||
|
|
||||||
call test_h_u0
|
! call test_h_u0
|
||||||
! call test_slater_tc_opt
|
! call test_slater_tc_opt
|
||||||
! call timing_tot
|
! call timing_tot
|
||||||
! call timing_diag
|
! call timing_diag
|
||||||
! call timing_single
|
! call timing_single
|
||||||
! call timing_double
|
! call timing_double
|
||||||
|
|
||||||
|
call test_no()
|
||||||
end
|
end
|
||||||
|
|
||||||
subroutine test_h_u0
|
subroutine test_h_u0
|
||||||
@ -252,3 +254,47 @@ subroutine timing_double
|
|||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
subroutine test_no()
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer :: i, j, k, l
|
||||||
|
double precision :: accu, contrib, new, ref, thr
|
||||||
|
|
||||||
|
print*, ' testing normal_two_body_bi_orth ...'
|
||||||
|
|
||||||
|
thr = 1d-8
|
||||||
|
|
||||||
|
PROVIDE normal_two_body_bi_orth_old
|
||||||
|
PROVIDE normal_two_body_bi_orth
|
||||||
|
|
||||||
|
accu = 0.d0
|
||||||
|
do i = 1, mo_num
|
||||||
|
do j = 1, mo_num
|
||||||
|
do k = 1, mo_num
|
||||||
|
do l = 1, mo_num
|
||||||
|
|
||||||
|
new = normal_two_body_bi_orth (l,k,j,i)
|
||||||
|
ref = normal_two_body_bi_orth_old(l,k,j,i)
|
||||||
|
contrib = dabs(new - ref)
|
||||||
|
accu += contrib
|
||||||
|
if(contrib .gt. thr) then
|
||||||
|
print*, ' problem on normal_two_body_bi_orth'
|
||||||
|
print*, l, k, j, i
|
||||||
|
print*, ref, new, contrib
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
print*, ' accu on normal_two_body_bi_orth = ', accu / dble(mo_num)**4
|
||||||
|
|
||||||
|
return
|
||||||
|
end
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user