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

Optimized cyclic 5idx

Browse Source
This commit is contained in:
Anthony Scemama 2023-06-02 10:34:05 +02:00
parent c4612318ae
commit 00bd8e2fcc
3 changed files with 152 additions and 202 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

65
src/bi_ort_ints/bi_ort_ints.irp.f
View File

@ -59,17 +59,18 @@ subroutine test_5idx
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'direct'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
! new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'direct'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_exch12_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_exch12_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
@ -81,6 +82,50 @@ subroutine test_5idx
! stop
! endif
! new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'cycle1'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
! new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'cycle2'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
! new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'exch23'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
new = three_e_5_idx_exch13_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
accu += contrib
if(contrib .gt. 1.d-10)then
print*,'exch13'
print*,i,k,j,l,m
print*,ref,new,contrib
stop
endif
enddo
enddo
enddo

286
src/bi_ort_ints/three_body_ijmkl.irp.f
View File

@ -1,4 +1,3 @@
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
@ -17,7 +16,6 @@
integer :: i, j, k, m, l
double precision :: wall1, wall0
integer :: ipoint
double precision :: weight
double precision, allocatable :: grad_mli(:,:,:), m2grad_r(:,:,:,:), m2grad_l(:,:,:,:)
double precision, allocatable :: tmp_mat(:,:,:,:), orb_mat(:,:,:)
allocate(m2grad_r(n_points_final_grid,3,mo_num,mo_num))
@ -45,17 +43,22 @@
do i=1,mo_num
do l=1,mo_num
do ipoint=1, n_points_final_grid
grad_mli(ipoint,l,i) = final_weight_at_r_vector(ipoint) * ( &
int2_grad1_u12_bimo_t(ipoint,1,m,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) + &
int2_grad1_u12_bimo_t(ipoint,2,m,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) + &
int2_grad1_u12_bimo_t(ipoint,3,m,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) )
orb_mat(ipoint,l,i) = mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,i)
m2grad_l(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) * final_weight_at_r_vector(ipoint)
m2grad_l(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) * final_weight_at_r_vector(ipoint)
m2grad_l(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) * final_weight_at_r_vector(ipoint)
m2grad_r(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i)
m2grad_r(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i)
m2grad_r(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i)
orb_mat(ipoint,l,i) = mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,i)
enddo
enddo
enddo
@ -108,239 +111,140 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = <mlk|-L|jim> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
double precision :: integral
integer :: i, j, k, m, l
double precision :: wall1, wall0
integer :: ipoint
double precision :: weight
double precision, allocatable :: grad_mli(:,:,:), m2grad_r(:,:,:,:), m2grad_l(:,:,:,:)
double precision, allocatable :: tmp_mat(:,:,:,:), orb_mat(:,:,:)
allocate(m2grad_r(n_points_final_grid,3,mo_num,mo_num))
allocate(m2grad_l(n_points_final_grid,3,mo_num,mo_num))
allocate(tmp_mat(mo_num,mo_num,mo_num,mo_num))
allocate(grad_mli(n_points_final_grid,mo_num,mo_num))
allocate(orb_mat(n_points_final_grid,mo_num,mo_num))
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t
print *, ' Providing the three_e_5_idx_cycle_1_bi_ort ...'
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_cycle_1_bi_ort)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, j, i, m, integral)
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_cycle_1_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE FIRST CYCLIC PERMUTATION TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = <mlk|-L|imj> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
double precision :: integral
integer :: i, j, k, m, l
double precision :: wall1, wall0
integer :: ipoint
double precision :: weight
double precision, allocatable :: grad_mli(:,:,:), m2grad_r(:,:,:,:), m2grad_l(:,:,:,:)
double precision, allocatable :: tmp_mat(:,:,:,:), orb_mat(:,:,:)
allocate(m2grad_r(n_points_final_grid,3,mo_num,mo_num))
allocate(m2grad_l(n_points_final_grid,3,mo_num,mo_num))
allocate(tmp_mat(mo_num,mo_num,mo_num,mo_num))
allocate(grad_mli(n_points_final_grid,mo_num,mo_num))
allocate(orb_mat(n_points_final_grid,mo_num,mo_num))
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t
print *, ' Providing the three_e_5_idx_cycle_2_bi_ort ...'
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_cycle_2_bi_ort)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
do l = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, i, m, j, integral)
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_cycle_2_bi_ort', wall1 - wall0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort , (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = <mlk|-L|jmi> ::: notice that i is the RIGHT MO and k is the LEFT MO
! three_e_5_idx_direct_bi_ort(m,l,j,k,i) = <mlk|-L|mji> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
double precision :: integral
integer :: i, j, k, m, l
double precision :: wall1, wall0
integer :: ipoint
double precision :: weight
double precision, allocatable :: grad_mli(:,:,:), m2grad_r(:,:,:,:), m2grad_l(:,:,:,:)
double precision, allocatable :: tmp_mat(:,:,:,:), orb_mat(:,:,:)
allocate(m2grad_r(n_points_final_grid,3,mo_num,mo_num))
allocate(m2grad_l(n_points_final_grid,3,mo_num,mo_num))
double precision, allocatable :: lk_grad_mi(:,:,:,:), rk_grad_im(:,:,:,:)
double precision, allocatable :: lm_grad_ik(:,:,:,:), rm_grad_ik(:,:,:,:)
double precision, allocatable :: tmp_mat(:,:,:,:)
allocate(lk_grad_mi(n_points_final_grid,3,mo_num,mo_num))
allocate(lm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
allocate(rk_grad_im(n_points_final_grid,3,mo_num,mo_num))
allocate(rm_grad_ik(n_points_final_grid,3,mo_num,mo_num))
allocate(tmp_mat(mo_num,mo_num,mo_num,mo_num))
allocate(grad_mli(n_points_final_grid,mo_num,mo_num))
allocate(orb_mat(n_points_final_grid,mo_num,mo_num))
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t
print *, ' Providing the three_e_5_idx_exch23_bi_ort ...'
print *, ' Providing the three_e_5_idx_cycle_bi_ort ...'
call wall_time(wall0)
do m = 1, mo_num
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch23_bi_ort)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
!$OMP PRIVATE (i,l,ipoint) &
!$OMP SHARED (m,mo_num,n_points_final_grid, &
!$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 rk_grad_im, rm_grad_ik, lk_grad_mi, lm_grad_ik, tmp_mat)
!$OMP DO COLLAPSE(2)
do i=1,mo_num
do k = 1, mo_num
do j = 1, mo_num
do l=1,mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, j, m, i, integral)
three_e_5_idx_exch23_bi_ort(m,l,j,k,i) = -1.d0 * integral
enddo
enddo
do ipoint=1, n_points_final_grid
lk_grad_mi(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,m,i) * final_weight_at_r_vector(ipoint)
lk_grad_mi(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,m,i) * final_weight_at_r_vector(ipoint)
lk_grad_mi(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,m,i) * final_weight_at_r_vector(ipoint)
lm_grad_ik(ipoint,1,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i) * final_weight_at_r_vector(ipoint)
lm_grad_ik(ipoint,2,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i) * final_weight_at_r_vector(ipoint)
lm_grad_ik(ipoint,3,l,i) = mos_l_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i) * final_weight_at_r_vector(ipoint)
rm_grad_ik(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,1,l,i)
rm_grad_ik(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,2,l,i)
rm_grad_ik(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,m) * int2_grad1_u12_bimo_t(ipoint,3,l,i)
rk_grad_im(ipoint,1,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,1,i,m)
rk_grad_im(ipoint,2,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,2,i,m)
rk_grad_im(ipoint,3,l,i) = mos_r_in_r_array_transp(ipoint,l) * int2_grad1_u12_bimo_t(ipoint,3,i,m)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_exch23_bi_ort', wall1 - wall0
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
lk_grad_mi, 3*n_points_final_grid, &
rm_grad_ik, 3*n_points_final_grid, 0.d0, &
tmp_mat, mo_num*mo_num)
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort, (mo_num, mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF DOUBLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = <mlk|-L|ijm> ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_3_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
END_DOC
implicit none
double precision :: integral
integer :: i, j, k, m, l
double precision :: wall1, wall0
integer :: ipoint
double precision :: weight
double precision, allocatable :: grad_mli(:,:,:), m2grad_r(:,:,:,:), m2grad_l(:,:,:,:)
double precision, allocatable :: tmp_mat(:,:,:,:), orb_mat(:,:,:)
allocate(m2grad_r(n_points_final_grid,3,mo_num,mo_num))
allocate(m2grad_l(n_points_final_grid,3,mo_num,mo_num))
allocate(tmp_mat(mo_num,mo_num,mo_num,mo_num))
allocate(grad_mli(n_points_final_grid,mo_num,mo_num))
allocate(orb_mat(n_points_final_grid,mo_num,mo_num))
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
PROVIDE mo_l_coef mo_r_coef int2_grad1_u12_bimo_t
print *, ' Providing the three_e_5_idx_exch13_bi_ort ...'
call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i,j,k,m,l,integral) &
!$OMP SHARED (mo_num,three_e_5_idx_exch13_bi_ort)
!$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
do m = 1, mo_num
call give_integrals_3_body_bi_ort(m, l, k, i, j, m, integral)
three_e_5_idx_exch13_bi_ort(m,l,j,k,i) = -1.d0 * integral
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = -tmp_mat(k,j,l,i)
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = -tmp_mat(l,i,k,j)
three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = -tmp_mat(l,j,k,i)
three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = -tmp_mat(k,i,l,j)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
lk_grad_mi, 3*n_points_final_grid, &
rk_grad_im, 3*n_points_final_grid, 0.d0, &
tmp_mat, mo_num*mo_num)
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(l,j,i,k)
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) - tmp_mat(k,i,j,l)
three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(k,j,i,l)
three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = three_e_5_idx_exch13_bi_ort (m,l,j,k,i) - tmp_mat(l,i,j,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
call dgemm('T','N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0, &
lm_grad_ik, 3*n_points_final_grid, &
rk_grad_im, 3*n_points_final_grid, 0.d0, &
tmp_mat, mo_num*mo_num)
!$OMP PARALLEL DO PRIVATE(i,j,k,l)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do l = 1, mo_num
three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i) - tmp_mat(l,i,j,k)
three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i) - tmp_mat(k,j,i,l)
three_e_5_idx_exch23_bi_ort (m,l,j,k,i) = three_e_5_idx_exch23_bi_ort (m,l,j,k,i) - tmp_mat(k,i,j,l)
three_e_5_idx_exch13_bi_ort (m,l,j,k,i) = three_e_5_idx_exch13_bi_ort (m,l,j,k,i) - tmp_mat(l,j,i,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
enddo
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print *, ' wall time for three_e_5_idx_exch13_bi_ort', wall1 - wall0
print *, ' wall time for three_e_5_idx_cycle_bi_ort', wall1 - wall0
END_PROVIDER
! ---

1
src/bi_ort_ints/three_body_ints_bi_ort.irp.f
View File

@ -85,6 +85,7 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral)
PROVIDE int2_grad1_u12_bimo_t
integral = 0.d0
! (n, l, k, m, j, i)
do ipoint = 1, n_points_final_grid
tmp = mos_l_in_r_array_transp(ipoint,k) * mos_r_in_r_array_transp(ipoint,i) &