LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-09 15:03:37 +01:00
Code Releases Activity
ee06ddf85e
qp2/src/bi_ort_ints/three_body_ijmk.irp.f
Abdallah Ammar a791a28523 working on memory footprint
2023-06-04 09:19:34 +02:00

485 lines
18 KiB
Fortran
Raw Blame History

! ---
BEGIN_PROVIDER [ double precision, three_e_4_idx_direct_bi_ort , (mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_4_idx_exch13_bi_ort , (mo_num, mo_num, mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_1_bi_ort, (mo_num, mo_num, mo_num, mo_num)]
!&BEGIN_PROVIDER [ double precision, three_e_4_idx_exch12_bi_ort , (mo_num, mo_num, mo_num, mo_num)]
!&BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_2_bi_ort, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_4_idx_direct_bi_ort (m,j,k,i) = < m j k | -L | m j i > ::: notice that i is the RIGHT MO and k is the LEFT MO
! three_e_4_idx_exch13_bi_ort (m,j,k,i) = < m j k | -L | i j m > ::: notice that i is the RIGHT MO and k is the LEFT MO
! three_e_4_idx_exch12_bi_ort (m,j,k,i) = < m j k | -L | m i j > ::: notice that i is the RIGHT MO and k is the LEFT MO
! = three_e_4_idx_exch13_bi_ort (j,m,k,i)
! three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = < m j k | -L | j i m > ::: notice that i is the RIGHT MO and k is the LEFT MO
! three_e_4_idx_cycle_2_bi_ort(m,j,k,i) = < m j k | -L | i m j > ::: notice that i is the RIGHT MO and k is the LEFT MO
! = three_e_4_idx_cycle_1_bi_ort(j,m,k,i)
!
! notice the -1 sign: in this way three_e_4_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
! three_e_4_idx_direct_bi_ort (m,j,k,i) : Lk Ri Imm Ijj + Lj Rj Imm Iki + Lm Rm Ijj Iki
! three_e_4_idx_exch13_bi_ort (m,j,k,i) : Lk Rm Imi Ijj + Lj Rj Imi Ikm + Lm Ri Ijj Ikm
! three_e_4_idx_cycle_1_bi_ort(m,j,k,i) : Lk Rm Imj Iji + Lj Ri Imj Ikm + Lm Rj Iji Ikm
!
END_DOC
implicit none
integer :: ipoint, i, j, k, l, m
double precision :: wall1, wall0
double precision, allocatable :: tmp1(:,:,:,:), tmp2(:,:,:,:), tmp3(:,:,:,:)
double precision, allocatable :: tmp_4d(:,:,:,:)
double precision, allocatable :: tmp4(:,:,:)
double precision, allocatable :: tmp5(:,:)
double precision, allocatable :: tmp_3d(:,:,:)
print *, ' Providing the three_e_4_idx_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
allocate(tmp_4d(mo_num,mo_num,mo_num,mo_num))
allocate(tmp1(n_points_final_grid,3,mo_num,mo_num))
allocate(tmp2(n_points_final_grid,3,mo_num,mo_num))
allocate(tmp3(n_points_final_grid,3,mo_num,mo_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, l, ipoint) &
!$OMP SHARED (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 tmp1, tmp2, tmp3)
!$OMP DO COLLAPSE(2)
do i = 1, mo_num
do l = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp1(ipoint,1,l,i) = int2_grad1_u12_bimo_t(ipoint,1,l,l) * mos_l_in_r_array_transp(ipoint,i) * final_weight_at_r_vector(ipoint)
tmp1(ipoint,2,l,i) = int2_grad1_u12_bimo_t(ipoint,2,l,l) * mos_l_in_r_array_transp(ipoint,i) * final_weight_at_r_vector(ipoint)
tmp1(ipoint,3,l,i) = int2_grad1_u12_bimo_t(ipoint,3,l,l) * mos_l_in_r_array_transp(ipoint,i) * final_weight_at_r_vector(ipoint)
tmp2(ipoint,1,l,i) = int2_grad1_u12_bimo_t(ipoint,1,l,l) * mos_r_in_r_array_transp(ipoint,i)
tmp2(ipoint,2,l,i) = int2_grad1_u12_bimo_t(ipoint,2,l,l) * mos_r_in_r_array_transp(ipoint,i)
tmp2(ipoint,3,l,i) = int2_grad1_u12_bimo_t(ipoint,3,l,l) * mos_r_in_r_array_transp(ipoint,i)
tmp3(ipoint,1,l,i) = int2_grad1_u12_bimo_t(ipoint,1,l,i) * mos_r_in_r_array_transp(ipoint,l)
tmp3(ipoint,2,l,i) = int2_grad1_u12_bimo_t(ipoint,2,l,i) * mos_r_in_r_array_transp(ipoint,l)
tmp3(ipoint,3,l,i) = int2_grad1_u12_bimo_t(ipoint,3,l,i) * mos_r_in_r_array_transp(ipoint,l)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( 'T', 'N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0 &
, tmp1, 3*n_points_final_grid, tmp2, 3*n_points_final_grid &
, 0.d0, tmp_4d, mo_num*mo_num)
!$OMP PARALLEL DO PRIVATE(i,j,k,m)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
three_e_4_idx_direct_bi_ort(m,j,k,i) = -tmp_4d(m,k,j,i)
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 &
, tmp3, 3*n_points_final_grid, tmp1, 3*n_points_final_grid &
, 0.d0, tmp_4d, mo_num*mo_num)
deallocate(tmp1)
!$OMP PARALLEL DO PRIVATE(i,j,k,m)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
three_e_4_idx_exch13_bi_ort(m,j,k,i) = -tmp_4d(m,i,j,k)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, l, ipoint) &
!$OMP SHARED (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 tmp1)
!$OMP DO COLLAPSE(2)
do i = 1, mo_num
do l = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp1(ipoint,1,l,i) = int2_grad1_u12_bimo_t(ipoint,1,i,l) * mos_l_in_r_array_transp(ipoint,l) * final_weight_at_r_vector(ipoint)
tmp1(ipoint,2,l,i) = int2_grad1_u12_bimo_t(ipoint,2,i,l) * mos_l_in_r_array_transp(ipoint,l) * final_weight_at_r_vector(ipoint)
tmp1(ipoint,3,l,i) = int2_grad1_u12_bimo_t(ipoint,3,i,l) * mos_l_in_r_array_transp(ipoint,l) * final_weight_at_r_vector(ipoint)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( 'T', 'N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0 &
, tmp1, 3*n_points_final_grid, tmp2, 3*n_points_final_grid &
, 0.d0, tmp_4d, mo_num*mo_num)
deallocate(tmp2)
!$OMP PARALLEL DO PRIVATE(i,j,k,m)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
three_e_4_idx_exch13_bi_ort(m,j,k,i) = three_e_4_idx_exch13_bi_ort(m,j,k,i) - tmp_4d(m,k,j,i)
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 &
, tmp1, 3*n_points_final_grid, tmp3, 3*n_points_final_grid &
, 0.d0, tmp_4d, mo_num*mo_num)
deallocate(tmp3)
!$OMP PARALLEL DO PRIVATE(i,j,k,m)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = -tmp_4d(m,k,j,i)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, l, ipoint) &
!$OMP SHARED (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 tmp1)
!$OMP DO COLLAPSE(2)
do i = 1, mo_num
do l = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp1(ipoint,1,l,i) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,l,l) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
tmp1(ipoint,2,l,i) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,l,l) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
tmp1(ipoint,3,l,i) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,l,l) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( 'T', 'N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0 &
, tmp1, 3*n_points_final_grid, int2_grad1_u12_bimo_t, 3*n_points_final_grid &
, 0.d0, tmp_4d, mo_num*mo_num)
deallocate(tmp1)
!$OMP PARALLEL DO PRIVATE(i,j,k,m)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
three_e_4_idx_direct_bi_ort(m,j,k,i) = three_e_4_idx_direct_bi_ort(m,j,k,i) - tmp_4d(m,j,k,i) - tmp_4d(j,m,k,i)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
deallocate(tmp_4d)
allocate(tmp_3d(mo_num,mo_num,mo_num))
allocate(tmp5(n_points_final_grid,mo_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, ipoint) &
!$OMP SHARED (mo_num, n_points_final_grid, &
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
!$OMP final_weight_at_r_vector, &
!$OMP tmp5)
!$OMP DO
do i = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp5(ipoint,i) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,i)
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
allocate(tmp4(n_points_final_grid,mo_num,mo_num))
do m = 1, mo_num
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, k, ipoint) &
!$OMP SHARED (mo_num, n_points_final_grid, m, &
!$OMP int2_grad1_u12_bimo_t, &
!$OMP tmp4)
!$OMP DO COLLAPSE(2)
do i = 1, mo_num
do k = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp4(ipoint,k,i) = int2_grad1_u12_bimo_t(ipoint,1,k,m) * int2_grad1_u12_bimo_t(ipoint,1,m,i) &
+ int2_grad1_u12_bimo_t(ipoint,2,k,m) * int2_grad1_u12_bimo_t(ipoint,2,m,i) &
+ int2_grad1_u12_bimo_t(ipoint,3,k,m) * int2_grad1_u12_bimo_t(ipoint,3,m,i)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( 'T', 'N', mo_num, mo_num*mo_num, n_points_final_grid, 1.d0 &
, tmp5, n_points_final_grid, tmp4, n_points_final_grid &
, 0.d0, tmp_3d, mo_num)
!$OMP PARALLEL DO PRIVATE(i,j,k)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
three_e_4_idx_exch13_bi_ort(m,j,k,i) = three_e_4_idx_exch13_bi_ort(m,j,k,i) - tmp_3d(j,k,i)
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j, k, ipoint) &
!$OMP SHARED (mo_num, n_points_final_grid, m, &
!$OMP mos_l_in_r_array_transp, &
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
!$OMP tmp4)
!$OMP DO COLLAPSE(2)
do k = 1, mo_num
do j = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp4(ipoint,j,k) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,j) &
* ( int2_grad1_u12_bimo_t(ipoint,1,m,j) * int2_grad1_u12_bimo_t(ipoint,1,k,m) &
+ int2_grad1_u12_bimo_t(ipoint,2,m,j) * int2_grad1_u12_bimo_t(ipoint,2,k,m) &
+ int2_grad1_u12_bimo_t(ipoint,3,m,j) * int2_grad1_u12_bimo_t(ipoint,3,k,m) )
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( 'T', 'N', mo_num*mo_num, mo_num, n_points_final_grid, 1.d0 &
, tmp4, n_points_final_grid, mos_r_in_r_array_transp, n_points_final_grid &
, 0.d0, tmp_3d, mo_num*mo_num)
!$OMP PARALLEL DO PRIVATE(i,j,k)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
three_e_4_idx_cycle_1_bi_ort(m,j,k,i) = three_e_4_idx_cycle_1_bi_ort(m,j,k,i) - tmp_3d(j,k,i)
enddo
enddo
enddo
!$OMP END PARALLEL DO
enddo
deallocate(tmp5)
deallocate(tmp_3d)
do i = 1, mo_num
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (m, j, ipoint) &
!$OMP SHARED (mo_num, n_points_final_grid, i, &
!$OMP mos_r_in_r_array_transp, &
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
!$OMP tmp4)
!$OMP DO COLLAPSE(2)
do j = 1, mo_num
do m = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp4(ipoint,m,j) = final_weight_at_r_vector(ipoint) * mos_r_in_r_array_transp(ipoint,m) &
* ( int2_grad1_u12_bimo_t(ipoint,1,m,j) * int2_grad1_u12_bimo_t(ipoint,1,j,i) &
+ int2_grad1_u12_bimo_t(ipoint,2,m,j) * int2_grad1_u12_bimo_t(ipoint,2,j,i) &
+ int2_grad1_u12_bimo_t(ipoint,3,m,j) * int2_grad1_u12_bimo_t(ipoint,3,j,i) )
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( 'T', 'N', mo_num*mo_num, mo_num, n_points_final_grid, -1.d0 &
, tmp4, n_points_final_grid, mos_l_in_r_array_transp, n_points_final_grid &
, 1.d0, three_e_4_idx_cycle_1_bi_ort(1,1,1,i), mo_num*mo_num)
enddo
deallocate(tmp4)
! !$OMP PARALLEL DO PRIVATE(i,j,k,m)
! do i = 1, mo_num
! do k = 1, mo_num
! do j = 1, mo_num
! do m = 1, mo_num
! three_e_4_idx_exch12_bi_ort (m,j,k,i) = three_e_4_idx_exch13_bi_ort (j,m,k,i)
! three_e_4_idx_cycle_2_bi_ort(m,j,k,i) = three_e_4_idx_cycle_1_bi_ort(j,m,k,i)
! enddo
! enddo
! enddo
! enddo
! !$OMP END PARALLEL DO
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_bi_ort', wall1 - wall0
call print_memory_usage()
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, three_e_4_idx_exch23_bi_ort , (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
!
! matrix element of the -L three-body operator FOR THE DIRECT TERMS OF SINGLE EXCITATIONS AND BI ORTHO MOs
!
! three_e_4_idx_exch23_bi_ort (m,j,k,i) = < m j k | -L | j m i > ::: notice that i is the RIGHT MO and k is the LEFT MO
!
! notice the -1 sign: in this way three_e_4_idx_direct_bi_ort can be directly used to compute Slater rules with a + sign
!
! three_e_4_idx_exch23_bi_ort (m,j,k,i) : Lk Ri Imj Ijm + Lj Rm Imj Iki + Lm Rj Ijm Iki
!
END_DOC
implicit none
integer :: i, j, k, l, m, ipoint
double precision :: wall1, wall0
double precision, allocatable :: tmp1(:,:,:,:), tmp_4d(:,:,:,:)
double precision, allocatable :: tmp5(:,:,:), tmp6(:,:,:)
print *, ' Providing the three_e_4_idx_exch23_bi_ort ...'
call wall_time(wall0)
provide mos_r_in_r_array_transp mos_l_in_r_array_transp
allocate(tmp5(n_points_final_grid,mo_num,mo_num))
allocate(tmp6(n_points_final_grid,mo_num,mo_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, l, ipoint) &
!$OMP SHARED (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 tmp5, tmp6)
!$OMP DO COLLAPSE(2)
do i = 1, mo_num
do l = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp5(ipoint,l,i) = int2_grad1_u12_bimo_t(ipoint,1,l,i) * int2_grad1_u12_bimo_t(ipoint,1,i,l) &
+ int2_grad1_u12_bimo_t(ipoint,2,l,i) * int2_grad1_u12_bimo_t(ipoint,2,i,l) &
+ int2_grad1_u12_bimo_t(ipoint,3,l,i) * int2_grad1_u12_bimo_t(ipoint,3,i,l)
tmp6(ipoint,l,i) = final_weight_at_r_vector(ipoint) * mos_l_in_r_array_transp(ipoint,l) * mos_r_in_r_array_transp(ipoint,i)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( 'T', 'N', mo_num*mo_num, mo_num*mo_num, n_points_final_grid, -1.d0 &
, tmp5, n_points_final_grid, tmp6, n_points_final_grid &
, 0.d0, three_e_4_idx_exch23_bi_ort, mo_num*mo_num)
deallocate(tmp5)
deallocate(tmp6)
allocate(tmp_4d(mo_num,mo_num,mo_num,mo_num))
allocate(tmp1(n_points_final_grid,3,mo_num,mo_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, l, ipoint) &
!$OMP SHARED (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 tmp1)
!$OMP DO COLLAPSE(2)
do i = 1, mo_num
do l = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp1(ipoint,1,l,i) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,l,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,l)
tmp1(ipoint,2,l,i) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,l,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,l)
tmp1(ipoint,3,l,i) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,l,i) * mos_l_in_r_array_transp(ipoint,i) * mos_r_in_r_array_transp(ipoint,l)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( 'T', 'N', mo_num*mo_num, mo_num*mo_num, 3*n_points_final_grid, 1.d0 &
, tmp1, 3*n_points_final_grid, int2_grad1_u12_bimo_t, 3*n_points_final_grid &
, 0.d0, tmp_4d, mo_num*mo_num)
deallocate(tmp1)
!$OMP PARALLEL DO PRIVATE(i,j,k,m)
do i = 1, mo_num
do k = 1, mo_num
do j = 1, mo_num
do m = 1, mo_num
three_e_4_idx_exch23_bi_ort(m,j,k,i) = three_e_4_idx_exch23_bi_ort(m,j,k,i) - tmp_4d(m,j,k,i) - tmp_4d(j,m,k,i)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
deallocate(tmp_4d)
call wall_time(wall1)
print *, ' wall time for three_e_4_idx_exch23_bi_ort', wall1 - wall0
call print_memory_usage()
END_PROVIDER
! ---
View Git Blame Copy Permalink